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)
\tau_{yx}(y)=\mu\frac{dv_x}{dy}.
Usando a solução analítica:
\frac{dv_x}{dy}=\frac{G}{2\mu}\,(H-2y),
logo
\tau_{yx}(y)=\mu\cdot\frac{G}{2\mu}\,(H-2y)=\frac{G}{2}\,(H-2y).
Note que \tau_{yx} é linear em y. Isso também nos dá uma referência para comparar com o resultado numérico.

Passo 4

Esquema numérico para v_x(y) — discretização por diferenças finitas (explicado passo a passo)
  • Dividimos o intervalo [0,H] em N subintervalos igualmente espaçados com passo \Delta y = H/N. Os nós são y_0=0,\;y_1=\Delta y,\;\dots,\;y_{N}=H.
  • Aproximação central para a segunda derivada no nó interno i ( 1\le i\le N-1):
\frac{d^2 v_x}{dy^2}\Big|_{y_i}\approx\frac{v_{i+1}-2v_i+v_{i-1}}{\Delta y^2},
onde v_i é a aproximação de v_x(y_i).
  • Substituindo na equação \dfrac{d^2 v_x}{dy^2} = -\dfrac{G}{\mu} obtemos para cada nó interno:
\frac{v_{i+1}-2v_i+v_{i-1}}{\Delta y^2} = -\frac{G}{\mu}.
  • Rearranjando:
- v_{i-1} + 2 v_i - v_{i+1} = \Delta y^2 \frac{G}{\mu}.
  • 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:
    \left.\frac{dv}{dy}\right|_{y_i} \approx \frac{v_{i+1}-v_{i-1}}{2\Delta y}.
    • 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).
  • Compute o erro relativo ou L2:
\text{L2 error}=\sqrt{\frac{\sum_{i=0}^{N}(v_i-v_{\text{ex}}(y_i))^2}{\sum_{i=0}^{N} v_{\text{ex}}(y_i)^2}}.
  • Para \tau comparar com a fórmula analítica \tau_{\text{ex}}(y)=\dfrac{G}{2}(H-2y).

Passo 7

Resumo do algoritmo que vou implementar em Python (passo a passo para você entender antes do código)
  1. Escolher parâmetros: H,\;N,\;\mu,\;G (ex.: H=1.0, N=100, \mu=1.0, G=1.0).
  2. Montar o grid: \Delta y=H/N, nós y_i.
  3. 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}.
  4. Resolver A\mathbf{v}=\mathbf{b} para \mathbf{v}=[v_1,\dots,v_{N-1}]^T. Colocar v_0=0 e v_N=0.
  5. Calcular as derivadas \dfrac{dv}{dy} com diferenças centrais para nós internos e diferenças unilaterais de segunda ordem nas bordas.
  6. Calcular \tau_i=\mu\,dv/dy.
  7. Fazer os gráficos: v vs y e \tau vs y. Traçar também as soluções analíticas para comparação.
  8. 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, matplotlib

import numpy as np
import matplotlib.pyplot as plt

def solve_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 in range(nint):
        A[i, i] = 2.0
        if i-1 >= 0:
            A[i, i-1] = -1.0
        if 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
    # internos
    for i in range(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

def analytical_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 is None:
        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()

Exercícios relacionados