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 [0,1/2] e [1/2,1] com condutividades ka e kb. Introduzimos r=ka/kb e temos os PDEs
partial^2 upartial x^2=partial upartial t, 0<= x<=12,
r,partial^2 upartial x^2=partial upartial t, 12<= x<=1,
com condições de contorno u(0,t)=1, u(1,t)=1, continuidade de fluxo em x=1/2, e condição inicial u(x,0)=0. Vamos usar método de diferenças finitas de segunda ordem no espaço e esquema de Crank–Nicolson no tempo para aproximar u(x,t), escrever um programa em Python e plotar u versus x para vários t e para r=1,0.1,0.01,0.001,0.

Passo 2

Discretizamos o domínio em N+1 pontos uniformes xi=iDe<a x, De<a x=1/N. Separamos índice m tal que xm=1/2. Para cada subdomínio usamos a diferença central de segunda ordem
partial^2 upartial x^2|iapproxui-1-2ui+ui+1De<a x^2
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
l(I-tfracDe<a t2Lr)u^n+1 = l(I+tfracDe<a t2Lr)u^n + b,
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

def solve_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-1
    for i in range(1, N):
        if i < m:
            coeff = 1.0
        elif 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 tempo
    while 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:
        if any(abs(t - ti) < 1e-6 for 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.

Exercícios de Livros Relacionados