Dissertacao Gabriel

7/21/2019 Dissertacao Gabriel

1/87

MINISTRIO DA DEFESAEXRCITO BRASILEIRO

DEPARTAMENTO DE CINCIA E TECNOLOGIAINSTITUTO MILITAR DE ENGENHARIA

(Real Academia de Artilharia, Fortificao e Desenho, 1792)

PROGRAMA DE PS-GRADUAO EM ENGENHARIA DE DEFESA

DISSERTAO DE MESTRADO

ESTUDO COMPARATIVO ENTRE MTODOS RUNGE-KUTTA EXPLCITO, IMPLCITOE IMEX COM PRESERVAO DE ESTABILIDADE NUMRICA NO-LINEAR

Aluno: GABRIEL DE MORAESNmero de matrcula: ED 11306Orientador: Leonardo Santos de Brito Alves Ph. D.rea de concentrao: Engenharia de DefesaLinha de pesquisa: Modelagem e Simulao em Sistemas de Defesa

Data da aprovao da dissertao da Seo de Ensino:Rio de Janeiro - RJ, 03/02/2014.


2/87

INSTITUTO MILITAR DE ENGENHARIA

GABRIEL DE MORAES

ESTUDO COMPARATIVO ENTRE MTODOS RUNGE-KUTTA

EXPLCITO, IMPLCITO E IMEX COM PRESERVAO DEESTABILIDADE NUMRICA NO-LINEAR

Dissertao de mestrado apresentada ao curso de Ps-graduao em Engenharia de Defesa (PGED) do InstitutoMilitar de Engenharia, como requisito parcial para a obten-o do ttulo de Mestre em Engenharia de Defesa.

Orientador: Leonardo S. de B. Alves, Ph.D.

Rio de Janeiro2014


3/87

c2014

INSTITUTO MILITAR DE ENGENHARIAPraa General Tibrcio, 80-Praia Vermelha

Rio de Janeiro-RJ CEP 22290-270

Este exemplar de propriedade do Instituto militar de Engenharia, que poder inclu-lo em basede dados, armazenar em computador, microfilmar ou adotar qualquer forma de arquivamento.

permitida a meno, reproduo parcial ou integral e a transmisso entre bibliotecas destetrabalho, sem modificao de seu texto, em qualquer meio que esteja ou venha a ser fixado,

para peqsquisa acadmica, comentrios e citaes, desde que sem finalidade comercial e queseja feita a referncia bibliogrfica completa.

Os conceitos expressos neste trabalho so de responsabilidade do autor e do orientador.

620.1 Moraes, Gabriel de5235a Estudo Comparativo entre Mtodos Runge-Kutta

Explcito,Implcito e IMEX com preservao de estabilidade

no linear / Gabriel de Moraes; orientado porLeonardo Santos de Brito Alves. Rio de Janeiro:Instituto Militar de Engenharia, 2014.

xxp.: il.

Dissertao (mestrado) Instituto Militar de Engenharia.Rio de Janeiro, 2014.

1. Engenharia de Defesa. 2. Estabilidade Numrica3. Strong Stability Preserving. 4. Runge Kutta.I. Alves, Leonardo Santos de Brito. II. Instituto Militarde Engenharia.

CDD 620.1

2


4/87

INSTITUTO MILITAR DE ENGENHARIA

Programa de Ps-graduao em Engenharia de Defesa

Dissertao de mestrado apresentada ao curso de Ps-graduao em Engenharia de Defesa(PGED) do Instituto Militar de Engenharia, como requisito parcial para a obteno do ttulo deMestre em Engenharia de Defesa.

Aprovada em 03 de Fevereiro de 2014 pela seguinte Banca Examinadora:

Leonardo S. de B. Alves, Ph.D. - UFF

Itamar Borges Jr., D.Sc. - IME

Paulo Fernando Ferreira Rosa, D.Sc. - IME

Luiz Eduardo Bittencourt Sampaio , D.Sc. - UFF

Rio de Janeiro2014


5/87

AGRADECIMENTOS

Ao orientador Leonardo Alves, pelos ensinamentos que renderam um grande crescimento

pessoal e profissional.

A todos do Laboratrio de Motores e Propulso (LMP), do Instituto Militar de Engenharia,em especial aos amigos Francisco Eduardo, Kledson Flvio, Ricardo Dias e Renan Teixeira sem

o qual esse trabalho no seria possvel.

A todos do Laboratrio de Mecnica Terica e Aplicada (LMTA), da Universidade Federal

Fluminense. Ao amigo Ricardo Figueiredo, pelos conselhos, apoio e confiana. Ao Instituto

Militar de Engenharia e aos professores Itamar Borges e Paulo Rosa pelo amplo suporte e

exemplar dedicao ao programa de Ps-Graduao em Engenharia de Defesa do IME.

Juliana Pepeu, pelo amor incondicional e por ser uma grande companheira compreen-

siva e presente em todos os momentos. minha famlia, pela incessante motivao, apoio e

confiana.

Gabriel de Moraes

4


6/87

SUMRIO

LISTA DE ILUSTRAES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7

LISTA DE TABELAS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10

LISTA DE SMBOLOS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11

1 INTRODUO . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14

1.1 Motivao Tecnolgica . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15

1.2 Motivao Cientfica . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18

1.3 Objetivo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21

2 FORMULAO NUMRICA . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22

2.1 Mtodos Strong Stability Preserving . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22

2.1.1 Formulao Shu-Osher dos mtodos de Runge-Kutta . . . . . . . . . . . . . . . . . . . . . . . . . 23

2.1.2 Formulao Shu-Osher modificada. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24

2.2 Discretizao Espacial: Limitadores de Fluxo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25

2.3 Verificao e validao do cdigo computacional . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26

2.4 Metodologia para avaliao de eficincia computacional . . . . . . . . . . . . . . . . . . . . . . 28

3 RESULTADOS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29

3.1 Caso teste 01: Equao de Buckley-Leverett. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31

3.2 Caso teste 02: Equao de Burgers no viscosa . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35

3.2.1 Verificao Computacional . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36

3.2.2 Estudo de preservao de estabilidade no Linear-Barreira SSP . . . . . . . . . . . . . . . . 38

3.2.3 Estudo comparativo de eficincia computacional . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45

3.3 Caso teste 03: Equao de Burgers viscoso . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49

3.3.1 Estudo comparativo de eficincia computacional . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50

3.4 Caso teste 04: Equao de Buckley-Leverett . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54

3.5 Caso teste 05:Equao de Buckley-Leverett viscosa . . . . . . . . . . . . . . . . . . . . . . . . . . 59

3.6 Caso teste 06:Tubo de Choque . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63

3.7 Caso teste 07:Propagao de ondas termoacsticas em uma cavidade . . . . . . . . . . . . 65

4 CONCLUSES E TRABALHOS FUTUROS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67

5 REFERNCIAS BIBLIOGRFICAS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69

5


7/87

6 APNDICES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72

6.1 Apndice 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72

6.2 Apndice 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77

6.3 Apndice 3 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83

6.4 Apndice 4 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84


8/87

LISTA DE ILUSTRAES

FIG.1.1 Bateria Antiarea composta por msseis Patriot FONTE:

http://goo.gl/w9MlB . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16

FIG.1.2 Foguete Ariane 5FONTE: http://goo.gl/Asd1L . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17

FIG.2.1 Curvas tpicas na comparao entre eficincia computacional entre m-

todos mtodos. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28

FIG.3.1 IC_01 : = 0, = 0Esquerda:Figura 02 JMP 2012, Direita: Soluo

com formulao SSP e semi-analtica MoC. . . . . . . . . . . . . . . . . . . . . . . . . . . . 32

FIG.3.2 IC_02 := 0, = 0Esquerda: Figura 06 JMP 2012, Direita: Soluo

com formulao SSP e semi-analtica MoC. . . . . . . . . . . . . . . . . . . . . . . . . . . . 33FIG.3.3 IC_01 : = 0, = 105 Esquerda: Figura 05 JMP 2012, Direita:

Soluo com formulao SSP. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33

FIG.3.4 IC_02 : = 0, = 105 Esquerda:Figura 09 JMP, Direita: Soluo

com formulao SSP . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34

FIG.3.5 Perfis de velocidade : burges no viscoso. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36

FIG.3.6 Convergncia de malha espacial, burges no viscoso. . . . . . . . . . . . . . . . . . . . . . . 37

FIG.3.7 Verificao de cdigo SSP(2,2) com Superbee. . . . . . . . . . . . . . . . . . . . . . . . . . . . 37

FIG.3.8 Verificao de cdigo SDIRK(1,2) com Superbee. . . . . . . . . . . . . . . . . . . . . . . . . 38

FIG.3.9 Eficincia computacional mtodos explicito x implcito 1,2s burgers no

viscoso. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45


viscoso. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46


viscoso. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47

FIG.3.12 Perfis de velocidade : burges viscoso. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49FIG.3.13 Comparativo: Eficincia de CPU Burgers Viscoso 0,1s. . . . . . . . . . . . . . . . . . . . . 50

FIG.3.14 Comparativo: Eficincia de CPU Burgers Viscoso 0,5s. . . . . . . . . . . . . . . . . . . . . 52

FIG.3.15 Comparativo: Eficincia de CPU Burgers Viscoso 1,0s. . . . . . . . . . . . . . . . . . . . . 53

FIG.3.16 Perfis de velocidade : Buckley-Leverett no viscoso. . . . . . . . . . . . . . . . . . . . . . . 54

FIG.3.17 Comparativo: Eficincia de CPU BL 1,0s. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55



FIG.3.20 Soluo Buckley-Leverett viscoso . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59

7


9/87

FIG.3.21 Comparativo: Eficincia de CPU . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60



FIG.3.24 Soluo para a massa especfica emt = 0.01sand C F L= 0.1atravs

da soluo analtica e numrica. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64

FIG.3.25 Soluo numrica para a massa especfica emt = 0.01sand CF L =

0.1atravs dos mtodos SSP (2,2) e RK2 com discretizao Upwind . . . . . . 64

FIG.3.26 Soluo numrica para a massa especfica em t =

[0.25c, 1.0c, 1.5c, 2.0c] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66

FIG.3.27 Soluo numrica para a velocidade t =[0.25c, 1.0c, 1.5c, 2.0c] . . . . . . . . . . 66

FIG.6.1 Reproduo JMP: IC01 := 0, = 105 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72








FIG.6.9 Reproduo JMP: IC02 := 0, = 108 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76FIG.6.10 Reproduo JMP: IC02 := 0, = 100 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76

FIG.6.11 Erro obtido no refinamento de malhas espaciais resolvido atravs do

mtodo explcito. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77

FIG.6.12 Ordem espacial mtodo explcito. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77

FIG.6.13 Erro obtido no refinamento de malhas temporais resolvido atravs do

mtodo explcito. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78

FIG.6.14 Ordem espacial mtodo explcito. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78


mtodo implcito. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79

FIG.6.16 Ordem espacial mtodo implcito. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79


mtodo implcito. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80

FIG.6.18 Ordem espacial mtodo implcito. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80


mtodo IMEX. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81


10/87

FIG.6.20 Ordem espacial mtodo IMEX. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81


mtodo IMEX. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82

FIG.6.22 Ordem espacial mtodo IMEX. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82


mtodo explcito. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83


mtodo explcito. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83

FIG.6.25 Soluo para a velocidade emt = 0.01sandCF L= 0.1. . . . . . . . . . . . . . . . . . . 84

FIG.6.26 Soluo para a presso emt= 0.01sandC F L= 0.1. . . . . . . . . . . . . . . . . . . . . 84

FIG.6.27 Soluo da presso para os tempos t =[0.25c, 1.0c, 1.5c, 2.0c] . . . . . . . . . . . . 85

FIG.6.28 Evoluo da presso na parede da direita da cavidade ao longo de 30s . . . . . . 85FIG.6.29 Soluo numrica para a temperatura para os tempos caractersticos

t=[0.25,1,1.5,2] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 86


11/87

LISTA DE TABELAS

TAB.3.1 Eficincia computacional no Burgers Viscoso para 1.2 s . . . . . . . . . . . . . . . . . . 46

TAB.3.2 Eficincia computacional Burgers no Viscoso para 2.4 s . . . . . . . . . . . . . . . . . . 46

TAB.3.3 Eficincia computacional Burgers no Viscoso para 4.8 s . . . . . . . . . . . . . . . . . . 47

TAB.3.4 Eficincia computacional em 0.1 s . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51



TAB.3.7 Eficincia computacional Buckley-Leverett no Viscoso para 1.0 s . . . . . . . . . . 56



TAB.3.10 Eficincia computacional em 1.0 s: Buckley-Leverett viscoso . . . . . . . . . . . . . . 60TAB.3.11 Eficincia computacional em 3.0 s: Buckley-Leverett viscoso . . . . . . . . . . . . . . 61

TAB.3.12 Eficincia computacional em 5.0 s: Buckley-Leverett viscoso . . . . . . . . . . . . . . 62

10


12/87

LISTA DE SMBOLOS

E Energia interna total por unidade de massa

h Coeficiente de transferncia de calor por conveco

L Comprimento do meio

P Presso

q Termo fonte

T Temperatura

TL Temperatura do fluido externo

k Condutividade trmica

K Nmero de onda modificado

ui Componentes do vetor velocidade Viscosidade dinmica

Parmetro de controle dos gradientes, e nmero de pontos na malha espacial

Massa especfica

Cp Capacidade trmica a presso constante

C Coeficiente no linear SSP

11


13/87

RESUMO

Diversos problemas complexos em mecnica dos fluidos requerem solues com altas ordens depreciso. Na tentativa de se aumentar a preciso dos esquemas numricos empregados acaba-se, geralmente, causando perda de estabilidade tanto linear quando no linear. Nesse contexto,os mtodos Strong Stability Preserving (SSP) so mtodos capazes de alcanar altas ordens depreciso temporal enquanto mantm a forte estabilidade do mtodo avanado de Euler que dorigem ao esquema numrico. Essa metodologia desenvolvida basicamente reescrevendo-seos mtodos multiestgio como uma combinao de passos do mtodo de Euler e limitando-se otamanho do passo no tempo empregado. Atualmente, estuda-se o desenvolvimento de mtodosSSP com maior eficincia, i.e. com maiores coeficientes. Apesar da grande popularidade queessa classe de mtodos tem conquistado, algumas questes permanecem sem resposta. A forteestabilidade linear dos mtodos implcitos possibilita o emprego de grandes passos no tempo, oque representa uma vantagem significante em comparao aos mtodos explcitos. Entretanto,a limitao no linear do passo no tempo imposta pela metodologia SSP pode tornar os mtodosimplcitos comparativamente ineficientes. Portanto, o presente estudo tem por objetivo princi-

pal avaliar comparativamente a eficincia computacional dos mtodos SSP implcitos, explci-tos e hbridos implcito-explcito na tentativa de descobrir qual metodologia apresenta maioreficincia numrica frente a limitao no linear de passo no tempo. Todos os trs esquemasimplementados so de segunda ordem e os termos espaciais foram derivados atravs da tcnicalimitadora de fluxo no linear de segunda ordem e foram testados em sete casos clssicos damecnica dos fluidos.

12


14/87

ABSTRACT

Several unsteady problems in transport phenomena require highly accurate solutions. Attemptsto increase accuracy order of most numerical schemes, however, often weakens their linearand/or nonlinear stability. On the other hand, Strong Stability Preserving (SSP) methods areable to increase their accuracy order in time while still maintaining the overall stability prop-erties of the original forward Euler method from which they were generated. This is achievedby, among other things, restricting the maximum allowed time step of these schemes. Hence,most recent studies have focused on developing optimal SSP schemes, i.e., minimally restric-tive. Under this context, a significant variety of explicit and implicit formulations for multi-stepand multi-stage marching schemes of different accuracy orders has been created. Despite theirpopularity, some open issues still remain. Linear stability of implicit schemes usually allowsvery large time steps, making them cost effective for many applications when compared to theirexplicit counterparts. Hence, the SSP time step restriction might render these schemes compar-atively inefficient. The present study evaluates different explicit, implicit and implicit explicitSSP time integration schemes for a series of test cases in an attempt to distinguish the most

efficient scheme according to an error / computer time analysis. All three schemes employed asecond order TVD flux limiter discretization for the spatial derivatives and were tested in sevenclassic case studies.

13


15/87

1 INTRODUO

O desenvolvimento industrial e tecnolgico de um pas comumente associado com a exis-

tncia de uma forte indstria de defesa (E.C., 2006). Em mbito nacional, iniciativas como a

publicao da Poltica de Defesa Nacional em 2005 (BRASIL, 2005) e a elaborao do Plano

Estratgico de Defesa Nacional em 2008 (BRASIL, 2008) traaram diretrizes para o desen-

volvimento desse setor. A partir de tais resolues, foram implementadas aes estratgicas

com o intuito de reestruturar o setor de defesa e acelerar o crescimento das reas de Cincia,

Tecnologia e Inovao.

Atuando de forma sinrgica com essa nova poltica de Estado, o Exrcito brasileiro (EB),atravs do Instituto Militar de Engenharia (IME), instituiu o programa de Ps-Graduao em

Engenharia de Defesa (PGED). De forma geral, pode-se definir a Engenharia de Defesa como a

rea da engenharia que trata de todos os ramos relacionados indstria de defesa e aos sistemas

de defesa. um empreendimento multi e interdisciplinar que se desenvolve em um ambiente

transdisciplinar, integra conhecimentos originrios de engenharias, fsica, qumica, biologia e

cincia dos materiais e se configura como uma rea complexa que engloba aspectos de anlise

e sntese relativos ao desenvolvimento, projeto, otimizao, integrao, certificao, avaliao,

operao e logstica de sistemas aplicados defesa (ADES et al., 2010).

Dessa forma, a Engenharia de Defesa abarca conhecimentos de diferentes segmentos da

Engenharia e das Cincias, focando na pesquisa pura e aplicada visando o desenvolvimento de

sistemas de defesa. Os sistemas de defesa, por sua vez, compreendem todos os aparatos tecno-

lgicos capazes de auxiliar em situaes crticas de conflito ou ameaa soberania nacional.

Fica claro, portanto, que o PGED assume, por princpio bsico, a pesquisa e o desenvol-

vimento de novas tecnologias, para que essas possam auxiliar tanto as Foras Armadas como

tambm a sociedade. O PGED concentra suas pesquisas e desenvolvimentos na rea de concen-trao (AC) de engenharia de defesa, sendo esta composta por trs linhas de pesquisa (LP):

LP1 - Comunicaes e Inteligncia em Sistemas de Defesa

LP2 - Mecatrnica e Sistemas de Armas

LP3 - Modelagem e Simulaes em Sistemas de Defesa

A linha de pesquisa trs (LP3) Modelagem e Simulao em Sistemas de Defesa, na qual

o presente trabalho se enquadra, lida com o emprego de mtodos numricos, modelos matem-

14


16/87

ticos, algoritmos computacionais avanados e computao de alto desempenho para investigar

aspectos dos fenmenos fsicos e qumicos relacionados aos sistemas de defesa (ADES et al.,

2010).

Em consonncia com um dos pr-requisitos bsicos do programa, a multidisciplinaridade,

importante ressaltar que o presente trabalho ser desenvolvido baseando-se em conceitos oriun-

dos das reas de engenharia mecnica atravs do estudo aprofundado de fenmentos relacio-

nados com a mecnica dos fluidos, matemtica aplicada, atravs da utilizao de modelos ma-

temticos e desenvolvimento de mtodos numricos, alm de interesses relacionados defesa

nacional, atravs da aplicao dos mtodos propostos em problemas tpicos aos encontrados

nos sistemas de defesa.

1.1 MOTIVAO TECNOLGICA

Dcadas de pesquisas e investimentos na rea de fluidodinmica computacional (CFD) re-

sultaram no desenvolvimento de inmeras metodologias numricas capazes de solucionar os

mais diversos problemas da engenharia moderna. Atualmente, os pacotes comerciais de CFD

empregados na indstria utilizam, em geral, mtodos de baixa ordem de preciso limitados,

muitas das vezes, a esquemas de segunda ordem (PREMASUTHAN, 2010). Dentre os prin-

cipais fatores que influenciaram na popularidade dos mtodos de baixa ordem em aplicaes

industriais aponta-se, em especial, a baixa preciso dos resultados requeridos em aplicaes deengenharia, ou ainda, a necessidade de lidar com problemas complexos, onde a estabilidade nu-

mrica passa a ser uma fonte de preocupao. Entretanto, em projetos de engenharia avanada

ou, simplesmente, quando um valor com alta preciso necessrio, os esquemas limitados

segunda ordem de acurcia so inadequados. Esse o caso de estudos relacionados s reas

de receptividade, aeroacstica e turbulncia (ALVES, 2010). Esquemas de baixa ordem no

so capazes de capturar fenmenos com geometrias complexas com preciso e, na tentativa

de utilizar malhas extremamente refinadas para avaliar tais fenmenos, o custo computacional

torna-se proibitivo. Torna-se evidente, portanto, a imperativa necessidade do emprego de m-

todos de alta ordem de preciso que sejam capazes de simular escoamentos complexos com

tempos computacionais ao menos equivalentes aos mtodos atuais de baixa ordem.

Nesse contexto, a pesquisa cientfica desenvolvida pelos integrantes da LP3 tem por obje-

tivo criar e aperfeioar mtodos numricos capazes de simular problemas fsicos relacionados

aos sistemas de defesa. Os fenmenos fluido-dinmicos presentes nos sistemas de defesa, apre-

sentam caractersticas fsicas complexas, tais como: a presena de fortes descontinuidades nas

propriedades, pontos de estagnao e elevados gradientes. O estudo de tais fenmenos est dire-

15


17/87

tamente relacionado aos projetos de sistemas de propulso e aerodinmica de msseis, foguetes

e aeronaves (MEDEIROS, 2012).

Em aplicaes industriais, a robustez de um cdigo computacional uma necessidade b-

sica. Os curtos prazos em que os projetos devem ser desenvolvidos e a limitao oramentria

exigem eficincia das simulaes realizadas, objetivando, com isso, a reduo do retrabalho ou

da necessidade de executar um grande nmero de prottipos para testes.

Outro fator crtico no emprego de cdigos computacionais para simulao de fenmenos

fsicos na indstria a incapacidade tcnica do usurio final de efetuar ajustes nos cdigos

durante eventuais mudanas nas condies fsicas do problema a ser resolvido. A seguir, so

apresentados dois casos de falhas em cdigos computacionais, no robustos, em aplicaes de

defesa, que causaram mortes e enormes prejuzos ao longo das ltimas dcadas.

Em 25 de fevereiro de 1991, durante a guerra do Golfo, um mssil Scud atravessou a bateriaanti-area composta por msseis Patriot em Dharan, Arbia Saudita, provocando a morte de 28

soldados americanos.Aps a investigao do fato, foi concludo que o mal funcionamento do

sistema antimsseis ocorreu devido perda de preciso dos clculos efetuados (SKEEL, 1992).

Tal sistema foi desenvolvido para um cenrio de guerra europeu, concebido para funcionar du-

rante, no mximo, algumas horas antes de ser transportado para uma nova localizao. No

Golfo, esse sistema estava em funcionamento contnuo h mais de 100 horas, provocando perda

de preciso dos clculos efetuados. Tal sistema realizava clculos para a localizao de msseis

com base nos valores de velocidade e de tempo, que eram considerados como nmeros intei-

ros. A preciso dos clculos estava limitada pelas converses desses inteiros em nmeros reais

computados com 24 bits. A propagao sistemtica de erros durante a operao do sistema

culminou com um atraso de trezentos e quarenta centsimos de segundo, tempo suficiente para

tornar o sistema ineficiente.

FIG. 1.1: Bateria Antiarea composta por msseis PatriotFONTE: http://goo.gl/w9MlB

16


18/87

No evento em questo, foi possvel observar que a modificao das condies de operao

para as quais o sistema foi projetado, colocou em evidncia um modo de avaria crtico que,

anteriormente, no havia sido detectado. Um dos pilares bsicos do desenvolvimento de um

cdigo robusto o controle do acmulo e da propagao de erros numricos.

Outro exemplo da importncia da robustez numrica em sistemas de defesa ocorreu durante

o projeto dos foguetes Ariane 5. Aps uma dcada de pesquisas, desenvolvimentos e investi-

mentos em um projeto de sete bilhes de dlares, a agncia espacial europeia produziu o foguete

Ariane 5, individualmente avaliado em 500 milhes de dlares. Em 4 de junho de 1996, apro-

ximadamente quarenta segundos aps seu lanamento e, acerca de 3700 metros de altitude, o

foguete Ariane 5 explodiu. Investigaes posteriores apontaram que o acidente ocorreu devido

FIG. 1.2: Foguete Ariane 5FONTE: http://goo.gl/Asd1L

a uma falha na converso de dados da velocidade da plataforma inercial de 64/16 bits causandoum erro de 36,7s no clculo da trajetria do foguete (ROBINSON, 1996).

Ambos os casos apresentados mostram a vital importncia da robustez em um sistema de

defesa. Portanto, em um contexto tecnolgico, um sistema robusto deve ser aquele capaz de

cumprir aos objetivos esperados mesmo que operando em condies adversas.

17


19/87

1.2 MOTIVAO CIENTFICA

O mais geral conjunto de leis de governo que descreve o movimento macroscpico de flui-

dos conhecido por equaes de Navier-Stokes. Tais equaes formam um sistema denomi-

nado misto hiperblico/parablico, quando os efeitos de compressibilidade so considerados.Nesse contexto, fica evidente a grande importncia de ser capaz de resolver com eficincia um

conjunto de equaes diferenciais parciais (EDPs) hiperblicas.

Os sistemas de EDPs hiperblicas so empregados rotineiramente para modelar fenmenos

fsicos nos quais o transporte advectivo de substncias importante (RHODES, 2001). Uma

das principais dificuldades observadas durante a soluo numrica de tais equaes a obteno

de perfis descontnuos. Tal comportamento tende a gerar oscilaes esprias e instabilidade

numrica (LEVEQUE, 2005).

A soluo exata da equao hiperblica que modela a adveco de uma substncia se pro-

pagando com velocidade mdia sem alterao em seu perfil apresenta a propriedade de variao

total (TV) constante no tempo. A soluo numrica dessa mesma equao, geralmente, no

apresentar variao total constante e, se o mtodo em questo introduz oscilaes soluo,

ento, a variao total dever aumentar com o tempo. Para evitar que tais oscilaes gerem

solues no fsicas, ou ainda, causem instabilidade numrica, mtodos com preservao de es-

tabilidade no-linear que apresentem a propriedade de variao total diminuda (TVD) devem

ser empregados.Atualmente, vrias abordagens tm sido empregadas no desenvolvimento de tais mtodos,

no-lineares, de ordem superior (PREMASUTHAN, 2010). A maioria dessas abordagens faz

uso de uma tcnica conhecida na literatura por mtodo das linhas que consiste em discretizar,

inicialmente, o sistema de PDEs no espao, para se obter um sistema de ODEs dependentes

somente do tempo. O sistema semi-discreto resultante, por sua vez, temporalmente integrado

utilizando, em geral, o mtodo de Runge-Kutta com forte estabilidade, capaz de preservar pro-

priedades no lineares. Na literatura, essa classe de mtodos recebe o nome de "Strong Stability

Preserving"(SSP).

Strong Stability Preserving methods so mtodos de discretizao temporal capazes de pre-

servar propriedades de estabilidade numrica no-linear como a monotonicidade e variao

total. Essa metodologia tem sido amplamente empregada para resolver escoamentos com-

pressveis[122TK], incompressveis[92TK] e bifsicos. Ao longo dos ltimos anos, essa tc-

nica foi combinada com uma srie de mtodos de discretizao espacial, como Galerkin des-

contnuo [23], ENO[13,28,1], WENO[4,19,115,30,77,2,125,91] alm de mtodos espectrais

[114,22,122,123].

18


20/87

A primeira metodologia desenvolvida para discretizao temporal com propriedades SSP

foi apresentada em (SHU, 1988) e (SHU e Osher, 1988b), baseados nos conceitos de variao

total diminuda. Os trabalhos em questo foram desenvolvidos atravs do emprego do mtodo

das linhas para tratamento dos termos espaciais e discretizao temporal atravs do mtodo de

Euler. Em (SHU, 1988), alm de um mtodo de Runge-Kutta de primeira ordem foram tambm

analisadas metodologias multipassos. Ainda em (SHU e Osher, 1988b), mtodos Runge-Kutta

SSP de segunda a quinta ordem foram implementados e testes realizados para comprovar a

manuteno de estabilidade numrica no linear.

Atualmente, diversos pesquisadores trabalham no intuito de desenvolver metodologias SSP

mais eficientes, isso tem sido feito, basicamente, atravs da busca por parmetros timos (GOT-

TLIEB et al., 2009a). Com essa nova classe de mtodos emergindo, um questionando cls-

sico de estabilidade numrica linear sobre a eficincia computacional entre os mtodos implci-tos/explcitos voltou a ser pertinente, no entanto, com um vis no-linear.

Em (GOTTLIEB e Shu, 1998), Shu e Gotilieb realizaram um estudo sistemtico de mto-

dos Runge-Kutta, onde apresentaram parmetros timos para os mtodos de segunda e terceira

ordem de dois e trs estgios respectivamente e, ainda, provaram que no possvel se obter

um mtodo de quatro estgios e quarta ordem com coeficientes no negativos.

O termo "Strong Stability Preserving"(SSP) foi utilizado pela primeira vez em (SHU et al.,

2001). No trabalho em questo, Gotilieb, Shu e Tadmor apresentaram uma reviso de mtodos

multipassos e Runge-Kutta com preservao de estabilidade no linear existentes na literatura.

Todos os artigos at ento, utilizavam, indiscriminadamente, o termo TVD para tcnicas de

integrao espacial e temporal.

Spiteri e Ruuth, em (SPITERI e Ruuth, 2002), desenvolveram uma nova classe de mtodos

Runge-Kutta SSP que empregam um nmero maior de estgios do que a respectiva ordem de

preciso, com isso, conseguiram utilizar maiores coeficientes e, ainda garantiar a propriedade

SSP e, Gotilieb, em (GOTTLIEB e Gottlieb, 2003), obteve parmetros timos para essa nova

metodologia.Em (SPITERI e Ruuth, 2002), novamente Spiteri e Ruuth provaram que no possvel a

obteno de um mtodo Runge Kutta SSP de quinta ordem com o emprego de coeficientes no

negativos. Nos trabalhos (SPITERI e Ruuth, 2002), (RUUTH. e Spiteri, 2004) e (RUUTH,

2006), os mesmos autores utilizaram a teoria e a otimizao numrica de Shu-Osher para de-

senvolver metodologias downwind com propriedades SSP.

Kraaijevanger, em (KRAAIJEVANGER, 1991), realizou importantes anlises sobre a li-

mitao de ordem nos mtodos Runge-Kutta SSP otimizados. Em (SPIJKER, 2007), Spijker

apresentou um estudo detalhado sobre a teoria SSP na aplicao de mtodos multiestgios para

19


21/87

soluo de equaes no lineares.

A equivalncia entre a teoria de Shu-Osher e da monotonicidade absoluta foram desenvolvi-

das independentemente por Ferracina e Spijker (FERRACINA e Spijker, 2004), (FERRACINA

e Spijker, 2005) e Higueras (HIGUERAS, 2005). A unificao dessas teorias promoveu a con-

solidao de um embasamento terico robusto para os estudos sequentes.

Esquemas essencialmente no-oscilatrios (ENO) e essencialmente no-oscilatrio ponde-

rado (WENO) utilizando diferenas finitas e volumes finitos apresentados em (SHU e Osher

(1988a), SHU e Osher (1989) e BALSARA e Shu. (2000)) so exemplos de esquemas espaciais

capazes de preservar propriedades no lineares. Tais metodologias so essenciais para a imple-

mentao de um esquema SSP. Em linhas gerais, essa metodologia realiza interpolaes para

obter um estncil numrico capaz de capturar, eficientemente, gradientes severos e descontinui-

dades. No entanto, o elevado nmero de avaliaes e ponderaes faz com que esses mtodosrequeiram grande demanda computacional, tornando-os pouco eficientes em problemas com

solues suaves.

Da teoria de anlise numrica linear consenso que os mtodos implcitos possuam uma

zona de estabilidade maior que os mtodos explcitos. Tal estabilidade numrica superior pos-

sibilita a realizao de simulaes com o emprego de um incremento (passo no tempo) maior,

essa prtica muito vantajosa para a obteno da solues de regime permanente ou quando

se deseja simular problemas transientes rgidos onde o passo no tempo demasiadamente limi-

tado para mtodos condicionalmente estveis. No entanto, o nmero de clculos por iterao

requeridos para resolver a matriz implcita superior ao requerido pelos mtodos explcitos.

Portanto, quando se deseja simular problemas em regime permanente onde o emprego do passo

no tempo pode ser o maior possvel os mtodos implcitos so os mais indicados. Por outro

lado, quando se deseja simular problemas em regime transiente onde a magnitude do erro tem-

poral influenciar na soluo desejada o passo mximo no tempo deve ser limitado e, portanto,

os mtodos explcitos so os mais indicados.

Recentemente, os mtodos hbridos Implcitos-explcitos (IMEX) tm ganhado grande des-taque por representar um meio termo entre os mtodos. Atravs do emprego dessa metodologia

possvel resolver parte de uma equao de forma explcita e parte implcita. Isso representa

uma grande vantagem quando deseja-se resolver equaes rgidas contendo termos no lineares.

No entanto, para a simulao de problemas no lineares com a presena de choques e/ou

descontinuidades na soluo, onde o emprego da metodologia SSP necessria, o mesmo enten-

dimento no pode ser replicado. Para manter a estabilidade numrica no linear, uma barreira

SSP deve ser respeitada, com isso, em ( GOTTLIEB et al. (2009b)) questionamentos foram

realizados sobre qual seria a metodologia tima para simular problemas no lineares com o

20


22/87

emprego da metodologia SSP. At o presente momento tal questionamento permanece sem res-

posta.

1.3 OBJETIVO

O presente trabalho tem por objetivo geral realizar um estudo numrico comparativo vi-

sando responder ao questionamento formulado em (GOTTLIEB et al. (2009b)) sobre a eficin-

cia computacional dos mtodos Strong Stability Preserving implcitos, explcitos e implcitos-

explcitos (IMEX) sob um contexto no-linear. Para alcanar este objetivo geral, os seguintes

objetivos especficos foram identificados:

Implementao / Avaliao de esquemas SSP timos para a discretizao temporal com

forte preservao de estabilidade numrica.

Implementao / Avaliao de esquemas de discretizao espacial no oscilatrio reque-

ridos para o xito do esquema SSP implementado.

Realizao de casos teste com fortes gradientes e variaes abruptas de propriedades para

avaliao das propriedades SSP.

Anlise comparativa de custo computacional das formulaes explicitas, implcitas e h-

bridas implcitas-explcitas (IMEX) em cada caso teste estabelecido.

21


23/87

2 FORMULAO NUMRICA

2.1 MTODOS STRONG STABILITY PRESERVING

Strong Stability Preserving (SSP) so mtodos de discretizao temporal de alta ordem ca-

pazes de preservar propriedades de estabilidade numrica no-linear. Essa metodologia am-

plamente empregada para resolver equaes diferenciais parciais hiperblicas, uma vez que os

perfis gerados pela soluo dessas equaes tendem a apresentar variaes abruptas de pro-

priedades e descontinuidades (SHU et al., 2011). Esse conceito pode ser melhor entendido

utilizando o mtodo das linhas, partindo de uma aproximao semi-discreta de uma lei de con-

versao hiperblica unidimensional

ut+f(u)x = 0 (2.1)

ondeu funo de xe t, e a derivadaf(u)x, deve ser discretizada atravs de um esquema es-

pacial que ser aqui representada porF(u). Dessa forma, obtm-se uma equao semidiscreta

ut = F(u) (2.2)

importante perceber que a escolha da tcnica de discretizao espacial aplicada na equa-

o 2.2 determina o tipo de estabilidade que o sistema semidiscreto apresentar. Em seguida,

espera-se discretizar temporalmente esse sistema de forma que as propriedades de estabilidade

no lineares existentes na verso semidiscreta seja mantida. Dessa forma, empregando o m-

todo de Euler avanado para discretizao da derivada temporal, obtm-se

un+1 =un + tF(un) (2.3)

considerando uma metodologia de discretizao espacial capaz de manter alguma propriedade

de estabilidade no linear possvel afirmar que o esquema totalmente discreto obtido em 2.3

tambm apresentar as mesmas propriedades de estabilidade no linear (SHU, 1988). Esse

esquema, no entanto, apresenta somente primeira ordem de preciso temporal. Nesse contexto,

os mtodos SSP so formulaes de altas ordens de preciso capazes de manter a estabilidade

numrica do esquema 2.3. Ou seja, assumindo que um sistema na forma 2.3 discretizado pelo

mtodo de Euler avanado estvel sob alguma norma, semi-norma ou funo convexa

||un+1|| ||un|| (2.4)

22


24/87

ento, um mtodo SSP capaz de manter a mesma estabilidade limitado a uma restrio adicio-

nal do passo no tempo alm daquela imposta pela condio de Courant-Friedrichs-Levy (CFL)

t CtFE. (2.5)

ondeC obtido atravs da diviso entre constantes do mtodo empregado.

2.1.1 FORMULAO SHU-OSHER DOS MTODOS DE RUNGE-KUTTA

O primeiro mtodo de Runge-Kutta explicito SSP foi apresentado em (SHU e Osher, 1988a).

No trabalho em questo, Shu e Osher empregaram uma discretizao espacial com variao

total diminuida (TVD) quando empregada em conjunto com o mtodo de Euler avanado. Por

isso, essa metodologia foi inicialmente conhecida como TVD. Entretanto, a grande contribuio

desse trabalho foi a descrio matemtica de como os mtodos de Runge-Kutta podiam ser

reescritos em combinaes de passos do mtodo de Euler. Cada estgio de um mtodo explcito

de Runge-Kutta escrito na forma de Butcher

u(i) =un + tmj=1

aijF(u(j)) (1 i m) (2.6)

u(n+1) =un + tm

j=1bjF(u

(j)) (2.7)

ondem o nmero de estgios. Essa formulao pode ser reescrita na forma

u(0) =un, (2.8)

u(i) =i1j=0

(i,ju(j) + ti,jF(u

(j))) (1 i m) (2.9)

un+1 =u(m) (2.10)

para mtodos Runge-Kutta, a consistncia exige quei1k=0i,k = 1. Dessa forma, se todos os

coeficientes forem no negativos, cada estgio do Runge-Kutta pode ser reescrito como com-

binaes de um mtodo de Euler avanado, no entanto, com um passo no tempo modificado,

onde,t substitudo por i,ki,k

t.

u(i) = i1j=0

i,ju

(j) + ti,jF(u(j)) i1j=0

i,j

u(j) + t i,ji,j

F(u(j)) (2.11)

A partir dessa observao, foi possvel elaborar o seguinte teorema (SHU, 1988):

23


25/87

Teorema 1: Se um esquema SSP for obtido ao se aplicar o mtodo de Euler avanado

eq.(2.1) sob uma restriot tFEe ainda sei,j, i,j 0, ento, a soluo obtida por um

mtodo de Runge-Kutta tambm ser SSP sob a restrio:

t C(, )tFE onde C(, ) =minijij

ij (2.12)

Essa abordagem pode ser generalizada para todos os mtodos Runge-Kutta. Portanto, exis-

tem condies suficientes para desenvolver uma metodologia RKSSP seja ela explcita, impl-

cita ou IMEX desde que a restrio do passo no tempo no seja violada (SHU et al., 2011).

2.1.2 FORMULAO SHU-OSHER MODIFICADA

At o presente momento foram consideradas somente abordagens explcitas, uma vez que

na equao 2.9, cada estgiou(i) funo deu(n) e de estgios anterioresu(1), u(2),...,u(i1).

Entretanto, pode-se reescrever essa formulao de forma mais geral e, com isso, estender as

futuras anlises tambm para formulaes implcitas. Essa representao foi apresentada na

literatura como formulao de Shu-Osher modificada (FERRACINA e Spijker (2005);HI-

GUERAS (2005)).

u(i) =iun +m

j=0(i,ju

(j) + ti,jF(u(j))) (1 i m+ 1) (2.13)

un+1 =u(m+1) (2.14)

Reagrupando os termos da expresso anterior na forma

u(i) =iun +mj=0

i,j(u(j) + t

i,ji,j

F(u(j))) (1 i m+ 1) (2.15)

e considerando que para que o mtodo seja consistente deve respeitar

iun +

mj=0

i,j = 1 (1 i m+ 1) (2.16)

obtm-se uma combinao convexa de passos de Euler avanado sempre que ij,ijeiforem

positivos. Com isso, possvel generalizar o teorema 1 tambm para esquemas implcitos.

Para melhor ilustrar a formulao de Shu-Osher modificada, reproduzido abaixo um m-

todo de Runge-Kutta trapezoidal de dois estgios:

u(1) =un (2.17)

u(2) =u(1) + tF(u(1)) (2.18)

24


26/87

un+1 =3

4u(1) +

1

4tF(u(1)) +

1

4u(2) +

1

2tF(u(2)) (2.19)

Que pode ser representado na forma

=

0 0

1 034

14

= 0 0

1 014

12

= 1

00

2.2 DISCRETIZAO ESPACIAL: LIMITADORES DE FLUXO

Conforme mencionado anteriormente, para que um esquema seja SSP, alm de tcnicas

especficas de discretizao temporal, uma tcnica no oscilatria deve ser empregada para a

discretizao dos termos espaciais. No presente trabalho foram utilizadas as tcnicas limita-

doras de fluxo com filtro no linear Superbee (TVD) de segunda ordem para a discretizaoespacial. As tcnicas limitadoras baseiam-se na decomposio dos fluxos de altas ordens na

soma de um fluxo de baixa ordem e um termo de ordem superior. De forma geral, possvel

arbitrar uma funo limitadora de fluxo tal que

hni+1/2 = hnLi+1/2

+ni[hHi+1/2 hLi+1/2] (2.20)

ondehLe hHso funes de fluxo numrico de um esquema conservativo de baixa ordem e de

um esquema de alta ordem respectivamente e a funoni limitante de fluxo. O objetivo dessa

metodologia desenvolver um esquema onde a funo de fluxo numrico tenha a capacidadede suavizao do esquema de baixa ordem e a preciso do esquema de alta ordem. O sucesso

dessa metodologia depende, portanto, da definio de um ni adequado de tal modo que esse

termo seja igual a 1 nas regies suaves e aproximadamente 0 nas regies de gradientes elevados

ou descontinuidades. Existem diversas maneiras de se definir a funo limitadora de fluxo, no

presente trabalho foi convencionadoni = (ni)onde

ni conhecido como um parmetro de

suavidade e pode ser definido como

n

i =

uni uni

1

uni+1 uni (2.21)A funopor sua vez dada por uma funo limitadora de fluxo no linear. A seguir po-

dem ser observadas algumas funes limitadoras de fluxo no lineares TVD de segunda ordem:

Minimod : mn() =max[0,min(1, )]

Osher :os() =max[0,min(, )], (1 2);lim os() =

Superbee:sb() =max[0,min(2, 1),min(2, )]

Van Leer:vl() = +||1+||

25


27/87

2.3 VERIFICAO E VALIDAO DO CDIGO COMPUTACIONAL

Em computao cientfica, a verificao de um cdigo tem por objetivo assegurar que este

represente fielmente o modelo matemtico proposto. A verificao uma etapa fundamental

durante o desenvolvimento de qualquer cdigo computacional. A preciso do resultado num-rico obtido durante a simulao de um problema fsico depender de diversos valores arbitrados

durante o desenvolvimento do cdigo, tais como o tipo de discretizao adotado, o grau de

refinamento da malha computacional, as tolerncias adotadas, entre outros. No entanto, in-

dependente dos parmetros escolhidos, possvel certificar-se de que o cdigo desenvolvido

est livre de erros de programao, atravs de procedimentos de verificao que utilizam os

resultados obtidos pelo prprio cdigo implementado. Por outro lado, a validao do cdigo

computacional e do modelo matemtico proposto analisa se estes so capazes de reproduzir de

forma aceitvel o comportamento do fenmeno fsico que se deseja simular (OBERKAMPF e

Roy, 2010).

A anlise da ordem de erro do mtodo programado um rigoroso critrio de verificao,

onde possvel analisar no somente a convergncia da soluo numrica, como tambm se o

erro de discretizao reduz, com a mesma taxa da ordem terica do mtodo, quando a malha

espacial e/ou temporal refinada. A ordem terica de um cdigo definida atravs de anlise

da ordem do erro de truncamento que surge durante a discretizao. J a ordem real do cdigo

programado pode ser calculada por um procedimento prtico, onde deve ser avaliado se a taxacom que o erro de discretizao decresce condizente com a ordem terica do mtodo. Atravs

da anlise de ordem, possvel identificar diversas fontes de erros que possam existir na pro-

gramao do cdigo implementado. Uma fonte comum de erro durante o desenvolvimento de

um cdigo computacional ocorre nos contornos do domnio, devido imposio de condies

equivocadas. Com a anlise de ordem, possvel verificar se as condies de contorno impostas

esto coerentes com o problema simulado, gerando, com isso, uma soluo de mesma ordem

em todo o domnio (SALARI, 2000).

De forma geral, pode-se dizer que a verificao da ordem de preciso de um cdigo analisa

se as equaes governantes esto sendo resolvidas de forma correta, enquanto que a valida-

o demonstra que as equaes a serem resolvidas representam, com fidelidade, o fenmeno

fsico em questo. Desta forma, fica claro que se deve realizar a verificao do cdigo antes

da validao, uma vez que um cdigo no verificado pode produzir resultados insatisfatrios,

devido a algum erro de programao ou, simplesmente, por terem sido utilizados parmetros

insuficientemente rigorosos durante a soluo numrica.

Segundo SALARI (2000), um mtodo frequentemente empregado para calcular a ordem de

26


28/87

preciso da soluo gerada por um cdigo computacional consiste em executar o cdigo em

questo para trs malhas consecutivamente refinadas. Selecionando-se duas solues, calcula-

se o erro entre as mesmas considerando a soluo obtida pela malha mais refinada como uma

soluo exata. Repete-se o procedimento realizado, utilizando uma nova soluo e uma das

duas solues j empregadas. Com isso, possvel calcular o erro de discretizao em funo

de h, sendo este parmetro a distncia entre os pontos na malha para o clculo da ordem espacial

ou o tamanho do passo no tempo para o clculo da ordem temporal.

Portanto, para qualquer mtodo discreto, sabe-se que a ordem do erro da soluo uma

funo deh, ou seja:

E E(h) (2.22)

Sabendo que o erro de truncamento em uma discretizao da ordem de hp, sendop a ordem

terica da discretizao empregada, para duas malhas consecutivamente refinadas tem-se:

EMalha1 EMalha2 O(hp) (2.23)

EMalha2 EMalha3 O

h

r

p (2.24)

sendo r o refinamento empregado na malha. O produto entre os erros descritos acima dado

por:

EMalha1 EMalha2EMalha2 EMalha3

=rp (2.25)

Com isso, a ordem real do cdigo implementado pode ser obtida conforme a equao (2.26).

Entretanto, cabe ressaltar que no existe uma maneira nica para clculo da ordem, podendo

ser usadas outras frmulas.

p=

log

EMalha1 EMalha2EMalha2 EMalha3

log(r) (2.26)Quando a ordem calculada menor que a ordem terica do mtodo, deve-se rever a pro-

gramao em busca de possveis erros no cdigo. Nos casos em que a ordem real calculada

igual ou aproximadamente igual ordem terica, aps alguns testes com diferentes malhas,

significa que a ordem do cdigo programado foi verificada com sucesso. Um cdigo com or-

dem real verificada deve ser, por fim, validado para que seja considerado confivel. Os casos

teste simulados no presente trabalho so problemas clssicos do estudo de mecnica dos flui-

dos e apresentam soluo analtica. Portanto, no ser uma preocupao a obteno de dados

experimentais para validao dos cdigos computacionais desenvolvidos.

27


29/87

2.4 METODOLOGIA PARA AVALIAO DE EFICINCIA COMPUTACIONAL

Para avaliar o impacto da limitao no linear do passo no tempo, foram desenvolvidos no

presente trabalho estudos sobre eficincia computacional de cada mtodo implementado. Esse

estudo consiste em realizar sucessivas rodadas de simulaes, sempre com o mesmo tempo f-sico final, no entanto, com diferentes valores de passo no tempo. Portanto, partindo do passo

mximo permitido pela barreira no linear, realiza-se sucessivos refinamentos na malha tempo-

ral. Para cada valor de passo no tempo, um par ordenado de erro/tempo computacional obtido.

Consequentemente, quanto menor o passo no tempo, maior ser o tempo computacional neces-

srio para alcanar o tempo final estipulado, no entanto, menor ser o erro agregado a essa

soluo. Abaixo possvel visualizar a essncia da anlise: Na figura 2.1 possvel observar as

CPU Time (s)

L

Error

10-1

100

101

10-7

10-6

10-5

10-4

10-3

10-2

10-1

A

B

C

FIG. 2.1: Curvas tpicas na comparao entre eficincia computacional entre mtodos mtodos.

curvas de custo computacional de dois mtodos hipotticos A, B e C. De forma geral, pode-seextrair a informao de qual mtodo apresenta melhor eficincia traando uma reta paralela ao

eixo das abscissas partindo de um valor de erro considerado aceitvel para a simulao. A curva

que estiver mais prxima do eixo das ordenadas ser cortada primeiro, ou seja, ser necessrio

um tempo computacional menor para realizar tal simulao. No exemplo dado na figura 2.1 no

ponto onde as curvas dos mtodos A e B se cortam os mtodos apresentam a mesma eficincia

computacional, ou seja, ambos os mtodos necessitam do mesmo tempo computacional para

gerar solues com o mesmo erro. No entanto, importante perceber que a escolha do mtodo

mais eficiente vai depender da tolerncia imposta pelo usurio.

28


30/87

3 RESULTADOS

Para alcanar os objetivos estabelecidos no captulo inicial do presente trabalho, foram rea-

lizados, basicamente, dois tipos de anlise. Inicialmente, foi realizada a reproduo de um caso

teste apresentado em MENG et al. (2012) visando demonstrar a grande importncia do emprego

de mtodos SSP na simulao de escoamentos com gradientes abruptos ou descontinuidades e,

posteriormente, foram realizados quatro casos-teste com anlises comparativas de eficincia en-

tre os mtodos implcitos, explcitos e IMEX uma vez que o passo no tempo deve ser limitado

pela barreira SSP. Para tal, foram considerados trs mtodos SSP com estabilidade linearL2 de

segunda ordem de preciso temporal.

Explcito de dois estgios (SHU et al., 2001),

=

0 0

0 0

0 12

=

0 0

1 0

0 12

=

1

1

12

SDIRK (GOTTLIEB et al., 2009a),

=

0 0

0 00 1

= 0 0

0 12

0 1

= 1

10

Implcito-Explcito (PARESCHI e Russo, 2005),

=

0 0

1 0

12

12

=

0

1 2

12

12

=

1

0

29


31/87

onde= 1 12

Para discretizao espacial, foram utilizados esquemas limitatores de fluxo com filtro Su-

perbee no-linear de segunda ordem. Com isso, os fluxos no viscosos de segunda ordem foramreescritos como uma soma de um fluxo de baixa ordem e um fluxo de ordem superior limitado

pelo filtro superbee. Para fluxos negativos:

ux =+3ui 4ui1+ui2

2x (3.1)

= + 1

x(ui ui1) +

sb12x

(ui ui1) sb22x

(ui1 ui2) (3.2)

para fluxos positivos:

u+x =3ui+ 4ui+1 ui+2

2x (3.3)

= 1

x(ui+1 ui)

sb12x

(ui+1 ui) + sb22x

(ui+2 ui+1) (3.4)

onde

sb() =max[0,min(2, 1),min(2, )] (3.5)

Todos os casos teste foram simulados em um computador pessoal Intel(R) Core(TM) i7-

2670QM CPU @2.20GHz, 8GB de memria RAM compilados em plataforma Linux-Ubuntu

12.05. O primeiro caso teste apresentado a seguir tem por objetivo demonstrar a grande impor-

tncia da utilizao de um esquema SSP para simulao de problemas com descontinuidades

e/ou variao abrupta de propriedades. Para tanto, o caso teste selecionado trata da reproduodos resultados apresentados em MENG et al. (2012) onde oscilaes numricas foram errone-

amente tratadas como um fenmeno fsico. Uma carta de correo foi enviada revista em

questo com uma detalhada explicao dos equvocos cometidos na publicao.

30


32/87

3.1 CASO TESTE 01: EQUAO DE BUCKLEY-LEVERETT

A equao de Buckley-Leverett (BL) usualmente empregada na modelagem de escoamen-

tos bifsicos em meios porosos (LEVEQUE (2007)). Em termos numricos, trata-se de uma

equao diferencial parcial hiperblica no linear. Devido sua natureza matemtica, a soluodessa equao tende a apresentar choques e descontinuidades, condies ideais para testar as

propriedades SSP.

Em sua forma mais simples, a equao de BL pode ser escrita na seguinte forma:

tu+x(f(u)) = 0, (3.6)

ondef(u) = u2

u2+M(1u)2 e M >0.

Adicionando efeitos viscosos sua formulao, obtm-se

tu+xf(, u) = xxu. (3.7)

Conforme descrito anteriormente, a equao de Buckley-Leverett usualmente empregada para

modelar escoamentos bifsicos em meios porosos e, em sua forma mais geral possvel con-

siderar tal escoamento com porosidade varivel. Para isso, um termo transiente adicionado

equao (3.7), com isso, chega-se na forma

tu+xf(, u) = xxu+xxtu (3.8)

ondef(, u) = u2

u2+M(u)2 ,0 u ,M >0.

Para avaliar a capacidade do mtodo SSP de resolver o escoamento sem o surgimento de

oscilaes esprias, foram realizadas diversas rodadas de testes com as condies iniciais ex-

tradas de MENG et al. (2012):

IC_01 :u(x, 0) = 1 0.4exp(10(x 5)2), 0 x 10, (3.9)

IC_02 :u(x, 0) = 1 0.9exp(10(x 5)2), 0 x 10, (3.10)

sendoM= 1.

Estudos realizados:

IC_01 := 0, = 105, 106, 107, 108 (3.11)

IC_01 := 105, = 105 (3.12)

(3.13)

31


33/87

IC_02 := 0, = 105, 106, 107, 108 (3.14)

IC_02 := 105, = 105 (3.15)

(3.16)

A seguir so apresentados dois grficos para cada condio inicial descrita acima, no en-

tanto, todos os grficos existentes em MENG et al. (2012) foram reproduzidos e podem ser

encontrados no APNDICE 2. Para cada figura apresentada a seguir pode-se perceber sempre o

grfico original a direita e na esquerda, os grficos desenvolvidos com a finalidade de comparar

as solues obtidas por cada metodologia. Para reproduzir os resultados exatamente como apre-

sentados em MENG et al. (2012), foi implementado o mesmo esquema numrico utilizado pelo

autor em questo, RK4 com uma discretizao centrada de quarta ordem no espao. Alm disso,

uma metodologia SSP explcita de segunda ordem em conjunto com uma tcnica limitadora defluxo para as derivadas espaciais foi implementada para demonstrar a grande importncia de se

utilizar a metodologia capaz de previvir oscilaes numricas.

Para ter certeza que o mtodo foi implementado corretamente e o esquema escolhido capaz

de reproduzir o fenmeno com preciso, a soluo semi-analtica da equao de BL foi obtida

atravs do mtodo das caractersticas. Com isso, possvel perceber na figura 3.1 e 3.2 que o

cdigo SSP foi capaz de reproduzir a soluo da equao de BL com preciso sem oscilaes

numricas. A soluo atravs do mtodo RK4 utilizado em MENG et al. (2012) apresentou

severa oscilao numrica causando perda de preciso na soluo para a segunda condio

inicial como pode ser observado na figura 3.2.

FIG. 3.1: IC_01 := 0, = 0Esquerda:Figura 02 JMP 2012, Direita: Soluo com formulao SSP e semi-analtica MoC.

32


34/87

FIG. 3.2: IC_02 := 0, = 0Esquerda: Figura 06 JMP 2012, Direita: Soluo com formulao SSP e semi-analtica MoC.

Nas figuras 3.3 e 3.4 os efeitos viscosos so considerados e, novamente, possvel perceber

que o mtodo SSP foi capaz de reproduzir o fenmeno fisico com preciso enquanto que o

mtodo no SSP utilizado em MENG et al. (2012) obteve perfis com severa oscilao numrica

e m convergncia para a segunda condio inicial como pode ser observado na figura 3.4.

FIG. 3.3: IC_01 := 0, = 105

Esquerda: Figura 05 JMP 2012, Direita: Soluo com formulao SSP.

Dos resultados apresentados nesse primeiro caso teste possvel perceber a grande im-

portncia de se empregar um esquema SSP para simular problemas com variao abrupta de

propriedades e/ou descontinuidades. Tendo em mente a grande importncia dessa classe de

mtodos, o restante deste captulo dedicado a investigao numrica com fim de comparar a

33


35/87

FIG. 3.4: IC_02 := 0, = 105

Esquerda:Figura 09 JMP, Direita: Soluo com formulao SSP

eficincia computacional entre mtodos implcitos, explcitos e IMEX considerando a restrio

no linear de passo no tempo.

34


36/87

3.2 CASO TESTE 02: EQUAO DE BURGERS NO VISCOSA

O segundo caso-teste apresentado consiste na soluo da equao no viscosa de Burges. A

equao de Burgers uma equao diferencial parcial no-linear do tipo conveco-difuso e

considerada como uma forma simplificada das equaes de Navier-Stokes para os casos emque o gradiente da presso pode ser desprezado. Para o caso no viscoso, escreve-se

ut+f(u)x = 0, (3.17)

ondef(u) = u2

2. Para o primeiro caso teste realizado foram adotadas condies de contorno

peridicas e, a condio inicial

u(x, 0) =1

2 1

4 sin(x). (3.18)

O termo advectivo da equao de Burgers no-linear d origem a uma onda de choque que

se propaga com sentido determinado pelo mdulo da velocidade. Portanto, esse modelo traz

diversas oportunidades para estudo das propriedades SSP. Outra implementao interessante

nesse caso teste foi o uso de condies de contorno peridicas, uma vez que o emprego de tais

condies torna possvel avaliar resultados com longos tempos de simulao sem a necessidade

do emprego de um domnio demasiadamente extenso.

As simulaes foram realizadas em um domnio adimensional [0,1], dividido em N=802pontos para os tempos finais TF1 = 1, 2s; TF2 = 2, 4s; TF3 = 4, 8s. Para avaliar a soluo

numrica obtida, a soluo exata da equao de Burgers no viscosa foi calculada atravs do

mtodo das caractersticas fazendo:

u(x, t)/u0 = u(x f(u)t, 0) (3.19)

A partir dos resultados apresentados na figura 3.5 possvel visualizar a formao da onda

de choque e a eficincia da metodologia SSP na captura de tal fenmeno sem a ocorrncia de

oscilaes numricas. EmTF3 possvel perceber a onda de choque se propagando ao longo

do domnio devido as condies de contorno peridicas.

35


37/87

x / L

u/umax

0 0.2 0.4 0.6 0.8 1

MoC

SSPRK(2,2)SDIRK(1,2)

t = 2.4 s

t = 1.2 s

t = 4.2 s

FIG. 3.5: Perfis de velocidade : burges no viscoso.

3.2.1 VERIFICAO COMPUTACIONAL

Uma vez ciente de que as rotinas programadas so capazes de reproduzir o fenmeno fsico

desejado, a prxima etapa do estudo consiste na avaliao da ordem de erro obtida por cada

cdigo. De forma geral, pode-se entender ordem de erro como a taxa com que o erro numrico

reduz mediante refinamentos sucessivos de malha. Todos os mtodos programados apresen-

tam segunda ordem temporal/espacial e, para verificao da ordem de erro real foi realizado o

procedimento de clculo apresentado na seo 2.3.

Para avaliao da ordem de preciso temporal de cada cdigo programado inicialmentedeve-se garantir que o erro dominante o temporal. Para isso, fixou-se o passo no tempo igual

a102 e sucessivos refinamentos de malha espacial foram realizados. Na figura 3.6 possvel

observar a estagnao do erro espacial a medida que a malha refinada.

Garantindo que o erro dominante o temporal, os cdigos implementados foram verificados

e podem ser visualizadas nas figuras 3.7 , 3.8.

36


38/87

dx

Error

0.01 0.02 0.03 0.04

6E-05

6.5E-05

7E-05

7.5E-05

8E-05

8.5E-05

9E-05

FIG. 3.6: Convergncia de malha espacial, burges no viscoso.

(a) Ordem Temporal

(b) Ordem Espacial

FIG. 3.7: Verificao de cdigo SSP(2,2) com Superbee.

37


39/87

(a) Ordem Temporal

(b) Ordem Espacial

FIG. 3.8: Verificao de cdigo SDIRK(1,2) com Superbee.

3.2.2 ESTUDO DE PRESERVAO DE ESTABILIDADE NO LINEAR-BARREIRA SSP

Como descrito no captulo anterior, para que o mtodo SSP seja capaz de preservar pro-

priedades no lineares, a restrio no linear do passo no tempo deve ser respeitada. Com o

intuito de verificar numericamente a barreira SSP de cada mtodo implementado, foram reali-

zadas anlises com diferentes passos no tempo para cada cdigo. Inicialmente, reescrevendo o

mtodo explicito RKSSP(2,2) na forma Shu-Osher, obtm-se:

k

(1)

=u

n

+ tF(u

n

), (3.20)un+1 =

1

2un +

1

2k(1) +

1

2tF(k(1)). (3.21)

atravs do Teorema 1 (SHU, 1988) sabe-se que o coeficiente SSP desse mtodo C = 1,

reescrevendo tambm o mtodo implcito na forma Shu-Osher, obtm-se:

k(1) =un +1

2tF(k(1)), (3.22)

un+1 =k(1) +12

tF(k(1)), (3.23)

38


40/87

procedendo de forma anloga, calcula-se o coeficiente SSP,C = 2.

Alm do mtodo escolhido para a marcha temporal, a tcnica empregada na discretizao

espacial outro fator determinante para a definio do passo mximo no tempo. Os termos

espaciais em todos os casos teste desse estudo foram discretizados atravs do esquema limitador

de fluxo com filtro superbee (segunda ordem no-linear). Como descrito no captulo anterior,

quando a soluo torna-se abrupta o limitador de fluxo reduz os valores dos multiplicadores de

ordens superiores at o limite em que a discretizao espacial torna-se um esquema upwind de

primeira ordem. Dessa forma, o termo espacial da equao de burgers no viscoso reduz-se a:

ux =uj uj1

x (3.24)

portanto, empregando o mtodo de Euler avanado para discretizao do termo temporal da

equao de Burgers, obtm-se:

un+1j =unj

t

x(unj u

nj1) (3.25)

atravs da anlise de estabilidade linear (von Neumann), obtm-se que a equao 3.25 estvel

para tx 1. Portanto, espera-se que o cdigo explcito seja capaz de realizar simulaes sem

a presena de oscilaes numricas para t = 1t, enquanto que para o mtodo implcito,

espera-se empregar o mtodo SSP com sucesso para um valor mximo de passo no tempo

t= 2t.Inicialmente, foram realizadas simulaes com tempo fsico finalTF = 3spara o mtodo

explcito com o passo no tempo limite, ou seja, na barreira SSP tericadt = dx.

Como pode ser observado nas figuras 3.9(a)-(b), para as simulaes utilizandodt mximo

(barreira SSP) o mtodo numrico preserva ambas ordem de erro e variao total diminuda

mostrando boa consistncia com a teoria SSP. A segunda rodada de testes foi realizada nova-

mente para um tempo final TF = 3s, no entanto, com um passo no tempo acima da barreira

SSP terica (dt= 1.2dx).

39


41/87

x/L

AverageOrder

0.25 0.5 0.75 1

1.25

1.5

1.75

2

2.25

2.5

2.75

3

(a) Ordem temporal preservada.

time(s)

TV

0.5 1 1.5 2 2.5 3

0.88

0.89

0.9

0.91

0.92

0.93

0.94

0.95

0.96

0.97

0.98

0.99

1

(b) Propriedade TVD preservada.

x/L

Ave

rageOrder

0.25 0.5 0.75 1

0.75

1

1.25

1.5

(c) Perda de Ordem temporal

time(s)

TV

0.5 1 1.5 2 2.5 3

0.88

0.89

0.9

0.91

0.92

0.93

0.94

0.95

0.96

0.97

0.98

0.99

1

(d) Propriedade TVD preservada.

40


42/87

Nas figuras 3.9(c)-(d), percebe-se que ao se ultrapassar a barreira em 20% de seu valor

terico, o mtodo explcito continua fornecendo uma soluo TVD, no entanto, a ordem de erro

da soluo j degradada. Essa perda de ordem evidencia que mesmo antes da percepo das

oscilaes numricas, o erro causado pelo aumento do passo no tempo causa a perda da preciso

original do mtodo. Nesse caso, a perda de ordem tem incio no choque, onde o maior valor de

erro ocorre e, ao longo do processo iterativo, esse erro se propaga para todo o domnio.

Por fim, uma terceira rodada de simulaes foi executada comdt = 1.4dx.

x/L

AverageOrder

0.25 0.5 0.75 1 1.25 1.5 1.75 2-17

-16

-15

-14

-13

-12

-11

-10

-9

-8

-7

-6

(e) Perda de preciso temporal

time(s)

TV

0.5 1 1.5 2 2.5 31

2

3

4

5

6

7

8

9

10

11

12

13

(f) Propriedade TVD violada

Conforme pode-se observar nas figuras 3.9(e)-(f), para um valor de passo no tempo 40%

acima da barreira SSP, o resultado perde a preciso temporal e a propriedade TVD. Com isso,

possvel observar o surgimento de oscilaes numricas (3.9(g)-(h)) que foram amplificadas ao

longo da simulao e causaram divergncia do cdigo emTF = 3, 5s.

Analogamente para o mtodo implcito foram realizados testes para avaliao da barreira

SSP terica. Inicialmente, foram realizados testes no limite da barreira comdt = 2.0dx.

41


43/87


44/87

A partir dos resultados apresentados nas figuras 3.9(i)-(j), possvel perceber que a ordem

de erro foi mantida como tambm a propriedade TVD, observando-se, com isso, boa consis-

tncia com a anlise numrica terica. Em seguida realizou-se uma nova rodada de experi-

mentos numricos acima da barreira SSP terica. Para essa nova rodada de testes utilizou-se

dt= 3.0dx.

x/L

Avera

geOrder

0 0.5 1 1.5 2

1.25

1.5

1.75

2

2.25

2.5

2.75

3

(k) Ordem temporal preservada

Time (s)

TV

1 2 3

0.92

0.93

0.94

0.95

0.96

0.97

0.98

0.99

1

1.01

1.02

1.03

dt = 3 dx

(l) Propriedade TVD violada

Como pode ser observado nas figuras 3.9(k)-(l), a ordem de erro manteve-se de segunda

ordem mesmo com um valor de passo no tempo 50%maior que a barreira SSP, no entanto, a

propriedade TVD foi violada e, portanto, oscilaes numricas se desenvolveram durante esse

experimento como pode ser observado nas figuras 3.9(m) e (n). possvel perceber, portanto,

que mesmo com um valor elevado de erro capaz de originar oscilaes esprias na regio do

choque, o mtodo implcito capaz de preservar a ordem ao longo do domnio, mostrando, com

isso, uma estabilidade superior a apresentada pelo mtodo explcito.

Da breve anlise apresentada possvel concluir que os cdigos programados esto deacordo com a teoria e que a barreira SSP numrica ligeiramente maior que a terica con-

forme j mencionado na literatura (GOTTLIEB et al. (2009a)).

43


45/87


46/87

3.2.3 ESTUDO COMPARATIVO DE EFICINCIA COMPUTACIONAL

Conforme descrito no captulo anterior, para cada metodologia implementada um estudo de

eficincia computacional foi realizado levando em considerao a restrio no linear do passo

no tempo. Dessa forma, para a equao de Burgers no viscosa, com contornos peridicos,foram realizadas avaliaes de eficincia computacional emt =1,2; 2,4 e 4,8 segundos.

Considerando a metodologia de anlise apresentada na seo 2.4, pode-se observar na figura

3.9 um desempenho superior do mtodo explcito para as simulaes com 1,2 segundos. Uma

exceo ocorre para o ponto onde t = 1, 2xcom o mtodo explcito et = 2, 4xcom

o mtodo implcito. Nesse ponto, ambos os mtodos apresentam eficincia equivalente. Os

valores de passo no tempo, erro obtido e tempo de CPU podem ser observados na tabela 3.1.

CPU Time (s)

L

Error

10-2

10-1

100

10-6

10-5

10-4

10-3

10-2

10-1

SDIRK (1,2)

RKSSP (2,2)

FIG. 3.9: Eficincia computacional mtodos explicito x implcito 1,2s burgers no viscoso.

Na segunda rodada de simulaes foram realizados testes para avaliar a eficincia com-

putacional em 2,4 s. Como pode ser observado na figura 3.10, o mtodo explcito apresenta

desempenho superior para todos os passos testados. Observa-se tambm que, a superioridade

do mtodo explcito tende a crescer a medida que o passo no tempo diminui.Isso se deve prin-

cipalmente ao fato do mtodo implcito ter um custo maior por iterao, portanto, quanto maior

o nmero de iteraes, maior ser a diferena entre os custos finais dos mtodos. Os valores de

tamanho de passo no tempo, erro obtido e tempo de CPU podem ser observados na tabela 3.2

45


47/87

Explcito Implcitot CPU(s) ERRO L CPU(s) ERRO L

4.8x nan nan 1.56E-02 1.59E-022.4x nan nan 4.67E-02 5.99E-031.2x 4.67E-02 5.70E-003 9.36E-02 1.95E-03

0.6x 9.36E-02 8.00E-004 0.15E-00 6.40E-040.3x 0.18E-00 1.98E-004 0.32E-00 2.18E-04

0.15x 0.35E-00 4.71E-005 0.68E-00 7.57E-050.075x 0.71E-00 9.43E-006 1.32E-00 2.46E-05

TAB. 3.1: Eficincia computacional no Burgers Viscoso para 1.2 s

CPU Time (s)

L

Error

10-1

100

101

10-4

10-3

10-2

10-1

SDIRK (1,2)

RKSSP (2,2)


Explcito Implcito

t CPU(s) ERRO L CPU(s) ERRO L4.8x nan nan 4.68E-02 0.24E-002.4x nan nan 9.36E-02 8.50E-021.2x 9.36E-02 2.74-02 0.20E-00 2.00E-020.6x 0.17E-00 8.12E-03 0.34E-00 7.10E-030.3x 0.35E-00 1.82E-03 0.67E-00 2.80E-03

0.15x 0.73E-00 4.07E-04 1.34E-00 1.10E-030.075x 1.46E-00 7.81E-05 2.62E-00 3.90E-04

TAB. 3.2: Eficincia computacional Burgers no Viscoso para 2.4 s

46


48/87

Por fim, simulaes foram realizadas para 4,8 segundos. Como pode ser observado na figura

3.11 o mtodo explcito apresentou-se mais eficiente que o mtodo implcito. Aps uma inves-

tigao numrica rigorosa, foi possvel concluir que tal superioridade se deve, principalmente,

por dois motivos: condies de contorno peridicas e devido tcnica de discretizao espacial

empregada.

CPU Time (s)

L

Error

10-1

100

101

10-4

10-3

10-2

10-1

SDIRK (1,2)

RKSSP (2,2)


Explcito Implcitot CPU(s) ERRO L CPU(s) ERRO L

4.8x nan nan 7.80E-02 0.15E-002.4x nan nan 0.17E-00 4.69E-02

1.2x 0.18E-00 1.84E-02 0.37E-00 1.39E-020.6x 0.37E-00 4.97E-03 0.67E-00 5.05E-030.3x 0.73E-00 1.16E-03 1.32E-00 2.01E-03

0.15x 1.46E-00 2.67E-04 2.63E-00 7.91E-040.075x 2.90E-00 5.21E-05 5.21E-00 2.52E-04

TAB. 3.3: Eficincia computacional Burgers no Viscoso para 4.8 s

Para resolver problemas peridicos atravs de uma metodologia implcita, deve-se inserir

tais condies na matriz implcita a ser resolvida. Com isso, a matriz esparsa que antes era

resolvida com o mtodo de Thomas para soluo de matrizes tridiagonais, agora deve ser resol-

47


49/87

vida por alguma tcnica especfica para sistemas cclicos pentadiagonais que chegam a deman-

dar at quatro vezes mais tempo computacional para serem resolvidas (NAVON (1987), PRESS

(2007)). J o emprego da tcnica de discretizao espacial limitadora de fluxo aumenta con-

sideravelmente o nmero de clculos realizados na matriz implcita causando, dessa forma,um

aumento no erro produzido pela soluo final da metodologia implcita em conjunto com essa

tcnica.

48


50/87

3.3 CASO TESTE 03: EQUAO DE BURGERS VISCOSO

O terceiro caso teste trata novamente da equao de Burgers, no entanto, efeitos viscosos

foram includos. Com isso, o termo responsvel por computar os efeitos dissipativos so adici-

onados equao (3.19).ut+f(u)x = g(u)xx, (3.26)

Na figura 3.12 pode-se observar solues da equao de burgers viscoso em diferentes tempos

considerando uma condio inicial senoidal

u(x, 0) =sin(x). (3.27)

alm de condies de contorno de primeiro tipo. As simulaes foram realizadas em um dom-

nio adimensional [0,1], dividido em N=801 pontos para os tempos finais TF1 = 0, 1s; TF2 =

0, 5s; TF3= 1, 0s.

x / L

u

/U

0 0.25 0.5 0.75 10

0.2

0.4

0.6

0.8

1

SSP (2,2)

SDIRK (1,2)IMEX

0,1s

0,5s

1,0s

FIG. 3.12: Perfis de velocidade : burges viscoso.

Como pode ser observado na figura 3.12, o fenmeno foi simulado tambm pelo mtodo h-

brido IMEX sendo resolvido o termo no linear de forma explcita e o termo difusivo de forma

49


51/87

implcita. Com isso, possvel obter alguns benefcios como uma maior estabilidade num-

rica alm de minimizar erros, uma vez que no necessrio realizar a linearizao do termo

advectivo. Considerando que o termo difusivo foi adicionado em todos os cdigos validados

anteriormente, foram realizadas, novamente, as validaes dos cdigos implcito, explcito e

IMEX. Todos os grficos podem ser observados no APNDICE 2.

3.3.1 ESTUDO COMPARATIVO DE EFICINCIA COMPUTACIONAL

Para a primeira rodada de simulaes emTF = 0, 1s, o choque ainda no est formado e,

como pode-se observar na figura 3.13, todos os mtodos tm um desempenho semelhante. No

entanto, percebe-se que o mtodo implcito perde eficincia perante o explcito e o IMEX, a

medida que o tamanho do passo no tempo reduzido. Novamente, isso se deve ao fato do m-

todo implcito realizar um maior nmero de operaes longas e tendo que resolver uma matriz

pentadiagonal a cada iterao. Outra informao relevante observada foi a grande estabilidade

numrica apresentada pelo mtodo hbrido que suportou um passo 64x maior que o mtodo

explcito e 16x maior que o mtodo implcito. Os valores de tamanho de passo no tempo em-

pregado, erro obtido e tempo de CPU podem ser observados na tabela 3.4

CPU Time (s)

L

Error

10-2

10-1

100

101

10-13

10-12

10-11

10-10

10-9

10-8

10-7

10-6

10-5

10-4

IMEX (2,2,2)

SDIRK (1,2)

RKSSP (2,2)

FIG. 3.13: Comparativo: Eficincia de CPU Burgers Viscoso 0,1s.

50


52/87

Explcito Implcito IMEXt CPU(s) ERRO L CPU(s) ERRO L CPU(s) ERRO L

x / 0.25 nan nan nan nan 3.89E-03 3.48E-04x / 4.00 nan nan 1.11E-02 1.39E-07 6.24E-02 1.29E-06x / 8.00 nan nan 3.12E-02 3.56E-08 - -

x / 16.0 3.12E-02 2.77E-08 4.68E-02 9.287E-09 9.36E-02 5.69E-08x / 32.0 6.24E-02 6.94E-09 0.10E-00 2.512E-09 0.20E-00 1.42E-08x / 64.0 0.15E-00 1.73E-09 0.24E-00 7.22E-10 0.37E-00 3.61E-09x / 128 0.32E-00 4.34E-10 0.48E-00 2.27E-10 0.78E-00 9.08E-10x / 256 0.60E-00 1.08E-10 0.96E-00 8.02E-11 1.57E-00 2.27E-10x / 512 1.18E-00 2.71E-11 1.93E-00 3.16E-11 3.12E-00 5.70E-11

x / 1024 2.37E-00 7.14E-12 3.85E-00 1.35E-11 6.19E-00 1.42E-11x / 2048 4.78E-00 3.21E-12 7.73E-00 6.06E-12 12.38E-00 3.55E-12x / 4096 9.53E-00 1.10E-12 15.41E-00 2.72E-12 25.33E-00 8.79E-13x / 8192 19.20E-00 2.83E-13 30.76E-00 1.14E-12 49.78E-00 2.96E-13

x / 16384 38.23E-00 9.41E-14 61.74E-00 3.80E-13 99.67E-00 4.52E-14

TAB. 3.4: Eficincia computacional em 0.1 s

Em 0,5s o choque comea a se formar e, com isso, possvel observar na figura 3.14 a maior

eficincia computacional dos mtodos implcitos para valores de erro entre [107, 1011]. Para

erros menores que 1011 os mtodos explcitos passam a ser mais eficintes que os mtodos

implcitos.Os valores de tamanho de passo no tempo empregado, erro obtido e tempo de CPU

podem ser observados na tabela 3.5.Explcito Implcito IMEX

t CPU(s) ERRO L CPU(s) ERRO L CPU(s) ERRO Lx / 0.25 nan nan nan nan 0.01 5.32E-03x / 4.00 nan nan 0.07 2.01E-08 0.12 6.66E-05x / 8.00 nan nan 0.17 6.03E-09 - -x / 16.0 0.20 1.24E-07 0.29 2.02E-09 0.48 4.81E-06x / 32.0 0.37 3.10E-08 0.59 7.29E-10 0.95 1.23E-06x / 64.0 0.76 7.76E-09 1.18 2.94E-10 1.91 3.13E-07

x / 128 1.52 1.94E-09 2.41 1.29E-10 3.83 7.88E-08

x / 256 3.04 4.85E-10 4.77 5.99E-11 7.65 1.97E-08x / 512 6.11 1.21E-10 9.56 2.86E-11 15.33 4.95E-09

x / 1024 12.19 2.92E-11 19.14 1.38E-11 30.60 1.23E-09x / 2048 24.21 1.11E-11 38.43 6.65E-12 61.19 3.08E-10x / 4096 48.53 4.11E-12 76.86 3.09E-12 122.32 7.63E-11x / 8192 95.89 1.02E-12 153.83 1.33E-12 244.74 1.81E-11

x / 16384 191.49 2.41E-13 307.35 4.67E-13 489.26 3.64E-12


51


53/87

CPU Time (s)

L

Error

10-2

10-1

100

101

102

10-11

10-9

10-7

10-5

10-3

IMEX (2,2,2)SDIRK (1,2)

RKSSP (2,2)


Em 1,0s, o choque j est completamente formado e como pode-se observar na figura 3.15,

os mtodos implcito e explcito necessitam de um passo no tempo consideravelmente menor

que o IMEX para convergir. Com isso, para erros maiores que 108 os mtodos hbridos so

a nica opo. Para solues com erro entre [108, 1010], os mtodos implcitos so os mais

eficientes. No entanto, devido ao seu custo elevado por iterao, para erros menores que 1010,

os mtodos explcitos tornam-se os mais indicados. Os valores de tamanho de passo no tempo

empregado, erro obtido e tempo de CPU podem ser observados na tabela 3.6.

52


54/87

CPU Time (s)

L

Error

10-1

100

101

102

10310

-14

10-12

10-10

10-8

10-6

10-4

IMEX (2,2,2)

SDIRK (1,2)

RKSSP (2,2)


Explcito Implcito IMEXt CPU(s) ERRO L CPU(s) ERRO L CPU(s) ERRO L

x / 0.25 nan nan nan nan 0.01 1.74E-03x / 4.00 nan nan 0.15 1.34E-08 0.23 2.32E-05x / 8.00 nan nan 0.29 2.50E-09 - -x / 16.0 0.37 1.13E-08 0.59 1.96E-10 0.96 1.68E-06x / 32.0 0.73 2.83E-09 1.18 1.64E-10 1.91 4.33E-07x / 64.0 1.48 7.09E-10 2.35 1.48E-10 3.85 1.09E-07

x / 128 2.97 1.77E-10 4.71 9.03E-11 7.72 2.76E-08x / 256 6.02 4.42E-11 9.453 4.90E-11 15.4 6.94E-09x / 512 11.90 1.10E-11 18.87 2.53E-11 30.88 1.74E-09

x / 1024 23.79 3.63E-12 37.75 1.27E-11 61.65 4.35E-10x / 2048 47.51 1.50E-12 75.53 6.21E-12 123.44 1.08E-10x / 4096 95.22 4.04E-13 151.07 2.90E-12 246.69 2.68E-11x / 8192 190.08 1.04E-13 302.17 1.24E-12 493.39 6.38E-12

x / 16384 380.81 3.26E-14 603.81 4.06E-13 987.40 1.27E-12


53


55/87

3.4 CASO TESTE 04: EQUAO DE BUCKLEY-LEVERETT

O quarto caso teste apresentado consiste na soluo da equao de Buckley-Leverett (BL).

Essa equao usualmente empregada na modelagem de escoamentos bifsicos em meios po-

rosos. Em termos numricos, alguma complexidade adicionada para resolver tal formulaoquando comparada com a equao de Burgers, uma vez que a funo no-linear empregada na

formulao de Burges agora substituda por f(u) = u2

u2+M(1u)2 . Portanto, foi resolvida

ut+f(u)x = 0, (3.28)

considerando, M = 0.25, condies de contorno prescritas e as condies iniciais

u(x, 0) = 1.0 ifx 0.10 ifx >0.1.

(3.29)

Na figura 3.16, so apresentados os perfis obtidos para trs tempos finais de simulao t =

1,0s ; 3,0s e 5,0s em um domnio adimensional [0,1] dividido em N=801 pontos. Com isso,

pode-se observar novamente a capacidade das metodologias SSP de reproduzir fenmenos com

elevados gradientes com acurcia e sem a ocorrncia de oscilaes numricas. Para verificar se

as solues numricas obtidas so as corretas foram geradas solues semi-analticas, a partir

do mtodo das caractersticas.

x / L

u

/U

0 0.2 0.4 0.6 0.8 1

0

0.2

0.4

0.6

0.8

1

EXPLICIT

IMPLICIT

M

Documents

Dissertacao Gabriel