Com as seguintes condições iniciais e de contorno: , para , Condição Inicial; , para , Condição de Contorno 1 ; , para , Condição de Contorno 2;
Calcular, utilizando o método das diferenças finitas, a propagação transiente de , com a dispersão dia e velocidade dia. O meio possui comprimento de 100 m , devendo ser dividido em, no mínimo, 20 partes iguais e o tempo variando de 0 a 300 dias. Esboçar o gráfico de . Apresentar os algoritmos elaborados, os desenvolvimentos teórico e de cálculo.
Passo 1
Beleza! O problema pede pra resolver a equação de advecção-dispersão em 1D
com
- CI: C(x,0)=0 para x>0
- CC1: C(0,t)=2
- CC2: C(∞,t)=0
onde
- D(x)=1{,}0+e^{-x} (m²/dia)
- U=0{,}5 (m/dia)
- 0≤x≤100 m dividido em ≥20 partes
- 0≤t≤300 dias
Conteúdo: equações em derivadas parciais, método de diferenças finitas, estabilidade de esquema numérico. Vamos precisar:
Discretizar x e t (dx, dt)
Montar esquema de diferenças finitas
Escolher dt segundo condição de estabilidade
Implementar passo a passo
Passo 2
Vamos fazer o esquema explícito com upwind para advecção (porque U>0) e central para difusão:
Discretização
x_i = i·Δx, i=0…N, com Δx = L/N, L=100 m, N=20 ⇒ Δx=5 m
t^n = n·Δt, n=0…Nt, com Δt que satisfaça estabilidade
Podemos escolher Δt=5 dias ⇒ Nt=300/5=60 iterações.
Passo 3
Algoritmo (pseudocódigo tipo Python):
# parâmetros
L=100; N=20; dx=L/N; dt=5; Nt=300/dt
U=0.5
# malha de x
x = [i*dx for i in range(N+1)]
# precomputar D em cada nó
D = [1.0 + exp(-xi) for xi in x]
# inicializar concentração: zeros
C = [[0.0]*(N+1) for n in range(Nt+1)]
# aplicar CC em todo t
for n in range(Nt+1):
C[n][0] = 2.0 # C(0,t)=2
C[n][N] = 0.0 # C(∞,t)=0 ~ C(100,t)=0
# iteração no tempo
for n in range(0, Nt):
for i in range(1, N): # nós internos
adv = U*(C[n][i] - C[n][i-1])/dx
dif = D[i]*(C[n][i+1] - 2*C[n][i] + C[n][i-1])/(dx*dx)
C[n+1][i] = C[n][i] + dt*(dif - adv)
Passo 4
Explicando a teoria:
Usamos forward-difference em (por isso “explícito”)
Upwind para ∂C/∂x porque garante melhor estabilidade quando o escoamento tem direção definida (U>0)
Central para ∂²C/∂x² pois é segunda ordem em espaço
A condição Δt≤Δx²/(2D) evita oscilações numéricas na parte difusiva
As CC (bordas) são impostas a cada iteração para manter C(0)=2 e C(100)=0
Passo 5
Para esboçar o gráfico de :
Depois de rodar o algoritmo, para instantes dias pegue os vetores C[n50][:], C[n100][:]…
Plote cada perfil vs. no mesmo gráfico
Você verá que, conforme cresce, a concentração se propaga para a direita e tende a um estado quase estacionário
Resposta
Segue o esquema final:
Defina Δx=5 m, Δt=5 dias, Nt=60 passos de tempo.
Use Ci^n+1 =Ci^n +Decom , .
Implemente o loop em tempo e espaço como no pseudocódigo.
Para o gráfico, plote perfis de em instantes representativos.