Este material tem como objetivo contribuir para o entendimento de Python (notebook jupyter) no uso em escoamento potencial, principalmente para analisar e visualizar na "prática" o escoamento. E assim expressar numericamente o campo de fluxo de um fluxo potencial e que possamos traçar esses fluxos e admirá-los.

Vamos recapitular. Nos primeiros módulos, nós descrevemos alguns conceitos teóricos como:

módulo 00 introdução e revisão - aprendemos que a mecânica dos fluidos é um tema excitante e fascinante, e que a função potencial de velocidade e de corrente são importantíssima na área do escoamento potencial. Vemos também que as linhas equipotenciais são ortogonais às linhas de fluxo.
módulo 01 soluções elementares ... , vemos que podemos modelar o fluxo usando a função potencial de velocidade e ou de corrente, encontradas a partir do campo velocidade. Definimos aqui também, essas funções para as três soluções elementares de escoamento no plano (uniforme, fonte ou sumidouro e vórtice irrotacional).
módulo superposição de soluções ... - este módulo é maravilhoso, porque usamos do princípio da superposição, construindo novas soluções adicionando soluções conhecidas do módulo anterior. E encontramos diversos modelos clássico de fluxo, e espero que você tenha divertido bastante. (Se não, pode confiar em mim, essa parte com o Python é tri legal!)

Mas qual é a coisa mais importante que queremos dessa parte da mecânica dos fluidos? Nós queremos fazer as coisas voarem, é claro! ... hehehehehe... É, para voar, deve haver uma força de sustentação aerodinâmica para contrabalançar o peso do objeto, e isso ainda não será discutido aqui, mas nesta seção aprenderemos como visualizar o fluxo de um escoamento potencial, no plano, com a linguagem Python (notebook jupyter).

Vamos começar!¶

Começamos importando as bibliotecas que usaremos neste material: a biblioteca de matrizes NumPy, a biblioteca álgebrica SymPy, a biblioteca de plotagem Matplotlib e as funções matemáticas no mathmódulo.

NumPy é uma biblioteca científica para criar e gerenciar matrizes e matrizes multidimensionais.
SymPy é uma biblioteca científica com ferramentas de álgebra computacional.
O Matplotlib é uma biblioteca de plotagem 2D que usaremos para visualizar nossos resultados.
o mathmódulo fornece as funções matemáticas definidas pelo padrão C.

Também incluímos o %matplotlib inlinecomando para configurar a biblioteca de plotagem para nos fornecer os gráficos incorporados em nosso bloco de anotações.

import math
import numpy
import sympy
from sympy.abc import x, y
from matplotlib import pyplot
# para plotar figuras do matplotlib no Notebook Jupyter
%matplotlib inline

O objetivo é visualizar as linhas de fluxo correspondentes a um escoamento. Para fazer isso, precisamos primeiro definir um conjunto de pontos onde as componentes da velocidade serão computadas.

1. Geração da malha¶

Vamos definir uma malha cartesiana de pontos uniformemente espaçada dentro de um domínio espacial que tem uma certa unidade de largura no eixo $x$ e também de certa unidade de largura no eixo $y$, isto é, os limites computacionais x_start, x_end, y_start e y_end. E que a variável Nserá o número de pontos que queremos em cada sentido $x$ e $y$.

Usamos a função NumPy linspace()para criar dois arrays 1D que contêm os valores uniformemente espaçados de $x$ e $y$ coordenadas, correspondendo aos nossos pontos da malha. A última linha do bloco de código abaixo chama a meshgrid()função, que gera matrizes contendo as coordenadas dos pontos onde a solução numérica será calculada.

N = 50                                	        # número de pontos em cada direção 

x_inicial, x_final = -2.0, 2.0            	    # limites na direção x 
y_inicial, y_final = -1.0, 1.0            	    # limites na direção y 

x_1D = numpy.linspace(x_inicial, x_final, N)    # cria um array 1D para a coordenada x
y_1D = numpy.linspace(y_inicial, y_final, N)    # cria um array 1D para a coordenada y

X, Y = numpy.meshgrid(x_1D, y_1D)               # gera uma malha

Agora que a grade de malha foi gerada, é hora de visualizá-la com o módulo pyplotda biblioteca matplotlibusando a função scatter().

Para obter os gráficos em nosso bloco de notas, usamos o comando %matplotlib inlineprimeiro. Este comando configura o módulo de plotagem para gerar gráficos em linha, como imagens PNG estáticas, e ele foi declarado no comando das bibliotecas.

# plotando a grade de pontos 
largura = 10.0
altura = (y_final - y_inicial) / (x_final - x_inicial) * largura
pyplot.figure(figsize=(largura, altura))                                # aqui definimos as dimensões da figura
pyplot.xlabel('x', fontsize=16); pyplot.ylabel('y', fontsize=16)        # aqui a nomeclatura dos eixo x e y                                
pyplot.xlim(x_inicial, x_final); pyplot.ylim(y_inicial, y_final)        # aqui definimos os limites de plotagem
pyplot.scatter(X, Y, s=5, color='#CD2305', marker='o')                  # com a função scatter() podemos visualizar a malha

<matplotlib.collections.PathCollection at 0x239fed072b0>

Em todos esses pontos bem ordenados, nós agora calcularemos o vetor de velocidade correspondente a um fluxo do escoamento em estudo. E a partir disso, podemos traçar as linhas de fluxo para esse escoamento. Vamos começar com as soluções elementares e depois entramos na superposição.

Ah, essas funções de fluxo, serão agora trabalhadas em coordenadas cilíndricas.

2. Fonte/Sumidouro na origem¶

Uma fonte/sumidouro é um ponto do qual imaginamos que o fluido está fluindo uniformemente. Assim, todos os streamlines irradiam/entram num único ponto como linhas retas e a velocidade radial diminui com a distância do ponto de origem. Vamos considerar primeiro o caso puramente bidimensional. Por causa da simetria radial, é conveniente usar um sistema de coordenadas cilíndricas, como foi comentado antes, $(r, \theta)$. O ângulo $\theta$ é $\tan^{-1} ( y/ x)$. Os componentes de velocidade (radial e tangencial) são: $$$$ $$ u_{r}(r,\theta) = \frac {Q/b}{2\pi r} \quad u_{\theta}(r,\theta) = 0 $$ $$$$ onde $(+Q/b)$ representa a força da fonte . O fato de a velocidade tangencial ser zero é óbvio em nossa exigência de que as linhas de fluxo irradiam linhas retas. Mas como conseguimos o componente radial da velocidade? Aplicar a condição de fluxo irrotacional, $ω=0$, em coordenadas cilíndricas, e você vai ter que a velocidade só pode ser uma função de $r$. Em seguida, aplique a equação de continuidade e você obterá o resultado. Reveja o módulo 01.

Lembre-se da função de fluxo para fonte, mas em coordenadas cilíndricas, então $\psi$ é obtido de: $$$$ $$ \frac{1}{r}\frac{\partial \psi}{\partial \theta} = u_{r} \quad -\frac{\partial \psi}{\partial r} = u_{\theta} $$ que integra para $$ \psi = \frac{Q/b}{2\pi}\theta + constante $$ $$$$ Em problemas práticos, estamos mais interessados nos componentes de velocidade que são obtidos pela diferenciação da função de fluxo, para que a constante possa ser descartada.

Em coordenadas cartesianas, o campo de velocidade $(u,v)$ na posição $(x,y)$ correspondendo a uma fonte de força $(Q/B)$ localizado em $(x_{fonte}, y_{fonte})$ é dada por: $$$$ $$ u = \frac{\partial \psi}{\partial y} = \frac{Q/b}{2\pi}\frac{x - x_{fonte}}{(x - x_{fonte})^{2} + (y - y_{fonte})^{2}} $$ e $$ v = - \frac{\partial \psi}{\partial x} = \frac{Q/b}{2\pi}\frac{y - y_{fonte}}{(x - x_{fonte})^{2} + (y - y_{fonte})^{2}} $$ $$$$ Vamos calcular o campo de velocidade usando a função Sympy lambdify() para a derivada primeira da função corrente em nossa malha de pontos. Adotando neste exemplo a fonte no local $(-1,0)$ e de força $(+Q/b)= 5$.

Em vez de escolher um ponto na grade e calcular sua velocidade (o que significa que teríamos que iterar sobre todas as posições [i,j]), calculamos diretamente as matrizes de velocidade (u_fonte, v_fonte) usando operadores aritméticos em matrizes. Sim, com a biblioteca Numpy, operadores aritméticos em array se aplicam elemento por elemento, e uma nova matriz é criada e preenchida com o resultado.

Vamos definir uma função chamada funcao_corrente_fonte()para obter a função corrente em coordenadas cilíndrica, e uma outra função, campo_velocidade(), que calcula os componentes de velocidade $(u,v)$, pela função Sympy `lambdify() nas coordenadas X e Y. No Python, uma função é definida com o comando def, seguido pelo nome escolhido para a função e quaisquer parâmetros dentro dos parênteses, e a linha termina com dois pontos:

intensidade_fonte = 5.0                         # intensidade da fonte
x_fonte, y_fonte = -1.0, 0.0                    # localização da fonte

def funcao_corrente_fonte(intensidade_fonte, x_fonte, y_fonte):
    """
    Calcula a função corrente de uma fonte.
    
    Parameters
    ----------
    intensidade_fonte: float
        Intensidade do fluxo da fonte.
    x_fonte: float
        Coordenada-x da posição da fonte.
    y_fonte: float
        Coordenada-y da posição da fonte.
    
    Returns
    -------
    (intensidade_fonte / (2*numpy.pi)) * theta: função simbólica
        função corrente em termos da matemática simbólica.
    """
    r = sympy.sqrt(x**2 + y**2)
    theta = sympy.atan2(y-y_fonte, x-x_fonte)
    return (intensidade_fonte / (2*numpy.pi)) * theta

def campo_velocidade(psi):
    """
    Calcula o campo velocidade.
    
    Parameters
    ----------
    psi: função simbólica.
    
    Returns
    -------
    u: função componete u-velocidade em termos da matemática simbólica,
        em funçao de (X,Y).
    v: função componete u-velocidade em termos da matemática simbólica,
        em funçao de (X,Y). 
    """
    u = sympy.lambdify((x, y), psi.diff(y), 'numpy')
    v = sympy.lambdify((x, y), -psi.diff(x), 'numpy')
    return u, v

# computando o campo velocidade na grade de malha 
psi = funcao_corrente_fonte(intensidade_fonte, x_fonte, y_fonte)
u, v = campo_velocidade(psi)

Vamos traçar as linhas de fluxo. Nós temos sorte que os contribuidores do OpenSource do Python adicionaram uma streamplot()função que faz tudo. Também usaremos a scatter()função para colocar um ponto vermelho na fonte.

# plot as linhas de fluxo da fonte
largura = 10.0
altura = (y_final - y_inicial) / (x_final - x_inicial) * largura
pyplot.figure(figsize=(largura, altura))
pyplot.xlabel('x', fontsize=16); pyplot.ylabel('y', fontsize=16)
pyplot.xlim(x_inicial, x_final); pyplot.ylim(y_inicial, y_final)
pyplot.streamplot(X, Y, u(X,Y), v(X,Y),
                  density=2, linewidth=1, arrowsize=2, arrowstyle='->')
pyplot.scatter(x_fonte, y_fonte,
               color='#CD2305', s=80, marker='o');

Legal, neh? É assim que uma fonte deveria parecer, apenas mais bonita hehehehe. Observe que adicionamos um ponto-e-vírgula na instrução mais recente da célula acima, para suprimir a exibição do feio <matplotlib ...>, que é apenas uma ID interna da imagem.

No fluxo de fonte, a intensidade foi escolhido para ser positivo. Uma fonte com força negativa é chamada de sumidouro. Em vez de irradiar a partir de um único ponto, as linhas de fluxo retas estão agora convergindo para um único ponto.

O campo de velocidade correspondente a um sumidouro é semelhante ao de uma fonte, exceto pela direção do fluxo. Assim, o código Python requer muito poucas modificações.

Vamos colocar o sumidouro no local $(1,0)$ e dar-lhe uma força igual à nossa fonte, mas negativa, claro.

intensidade_sumidouro = -5.0                            # intensidade do sumidouro
x_sumidouro, y_sumidouro = 1.0, 0.0                     # localização do sumidouro

def funcao_corrente_sumidouro(intensidade_sumidouro, x_sumidouro, y_sumidouro):
    """
    Calcula a função corrente de um sumidouro.
    
    Parameters
    ----------
    intensidade_fonte: float
        Intensidade do fluxo do sumidouro.
    x_fonte: float
        Coordenada-x da posição do sumidouro.
    y_fonte: float
        Coordenada-y da posição do sumidouro.
    
    Returns
    -------
    (intensidade_fonte / (2*numpy.pi)) * theta: função simbólica
        função corrente em termos da matemática simbólica.
    """
    r = sympy.sqrt(x**2 + y**2)
    theta = sympy.atan2(y-y_sumidouro, x-x_sumidouro)
    return (intensidade_sumidouro / (2*numpy.pi)) * theta

def campo_velocidade(psi):
    """
    Calcula o campo velocidade.
    
    Parameters
    ----------
    psi: função simbólica.
    
    Returns
    -------
    u: função componete u-velocidade em termos da matemática simbólica,
        em funçao de (X,Y).
    v: função componete u-velocidade em termos da matemática simbólica,
        em funçao de (X,Y). 
    """
    u = sympy.lambdify((x, y), psi.diff(y), 'numpy')
    v = sympy.lambdify((x, y), -psi.diff(x), 'numpy')
    return u, v

# computando o campo velocidade na grade de malha 
psi = funcao_corrente_fonte(intensidade_sumidouro, x_sumidouro, y_sumidouro)
u, v = campo_velocidade(psi)

# plotando as linhas de fluxo do sumidouro
largura = 10.0
altura = (y_final - y_inicial) / (x_final - x_inicial) * largura
pyplot.figure(figsize=(largura, altura))
pyplot.xlabel('x', fontsize=16); pyplot.ylabel('y', fontsize=16)
pyplot.xlim(x_inicial, x_final); pyplot.ylim(y_inicial, y_final)
pyplot.streamplot(X, Y, u(X,Y), v(X,Y),
                  density=2, linewidth=1, arrowsize=2, arrowstyle='->')
pyplot.scatter(x_sumidouro, y_sumidouro,
               color='#CD2305', s=80, marker='o');

3. Escoamento uniforme¶

No tópico anterior, já calculamos o campo de velocidade, e o seu fluxo, de uma fonte e do sumidouro. As linhas de fluxo de um escoamento uniforme, de velocidade $U$ e ângulo de ataque $α$ são dadas por: $$ \psi = U \> r \sin(\theta) $$ Pense nisso: as linhas de fluxo são linhas retas e paralelas que formam um ângulo $α$ com o eixo $x$. Se o fluxo estiver completamente horizontal, $\psi=Uy$. Integrar e você obtém esse $u=U_{∞}$ e $v=0$.

Vamos escrever algum código que preencha os arrays contendo o $u$-velocidade, o $v$-velocidade e a função de fluxo de um escoamento horizontal uniforme ($U_{∞}, α = 0$), em todos os pontos da nossa grade. Observe as úteis funções NumPy ones(), que cria uma nova matriz e a preenche com o valor 1 em todos os lugares e zeros(), o que cria uma matriz preenchida com zeros.

u_inf = 1.0        # velocidade do escoamento no infinito

# calcula o campo de velocidade do escoamento
u_uniforme = u_inf * numpy.ones((N, N), dtype=float)
v_uniforme = numpy.zeros((N, N), dtype=float)

# compute the stream-function
psi_uniforme = u_inf * Y

Observe que calculamos todos os valores de psi_uniformeuma só vez. Não há necessidade desses complicados loops duplos e, quando os arrays são grandes, é muito mais rápido computar! Obrigado Numpy!

Vamos traçar as linhas de fluxo. Usaremos também a função streamplot()também aqui que faz tudo.

# plotando as linhas de fluxo uniforme
largura = 10.0
altura = (y_final - y_inicial) / (x_final - x_inicial) * largura
pyplot.figure(figsize=(largura, altura))
pyplot.xlabel('x', fontsize=16); pyplot.ylabel('y', fontsize=16)
pyplot.xlim(x_inicial, x_final); pyplot.ylim(y_inicial, y_final)
pyplot.streamplot(X, Y, u_uniforme, v_uniforme,
                  density=0.5, linewidth=1, arrowsize=3, arrowstyle='->');

4. Vórtice¶

O set-up é o mesmo de antes: carregamos nossas bibliotecas favoritas e criamos uma grade de pontos para avaliar o campo de velocidade.

Na sequência damos a intensidade e posição do vortex, neste exemplo: $Γ=5$ e de posição no centro do seu domínio.

gamma = 5.0                      # intensidade do vortex
x_vortex, y_vortex = 0.0, 0.0    # posição do vortex

Definimos agoras as duas funções que criamos,

função_campo_velocidade(),
função_corrente_vortex(), e

para calcular a função corrente vortex modificamos apenas a linha de return, e em seguida, usaremos nossa função personalizada para calcular os componentes de velocidade, sendo a mesma função dos itens anteriores.

def funcao_corrente_vortex(intensidade_vortex, x_vortex, y_vortex):
    """
    Calcula a função corrente de uma fonte.
    
    Parameters
    ----------
    intensidade_fonte: float
        Intensidade do fluxo do vortex.
    x_fonte: float
        Coordenada-x da posição do vortex.
    y_fonte: float
        Coordenada-y da posição do vortex.
    
    Returns
    -------
    (intensidade_fonte / (2*numpy.pi)) * theta: função simbólica
        função corrente em termos da matemática simbólica.
    """
    r = sympy.sqrt(x**2 + y**2)
    theta = sympy.atan2(y-y_vortex, x-x_vortex)
    return (- intensidade_vortex / (2*numpy.pi)) * sympy.ln(r)

Agora, chame as funções com a força e posição do vórtice, mais as coordenadas da grade de avaliação, para obter a velocidade e a função de fluxo do vórtice.

intensidade_vortex = -5.0                         # intensidade da fonte
x_vortex, y_vortex = 0.0, 0.0                    # localização da fonte

# computando o campo velocidade na grade de malha 
psi = funcao_corrente_vortex(intensidade_vortex, x_vortex, y_vortex)
u_vortex, v_vortex = campo_velocidade(psi)

Agora podemos visualizar as linhas de fluxo de um vórtice e eles se parecem com círculos concêntricos ao redor do centro do vórtice, como esperado.

# plot as linhas de fluxo da fonte
largura = 10.0
altura = (y_final - y_inicial) / (x_final - x_inicial) * largura
pyplot.figure(figsize=(largura, altura))
pyplot.xlabel('x', fontsize=16); pyplot.ylabel('y', fontsize=16)
pyplot.xlim(x_inicial, x_final); pyplot.ylim(y_inicial, y_final)
pyplot.streamplot(X, Y, u_vortex(X,Y), v_vortex(X,Y),
                  density=2, linewidth=1, arrowsize=2, arrowstyle='->');

Batystuta Rocha

Vamos começar!¶

1. Geração da malha¶

2. Fonte/Sumidouro na origem¶

3. Escoamento uniforme¶

4. Vórtice¶