Equações Diferenciais Parciais 1
Equações Diferenciais Parciais Modelagem de Problemas Físicos O objetivo da modelagem de problemas físicos é encontrar um conjunto de equações matemáticas que descrevam adequadamente um fenômeno físico e possibilitem encontrar uma solução exata ou solução aproximada. A solução exata usualmente é fruto de um método de solução analítica encontrado através de métodos algébricos e diferenciais; enquanto que a solução aproximada é resultante da aproximação da solução analítica empregando-se métodos numéricos, que usualmente baseiam-se em operações aritméticas elementares. Classificação de Problemas Físicos A maioria dos problemas físicos pode ser dividida em sistemas discretos e contínuos. Um sistema discreto consiste de um número finito de elementos interconectados, enquanto que um sistema contínuo envolve um fenômeno que ocorre sobre uma região contínua. Exemplo do primeiro é um conjunto de massas conectadas entre si por um sistema de molas. Um exemplo para o segundo tipo de sistema é a condução de calor numa chapa. Problemas de equilíbrio: são aqueles em que o sistema permanece constante em relação ao tempo. São também conhecidos como problemas em regime permanente ou no estado estacionário. Exemplos são a estática de estruturas, o escoamento compressível em regime permanente e a distribuição de campos eletrostático e magnético estacionários. Geralmente são descritos por um sistema de equações lineares e, em casos particulares, por sistema de equações não lineares ou por equações diferenciais ordinárias com condições de contorno fechado. Problemas de autovalores: são considerados como extensão dos problemas de equilíbrio nos quais valores críticos ou específicos de certos parâmetros adicionais devem ser determinados, além daqueles correspondentes ao estado estacionário. Exemplos desta categoria de problemas físicos são a flambagem e a estabilidade de estruturas, problemas de freqüência natural em sistemas mecânicos e a determinação de ressonância em circuitos elétricos. A solução deste tipo de problema envolve tanto a solução de sistema de equações lineares com matriz de coeficientes singular ou por equações diferenciais com as condições de contorno fechado. Problemas de propagação: incluem os fenômenos transitórios e de regime não-permanente e são aqueles nos quais os estados subseqüentes de um sistema devem ser relacionados com um estado conhecido inicialmente. São exemplos deste tipo de problema a propagação de ondas em meios elásticos contínuos, vibrações auto-excitadas e a condução térmica em regime transiente. A tabela I esquematiza as relações entre os principais tipos de problemas físicos e as correspondentes equações matemáticas que modelam os fenômenos físicos. TABELA I - Relação entre os tipos de problemas e o correspondente conjunto de equações matemáticas que governam o fenômeno.
Classificação do Equações que governam o fenômeno físico
Equações Diferenciais Parciais 2
problema físico Discreto Continuo
Equilíbrio Sistema de equações algébricas simultâneas
Equações diferenciais ordinárias ou parciais com condições de contorno fechadas
Autovalores Sistema de equações algébricas simultâneas ou equações diferenciais ordinárias redutíveis à forma algébrica
Equações diferenciais ordinárias ou parciais com condições de contorno fechadas
Propagação Sistema de equações diferenciais ordinárias simultâneas com condições iniciais conhecidas
Equações diferenciais parciais com condições iniciais prescritas e condições de contorno abertas
Neste capítulo, trataremos dos problemas descritos por equações diferenciais parciais de 2ª ordem lineares, utilizando o método das diferenças finitas para obter a solução numérica. Classificação das Equações Diferenciais Parciais de 2a Ordem As equações diferenciais parciais (EDP) de 2ª ordem podem ser classificadas em três tipos: elípticas, parabólicas e hiperbólicas. A classificação das EDPs objetiva especificar os métodos numéricos mais apropriados em sua resolução. Geralmente, os problemas de fronteira fechada são descritos por equações elípticas, enquanto que os problemas de fronteira aberta são descritos por equações parabólicas ou hiperbólicas. Se a solução de um problema for descrito pela variável u = u(x,y), a EDP que expressa a relação entre u e as variáveis independentes x e y pode ser escrita, genericamente, como:
gfueuduucubua yxyyxyxx =+++++ 2 (1)
na qual a, b, c, d, f e g são constantes ou funções das variáveis independentes x e y. Os coeficientes a, b e c são tais que:
0222 ≠++ cba (2) Se os coeficientes a, b e c forem constantes, pode-se classificar as EDPs lineares de forma análoga às curvas cônicas tridimensionais (Fig. 1) através das seguintes relações:
• EDPs hiperbólicas: b2 – ac > 0, raízes reais e distintas • EDPs parabólicas: b2 – ac = 0, raízes reais e idênticas • EDPs elípticas: b2 – ac < 0, raízes conjugadas complexas
Equações Diferenciais Parciais 3
x
y
z
Hipérbole
Elipse
Parábola
Fig. 1. Curvas cônicas representativas das EDPs lineares. A classificação das EDPs em três grupos representativos tem importância na sua análise teórica, na descrição de métodos numéricos e nas aplicações. A tabela II apresenta os principais tipos de equações diferenciais parciais.
TABELA II - Tipos de equações diferenciais parciais
NOME TIPO EQUAÇÃO
Equação de Laplace elíptica ∇ =+ + =
2 0
0
u
u u uxx yy zz
Equação de Poisson elíptica
cuuu
cu
zzyyxx =++=∇2
Equação de Fourier parabólica
0)(
2
=−++α∂∂=∇α
tzzyyxx uuuu
t
uu
Equação da onda hiperbólica
0)(2
2
222
=−++
∂∂=∇
ttzzyyxx uuuuc
t
uuc
Equação de Lorentz hiperbólica
fuuuuc
t,z,y,xft
u
cu
ttzzyyxx −=−++
−=∂∂−∇
)(
)(1
2
2
2
22
Equações Diferenciais Parciais 4
Método das Diferenças Finitas O Método das Diferenças Finitas (MDF) é um método numérico de cálculo de problemas de valor de fronteira bastante popular por causa da sua simplicidade e facilidade de implementação computacional. Ele baseia-se na aproximação das derivadas de primeira e de segunda ordem da função u(x,y) pelas respectivas equações de diferenças divididas de primeira e de segunda ordem em x e y: Diferença central:
l2
)()(
2
)()(
11
11
11
11
11
11
−+
−+
−+
−+
−+
−+
−=
−−
≈∂∂
−=
−−
≈∂∂
j,ij,i
ii
jiji
j,ij,i
ii
jiji
uu
yy
y,xuy,xu
y
u
h
uu
xx
y,xuy,xu
x
u
, (3)
Diferença de 2ª ordem:
( )
( ) 2
11
211
11
2
2
2
11
211
11
2
2
2)()(2)(
2)()(2)(
l
+−
−+
+−
+−
−+
+−
+−=
−+−
≈∂∂
+−=
−+−
≈∂∂
j,ij,ij,i
jj
jijiji
j,ij,iji
ii
jijiji
uuu
yy
y,xuy,xuy,xu
y
u
h
uu,u
xx
y,xuy,xuy,xu
x
u
(4)
Equação de Laplace em coordenadas retangulares
A equação de Laplace 02 =∇ u é uma típica EDP elíptica que pode ser representada pelo MDF como:
022
00
2
11
2
11
2
2
2
2
=+−
++−
⇒=∂∂+
∂∂
⇒=+
+−+−
l
j,ij,ij,ij,ij,ij,i
yyxx
uuu
h
uuu
y
u
x
uuu
(5)
Se adotarmos o incremento h = ℓ, a equação (5) pode ser re-escrita como:
{ } 041
11112=++−+ +−+− j,ij,ij,ij,ij,i uuuuu
h (6)
Se escrevermos o valor da função u(x,y) no ponto central (xi, yj), vem que:
Equações Diferenciais Parciais 5
41111 +−+− +++
= j,ij,ij,ij,ij,i
uuuuu (7)
A equação (7) mostra que a equação de Laplace calcula o valor da função u(xi,yj) pela média aritmética dos quatro valores vizinhos ao ponto de coordenadas (xi,yj). Observar que cinco pontos são envolvidos no cálculo do valor central ui j na equação (7). O operador laplaciano pode ser representado pictorialmente pelos quatro pontos que circundam o ponto central na forma:
−=∇1
141
11
22
hu ji (8)
A representação pictorial em cinco pontos do operador laplaciano possui erro proporcional à h2. Em adição à fórmula de cinco pontos, pode-se calcular a fórmula do operador laplaciano com nove pontos, com erro proporcional à h6:
−=∇141
4204
141
6
12
2
hu ji (9)
Exemplo: Problema de transferência de calor por condução em regime permanente. É um típico problema bidimensional (2D) que representa a aplicação prática de uma EDP elíptica. O exemplo ilustra como o método das diferenças finitas transforma uma equação diferencial em um sistema de equações algébricas para calcular o potencial (temperatura) em pontos de uma malha quadrada (h = ℓ).
02
2
2
2
=∂∂+
∂∂
y
T
x
T
Condições de contorno: 0)0( =,xT
0)10( =,xT 0)0( =y,T 100)20( =y,T
Equações Diferenciais Parciais 6
x
y
T = 0
T = 0 T = 100
T = 0
10
20
Fig. 2. Placa retangular com temperaturas definidas na fronteira do problema.
[ ] 041
11112=−+++ −+−+ j,ij,ij,ij,ij,i TTTTT
h
−=∇1
141
11
22
hT ji
Se o domínio do problema for discretizado com uma malha de diferenças finitas com h = 5 nas duas direções, obtém-se a geometria mostrada na Fig. 3.
0 Co
0 Co
0 Co0 Co
0 Co0 Co0 Co
100 CoT1 T3T2
Fig. 3. Discretização do domínio do problema de transferência de calor em uma placa retangular (h = 5). O sistema de equações é descrito por:
( ) 040005
1122
=−+++ TT
( ) 04005
12312
=−+++ TTT
Equações Diferenciais Parciais 7
( ) 04010005
1322
=−+++ TT
Na forma matricial, o sistema de equações pode ser descrito como:
−=
⋅
−−
−
100
0
0
410
141
014
3
2
1
T
T
T
A solução deste sistema fornece os valores: T1 = 1,79oC, T2 = 7,14 oC e T3 = 26,79 oC.
0 Co
0 Co
0 Co0 Co
0 Co0 Co0 Co
100 Co1,79oC 7,14oC 26,79oC
Fig. 4. Solução do problema de transferência de calor em uma placa retangular. A solução da equação de Laplace pode ser refinada aumentando o número de linhas e colunas, como mostra a Fig. 5.
T3T1 T2 T4 T7T6T50 Co
0 Co
0 Co
0 Co 0 Co0 Co0 Co0 Co 0 Co0 Co
0 Co 0 Co0 Co0 Co0 Co 0 Co0 Co
100 Co
100 Co
100 Co
T10T8 T9 T14T13T12T11
T17 T21T20T19T18T16T15
Fig. 5. Discretização do domínio do problema de transferência de calor em uma placa retangular. Fazendo-se n = 7 colunas e m = 3 linhas, resulta a matriz A com valores -4 na diagonal principal e 1 nas bandas adjacentes à diagonal. A Fig. 6 mostra a matriz banda tridiagonal A, juntamente com os vetores incógnitas Tj e o vetor contendo as condições de contorno.
Equações Diferenciais Parciais 8
Observar que o padrão da matriz esparsa apresenta a seqüência de elementos [1 -4 1 ] em torno da diagonal e um ou dois valores iguais a 1 fora da banda tridiagonal. A exceção a seqüência são a sétima, oitava, a 14ª e 15ª linhas, por causa do esquema de numeração dos nós (vide Fig. 5). A Fig. 7 mostra o diagrama esquemático da matriz banda A obtido pela função Matlab spy(A) .
−
−
−
=
−−
−−
−−
−−
−−
−−
−−
−−
−−
−−
−
100
0
0
0
0
0
0
100
0
0
0
0
0
0
100
0
0
0
0
0
0
411
1411
1411
1411
1411
1411
141
1411
11411
11411
11411
11411
11411
1141
141
1141
1141
1141
1141
1141
114
21
20
19
18
17
16
15
14
13
12
11
10
9
8
7
6
5
4
3
2
1
T
T
T
T
T
T
T
T
T
T
T
T
T
T
T
T
T
T
T
T
T
Fig. 6. Sistema de equações na forma matricial.
0 2 4 6 8 10 12 14 16 18 20 22
0
2
4
6
8
10
12
14
16
18
20
22
nz = 85
Fig. 7. Representação esquemática da matriz tridiagonal A mostrando os coeficientes não nulos na forma de pontos.
Equações Diferenciais Parciais 9
A tabela III apresenta o perfil 2D de temperatura calculado pelo programa Matlab laplace.m .
TABELA III - Perfil de temperatura na placa retangular 0 0 0 0 0 0 0 0 50.0000 0 0.3530 0.9132 2.0103 4.2957 9.15 32 19.6632 43.2101 100.0000 0 0.4989 1.2894 2.8324 6.0194 12.65 38 26.2894 53.1774 100.0000 0 0.3530 0.9132 2.0103 4.2957 9.15 32 19.6632 43.2101 100.0000 0 0 0 0 0 0 0 0 50.0000
Programa Matlab, para o cálculo pelo MDF do problema de transferência de calor em regime permanente para uma placa retangular. % LAPLACE.M Programa para o calculo do laplaciano sobre uma regiao % retangular com n pontos na direcao x e m pontos na direcao y. % Metodo das diferencas finitas clear; % Numero de linhas e colunas m = input('Digite o numero de linhas '); n = input('Digite o numero de colunas '); % Define a matriz de coeficientes A = zeros(n,m); nl1 = n - 1; np1 = n + 1; msize = n*m; msln = msize - n; % Colocando o valor -4 na diagonal principal for i = 1:msize A(i,i) = -4.0; end % Colocando o valor 1 acima e abaixo da diagonal pr incipal for i = 2:msize A(i-1,i) = 1.0; A(i,i-1) = 1.0; end % Substituindo alguns valores 1 por 0 for i = n:n:msln A(i,i+1) = 0.0; A(i+1,i) = 0.0; end for i = np1:msize A(i,i-n) = 1.0; A(i-n,i) = 1.0; end % Montando a matriz de termos independentes (valore s de fronteira) for i = 1:msize k = mod(i,n); if k == 0 b(i,1) = -100; else b(i,1) = 0;
Equações Diferenciais Parciais 10
end end % Solucao da EDP u = A\b % Ordenamento matricial dos valores de temperatura k = 0; for i = 1:m for j = 1:n k = k + 1; Temp(i,j) = u(k); end end Temp
Condições de contorno As condições de contorno que devem ser especificadas para o cálculo de equações diferenciais podem ser de três tipos: 1) Condição de contorno de Dirichlet, na qual os valores da função potencial u são conhecidos
na fronteira do problema; 2) Condição de contorno de Neumann ou condição de contorno natural, na qual são conhecidos
os valores da derivada primeira de u na fronteira do problema; 3) Condição de contorno mista, na qual é constituído pelos valores da função u e das suas
derivadas. A condição de contorno de Neumann pode ser calculada aplicando-se o esquema de diferença dividida de 1ª ordem regressiva ou progressiva, conforme a fronteira do problema:
ll
j,ij,ij,ij,i
j,ij,ij,ij,i
uuuu
y
u
h
uu
h
uu
x
u
−=
−≈
∂∂
−=
−≈
∂∂
+−
+−
11
11
(10)
Exemplo: Vamos utilizar o problema de transferência de calor em regime permanente na placa retangular do exemplo anterior, considerando as seguintes condições de contorno:
0)0( =,xT ; 0)10( =,xT ; 0)0( =y,T ; 5020
=∂∂
= y,xx
T
Observar que no lado direito da placa, o gradiente sendo positivo indica que o fluxo de calor se dá no sentido de fora para dentro da placa:
hTTh
TT
x
Tj,ij,i
j,ij,i 5050 11 +=⇒=
−≈
∂∂
−− (11)
Equações Diferenciais Parciais 11
ou seja, a temperatura na fronteira é maior do que a temperatura no interior do domínio do problema.
x
y
T = 0
T = 0dT/dx =
50
T = 0
10
20
Fig. 8. Placa retangular com temperaturas definidas na fronteira do problema.
0 Co
0 Co
0 Co0 Co
0 Co0 Co0 Co
T1 T3T2 T4
dT/dx = 50
Fig. 9. Discretização do domínio do problema de transferência de calor em uma placa retangular (h = 5). Considerando a discretização da Fig. 9, o valor da fronteira na face direita é uma incógnita, sendo que o seu valor será calculado pela equação (11). Desta forma, o sistema de equações será descrito pelas equações:
( ) 040005
1122
=−+++ TT
( ) 04005
12312
=−+++ TTT
( ) 04005
13422
=−+++ TTT
A condição de contorno de Neumann, descrita pela equação (11) fornece a relação:
25050 334 +=+= ThTT
Equações Diferenciais Parciais 12
Na forma matricial, o sistema de equações pode ser descrito como:
−=
⋅
−−
−−
250
100
0
0
1100
1410
0141
0014
4
3
2
1
T
T
T
T
A solução deste sistema fornece os valores: T1 = 8,54oC, T2 = 34,15 oC, T3 = 128,05 oC e T4 = 378,05 oC.
0 Co
0 Co
0 Co0 Co
0 Co0 Co0 Co
8,54oC
dT/dx = 50
34,15oC 128,05oC 378,05oC
Fig. 10. Solução do problema de transferência de calor em uma placa retangular com condição de contorno de Neumann. Equação de Poisson Os problemas em regime permanente que envolve a presença de fonte e/ou dreno de potencial, como no caso da transferência de calor com geração de calor interna ao domínio do problema, são descritos pela equação de Poisson:
02 =+∇ cu (12) na qual c é uma constante relacionada com a variação do potencial u. Outro tipo de problema descrito pela equação de Poisson é a torção de barras retangulares, na qual a função u é a função torção e a constante c é o produto do ângulo de torção por unidade de comprimento pelo módulo de rigidez do material. Aplicando-se as equações de diferenças finitas de 2ª ordem (4) em (12) resulta:
0
1
141
11
22 =+
−=+∇ cuh
cu ji (13)
A equação (13) permite obter o sistema de equações lineares de onde se obtém a solução do problema físico descrito pela Equação de Poisson de maneira análoga à obtida no exemplo da equação de Laplace acrescido da constante c que será adicionada aos termos da matriz de termos independentes.
Equações Diferenciais Parciais 13
Métodos numéricos para equações parabólicas A forma geral das EDPs parabólicas 2D pode ser descrita pela equação:
( ) ( ) ucububuauau yxyyxxt 12121 −+++= (14)
na qual a1 > 0, a2 > 0, c1 ≥ 0, b1 e b2 são funções das variáveis independentes x, y e t, sendo x e y variáveis espaciais e t a variável temporal. Na substituição das derivadas de (14) pelo método das diferenças finitas, podem ser aplicados dois esquemas distintos denominados método explícito, no qual a substituição é direta e método implícito de Crank-Nicolson. Método explícito Para ilustrar o procedimento de aplicação das equações de diferenças finitas pelo método implícito será realizado sobre a equação de difusão 1D:
t
u
Dx
u
∂∂=
∂∂ 1
2
2
(15)
na qual D é o coeficiente de difusão. Aplicando a equações de diferenças finitas (3) e (4) na equação (15), vem que:
t
uu
Dh
uuu ki
ki
ki
ki
ki
∆−=+− +
+−1
211 12
(16)
Os índices i na variável u na equação (16) se referem à coordenada xi, enquanto que os expoentes k se referem à variável temporal tk. O intervalo h = ∆x e ∆t é o incremento de tempo. Resolvendo a equação (16) em termos de 1+j
iu :
( ) ki
ki
ki
ki u
h
tDuu
h
tDu ⋅
∆⋅−++∆⋅= +−+
21121 2
1 (17)
Para resolução numérica da equação (17), determina-se o valor futuro 1+j
iu a partir dos
valores conhecidos no tempo tk. Desta forma, a solução de EDPs parabólicas requer o conhecimento não somente das condições de contorno em x mas também dos valores iniciais em t = t0. O tamanho relativo dos passos de integração h e ∆t afeta a equação (17). Se a razão ∆t / h2 for escolhida, tal que:
2
12 =∆⋅ h/tD (18)
a equação (17) será simplificada de tal forma que o último termo desaparecerá, sendo que o cálculo do valor futuro 1+j
iu será determinado pela média dos valores presentes jiu 1− e j
iu 1+ :
Equações Diferenciais Parciais 14
( )ki
ki
ki uuu 11
1
2
1+−
+ += (19)
Exemplo Considere uma placa retangular de aço com 2 cm de espessura e largura e comprimento muito grandes (Fig. 11). Se a temperatura inicial da placa for função da distância da extremidade da placa conforme as equações:
T = 100 x, 10 ≤≤ x T = 100 (2 – x), 21 ≤≤ x
Calcular as temperaturas como função de x e t se ambas as extremidades forem mantidas a T = 0oC.
x
T = 0oC T = 0oC
x = 0 x = 2x (cm)
T (oC)
0 1 2
100
50
25
75
0
2 cm
(a) (b)
Fig. 11. (a) Placa retangular de longo comprimento, com as extremidades mantidas a 0oC. (b) Perfil de temperatura inicial. O problema consiste na obtenção dos valores de temperatura T = T(x, t). Como o comprimento e a largura da placa são grandes o suficiente para considerarmos o problema geométrico 1D, podemos desprezar o fluxo térmico na direção y e z e utilizamos a equação (15) para descrever o problema físico. Tomando-se a equação (17) para resolução da EDP parabólica pelo MDF:
( ) ki
ki
ki
ki T
h
tDTT
h
tDT ⋅
∆⋅−++∆⋅= +−+
21121 2
1
O coeficiente de difusão térmico é expresso por D = k/ρCp., no qual k é a condutividade térmica, ρ a densidade e Cp o calor específico do aço. Assim, sendo os dados físicos para o aço: k = 0,13 cal/s.cm.oC, Cp = 0,11 cal/g.oC e ρ = 7,8 g/cm3, o coeficiente de difusão D = 0,1515 cm2/s. Vamos subdividir a espessura da placa em intervalos h = ∆x = 0,25 cm, de modo a ter oito sub-intervalos. Adotando o critério de convergência (18), temos que o intervalo ∆t será dado por:
Equações Diferenciais Parciais 15
( ) s,,
,
D
ht,h/tD 206250
151502
250
22
1 222 =
⋅==∆=∆⋅
A equação de diferenças finitas (19) será expressa como:
( ) ( ) ( )[ ]kikiki t,xTt,xTt,xT 111 2
1+−+ +=
As condições de contorno do problema são descritas por:
T(0, t) = 0, T(2, t) = 0 As condições iniciais são descritas por:
T(x, 0) = 100 x, 10 ≤≤ x T(x, 0) = 100 (2 – x), 21 ≤≤ x
Os resultados de cálculo pelo MDF estão apresentados na tabela IV. Para comparação, este problema possui solução analítica, cuja expressão é dada por:
( ) ( )( ) ( ) ( ) tn,
n
exn
cosn
t,xT21237380
0
22 2
112
12
1800 +−
∞
=
⋅−+⋅+
= ∑ ππ
(20)
Tabela IV – Resultado do cálculo da EDP parabólica pelo MDF.
t (s) x = 0 x = 0,25 x = 0,50 x = 0,75 x = 1,00 x = 1,25 x = 1,50 x = 1,75 x = 2,00 0 0 25 50 75 100 75 50 25 0
0.20625 0 25 50 75 75 75 50 25 0 0.4125 0 25 50 62.5 75 62.5 50 25 0 0.61875 0 25 43.75 62.5 62.5 62.5 43.75 25 0 0.825 0 21.875 43.75 53.125 62.5 53.125 43.75 21.875 0
1.0313 0 21.875 37.5 53.125 53.125 53.125 37.5 21.875 0 1.2375 0 18.75 37.5 45.313 53.125 45.313 37.5 18.75 0 1.4438 0 18.75 32.031 45.313 45.313 45.313 32.031 18.75 0 1.65 0 16.016 32.031 38.672 45.313 38.672 32.031 16.016 0
1.8563 0 16.016 27.344 38.672 38.672 38.672 27.344 16.016 0 2.0625 0 13.672 27.344 33.008 38.672 33.008 27.344 13.672 0
Observar que os valores de temperaturas são simétricos em relação ao valor médio da espessura da placa (x = 1,0 cm) mostrando que o problema possui simetria em relação ao eixo médio da placa. A figura 12 apresenta a comparação entre as curvas de temperatura em função do tempo para duas posições diferentes dentro da placa calculadas pela solução numérica pelo MDF e pela solução analítica (20). A figura 13 mostra duas curvas da temperatura em função da posição dentro da placa tendo como parâmetro o tempo, calculadas pelo MDF (símbolo) e pela solução analítica (linha contínua). Observar que o perfil de temperatura é simétrico em relação ao centro da placa, conforme mostrado nos valores numéricos da tabela IV.
Equações Diferenciais Parciais 16
Fig. 12. Variação da temperatura em função do tempo para duas posições internas à chapa.
Fig. 13. Variação da temperatura em função da posição para t = 0, 0,825 e 1,65 s, exibindo perfil simétrico em relação ao centro. Parâmetros de cálculo: h = 0,25 cm, ∆t = 0,2065 s e
=∆⋅ 2h/tD 0,5. Comparando os resultados nas figuras 12 e 13, observa-se que o resultado numérico apresenta uma pequena oscilação com relação ao resultado analítico. Em geral, o erro é menor do que 4%. Quando o valor da relação rh/tD =∆⋅ 2 em (18) for diferente de ½, a equação de diferenças finitas (17) é ligeiramente mais complexa, embora do ponto de vista de esforço
Equações Diferenciais Parciais 17
computacional não seja mais lenta do que a relação (19). Se r > ½ mantendo-se constante o valor de h haverá aumento no passo temporal ∆t, com conseqüente redução no número de cálculo. Entretanto, ocorrerá instabilidade numérica como pode ser observado nos gráficos da figura 14.
Fig. 14. Variação da temperatura em função da posição para t = 0, 1,0 e 2,0 s, exibindo instabilidade numérica. Parâmetros de cálculo: h = 0,25 cm, ∆t = 0,25 s e =∆⋅ 2h/tD 0,606. A seguir é apresentado o programa Matlab utilizado no cálculo da distribuição de temperatura em função do tempo e da posição dentro da placa, utilizando o critério de convergência (18) para o cálculo do incremento de tempo ∆t. Programa Matlab % PARABOLICA.M Calculo da equacao de difusao de calor % Placa de espessura 2 cm: geometri a 1D % Condicoes iniciais: T(x,t=0) = 10 0x, x <= 1 % T(x,t=0) = 10 0(2 - x), x > 1 % Condicoes de contorno: T(x=0,t) = 0, T(x=2,t) = 0 clear; % Propriedades fisicas do aco k = 0.13; Cp = 0.11; dens = 7.8; D = k/dens/Cp; % Difusividade termica h = input('Incremento h = '); % Passo de x nx = round(2/h); nx = nx + 1; % Numero de pas sos de x dt = h^2/(2*D); % Calculo do passo de tempo (c riterio de convergência) nt = 11; % Numero de passos de tempo % Condicoes iniciais for i = 1:nx x(i) = (i-1)*h; if x(i) <= 1; T(i,1) = 100*x(i); elseif x(i) <= 2; T(i,1) = 100*(2 - x(i));
Equações Diferenciais Parciais 18
end end % Condicoes de contorno for k = 1:nt T(1,k) = 0; T(nx,k) = 0; end % Calculo da EDP t(1) = 0; for k = 1:nt-1 t(k+1) = t(k) + dt; for i = 2:nx-1 T(i,k+1) = 1/2*(T(i-1,k) + T(i+1,k)); end end % Transposicao da matriz T e do vetor t T = T'; t = t'; % Solucao analitica for k = 1:nt for i = 1:nx Te(i,k) = 0; for n = 0:1000 Te(i,k) = Te(i,k) + ... 1/(pi^2*(2*n+1).^2)*cos((2*n+1)*pi* (x(i)-1)/2)*... exp(-0.3738*(2*n+1).^2*t(k)); end end end % Transposicao da matriz Te Te = 800*Te';
Método implícito de Crank-Nicolson O método implícito de Crank-Nicolson objetiva reduzir o problema de instabilidade numérica observado no método explícito através da substituição do esquema de diferença de primeira ordem progressiva no tempo pelo esquema de diferença de primeira ordem central:
∆−⋅=
+−++− +++
++−+−
t
uu
Dh
uuu
h
uuu ki
ki
ki
ki
ki
ki
ki
ki
1
2
11
111
211 122
2
1 (20)
Diferença central Diferença central Diferença central em tk em tk+1 em tk+1/2
Definindo 2h
tDr
∆⋅= , substituindo e rearranjando na equação (20), vem que:
( ) ( ) k
iki
ki
ki
ki
ki urururururur 11
11
111 2222 +−
++
++− +−+=−++− (21)
Fazendo r = 1, obtém-se a simplificação da equação (21):
ki
ki
ki
ki
ki uuuuu 11
11
111 4 +−
++
++− +=−+− (22)
Equações Diferenciais Parciais 19
A grande vantagem do método implícito de Crank-Nicolson é a sua estabilidade para qualquer valor de r, embora maior precisão no cálculo seja obtida para pequenos valores de r. Deste modo, a equação (21) é a forma usual do método de Crank-Nicolson. Exemplo Retomaremos o exemplo numérico utilizado para exemplificar o método explícito. Aplicando a equação de diferenças finitas pelo método implícito de Crank-Nicolson, com razão r = 1 (equação 22) ao problema de transferência de calor em regime transiente em uma placa de aço com 2,0 cm de espessura, obtém-se:
ki
ki
ki
ki
ki TTTTT 11
11
111 4 +−
++
++− +=−+−
Adotando-se h = 0,25 cm, o valor de ∆t será calculado a partir de =∆⋅= 2h/tDr 1, ou seja,
∆t = s,,/,D/h 4125015150250 22 == Assim, teremos nx = 8 intervalos em x e calcularemos a temperatura em nt = 10 passos de tempo. Para t1 = 0,4125 s, ao final do primeiro passo de tempo, as temperaturas na placa são incógnitas que serão calculadas a partir das temperaturas iniciais em t0 = 0:
08
06
18
17
16
07
05
17
16
15
06
04
16
15
14
05
03
15
14
13
04
02
14
13
12
03
01
13
12
11
02
00
12
11
10
4
4
4
4
4
4
4
TTTTT
TTTTT
TTTTT
TTTTT
TTTTT
TTTTT
TTTTT
+=−+−
+=−+−
+=−+−
+=−+−
+=−+−
+=−+−
+=−+−
Como os valores de temperatura ( )00
0 === t,xTT ik
i (i = 0, 1, ..., 8) no lado direito do
sistema de equações acima são fornecidos na condição inicial e nas condições de contorno, o sistema de equações fornecerá os valores de temperatura no passo de tempo seguinte,
( )4125011 ,t,xTT i
ki === . Uma vez obtido esses valores, no próximo passo de tempo, em
t2 = 0,825 s, eles serão os valores que serão alocados no lado direito do sistema de equações. Esse procedimento será realizado da mesma forma para cada intervalo de tempo. Substituindo os valores de temperatura inicial e os valores de temperatura do contorno do problema, o sistema de equações geral, para um intervalo de tempo tk+1 será dado pelo seguinte sistema de equações:
Equações Diferenciais Parciais 20
004
4
4
4
4
4
040
61
71
6
751
71
61
5
641
61
51
4
531
51
41
3
421
41
31
2
311
31
21
1
21
21
1
+=−+−
+=−+−
+=−+−
+=−+−
+=−+−
+=−+−
+=−+−
++
+++
+++
+++
+++
+++
++
kkk
kkkkk
kkkkk
kkkkk
kkkkk
kkkkk
kkk
TTT
TTTTT
TTTTT
TTTTT
TTTTT
TTTTT
TTT
que pode ser re-escrito na forma matricial:
++++++
+
=
⋅
−−−
−−−−
−−−−
−
+
+
+
+
+
+
+
0
0
4100000
1410000
0141000
0014100
0001410
0000141
0000014
6
75
64
53
42
31
2
17
16
15
14
13
12
11
k
kk
kk
kk
kk
kk
k
k
k
k
k
k
k
k
T
TT
TT
TT
TT
TT
T
T
T
T
T
T
T
T
Para cada intervalo de tempo, a matriz de coeficientes A do sistema acima será a mesma, de modo que um método de solução do sistema de equações lineares apropriado para este tipo de problema é o método LU, no qual a matriz de coeficientes A é decomposta em duas matrizes L e U, de modo que o produto LU = A. Para o sistema de equações A x = b, a substituição A = LU nos leva a:
A x = b → LU x = b Se fizermos U x = y, a solução do sistema será obtida a partir de Ly = b, que por substituição reversa fornecerá os valores de y na equação U x = y. O cálculo do problema térmico está mostrado na figura 15 na qual pode-se observar que não ocorre instabilização numérica.
Fig. 15. Cálculo da variação da temperatura em função da posição e do tempo através do método implícito.
Equações Diferenciais Parciais 21
Programas Matlab A função implicito resolve o problema da EDP parabólica utilizando o método das diferenças finitas implícito e a função trid para o cálculo do sistema de equações tridiagonal (22) pelo método LU. function [u,x,t] = implicito(D,xf,tf,it0,bx0,bxf,M, N) % EXPLICITO.M Solucao da EDP parabolica 1D: u_xx = u_t para 0 <= x <= xf, 0 <= t <= tf % Metodo de diferencas finitas explici to % Condicao inicial: u(x,0) = it0(x) % Condicao de contorno: u(0,t) = bx0(t ), u(xf,t) = bxf(t) % M: no. de subintervalos na direcao x % N: no. de subintervalos de tempo dx = xf/M; x = [0:M]'*dx; dt = tf/N; t = [0:N]*dt; % Condicao inicial for i = 1:M+1 u(i,1) = it0(x(i)); end % Condicao de contorno for k = 1:N+1 u([1 M+1],k) = [bx0(t(k)); bxf(t(k))]; end r = D*dt/dx/dx; r2 = 1+2*r; for i = 1:M-1 A(i,i) = r2; if i > 1, A(i-1,i) = -r; A(i,i-1) = -r; end end for k = 2:N+1 b = [r*u(1,k); zeros(M-3,1); r*u(M+1,k)] + u(2:M ,k-1); u(2:M,k) = trid(A,b); end
function x = trid(A,b) % TRID.M Resolucao de sistema de equacoes tridia gonal % Triangularizacao superior N = size(A,2); row=0; if size(b,1) == 1, b = b(:); row = 1; end for m = 2:N tmp = A(m,m-1)/A(m-1,m-1); A(m,m) = A(m,m) - A(m-1,m)*tmp; A(m,m-1) = 0; b(m,:) = b(m,:) - b(m-1,:)*tmp; end % Substituicao reversa x(N,:) = b(N,:)/A(N,N); for m = N-1: -1: 1 x(m,:) = (b(m,:) - A(m,m+1)*x(m+1))/A(m,m); end if row == 1, x = x'; end
Equações Diferenciais Parciais 22
Métodos numéricos para equações hiperbólicas Uma EDP hiperbólica típica e de grande interesse prático é a equação de onda 2D:
2
2
2
2
2
2
2
1
t
u
y
u
x
u
v ∂
∂=
∂
∂+∂
∂ (23)
na qual v é a velocidade de propagação da onda. Substituindo cada derivada pela respectiva diferença finita central e admitindo-se que h = ∆x = ∆y, vem que:
2
11
2
1111
2 )(
241
t
uuu
h
uuuuu
v
kj,i
kj,i
kj,i
kj,i
kj,i
kj,i
kj,i
kj,i
∆+−
=
++−+⋅
−++−+− (24)
Resolvendo a equação para o tempo tk+1:
( ) kj,i
kj,i
kj,i
kj,i
kj,i
kj,i
kj,i u
hv
tuuuuu
hv
tu
∆−+−+++∆= −+−+−
+22
21
111122
21 )(
42)(
(25)
Nas equações acima, os expoentes k, k-1 e k+1 representam os valores de u nos tempos tk, tk-1 e tk+1, respectivamente. Se escolhermos (∆t)2/v2 h2 = ½, o último termo da equação (25) desaparece, de modo que ela pode ser re-escrita como:
( ) 1111122
21 )( −
+−+−+ −+++∆= k
j,ik
j,ik
j,ik
j,ik
j,ik
j,i uuuuuhv
tu (26)
Referências Bibliográficas W. E. Williams, PARTIAL DIFFERENTIAL EQUATIONS - Oxford Applied Mathematics and Computing Science Series, New York: Oxford University Press, 1980. H. Pina, MÉTODOS NUMÉRICOS, Lisboa: McGraw-Hill, 1995. M. C. Cunha, MÉTODOS NUMÉRICOS, Campinas: Editora da Unicamp, 2000. D. H. Norrie, G. De Vries, THE FINITE ELEMENT METHOD, Fundamentals and Applications, New York: Academic Press, 1973. C. F. Gerald, P. O. Wheatley, APPLIED NUMERICAL ANALYSIS, Reading, MA: Addison-Wesley, 1989. T. J. Akai, APPLIED NUMERICAL METHODS FOR ENGINEERS, New York: John Wiley, 1994.
Recommended
