Aplicação de Métodos Espectrais ao Cálculo Numérico de
Escoamentos Compressíveis
André Rosalen
Instituto Tecnológico de Aeronáutica
Praça Marechal Eduardo Gomes, 50 - Vila das Acácias – São José dos Campos – São Paulo. CEP: 12228-461
Bolsista PIBIC – CNPq
[email protected]
....................................................................................................................................................................................................
O objetivo do presente estudo é entender e aplicar os métodos espectrais como ferramenta na solução, por
cálculo numérico, de equações diferenciais parciais que modelam escoamento compressíveis.
A metodologia empregada foi, em uma primeira aproximação, conhecer o método através de referências
bibliográficas e posterior aplicação numérica em equações diferenciais, para verificação das características que esse
moderno método traz consigo.
Conforme será mostrado adiante, os resultados obtidos mostraram traços da potencialidade do método, não em
sua plenitude), no entanto é possível identificar as características do método.
1. Introdução
A proposta da bolsa de iniciação científica é o desenvolvimento de um algoritmo numérico para o estudo
de escoamento através de bocais transônicos, assim como aerofólios, através dos métodos espectrais. Como o
estudo realizado para tal são através de equações diferenciais parcias resolvidas numericamente, o intuito é
estudar um dos métodos mais modernos que existe atualmente em termos de convergência de resultados e
mostrar suas vantagens quanto a tempo e requisito de capacidade de processamento que são exigidos para se
obter os mesmos resultados (conforme veremos, melhores ainda) para o estudo de uma equação diferencial
qualquer, que venha a modelar um dos objetos descritos acima. Como sabemos, o cálculo numérico para
resolução de equações diferenciais parciais ocorre fundamentalmente devido às aproximações das derivadas
nas equações por discretizações, e nesse aspecto há inúmeras formas de se realizar tal feito, os quais se
aprimoram cada vez mais. Uma grande característica da maioria das formas de discretização das derivadas
existentes para o método espectral é que aquela , para um determinado ponto numa malha discretizada, sua
derivada depende de como a função do problema se comporta em sua vizinhança, enquanto que no método
espectral, para qualquer ponto estudado em questão, sua derivada tomará a informação da malha inteira; aqui
vemos uma diferença muito importante, principalmente na acuracidade que o método traz consigo, pois esse
utiliza mais informações do domínio da função estudada que os demais métodos.
2. Apresentação e desenvolvimento de algumas ferramentas do método espectral
Para mostrarmos o estudo, vejamos primeiramente um caso oscilatório. Para tal, precisaremos dos
conceitos de expansão em série de Fourier, para posteriormente chegarmos ao conceitos de polinômio
interpolador e matriz de derivação de uma função dada, para uma determinada discretização no domínio da
mesma.
Considerando-se uma função escalar de duas variáveis,
determindo instante de tempo, com dominío
u ( x, t ) , x denotando uma distância e t um
x ∈ [ 0, l ] , período l e t ∈ [ 0, ∞ ] , então a função pode ser
reprensentada por uma expansão em série de Fourier, e em sua forma complexa é dada por (suponha que L =
semi período da função):
u ( x, t ) =
∑ a ( t ).e
 nπ x 
i .

 L 
n
,
(1)
n ≤∞
sendo:
 nπ x 
L 
2L
− i .
1
. ∫ u ( x, t ).e 
an ( t ) =
2L 0
Tomando-se
lim ∑ an ( t ).e
dx
 nπ x 
i .

 L 
n →∞
(2)
teremos a representação de nossa função em série de Fourier.
Para chegarmos ao nosso polinômio interpolador da função, primeiramente vamos discretizar o domínio
em
x da função. Seja então:
x j = j.∆x = j.
2π
, j ∈ [ 0,..., N ]
N +1
Com isso, para
an ( t ) (realizando-se a integral pela quadratura trapezoidal), teremos que:
 nπ x j 

L 
− i .
1 N −1
an ( t ) =
.∑ u ( x j , t ) .e 
N .cn j = 0
~
com
(3)
cn = 2 se n =
,
(4)
N
N
, ou cn = 1 se n <
.
2
2
Por essa aproximação, chegamos então ao polinômio interpolador de
u ( x, t ) , apenas no domínio
espacial (onde a função possuí a característica ondulatória):
In (u ) =
∑
an ( t ) .e
 nπ x 
i .

 L 
(5)
N
n≤
2
Escrevendo a equação acima em termos de cossenos e senos, e realizando-se as devidas somas, é
demonstrável que o polinômio interpolador pode ser dado por:
N −1
I n ( u ) = ∑ u ( x j , t ) .g j ( x ) ,
(6)
j =0
onde:
g j ( x) =

( x − x j )  .cot  ( x − x j ) 
1
.s en ( 2 N − 1) .
2N
2 

 2 

(7)
Com esse resultado, podemos agora determinar a derivada com respeito a
x da função (note que a
~
informação do tempo está contida nos coeficientes
an ( t ) ),
para posterior aplicação nas equações
diferenciais parciais.
No entanto, o elemento
g j ( x ) é a aproximação para a derivada de primeira ordem. É possível
demonstrar (a demonstração encontra-se em Jan S. Hesthaven, Sigal Gottlieb, David Gottlieb) que para
derivadas de ordem superior, a aproximação se dá na seguinte forma segundo os elementos xi das
discretizações da malha espacial:

ξ 

sen ( 2 N − 1)  

1 d
2



gi , j ( p ) ( xi ) =
, j≠i
p
2 N dξ 
ξ 

sen  


2
ξ = xi − x j
p
(8)
0, i = j , p = ímpar
gi , j
( p)
( xi )
( −1)
=
p 2
2N
(9)
N
∑k
p
, i = j , p = par
(10)
k =− N
Para o nosso caso, como não há derivada, chegamos ao resultado de
g j ( x) .
Diferenciando-se e aplicando-se a derivada a um determinado ponto da malha (seja
xi , por exemplo),
teremos então que a derivada será dada por:
N −1
∂g j ( x )
∂I n ( u )
= ∑ u ( x j , t ).
∂x x = xi j =0
∂x x = xi
Di , j
( −1)
=
i+ j
2
Di , j = 0, i = j
.cot
(x − x ) , j ≠ i
i
j
2
(11)
(12)
(13)
Como também é usual obtermos derivadas de segunda ordem no modelamento de problemas, segue
∂ 2 gi, j ( x )
abaixo a forma de
Di , j =
( −1)
i + j +1
2
∂x
2
= D 2i , j :


1
.sen 

2  xi − x j
 sen 
2




,i ≠ j



(14)
2N 2 +1
,i = j
Di , j = −
6
(15)
Possuímos então capacidade para aplicarmos a discretização de derivadas de qualquer ordem, tomando-se
informação de todos os pontos na malha do domínio, que no caso, entra no ramo dos métodos espectrais.
Apesar da formulação espectral ser relativamente recente, há alguns métodos para resolução numérica,
como o método de “Galerkin” e o da “Colocação”.
Para se traduzir matematicamente é relativamente simples, o importante é a filosofia associada a cada um
dos métodos para se obter as equações necessárias para a obtenção dos resultados da função nos pontos
desejados da malha.
• Método de Galerkin: Seleciona-se uma base de funções e impõe-se ortogonalidade com o resíduo da
equação diferencial. Aqui, o intuito é minimizar o resíduo, e nesse processo levantam-se um
conjunto de equações para determinar-se os coeficientes lineares que acompanham cada um dos
elementos da base de função, obtendo-se assim uma solução aproximada pelo método espectral;
• Método da Colocação: Aqui impõe-se resíduo nulo em um conjunto de pontos colocados no
domínio. Para simplificações, geralmente utilizam-se os próprios pontos escolhidos como malha da
equação.
2.1. Características de convergência e erros de truncamentos na aplicação do método espectral
Antes de prosseguirmos, é importante mostrar algumas características da convergência espectral , por
exemplo, apresenta comportamento exponencial se a função a ser discretizada na dada variável segue algumas
características.
1 - Seja a identidade de Parseval:
u − P2 N u
Se
2
L2 [ 0,2π ]
= 2π . ∑ uˆn
2
(16)
n >N
u ( x ) ∈ H r p [ 0, 2π ] , segundo a norma-r de Sobolev, então ∃ C tal que:
u − P2 N u
L2 [ 0,2π ]
≤ C.N − q . u ( q )
(17)
L2 [ 0,2π ]
Ainda, se u(x) é analítica, então
u(q)
L2 [ 0,2π ]
q ! ∼ q q .e− q , e assumindo q ∝ N , então segue que:
≤ C.q !. u
L2 [0,2π ]
, e segundo a fórmula de Stirling:
u − P2 N u
( N ) .e
≤ ∼ C. q
L2 [ 0,2π ]
q
−q
.u
L2 [ 0,2π ]
∼ K .e − C . N . u
L2 [ 0,2π ]
(18)
Vemos aqui o fator exponencial indicando a característica da convergência.
2 – Em Jan S. Hesthaven, Sigal Gottlieb, David Gottlieb, é possível demonstrar as seguintes
características para a acuracidade das derivadas expandidas em série de Fourier:
a) Se
u ( x ) ∈ W p r [ 0, 2π ] , então ∃ C tal que:
u − P2 N u
Wq p [ 0,2π ]
≤
C
.u
N r −q
(19)
W p r [0,2π ]
Notemos aqui a convergência de ordem r-q em N para a derivada expandidade em série de Fourier.
b) Se
u ( x ) ∈ C p q [ 0, 2π ] , então ∃ C tal que:
u − P2 N u
L∞
1
≤ C.
N
q− 1
. u(
2
q)
L2 [ 0,2π ]
,q > 1
2
(20)
Novamente, para esse caso vemos que a convergência ocorre em ordem q-1/2 em N.Por último,
verifiquemos o erro de truncamento para o operador diferencial:
3 – Para um operador
L com coeficientes constantes tal que:
d ju
Lu = ∑ a j . j ,
dx
j =1
então para ∀q, r com 0 ≤ q + s ≤ r , ∃ C tal que:
(21)
≤ C .N − ( r − q − s ) . u
(22)
s
Lu − LPN u
W p q [ 0,2π ]
W p r [ 0,2π ]
Novamente, notemos a convergência do operador diferencial, para esse caso, de ordem r-q-s em N.
Apesar das características expostas específicas, a interpretação a se retirar é que as dependências de
convergências estão fortemente ligadas a ordens em N (número de pontos da malha discretizada na variável
escolhida), seja exponencialmente ou por ordem de decaimento. Esse tipo de aspecto é de grande valia,
sobretudo na acuracidade dos cálculos e na capacidade de processamento computacional requerida, pois N se
reduz substancialmente para se obter os resultados desejados segundo um determinado nível de qualidade de
dados, quando comparados a outros métodos utilizados até então (Runge-Kutta, por exemplo).
2.2. Exemplo de aplicação do método espectral, pelo método da Colocação
Como forma de demonstrar a aplicação do conceito, estudou-se problemas de natureza oscilatória. Segue
abaixo um exemplo:
∂u ( x, t )
∂u ( x, t )
=−
, x ∈ [ 0, 2π ] , t ≥ 0
,
∂t
∂x
u ( x, 0 ) = sen ( x )
com fronteiras periódicas.
(23)
Aqui foi dado o formato inicial da onda, pois assim poderemos comparar os resultados numéricos com o
analítico para o problema proposto acima. Sabemos que esse problema se trata de uma onda viajando da
esquerda para direita com função
domínio em
u ( x, t ) = sen ( x − t ) . Para tal, tomou-se uma malha com 16 pontos no
x , discretizou-se a parcela no tempo por Runge-Kutta de segunda ordem e utilizou-se a razão
∆t
= 0, 025465 ( ∆t = 0, 01 s); no domínio em x, utilizou-se a derivada de primeira ordem deduzida
∆x
acima. Em contraste com a solução analítica, foram obtidos os seguintes resultados:
Tabela 4.1. Dados relativos a discretização da malha por uma uniforme, de espaçamento
Número de
pontos na
malha
16
Pontos na
Índice dos
malha
associados no
pontos
observados[i] intervalo 2*pi
0
0,00
1
0,39
2
0,79
3
1,18
4
1,57
5
1,96
6
2,36
7
2,75
8
3,14
9
3,53
10
3,93
11
4,32
12
4,71
13
5,11
14
5,50
15
5,89
∆x =
2π
16 + 1
Tabela 4.2. Dados obtidos pelo método numérico em contraste com a solução analítica, para instante inicial
em t = 0 s
Tempos
tomados
(discretização)
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
Condição de
partida para os
valores em j da
função (tempo t)
0,000000
0,382683
0,707107
0,923880
1,000000
0,923880
0,707107
0,382683
0,000000
-0,382683
-0,707107
-0,923880
-1,000000
-0,923880
-0,707107
-0,382683
Derivada
temporal da
função no
ponto
determinado
Valor da
segundo o
função para o
tempo dado
tempo t + Dt
-2,000000
-0,020000
-1,847759
0,364206
-1,414214
0,692965
-0,765367
0,916226
0,000000
1,000000
0,765367
0,931533
1,414214
0,721249
1,847759
0,401161
2,000000
0,020000
1,847759
-0,364206
1,414214
-0,692965
0,765367
-0,916226
0,000000
-1,000000
-0,765367
-0,931533
-1,414214
-0,721249
-1,847759
-0,401161
Valor da
função exata
para (x,t)
dados
-0,010000
0,373426
0,700000
0,920007
0,999950
0,927660
0,714142
0,391903
0,010000
-0,373426
-0,700000
-0,920007
-0,999950
-0,927660
-0,714142
-0,391903
1,500000
1,000000
0,500000
Solução analítica
0,000000
1 2
3 4 5
6 7 8
9 10 11 12 13 14 15 16
Solução numérica
-0,500000
-1,000000
-1,500000
Figura 4.1. Apresentação gráfica das soluções obtidas para o tempo t = 0 s
Tabela 4.3. Dados obtidos pelo método numérico em contraste com a solução analítica, para instante inicial
em t = 0,01 s
Condição de
partida para os
valores em j da
função (tempo t)
-0,020000
0,364206
0,692965
0,916226
1,000000
0,931533
0,721249
0,401161
0,020000
-0,364206
-0,692965
-0,916226
-1,000000
-0,931533
-0,721249
-0,401161
Tempos
tomados
(discretização)
0,01
0,01
0,01
0,01
0,01
0,01
0,01
0,01
0,01
0,01
0,01
0,01
0,01
0,01
0,01
0,01
Derivada
temporal da
função no
ponto
determinado
Valor da
segundo o
função para o
tempo dado
tempo t + Dt
-2,000000
-0,040000
-1,863067
0,345575
-1,442498
0,678539
-0,802322
0,908203
-0,040000
0,999600
0,728411
0,938817
1,385929
0,735108
1,832452
0,419486
2,000000
0,040000
1,863067
-0,345575
1,442498
-0,678539
0,802322
-0,908203
0,040000
-0,999600
-0,728411
-0,938817
-1,385929
-0,735108
-1,832452
-0,419486
Valor da
função exata
para (x,t)
dados
-0,019999
0,364130
0,692824
0,916042
0,999800
0,931348
0,721107
0,401083
0,019999
-0,364130
-0,692824
-0,916042
-0,999800
-0,931348
-0,721107
-0,401083
1,500000
1,000000
0,500000
Solução analítica
0,000000
1
2
3
4
5
6
7
8
9 10 11 12 13 14 15 16
Solução numérica
-0,500000
-1,000000
-1,500000
Figura 4.2. Apresentação gráfica das soluções obtidas para o tempo t = 0,01 s
Tabela 4.4 - Dados obtidos pelo método numérico em contraste com a solução analítica, para instante inicial
em t = 0,02 s
Condição de partida
para os valores em j
da função (tempo t)
-0,040000
0,345575
0,678539
0,908203
0,999600
0,938817
0,735108
0,419486
0,040000
-0,345575
-0,678539
-0,908203
-0,999600
-0,938817
-0,735108
-0,419486
Tempos
tomados
(discretização)
0,02
0,02
0,02
0,02
0,02
0,02
0,02
0,02
0,02
0,02
0,02
0,02
0,02
0,02
0,02
0,02
Derivada
temporal da
função no ponto
determinado
segundo o
tempo dado
-1,999200
-1,877635
-1,470217
-0,838972
-0,080001
0,691150
1,357079
1,816405
1,999200
1,877635
1,470217
0,838972
0,080001
-0,691150
-1,357079
-1,816405
Valor da
Valor da
função para o
função exata
tempo t + Dt para (x,t) dados
-0,059992
-0,029996
0,326798
0,354799
0,663837
0,685578
0,899813
0,911985
0,998800
0,999550
0,945729
0,934943
0,748679
0,727999
0,437650
0,410224
0,059992
0,029996
-0,326798
-0,354799
-0,663837
-0,685578
-0,899813
-0,911985
-0,998800
-0,999550
-0,945729
-0,934943
-0,748679
-0,727999
-0,437650
-0,410224
1,500000
1,000000
0,500000
Solução analítica
0,000000
1
2
3
4
5
6
7
8
9 10 11 12 13 14 15 16
Solução numérica
-0,500000
-1,000000
-1,500000
Figura 4.3. Apresentação gráfica das soluções obtidas para o tempo t = 0,02 s
Os resultados obtidos acima ainda não possuem a característica precisão que o método espectral propõe,
pois:
1. Há discretização de segunda ordem na parcela temporal, o que diminui a precisão do método;
2.
A razão
∆t
= 0, 025465 , apesar da discretização optada para a parcela temporal, ainda é
∆x
grosseira;
Esses resultados foram apresentados para demonstrar a potencialidade do método espectral. Note que a
convergência ocorre até no máximo na segunda casa decimal, o que ainda é comum encontrar em
procedimentos pouco sofisticados de discretização. Alterando
∆t
para 10-3 (portanto, alterando apenas
∆x
∆t para 0,000393s), vejamos os resultados para o instante em t = 0,001s:
Tabela 4.5 - Dados obtidos pelo método numérico em contraste com a solução analítica, para instante inicial
em t = 0,001 s
Tempos
tomados
(discretização)
0,001
0,001
0,001
0,001
0,001
0,001
0,001
0,001
0,001
0,001
0,001
0,001
0,001
0,001
0,001
0,001
Condição de
partida para os
valores em j da
função (tempo t)
-0,001571
0,381232
0,705996
0,923278
0,999999
0,924480
0,708217
0,384134
0,001571
-0,381232
-0,705996
-0,923278
-0,999999
-0,924480
-0,708217
-0,384134
Derivada
temporal da
função no
ponto
determinado
Valor da
segundo o
função para o
tempo dado
tempo t + Dt
-1,999999
-0,002356
-1,848960
0,380506
-1,416434
0,705439
-0,768269
0,922976
-0,003142
0,999998
0,762464
0,924779
1,411991
0,708772
1,846556
0,384860
1,999999
0,002356
1,848960
-0,380506
1,416434
-0,705439
0,768269
-0,922976
0,003142
-0,999998
-0,762464
-0,924779
-1,411991
-0,708772
-1,846556
-0,384860
Valor da
função exata
para (x,t)
dados
-0,001178
0,381595
0,706273
0,923428
0,999999
0,924330
0,707939
0,383772
0,001178
-0,381595
-0,706273
-0,923428
-0,999999
-0,924330
-0,707939
-0,383772
Aqui já notamos que a convergência ocorre na terceira casa decimal, com boa aproximação na quarta
casa. Ainda, a precisão dos resultados é diminuída devido a parcela temporal ser derivada por um método de
segunda ordem.
Esse exemplo foi para ilustrar que o método espectral, dependendo da forma como é aplicado, pode
possuir precisão acima daquelas obtidas por outros métodos. Por exemplo, se tivéssemos aplicado
discretização de segunda ordem no domínio em x , conseguiríamos no máximo precisão na segunda casa
decimal, com a razão
∆t
utilizada.
∆x
3. Agradecimentos
Os agradecimentos seguem ao meu orientador Marcos Aurélio Ortega, pela disposição e ajuda realizadas
durante os estudos do trabalho, ao ITA, pela oportunidade oferecida de ingressar em um estudo de alta
qualidade em uma área que foi interesse do orientado, e ao CNPq, pela oportunidade de um contrato para o
estudo em questão.
4. Referências
Canuto, C., Hussaini, M. Y., Quarteroni, Z., Zang, T. A., Spectral methods in fluid dynamics, SpringerVerlag, (1988).
Jan S. Hesthaven, Sigal Gottlieb, David Gottlieb, Spectral Methods for Time-Dependent Problems,
9780521792110, Cambridge University Press.
MIT – Opencourse ware: http://ocw.mit.edu/OcwWeb/web/home/home/index.htm.