32.13 Uma barra composta isolada é formada por duas partes encostadas nas extremidades e ambas as partes têm comprimentos iguais. A parte possui condutividade térmica , para , e a parte tem condutividade térmica , para . As equações de condução de calor transientes adimensionais que descrevem a temperatura ao longo do comprimento da barra composta são
onde é a temperatura, é a coordenada axial, é o tempo e . As condições de contorno e iniciais são
Condições de contorno
Condições iniciais
Resolva esse conjunto de equações para a distribuição de temperaturas como função do tempo. Use diferenças finitas com acurácia de segunda ordem para as derivadas com uma formulação de Crank-Nicolson para integrar no tempo. Escreva um programa computacional para a solução e escolha valores de e para uma boa acurácia. Trace a temperatura versus o comprimento para diversos valores de tempo . Gere uma curva separada para cada um dos seguintes valores do parâmetro e 0 .
Passo 1
O enunciado trata de condução de calor transiente unidimensional em uma barra composta de comprimento adimensional unitário dividida em duas metades e com condutividades e . Introduzimos e temos os PDEs
com condições de contorno , , continuidade de fluxo em , e condição inicial . Vamos usar método de diferenças finitas de segunda ordem no espaço e esquema de Crank–Nicolson no tempo para aproximar , escrever um programa em Python e plotar versus para vários e para .
Passo 2
Discretizamos o domínio em pontos uniformes , . Separamos índice tal que . Para cada subdomínio usamos a diferença central de segunda ordem
No trecho esquerdo (i=1,dots,m-1) atribuímos coeficiente 1; no direito (i=m+1,dots,N-1) coeficiente r. Em i=m (interface) impomos continuidade de derivada pela relação
um-um-1De<a x = um+1-umDe<a x ;Longrightarrow; um-1-2um+um+1=0
Passo 3
No esquema de Crank–Nicolson integramos de t^n a t^n+1=t^n+De<a t obtendo o sistema linear
onde L é a matriz tridiagonal resultante de partial^2/partial x^2 com pesos 1 ou r, e b incorpora condições de contorno (vetor com b1=bN+1=De<a t/(2De<a x^2)). Montamos duas matrizes tridiagonais A=I-tfracDe<a t2L e B=I+tfracDe<a t2L.
Passo 4
Implementamos em Python:
import numpy as np
import matplotlib.pyplot as plt
from scipy.linalg import solve_banded
defsolve_bar(r, N=100, dt=1e-4, t_max=0.1):
dx = 1.0 / N
m = int(0.5 / dx)
# Monta L em formato banda para solve_banded: (lower, upper, shapeN)
lower = np.zeros(N - 1); main = np.zeros(N + 1); upper = np.zeros(N - 1)
# pontos internos i=1..N-1for i inrange(1, N):
if i < m:
coeff = 1.0elif i > m:
coeff = r
else: # interface
coeff = 0.0
main[i] = -2 * coeff
if i > 1: lower[i - 2] = coeff
if i < N - 1: upper[i] = coeff
# Bordas
main[0] = main[N] = -1.0# para impor u=1, effectively diag=1 no sistema# Monta A e B
L_bands = np.vstack([np.hstack(([0], lower)), main, np.hstack((upper, [0]))])
A = np.eye(N + 1) - dt / (2 * dx * dx) * L_bands
B = np.eye(N + 1) + dt / (2 * dx * dx) * L_bands
# Inicial e BC
u = np.zeros(N + 1)
u[0] = u[N] = 1.0
t = 0.0
sol = [(t, u.copy())]
# Loop tempowhile t < t_max:
# monta RHS e corrige BC
rhs = B.dot(u)
rhs[0] = 1.0; rhs[N] = 1.0# resolve A u^{n+1} = rhs; aqui podemos usar solve_banded ou solver direto
u = np.linalg.solve(A, rhs)
t += dt
sol.append((t, u.copy()))
return sol, np.linspace(0, 1, N + 1)
# Executa para vários r e plota
rs = [1, 0.1, 0.01, 0.001, 0]
times = [0.001, 0.01, 0.05, 0.1]
for r in rs:
sol, x = solve_bar(r)
plt.figure()
for (t, u) in sol:
ifany(abs(t - ti) < 1e-6for ti in times):
plt.plot(x, u, label=f"t={t:.3f}")
plt.title(f"r={r}")
plt.xlabel("x"); plt.ylabel("u");
plt.legend(); plt.show()
Passo 5
Escolhemos N=100 (De<a x=0.01) e De<a t=10^-4 garantindo estabilidade do esquema implícito (CN é incondicional) e boa resolução temporal. A cada tempo selecionado, extraímos e plottamos u(x,t).
Passo 6
Ao variar r, observa-se que para rto0 a segunda metade praticamente não conduz calor e o perfil de u atinge 1 somente próximo a x=0 e x=1, enquanto para r=1 a barra é homogênea com perfil linear transiente.
Resposta
Segue o programa Python acima que implementa diferenças finitas de segunda ordem e Crank–Nicolson, e gera os gráficos de u(x,t) para r=1,0,1,0,01,0,001,0 nos tempos indicados.
Cálculo Numérico - Steven C. Chapra, Raymond P. Canale - Métodos Numéricos para Engenharia - 5ª Edição - Capítulo 32.Problemas - Ex. 32.13