Escreva a implementação em Python para encontrar o perfil de velocidade, também desenhando v_{x} vs y ( 2 pontos).
Desenvolva usando um esquema de solução numérica para encontrar o pertil de cisalhamento ( 2 pontos).
Escreva a implementação em Python para encontrar o pertil de cisalhamento, também desenhando \tau vs y (2 pontos).
Passo 1
O enunciado pede três coisas (vou listar bem certinho para não nos perder):
Escrever uma implementação em Python para encontrar o perfil de velocidade e também desenhar v_x versus y.
Desenvolver (ou justificar) um esquema numérico para encontrar o perfil de cisalhamento (isto é, a distribuição de \tau(y)).
Escrever a implementação em Python para encontrar o perfil de cisalhamento e desenhar \tau versus y.
Informações implícitas que vou usar sem te perguntar (porque você não forneceu um problema físico específico): vou considerar o caso clássico de escoamento unidimensional, estacionário e plenamente desenvolvido entre duas placas paralelas (plano de Poiseuille), ou seja, velocidade apenas na direção x que varia só com y. Isso é suficiente para demonstrar o esquema numérico e os gráficos pedidos e é exatamente o tipo de questão que costuma cair em prova. Vou usar condições de contorno do tipo Dirichlet (placas estacionárias: v_x(0)=0 e v_x(H)=0). Vou apresentar também a solução analítica para que você possa comparar a solução numérica.
Passo 2
Equação de partida (versão 1D de Navier–Stokes para escoamento totalmente desenvolvido e estacionário):
0 = -\frac{dp}{dx} + \mu \frac{d^2 v_x}{dy^2},
onde:
\frac{dp}{dx} é o gradiente de pressão (assumo constante ao longo de x);
\mu é a viscosidade dinâmica (constante);
v_x(y) é a componente de velocidade na direção x que depende apenas de y.
Rearranjando:
\mu \frac{d^2 v_x}{dy^2} = \frac{dp}{dx}.
Convém definir G = -\dfrac{dp}{dx}. Com isso a equação fica
\mu \frac{d^2 v_x}{dy^2} = -G,
ou
\frac{d^2 v_x}{dy^2} = -\frac{G}{\mu}.
As condições de contorno que usarei são
v_x(0)=0,\qquad v_x(H)=0,
com H a separação entre as placas.
Passo 3
Solução analítica (útil para checar numericamente). Integrando duas vezes:
\frac{d^2 v_x}{dy^2} = -\frac{G}{\mu}
integra uma vez:
\frac{d v_x}{dy} = -\frac{G}{\mu} y + C_1,
integra de novo:
v_x(y) = -\frac{G}{2\mu} y^2 + C_1 y + C_2.
Aplicando v_x(0)=0 e v_x(H)=0 obtém-se C_2=0 e C_1=\frac{G}{2\mu}H. Portanto a forma usual:
v_x(y)=\frac{G}{2\mu}\,y(H-y).
Observação de sinal: se \dfrac{dp}{dx}<0 então G>0 e v_x será positiva no meio do canal, como esperado.
A tensão de cisalhamento (definição clássica, para fluido newtoniano)
Para todos os i=1,\dots,N-1 teremos um sistema tridiagonal linear em v_1,\dots,v_{N-1}. Os valores v_0 e v_N são conhecidos (condições de contorno iguais a zero).
Em notação matricial: A\mathbf{v} = \mathbf{b}, onde A tem diagonais [ -1,\;2,\;-1 ] (multiplicadas conforme), e \mathbf{b} tem todos os termos iguais a \Delta y^2 \dfrac{G}{\mu}.
Observação: o sistema é simples e simétrico; pode-se resolver com um solver direto (por exemplo, usando NumPy) ou com o algoritmo de Thomas (algoritmo linear para tridiagonais). Para provas, basta escrever a montagem e usar uma resolução direta — o Thomas é um extra eficiente, mas não obrigatório aqui.
Passo 5
Cálculo numérico do cisalhamento \tau(y)
Após obter \{v_i\}, a derivada dv/dy pode ser obtida por diferenças finitas:
Para nós internos (1\le i \le N-1), uso a diferença central:
Para as fronteiras uso diferenças unilaterais (primeiro e último nó):
\left.\frac{dv}{dy}\right|_{y_0} \approx \frac{-3v_0+4v_1-v_2}{2\Delta y},\qquad \left.\frac{dv}{dy}\right|_{y_N} \approx \frac{3v_N-4v_{N-1}+v_{N-2}}{2\Delta y},que são fórmulas de alta precisão (diferença progressiva/degressiva de segunda ordem).
Depois calcula-se \tau_{yx,i}=\mu\left.\dfrac{dv}{dy}\right|_{y_i}.
Observação: também poderia obter \tau(y) integrando a equação de equilíbrio do momento escrita como
\frac{d\tau}{dy} = -\frac{dp}{dx} = G,
o que dá \tau(y)=\tau(0)+G y. Com \tau known at boundary we can set \tau consistent with BCs; porém o método por derivada de v é o pedido no enunciado (esquema numérico), e é o que vou implementar.
Passo 6
Checagem e controle de erro (boa prática e útil na prova)
Compare a solução numérica v_i com a solução analítica v_{\text{ex}}(y_i)=\dfrac{G}{2\mu}y_i(H-y_i).
Montar a matriz tridiagonal A (tamanho (N-1)\times(N-1)) com -1 na sub- e super-diagonal e 2 na diagonal principal; montar o vetor \mathbf{b} com entradas \Delta y^2\frac{G}{\mu}.
Resolver A\mathbf{v}=\mathbf{b} para \mathbf{v}=[v_1,\dots,v_{N-1}]^T. Colocar v_0=0 e v_N=0.
Calcular as derivadas \dfrac{dv}{dy} com diferenças centrais para nós internos e diferenças unilaterais de segunda ordem nas bordas.
Calcular \tau_i=\mu\,dv/dy.
Fazer os gráficos: v vs y e \tau vs y. Traçar também as soluções analíticas para comparação.
Calcular e imprimir o erro L2.
Agora que a teoria, esquema e passos estão claros, no próximo bloco coloco o código Python que implementa tudo isto e gera os gráficos pedidos.
Resposta
# Código Python: resolve o perfil de velocidade v_x(y) por diferenças finitas# e calcula o perfil de cisalhamento tau(y). Plota v_x vs y e tau vs y.# Dependências: numpy, matplotlibimport numpy as np
import matplotlib.pyplot as plt
defsolve_poiseuille(H=1.0, N=100, mu=1.0, G=1.0):
"""
Resolve mu * d2v/dy2 = -G em 0<=y<=H com v(0)=v(H)=0 por diferenças finitas.
Retorna y, v (com N+1 pontos), dvdy, tau.
Parâmetros:
H: distância entre placas
N: número de subintervalos (nós = N+1)
mu: viscosidade dinâmica
G: -dp/dx (use G>0 para dp/dx<0)
"""# malha
dy = H / N
y = np.linspace(0.0, H, N+1)
# montar sistema para v_1 ... v_{N-1}# A é (N-1)x(N-1) tridiagonal com [-1, 2, -1]
nint = N - 1
A = np.zeros((nint, nint))
b = np.ones(nint) * (dy**2) * (G / mu) # RHS = dy^2 * G/mu (veja derivação)for i inrange(nint):
A[i, i] = 2.0if i-1 >= 0:
A[i, i-1] = -1.0if i+1 < nint:
A[i, i+1] = -1.0# resolver sistema linear
v_interior = np.linalg.solve(A, b)
# construir vetor completo com condições de contorno v0 = vN = 0
v = np.zeros(N+1)
v[1:N] = v_interior # v[0] e v[N] permanecem 0# calcular dv/dy: usar diferença central interna e unilaterais de segunda ordem nas bordas
dvdy = np.zeros_like(v)
# fronteira esquerda (y0)if N >= 2:
dvdy[0] = (-3*v[0] + 4*v[1] - v[2]) / (2*dy)
else: # caso degenerado N=1
dvdy[0] = (v[1] - v[0]) / dy
# internosfor i inrange(1, N):
dvdy[i] = (v[i+1] - v[i-1]) / (2*dy)
# fronteira direita (yN)if N >= 2:
dvdy[N] = (3*v[N] - 4*v[N-1] + v[N-2]) / (2*dy)
else:
dvdy[N] = (v[N] - v[N-1]) / dy
# tau = mu * dv/dy
tau = mu * dvdy
return y, v, dvdy, tau
defanalytical_solution(H=1.0, mu=1.0, G=1.0, y=None):
"""
Solução analítica v(y) = (G/(2 mu)) y (H - y)
tau(y) = mu dv/dy = (G/2) (H - 2y)
"""if y isNone:
raise ValueError("Forneça o vetor y para avaliar")
v_exact = (G / (2.0 * mu)) * y * (H - y)
tau_exact = (G / 2.0) * (H - 2.0 * y)
return v_exact, tau_exact
if __name__ == "__main__":
# parâmetros (você pode mudar para testar)
H = 1.0# distância entre placas
N = 200# subdivisões da malha (quanto maior, mais preciso)
mu = 1.0# viscosidade
G = 1.0# G = -dp/dx (colocar positivo para dp/dx negativo)# resolver numericamente
y, v_num, dvdy_num, tau_num = solve_poiseuille(H=H, N=N, mu=mu, G=G)
# solução analítica para comparação
v_ex, tau_ex = analytical_solution(H=H, mu=mu, G=G, y=y)
# calcular erro relativo L2 para v e tau
l2_v = np.sqrt(np.sum((v_num - v_ex)**2) / np.sum(v_ex**2))
l2_tau = np.sqrt(np.sum((tau_num - tau_ex)**2) / np.sum(tau_ex**2))
print(f"L2 relative error in v: {l2_v:.6e}")
print(f"L2 relative error in tau: {l2_tau:.6e}")
# Plot v_x vs y
plt.figure(figsize=(6,4))
plt.plot(v_num, y, 'b-', label='v (numérico)')
plt.plot(v_ex, y, 'k--', label='v (analítico)')
plt.xlabel(r'\(v_x\)')
plt.ylabel(r'\(y\)')
plt.title(r'Perfil de velocidade \(v_x(y)\)')
plt.legend()
plt.grid(True)
plt.tight_layout()
plt.gca().invert_yaxis() # opção estética: mostrar y=0 no topo (opcional)
plt.show()
# Plot tau vs y
plt.figure(figsize=(6,4))
plt.plot(tau_num, y, 'r-', label=r'\(\tau_{yx}\) (numérico)')
plt.plot(tau_ex, y, 'k--', label=r'\(\tau_{yx}\) (analítico)')
plt.xlabel(r'\(\tau_{yx}\)')
plt.ylabel(r'\(y\)')
plt.title(r'Perfil de cisalhamento \(\tau_{yx}(y)\)')
plt.legend()
plt.grid(True)
plt.tight_layout()
plt.gca().invert_yaxis()
plt.show()