View
220
Download
0
Category
Preview:
Citation preview
7/26/2019 Dissertao Porosos 229
1/88
Luccas Cassimiro Campos
Modelagem do escoamento de fluidos em meios
porosos utilizando a estrutura de dadosAutonomous
Leaves Graph
Belo Horizonte - MG, Brasil
Dezembro de 2013
7/26/2019 Dissertao Porosos 229
2/88
Luccas Cassimiro Campos
Modelagem do escoamento de fluidos em meios
porosos utilizando a estrutura de dadosAutonomous
Leaves Graph
Dissertao apresentada para obteno do
Grau de Mestre em Matemtica pela
Universidade Federal de Minas Gerais
Orientador:
Rodney Josu Biezuner
Co-orientadora:
Denise Burgarelli Duczmal
DEPARTAMENTO DE MATEMTICA
INSTITUTO DE CINCIAS EXATAS
UNIVERSIDADE FEDERAL DE MINAS GERAIS
Belo Horizonte - MG, Brasil
Dezembro de 2013
7/26/2019 Dissertao Porosos 229
3/88
Dissertao de Mestrado sob o ttulo Modelagem do escoamento de fluidos em meios
porosos utilizando a estrutura de dadosAutonomous Leaves Graph, defendida por Luccas
Cassimiro Campos e aprovada em 09/12/2013, em Belo Horizonte, Estado de Minas
Gerais, pela banca examinadora constituda pelos professores:
Prof. Rodney Josu BiezunerOrientador
Profa. Denise Burgarelli Duczmal
Co-orientadora
Prof. Ricardo Hiroshi Caldeira TakahashiUniversidade Federal de Minas Gerais
Prof. Eduardo Cardoso de AbreuUniversidade Estadual de Campinas
7/26/2019 Dissertao Porosos 229
4/88
Agradecimentos
Agradeo aos meus pais e minha famlia pelo incentivo incondicional aos meus estudos.
Aos meus orientadores, Rodney e Denise, por acreditarem no meu trabalho e me empres-tarem sua experincia durante os ltimos anos.
Ao CNPq e CAPES, pela criao do PICME e pela concesso da bolsa durante o perodo
do mestrado.
Sylvie por levar adiante o projeto do PICME na UFMG e permitir a realizao desse
trabalho.
Aos professores e funcionrios do Departamento de Matemtica, por todo o apoio pres-
tado.
Por fim, agradeo Universidade Federal de Minas Gerais por ter me acolhido como aluno
desde o bacharelado.
4
7/26/2019 Dissertao Porosos 229
5/88
Foi o tempo que perdeste com
tua rosa que fez tua rosa to
importante
(Antoine de Saint-Exupry)
7/26/2019 Dissertao Porosos 229
6/88
Resumo
A simulao do escoamento multicomponente de fluidos em meios porosos requer o trata-
mento de fenmenos localizados, tais como frentes de concentrao e gradientes de presso,
assim como a considerao da geometria do meio. O fato de serem localizados nos sugere
a utilizao de uma malha adaptativa, altamente refinada nas regies onde necessriauma melhor resoluo, e pouco refinada onde o fluxo de massa menor. O presente traba-
lho tem como objetivo a utilizao da metodologia conhecida como Autonomous Leaves
Graph (ALG) na simulao do escoamento de fluidos em meios porosos. A estrutura ALG
tem se mostrado eficiente em problemas que envolvem refinamento adaptativo da malha,
permitindo um refinamento local e a coexistncia de regies vizinhas com nveis de refi-
namento arbitrariamente diferentes. Utilizaremos o mtodo misto-hbrido de elementos
finitos para a soluo da Equao Diferencial Parcial elptica dada pela Lei de Darcy, e o
mtododonor cell upwind schemede volumes finitos para avanar no tempo a concentra-
o dotracer(marcador) em cada clula. Aplicamos a tcnica de refinamento adaptativo
simulao do escoamento monofsico bicomponente em um meio poroso. So discutidos
os resultados numricos obtidos na simulao do escoamento bidimensional.
Palavras-chave: Escoamento em meios porosos. Refinamento adaptativo da malha
(AMR). Autonomous Leaves Graphs(ALG). Mtodo Misto-Hbrido. Donor cell upwind
scheme.
6
7/26/2019 Dissertao Porosos 229
7/88
Abstract
Multicomponent flow through porous media simulation involves the treatment of localized
phenomena, such as concentration fronts, pressure gradients or geometry of the media.
This suggests the use of an adaptive mesh, that is highly refined in the regions on which
a better resolution is needed and less refined where the mass flow is small. The presentwork intends to use the methodology known as Autonomous Leaves Graph (ALG) to
simulate fluid flow through porous media. The ALG structure have been showing being
efficient on problems that involve adaptive mesh refinement, allowing local refining and
the coexistence of neighbor regions with arbitrarily different refinement levels. We use
a mixed-hybrid finite element method to solve the elliptic Partial Differential Equation
given by Darcys Law, and the donor cell upwind schemefinite volume method to update
tracer concentrations within each cell. The adaptive refinement technique is applied to
monofasic bicomponent flow through a porous medium and numerical results of the bidi-
mensional flow simulation are discussed.
Keywords: Flow through porous media. Adaptive mesh refinement (AMR). Autono-
mous Leaves Graphs(ALG). Hybrid mixed method. Donor cell upwind scheme.
7
7/26/2019 Dissertao Porosos 229
8/88
Contedo
Lista de Figuras 10
Lista de Tabelas 12
Lista de Smbolos 13
Introduo 13
1 Mecnica dos Fluidos 16
1.1 Viscosidade . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
1.2 Hiptese do contnuo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17
1.3 Meios Porosos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17
1.4 Lei de Darcy. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 181.5 Equao da Continuidade . . . . . . . . . . . . . . . . . . . . . . . . . . . 19
1.6 Fluxo Mssico . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19
1.7 Conservao do Volume . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20
2 Mtodos Numricos para a resoluo de EDPs 21
2.1 O Mtodo de Galerkin . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
2.1.1 Propriedades da aproximao de Galerkin . . . . . . . . . . . . . . 25
2.2 O Mtodo dos Elementos Finitos (FEM) . . . . . . . . . . . . . . . . . . . 26
2.3 O Mtodo dos Volumes Finitos (FVM) . . . . . . . . . . . . . . . . . . . . 282.4 Mtodos numricos para a resoluo de Sistemas Lineares. . . . . . . . . . 28
2.4.1 O Mtodo do Gradiente Conjugado . . . . . . . . . . . . . . . . . . 28
2.4.2 O Mtodo do Gradiente Biconjugado Estabilizado . . . . . . . . . . 32
2.4.3 A fatorao LU incompleta . . . . . . . . . . . . . . . . . . . . . . 33
3 Discretizao da Lei de Darcy 34
3.1 Formulaes do problema elptico . . . . . . . . . . . . . . . . . . . . . . . 35
3.2 A formulao mista . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36
3.2.1 Existncia, unicidade e convergncia . . . . . . . . . . . . . . . . . 37
3.3 A formulao mista-hbrida . . . . . . . . . . . . . . . . . . . . . . . . . . 44
8
7/26/2019 Dissertao Porosos 229
9/88
3.4 A formulao mista-hbrida e o escoamento em meios porosos. . . . . . . . 47
3.5 Construo do sistema linear. . . . . . . . . . . . . . . . . . . . . . . . . . 48
3.6 Positividade definida do sistema linear . . . . . . . . . . . . . . . . . . . . 49
3.7 Implementao numrica do mtodo hbrido-misto . . . . . . . . . . . . . . 50
4 Discretizao das equaes convectivo-difusivas 53
4.1 Fluxo convectivo em 1D . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54
4.2 Fluxo convectivo em 2D . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56
4.3 Fluxo difusivo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59
5 Refinamento Adaptativo da Malha 60
5.1 Malhas Adaptativas. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60
5.2 A estrutura de dadosAutonomous Leaves Graph(ALG). . . . . . . . . . . 625.2.1 A estrutura da malha. . . . . . . . . . . . . . . . . . . . . . . . . . 62
5.3 Refinamento da malha . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63
5.4 Desrefino da malha . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65
5.5 Ordenao Total da Malha . . . . . . . . . . . . . . . . . . . . . . . . . . . 66
5.6 Refinando e Desrefinando. . . . . . . . . . . . . . . . . . . . . . . . . . . . 67
6 Implementao e resultados numricos 70
6.1 Clculo do campo de presses e de velocidades . . . . . . . . . . . . . . . . 70
6.1.1 Exemplo Numrico 1 . . . . . . . . . . . . . . . . . . . . . . . . . . 716.1.2 Exemplo Numrico 2 . . . . . . . . . . . . . . . . . . . . . . . . . . 71
6.1.3 Exemplo Numrico 3 . . . . . . . . . . . . . . . . . . . . . . . . . . 74
6.2 Transporte de massa e refinamento adaptativo . . . . . . . . . . . . . . . . 76
6.2.1 Exemplo Numrico 1 . . . . . . . . . . . . . . . . . . . . . . . . . . 76
6.2.2 Exemplo Numrico 2 . . . . . . . . . . . . . . . . . . . . . . . . . . 78
6.2.3 Exemplo Numrico 3 . . . . . . . . . . . . . . . . . . . . . . . . . . 79
7 Concluses 817.1 Perspectivas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81
7.1.1 Permeabilidade Heterognea . . . . . . . . . . . . . . . . . . . . . . 81
7.1.2 Poos (wells) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82
7.1.3 Escoamento Multifsico . . . . . . . . . . . . . . . . . . . . . . . . 82
7.1.4 Escoamento Tridimensional . . . . . . . . . . . . . . . . . . . . . . 83
7.2 Principais Contribuies . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83
9
7/26/2019 Dissertao Porosos 229
10/88
7/26/2019 Dissertao Porosos 229
11/88
6.4 Presena de clulas impermeveis . . . . . . . . . . . . . . . . . . . . . . . 75
6.5 Campo de concentraes em diferentes intervalos de tempo . . . . . . . . . 77
6.6 Campo de concentraes em diferentes intervalos de tempo . . . . . . . . . 78
6.7 Campo de concentraes em diferentes intervalos de tempo . . . . . . . . . 80
11
7/26/2019 Dissertao Porosos 229
12/88
Lista de Tabelas
6.1 Tempo gasto para atualizar as concentraes . . . . . . . . . . . . . . . . . 76
6.2 Tempo gasto para atualizar as concentraes . . . . . . . . . . . . . . . . . 79
12
7/26/2019 Dissertao Porosos 229
13/88
Lista de Smbolos
Principais caracteres romanos
p presso em atm =1.01325 105 PaK tensor permeabilidade, com cada entrada em milidarcy (mD) =9.869 1016m2v velocidade de Darcy, dada por v= u, em que u a velocidade do fluidoc frao volumtrica, em geral, dotracer
D parte do domnio na qual a presso especificadaN parte do domnio na qual a velocidade especificada
Principais caracteres gregos
densidade em g cm3
viscosidade em cP =1 102 g cm1 s1 porosidade do meio (razo entre o volume dos poros e o volume total) domnio computacional do problema
ij clula da malha, em queirepresenta sua posioxe jsua posio y
13
7/26/2019 Dissertao Porosos 229
14/88
Introduo
O processo de escoamento em meios porosos de interesse de uma grande variedade de
engenheiros e cientistas, bem como de polticos e economistas que reconhecem a impor-tncia dos lencis freticos e de uma variedade de processos de extrao de petrleo. A
lei unidimensional descoberta empiricamente em 1856 por Darcy serviu como ponto de
partida para numerosas aplicaes prticas e como um desafio constante para os tericos.
Enquanto as condies originais estudadas por Darcy so encontradas em vrias situaes
prticas, suas extenses para casos mais gerais necessitam de uma anlise terica especial,
por se tratarem de situaes nas quais experimentos so difceis de se realizar [55].
Alm disso, o escoamento multicomponente em meios porosos envolve fenmenos locali-
zados, tais como frentes de concentrao, gradientes de presso e mesmo a geometria do
meio. J as frentes de concentrao podem se dispersar por meio de transporte convectivo
ou difuso na matriz slida. Em[52], vemos uma abordagem que utiliza o refinamento
adaptativo da malha (AMR) para estudar o escoamento em meios porosos em um meio
de permeabilidade (possivelmente) heterognea, levando em conta os efeitos convectivos
e difusivos existentes. Nesse artigo, mtodos de alta resoluo para as frentes de concen-
trao so utlizados, de forma a tratar os fenmenos localizados, com grande nfase nas
tcnicas de lgebra linear utilizadas.
Utilizando ideias de [52], temos como objetivo neste trabalho a utilizao da estrutura
de dados Autonomous Leaves Graphs (ALG) na simulao do escoamento de fluidos em
meios porosos. Essa estrutura foi proposta em [11], e pode ser utilizada em problemas
que envolvem refinamento adaptativo da malha. Devido eficiencia e flexibilidade do
ALG no tocante ao refinamento e desrefinamento, o resultado esperado uma melhora
significativa no esforo computacional para o clculo da conveco.
A fim de estabelecer o backgroundterico, relembraremos no Captulo 1 alguns conceitosfundamentais da mecnica dos fluidos, como a viscosidade de um fluido e a permeabilidade
de um meio, e introduziremos as equaes que regem o escoamento em meios porosos (Lei
14
7/26/2019 Dissertao Porosos 229
15/88
de Darcy e Equao da Continuidade).
O Captulo 2 voltado para os conceitos de mtodos numricos, especialmente os Mtodos
dos Elementos Finitos e dos Volumes Finitos, bem como a parte de lgebra linear esparsa,que sero necessrios para as simulaes.
Reservamos o Captulo 3 para a apresentao e discusso das formulaes mista e hbrida-
mista do sistema elptico que surge da Lei de Darcy e da Equao da Continuidade,
incluindo teoremas de existncia, unicidade e aproximao.
A discretizao da Equao da Continuidade para o transporte de massa convectivo ser
vista no Captulo 4. L, apresentamos a forma discreta dos fluxos difusivo e convectivo,
levando tambm em conta a estabilidade do esquema.
O Captulo 5 apresenta a estrutura de dados Autonomous Leaves Graph (ALG), como
proposto em[11]. Descrevemos como se d o refinamento e o desrefinamento, bem como
a ordenao das clulas da malha.
Os resultados dos testes numricos so apresentados no Captulo 6, e a concluso se d
no Captulo 7, juntamente com as perspectivas futuras para o trabalho.
15
7/26/2019 Dissertao Porosos 229
16/88
Captulo 1
Mecnica dos Fluidos
Um fluido pode ser definido, informalmente, como uma substncia que no possui forma
definida, e que assume o formato do recipiente que o contm. Do ponto de vista mecnico,
um fluido uma substncia que se deforma continuamente sob a ao de uma tenso
cisalhante (i.e., o efeito de uma fora paralela superfcie). Assim, os fluidos compreendem
as fases lquida e gasosa (ou vapor) nas quais a matria existe. A diferena principal
entre um slido e um fluido a seguinte: um slido deforma-se quando uma tenso de
cisalhamento lhe aplicada, mas no continuamente [27]. Relembraremos aqui alguns
conceitos e equaes fundamentais no estudo do escoamento.
1.1 Viscosidade
Um fluido cuja taxa de deformao varia linearmente com as tenses de cisalhamento
chamadoFluido Newtoniano. Em termos um pouco mais precisos: tomando coordenadas
cartesianas, se ij representa a componente da tenso cisalhante perpendicular a ei e na
direo ej, um fluido newtoniano incompressvel satisfaz, para i=j [9]:
ij =
vixj
+vjxi
A constante de proporcionalidade chamada de viscosidade, e assume-se que variaapenas com a temperatura do fluido.
16
7/26/2019 Dissertao Porosos 229
17/88
1.2 Hiptese do contnuo
Fisicamente falando, os fluidos so compostos por molculas em constante movimento.
Porm, na maioria das aplicaes, em especial na Engenharia, estamos interessados apenas
nas propriedades macroscpicas de muitas molculas. Estamos interessados, portanto, em
propriedades mdias: so essas propriedades que geralmente podemos medir. Tratamos,
portanto, um fluido como algo arbitrariamente divisvel, algo contnuo, e deixamos de
lado o comportamento das molculas individuais.
Essa a chamada Hiptese do Contnuo. a base da mecnica dos fluidos clssica, e
ser usada nesse texto, implicitamente, em todas as equaes relacionadas ao movimento
de fluidos. Porm, devemos estar atentos s condies nas quais a Hiptese vlida.S podemos assumi-la quando a ordem de grandeza da menor dimenso significativa do
problema superior trajetria mdia livre das mleculas (cerca de 6 x 108 m nas
condies normais de temperatura e presso (CNTP)[27]). Caso contrrio, mtodos que
estudam fluidos a nvel molecular devem ser utilizados, rea conhecida como Mecnica
Estatsitica.
1.3 Meios Porosos
Em nossas simulaes, assumiremos a existncia de um meio poroso e rgido, de tal forma
que o escoamento se d apenas por meio dos poros. Definimos, ento, a porosidade
como a razo entre o volume total dos poros e o volume total do meio.
Outro conceito importante apermeabilidadeKdo meio poroso. Essa grandeza tensorial
corresponde resistncia oferecida pelo meio passagem do fluido. No caso de escoamentounidirecional em um meio homogneo, a permeabilidade uma grandeza escalar e podemos
escrever
u=K
dP
dx,
em queu a velocidade do fluido e Pa presso do fluido. Essa equao conhecida como
Lei de Darcy, e nos diz que afora motrizdo escoamento em meios porosos o gradientede presso.
17
7/26/2019 Dissertao Porosos 229
18/88
Quando se trata de meios porosos, normalmente se utiliza a chamada velocidade de Darcy
no lugar da velocidade usual. A relao entre as duas velocidades dada pela equao:
v= u, (1.1)
ondev a velocidade de Darcy e u a velocidade de fluido. De agora em diante, a menos
que dito expressamente o contrrio, trabalharemos apenas com a velocidade de Darcy.
1.4 Lei de Darcy
No caso mais geral, K uma matriz simtrica positivo-definida, chamada matriz de
permeabilidades, e a equao que rege o escoamento pode ser escrita como:
K1v dV =
pn ds. (1.2)
importante ressaltar que essa equao s vlida para o escoamento laminar; para
velocidades muito altas do fluido, devem ser usadas as equaes de Navier-Stokes. De fato,
a Lei de Darcy derivada das equaes de Navier-Stokes sob as hipteses de escoamento
laminar, incompressvel, permanente e de uma relao linear entre as foras viscosas e a
velocidade do fluido. Os interessados encontraro a deduo da Lei de Darcy a partir de
Navier-Stokes em [55].
Adicionalmente, convm notar que a expresso a permeabilidade em um determinado
ponto no possui sentido fsico: K s tem sentido em regies mensurveis do domnioque sejam grandes comparadas com o volume dos poros, e possivelmente menores que
o domnio total do escoamento (de maneira anloga Hiptese do Contnuo). Porm,
ignorando essas observaes, comum usar o teorema da divergncia e escrever a equao
(1.2) na forma diferencial:
v=K
p.
18
7/26/2019 Dissertao Porosos 229
19/88
1.5 Equao da Continuidade
Como estamos trabalhando com um modelo do mundo fsico, a massa deve ser conservada.
Essa conservao d origem a uma equao conhecida como Equao da Continuidade.
Podemos escrev-la como:
d
dt
cdV +
F n ds=
cw dV (1.3)
Na equao acima, representa a densidade do componente em questo, ca sua concen-
trao volumtrica, F o fluxo de massa pela fronteira do domnio , n o vetor normalexterior unitrio e w o termo de fonte (ou sorvedouro) de massa.
Consideraremos o escoamento incompressvel de uma mistura de gua e um soluto mis-
cvel, que chamaremos tracer. Esses componentes sero indicados pelos subndices w e
t, respectivamente. Dessa forma, as equaes da continuidade para os componentes em
uma regio so:
ddt
tct dV +
Ft n ds =
tctwt dV
d
dt
wcw dV +
Fw n ds =
wcwwwdV
(1.4)
Devido incompressibilidade do escoamento,t e w so considerados constantes.
1.6 Fluxo Mssico
Os fluxos de massa envolvem, principalmente, dois efeitos: o transporte convectivo pelo
fluido e a difuso. Porm, em muitas aplicaes, a contribuio difusiva desprezvel com
relao convectiva. Desse modo, desprezando-se a difuso, os fluxos de massa F podem
ser escritos como:
Os fluxos de massa Ft envolvem trs efeitos: o transporte convectivo pelo fluido, a dis-
perso devida s conexes aleatrias entre poros e a difuso molecular. Se n o vetornormal unitrio a, assumiremos que os vetores de fluxos de massaFte Fwsatisfazem:
19
7/26/2019 Dissertao Porosos 229
20/88
Ftnds=
ctv n v n
v
(l t) v ct v t+
n ctds
Fwnds=
cwv n v nv (l t) v cwv w+
n cwds (1.5)
Aqui, l e t so constantes fsicas, a saber: os coeficientes de disperso longitudinal e
transversal, parte do modelo de disperso Fickiana. A difuso molecular representada
pela difusividade e tortuosidade .
1.7 Conservao do Volume
Por definio, as soma das concentraes dos componentes um:
ct+ cw = 1. (1.6)
Podemos substituir os fluxos de massa em (1.5) nas equaes (1.4), dividir pelas respec-
tivas densidades (t e w) e somar sobre os dois componentes para obter:
v n ds=
w dV (1.7)
Uma vez que vn representa a componente normal do fluxo de fluido, essa equao
representa a conservao do volume.
Comumente, fazemos o volume de tender a zero e, utilizando o teorema da divergncia,
obtemos:
v= w (1.8)
Longe de poos (fontes ou sorvedouros de fluido), essa equao diz que a velocidade deDarcy v possui divergente nulo. Em 1D, isso significa que o campo de velocidades
constante entre poos.
20
7/26/2019 Dissertao Porosos 229
21/88
Captulo 2
Mtodos Numricos para a resoluo
de EDPs
Vimos, no Captulo1, as equaes que descrevem o comportamento de um fluido num
meio poroso. Agora precisamos de ferramentas para realizar a simulao numrica. Temos
dois problemas aqui envolvidos: o primeiro deles envolve a determinao do campo de
velocidades, pela Lei de Darcy (Eq. 1.2) e pela continuidade do volume (Eq. 1.7).Veremos mais adiante (Captulo3) que essa uma EDP elptica. Uma vez estabelecido o
escoamento do fluido, o segundo problema computar o transporte convectivo do tracer,
pela Equao da Continuidade (Eq. 1.4). Problemas de conveco so bastante conhecidos
por seu carter hiperblico.
Dessa forma, precisamos de mtodos tanto para equaes elpticas quanto para hiper-
blicas. Seguindo o procedimento descrito em [52], utilizaremos o mtodo de Elementos
Finitos (que baseado no mtodo de Galerkin) para o problema elptico, e Volumes
Finitospara o problema hiperblico. Vejamos como funciona cada um desses mtodos.
2.1 O Mtodo de Galerkin
Dada uma EDP na sua formulao variacional, introduziremos um mtodo para encontrar
uma soluo aproximada, que posteriormente ser utilizado no mtodo dos Elementos
Finitos.
21
7/26/2019 Dissertao Porosos 229
22/88
Dado um espao de Hilbert V, uma forma bilinear B(, ) : V V R(procedente daformulao variacional da EDP) elV (lrepresentando o termo independente da EDPe V o dual de V), queremos encontrar uVtal que:
B(u, v) =l(v) vV . (2.1)
Assumiremos tambm que a forma B(, ) limitada e coerciva. Ou seja, existem cons-tantes ke tais que
|B(u, v)| kuVvV vV (2.2)
e
B(v, v)v2V vV. (2.3)
A ideia do mtodo de Galerkin muito simples. A dificuldade de resolver (2.1) vem do
fato que V um espao muito grande (em geral, tem dimenso infinita!), impedindo a
criao de um procedimento para calcular a soluo exata. O mtodo de Galerkin se
baseia, ento, em uma famlia enumervel de subespaos de dimenso finita{Vn}n=1V,
tais quen=1
Vn=V (2.4)
e que satisfazem Vn Vn+1 e dim(Vn) = Nn
7/26/2019 Dissertao Porosos 229
23/88
N.
Demonstrao. A restrio da forma B(,
) a Vn
Vn elptica e coerciva. Da mesma
forma, a restrio do funcional l a Vn continua sendo linear e limitada, portanto l Vn.Assim, as hipteses do teorema de Lax-Milgram so satisfeitas, mostrando que a soluo
de (2.5) existe e nica.
A essncia do mtodo de Galerkin se deve ao fato de que a soluo un Vn pode serescrita como uma combinao linear finita de elementos da base deVn. Portanto, existem
escalares aj e, para qualquer vVn, escalaresbique satisfazem
un =Nnj=1
ajj, v=Nni=1
bii. (2.6)
Substituindo (2.6) em (2.5), e usando a linearidade de B(, ), obtemos
Nn
i=1
biNn
j=1
B(j, i)aj
l(i)= 0, (2.7)
ou, de forma mais clara:
Nni=1
bi
Nnj=1
Kijaj Fi= 0, (2.8)
em que
Kij :=B(j, i) e Fi:=l(i) (2.9)
so, respectivamente, uma matriz Nn Nn(chamada dematriz de rigidez) e um vetor dedimensoNn(chamado devetor de carga). Note que podem-se calcular explicitamenteKijeFi, posto que as funesiso conhecidas e tantoB quantol so dados explicitamente.
Uma vez que os coeficientes bi so arbitrrios, segue que o termo entre parnteses da
equao (2.8) deve ser zero, o que reduz a equao a um sistema linear:
23
7/26/2019 Dissertao Porosos 229
24/88
Ka= F. (2.10)
Uma vez resolvido esse sistema de equaes, pode-se escrever a soluo aproximada unpela equao (2.6).
Notemos, por curiosidade, uma propriedade inerente aos un. Caso B seja uma forma
simtrica (e coerciva, por hiptese), podemos definir um produto interno (, )B como(u, v)B =B (u, v). Esse produto interno induz uma norma B, chamada de norma daenergia. Subtraindo as equaes (2.1) e (2.5), obtemos
B(u un, v) = 0 vVn. (2.11)
Assim, no produto interno definido por B, uun ortogonal aVn, o que nos leva a concluirqueun uma espcie de projeo ortogonal deuemVn. De fato, temos o seguinte teorema:
Teorema 1.
u unB = infvVn
u vB.
Em outras palavras, na norma da energia, un a melhor aproximao deu emVn.
Demonstrao. Temos:
B(u
v, u
v) =B(u
un+ un
v, u
un+ un
v)
=B(u un, u un) + 2B(u un, un v) + B(un v, un v)
Como B(u un, un v) = 0, por (2.11), temos
u v2B =u un2B+ un v2B.
Fixando ue un, vemos queu vB assume um mnimo (global) quando un = v, isto
24
7/26/2019 Dissertao Porosos 229
25/88
u unB = infvVn
u vB
como queramos demonstrar.
Procuraremos, agora, cotas para uun usando a norma de V. Isso nos dar respaldopara concluir sobre a convergncia da sequncia{un}.
2.1.1 Propriedades da aproximao de Galerkin
Provaremos nesta seo os teoremas que faltam para estabelecer a convergncia de{un}para a soluo exata u. Comearemos pelo:
Teorema 2 (Lema de Ca). SeB satisfaz (2.2) e (2.3), ento
u
un
V
k
infvVn
u
v
V
Demonstrao. Usando (2.11), obtemos:
B(u un, u un) =B(u un, u v) + B(u un, un v)=B(u un, u v).
Ora, temos
B(u un, u un)u un2V
e
B(u un, u un)ku unVu vV vVn
donde
25
7/26/2019 Dissertao Porosos 229
26/88
u unV k
u vV vVn
o que completa a prova.
Estamos, finalmente, prontos para provar a convergncia das solues un:
Teorema 3. SejaVum espao de Hilbert eV1V2V3 uma famlia enumervel deespaos de dimenso finita que cumpra
n=1
Vn=V . (2.12)
Adicionalmente, suponhamos queB(, ) :V V Rseja uma forma bilinear limitada ecoerciva e que lV. Ento
limn
u unV = 0 (2.13)
Ou seja, o mtodo de Galerkin para o problema (2.1) convergente.
Demonstrao. Dada uma soluoude (2.1), por (2.12) possvel encontrar uma sequn-
cia tal que{vn} Vn para todo n Ne
limn
u vnV = 0
O resultado segue pelo Lema de Ca.
2.2 O Mtodo dos Elementos Finitos (FEM)
Uma vez compreendida a essncia do mtodo de aproximao de Galerkin, precisamos
encontrar os espaos Vn apropriados para que se ponha em prtica a simulao. Um
26
7/26/2019 Dissertao Porosos 229
27/88
procedimento simples e sistemtico para tal o Mtodo dos Elementos Finitos. Este
consiste em dividir o domnioem uma coleo de polgonosk(chamados deelementos).
A coleo desses elementos chamada de malha. Definimos, ento, a base do espao Vn
usando funes suaves (em geral, polinmios) em cada k. Esses elementos devem cobrirtodo o espao , e a interseo de dois elementos deve ser um vrtice, uma aresta ou o
conjunto vazio. Alm disso, a base de Vndefinida dessa maneira deve pertencer ao espao
V (pois VnV).
A ltima afirmativa causa certa dificuldade em aplicar o FEM: em geral, as funes de
Vn devem ser contnuas em . Como um vrtice pode pertencer a mais de um elemento,
seVnsatisfaz(x0)= 0para algumx0i j, ento as restries deao interiorde
i e de
jno podem ser identicamente nulas. Dessa forma, faz sentido definir as
funes-base por meio dos seus valores nos vrtices, com cada funo-base possuindo o
valor 1 em um vrtice e 0 em todos os outros. Tambm podem ser usados pontos das
arestas, faces ou clulas [18]. Uma tpica funo base linear por partes em uma malha
triangular mostrada na Fig. 2.1.
Figura 2.1: Funo base.Fonte: GIACCHINI (2012, p. 13)[28]
A restrio quanto continuidade das funes-base acaba por introduzir um acoplamento
entre as clulas. Retomando o mtodo de Galerkin, precisamos calcular os valores de
Kij =B(j, i).
Em geral, a forma bilinear B dada por uma integral envolvendoi, j e suas derivadas.
Assim, se i e j no so identicamente nulas em uma determinada clula, o termo Kijno ser necessariamente igual a zero. De modo a evitar essa dificuldade, gostaramos de
trabalhar em espaos de funes possivelmente descontnuas, que nos permitam zerar o
mximo de elementos da matriz K. Um dos trunfos da FormulaoMista-Hbridaparaa soluo da Lei de Darcy, que ser apresentada no Captulo 3, ser a possibilidade de
trabalhar em espaos descontnuos (no sendo, porm, o nico mtodo com essa carac-
27
7/26/2019 Dissertao Porosos 229
28/88
terstica: h tambm os mtodos de Galerkin descontnuos [30]). Mtodos de Elementos
Finitos que trabalham com funes-base descontnuas so chamados de no-conformes.
2.3 O Mtodo dos Volumes Finitos (FVM)
O mtodo dos Volumes Finitos um mtodo de discretizao utilizado para a simulao
numrica de leis de conservao (sejam elas elpticas, hiperblicas ou parablicas). Tem
sido usado extensivamente em vrios campos da engenharia, como a mecnica dos fluidos,
a transferncia de calor e massa e a engenharia de petrleo[25]. Alguns aspectos impor-
tantes do mtodo dos Volumes Finitos so similares aos dos Elementos Finitos [40]: podeser usado em geometrias arbitrrias e gera esquemas robustos. Um aspecto adicional a
conservao local dos fluxos numricos, i.e, o fluxo numrico conservado de uma clula
da discretizao para outra. Essa ltima caracterstica torna o FVM bastante atraente
na modelagem de problemas para os quais o fluxo importante, como na transferncia de
calor, massa e na mecnica dos fluidos. O FVM localmente conservativo por ser baseado
numa abordagem debalanos: um balano local escrito em cada clula da discretizao,
normalmente chamada de volume de controle; pelo teorema da divergncia, uma formu-
lao integral dos fluxos na fronteira do volume de controle ento obtida. Para uma
abordagem completa, veja por exemplo [37], [51],[49], [25] ou [26].
2.4 Mtodos numricos para a resoluo de Sistemas
Lineares
Vimos que, em nossas simulaes, precisaremos resolver sistemas lineares. Veremos tam-bm que esses sistemas sero de grande porte e esparsos. Para tal, mtodos convencionais
de resoluo de sistemas densos no so efetivos, pois o grande porte do sistema causa
uma vasta utilizao de memria. Alm disso, o extenso nmero de clculos feitos implica
na necessidade de uma grande capacidade de processamento. Surgem, ento, os mtodos
iterativos, adaptados aos sistemas esparsos e de grande porte.
2.4.1 O Mtodo do Gradiente Conjugado
O mtodo do Gradiente Conjugado (CG) efetivo na resoluo de sistemas da forma
28
7/26/2019 Dissertao Porosos 229
29/88
7/26/2019 Dissertao Porosos 229
30/88
O Mtodo da Mxima Descida
Vimos que a soluo do sistema (2.14) pode ser visto como um problema de minimizao de
uma forma bilinear. Para resov-lo, comeamos com uma estimativa inicial x0 e criamosuma sequncia{xn} que esperamos que convirja para a soluo. Para tal, em cadaiterao, escolhemos a direo na qual fdecresce mais rapidamente. Essa direo a
oposta af(xn). De acordo com a equao (2.15), a direo desejada f(xn) =bAxn.Se denotarmos oresduona ensima iterao por rn = bAxn, temos que o vetor buscadopara a prxima iterao tem a forma xn+1=xn+ rn.
Falta-nos determinar . Isso feito de forma a minimizar o valor de f(xn+1) = f(xn+
nrn). Um clculo simples nos mostra que
n= rtnrnrtnArn
. (2.17)
Para o mtodo da mxima descida, pode-se mostrar (vide [48]) que, se x a soluo de
(2.14), vale:
f(xn) f(x)f(x0) f(x)
k 1k+ 1
2n(2.18)
em que k=A A1 o nmero de condio da matriz A.
Vemos que o mtodo da mxima descida convergente, embora para nmeros de condio
altos essa convergncia possa ser lenta. Gostaramos, portanto, de alterar um pouco o
mtodo de forma a melhorar a convergncia.
Trabalharemos comdirees de buscapn, de forma que xn+1= xn+ npn. No mtodo da
mxima descida, tnhamos pn =rn. Generalizaremos esse mtodo, construindo os pn de
tal forma que xn+1seja o mnimo de frestrita ao espao gerado pelas direes de busca
j tomadas{p1, p2, , pn}. Assim, caso trabalhemos com aritmtica de preciso infinita(ou seja: apenas teoricamente), aps no mximon iteraes, atinge-se o mnimo de fem
todo o espao.
Esse o Mtodo do Gradiente Conjugado, proposto originalmente por Hestenes e Stiefel
30
7/26/2019 Dissertao Porosos 229
31/88
[33]. So construdas sequncias de vetores
rn+1=rn nApnpn+1=rn+ kpn
que satisfazem as condies de ortogonalidade e conjugao
rtirj = 0
ptiApj = 0
rtipj = 0
para todo i=j .
Os escalaresk e k so dados por
k = rtk1rk1
ptkApk(2.19)
k=rtk1rk1rtk2rk2
. (2.20)
Assim, como proposto por Reid [42], o algoritmo fica na forma
1: x0estimativa inicial2: tolerncia
3: r0=b Ax04: k= 1
5: p1 = r0
6: 1= (rt0r0) \ (pt1Ap1)
7: x1=x0+ 1r1
8: r1=r0 1Ap19: while
rk
> do
10: k= k + 1
11: k= (rtk1rk1)\(rtk2rk2)
31
7/26/2019 Dissertao Porosos 229
32/88
12: pk =rk1+ kpk1
13: k = (rtk1rk1)\(ptkApk)
14: xk =xk1+ kpk
15: rk=rk1 kApk16: end while
17: x= xk
Notemos que esse mtodo requer apenas uma multiplicao matriz-vetor por iterao.
Note tambm que somente 4 vetores so armazenados: x,r,pe Ap.
Analisemos, agora, a convergncia. Caso no haja erros de arredondamento, o mtodo se
encerra em exatamenteniteraes, sendona ordem da matrizA. Porm, em simulaes,esse nunca o caso. DefinindoxA=
xtAx, pode-se mostrar (vide [48]) que
xn xA2
k 1k+ 1
nx0 xA . (2.21)
Assim, o Gradiente Conjugado mais eficiente que o da mxima descida. Porm, em
qualquer um dos dois mtodos, a convergncia depende do nmero de condio da matriz
A.
2.4.2 O Mtodo do Gradiente Biconjugado Estabilizado
Uma generalizao do Mtodo do Gradiente Conjugado aplicvel a matrizes no-simtricas
foi proposta em 1992 por H. A. Van der Vorst [54]. chamado de Gradiente Biconjugado
Estabilizado (normalmente abreviado como BiCGSTAB), e pertence classe de mtodos
numricos conhecidos como Metdos de Krylov.
Dado um precondicionador K=K1K2A, o mtodo do BiCGSTAB descrito abaixo:1: x0estimativa inicial
2: > 0 tolerncia
3: r0=b Ax04: r0=r0
5: 0 = 0=0 = 1
6: p0 = v0 = 0
7: k= 1
32
7/26/2019 Dissertao Porosos 229
33/88
8: whilerk> do9: k= k + 1
10: k = rt0rk1
11: k= (k/k1)(k1/k1)12: pk =rk1+ k(pk1 k1vk1)13: Resolver Ky=pk14: vk =Ay
15: k =k/(rt0vi)
16: sk =rk1 kvk17: Resolver Kzk =sk18: tk=Azk
19: k
= (K11
tk
)t(K11
sk
)/ K11
tk2
20: xk =xk1+ kyk+ kzk
21: rk=sk ktk22: end while
23: x= xk
Como esperado, esse mtodo tem um maior custo computacional que o CG.
2.4.3 A fatorao LU incompleta
A fatorao LU incompleta (abreviada como iLU) uma alternativa dentre os vrios
possveis precondicionadores de um sistema. Consiste em realizar uma eliminao gaus-
ssiana aproximada na matriz, de forma que apenas a diagonal e termos significativos
mantenham-se no-nulos. Dessa forma, se a matriz esparsa, podemos tomar sua fa-
torao iLU (composta de uma matriz triangular inferior L e uma triangular superior
U) composta de matrizes esparsas, mantendo assim a eficincia na resoluo de sistemasesparsos e de grande porte.
Para um estudo mais completo sobre a escolha de precondicionadores ou sobre fatoraes
de matrizes aplicadas soluo de equaes, pode-se consultar, por exemplo, [7] e [50].
33
7/26/2019 Dissertao Porosos 229
34/88
Captulo 3
Discretizao da Lei de Darcy
Este captulo dedicado ao estudo do sistema associado Lei de Darcy e equao da
conservao do volume. Relembrando o captulo 1, essas equaes podem ser descritas,
na forma forte, como
v=
K
p
v= 0p= pD em D
v n= vN n em N
(3.1)
Aqui, D a parte da fronteira onde especificada a presso pD, e N a parte onde
especificada a componente normal da velocidade vN n. Essas equaes so, na verdade,um caso particular de equaes da forma:
T1v + p= g v= wp= pD(x) em D
v n= vN n em N
(3.2)
na qual T:= K simtrica e positivo-definida.
34
7/26/2019 Dissertao Porosos 229
35/88
Fisicamente, podemos atribuir o termo g ao da gravidade, por exemplo, e o termo
w a uma fonte ou sorvedouro de fluido. Para um tratamento completo, mostraremos os
resultados vlidos no caso geral (3.2). Veremos agora, maneiras de tratar esse problema,
que culminaro no mtodo hbrido-misto, o qual ser utilizado nas simulaes.
3.1 Formulaes do problema elptico
primeira vista, parece interessante substituir a primeira das equaes (3.2) na segunda,
obtendo uma equao somente emp:
(Tp) =w (Tg)
p= pD em D
(Tp) n= (vN Tg) n em N(3.3)
ou simplesmente:
(Tp) =fp= pD em D
(Tp) n= vN n em N(3.4)
Essa formulao a famosa forma divergente de um problema elptico, e uma variao
da conhecida Equao de Laplace (que corresponde ao caso T = I). Sabemos que esse
problema tem soluo nica (supondo a regularidade necessria, vide [29], [24]), e quepode ser visto como um problema de minimizao de um certo funcional. Porm, estamos
interessados no campo develocidades, obtido por meio do gradiente de p. Computacional-
mente, essa abordagem pouco interessante, pois derivao numrica implica em perda
de preciso. Dessa forma, precisamos de uma maneira de obter o campo de velocidades
diretamente. Surge, ento, o mtodo mistopara a soluo do problema.
35
7/26/2019 Dissertao Porosos 229
36/88
3.2 A formulao mista
Podemos tentar resolver diretamente o sistema (3.2). Porm, o tratamento das equaes
pode ser mais fcil em espaos de funes que requerem menor regularidade. Faz sentido,
ento, tentar reduzir a exigncia de suavidade das funes e buscar solues num espao
mais geral. Essa a essncia da formulao fracadas EDPs.
Consideraremos, por ora, vN n 0. Multipliquemos a primeira das equaes de (3.2)por uma funo arbitrria uUe a segunda por q Q(em que U e Qsero definidosposteriormente), e integremos em. Ento, se (u, p) satisfazem (3.2), tambm satisfazem
u T1v +
u p=
u g uU
q v=
q w qQ
p= pD em D
v n= 0 em N
(3.5)
Se tomarmos U =H0,N(div, ) :={u
L2() :
u
L2()e un|N= 0
}, podemos
aplicar o teorema da divergncia funo F= uppara obter
u T1v
( u)p=
u g
(u n)pD uH0,N(div, )
q v=
q w qL2()
pL2()v
H0,N(div, )
(3.6)
.
Note que agora podemos buscar uma soluo(u, p)H0,N(div, ) L2(). E mais: sepfor suficientemente suave, pelo teorema de Hahn-Banach, o problema (3.6) equivalente
ao problema (3.2) em H0,N(div, ) L2()com condies de Neumann homogneas emN.
A vantagem desse mtodo obter o campo de velocidades vdiretamente, resultando numa
maior preciso numrica. Resta, portanto, mostrar a existncia e unicidade da soluodo problema (3.6), obtendo com eles os teoremas de convergncia usados para estimar o
erro da simulao.
36
7/26/2019 Dissertao Porosos 229
37/88
Antes disso, porm, retiremos a homogeneidade da condio de Neumann. Seja v qualquer
funo que satisfaa v n= vNn. Isso pode ser feito, por exemplo, tomando uma soluoclssica (veja[29] ou [24]) de
(Tp) = 0
p= 0 em D
(Tp) n= vN n em N
(3.7)
Tomando a soluo (v0, p)do problema
u T1v +
( u)p=
u g
(u n)pD
u T1v uH0,N(div, )
q v=
q w
q v qH0D()
pL2()vH0,N(div, )
(3.8)
e definindo v = v+ v0, vemos que (v, p) a soluo procurada. Isso nos diz que, em
N, condies de contorno no homogneas significam, simplesmente, uma mudana no
lado direito das equaes (3.5). Vejamos agora um tratamento mais formal da formulao
mista.
3.2.1 Existncia, unicidade e convergncia
Trabalharemos, nesta seo, com um problema mais abstrato. Seja V um espao de
Hilbert com a normaV, produto interno(, )Ve dual V. Consideraremos uma formalinear limitada emV V
|a(u, v)| a uVvV . (3.9)
Essa forma define um operador contnuo A : VV
37
7/26/2019 Dissertao Porosos 229
38/88
(Au,v)VV =a(u, v)u, vV. (3.10)
Seja outro espao de Hilbert Qcom normaQ, produto interno (, )Q e dual Q. Sejatambm uma forma bilinear contnua em V Q
|b(u, q)| b uVqQ . (3.11)
Novamente, introduzimos um operador linearB : V Q e o operador transposto Bt :Q
V, definidos por
(Bu,q)QQ= (u, Btq)VV =b(u, q)uV, qQ. (3.12)
Sejam f V, gQ. Queremos encontrar uV, pQ tais que
a(v, u) + b(u, p) =f(u) uV
b(v, q) =g(q) qQ(3.13)
ou
Av+ Btp= f em V
Bv = g em Q(3.14)
Por exemplo, para o problema (3.6), temos as seguintes formas bilineares:
a:H0,N(div, ) H0,N(div, ) R(v, u)
T1v u (3.15)
b: H0,N(div, ) L2()R
(u, p) ( u)p (3.16)
.
38
7/26/2019 Dissertao Porosos 229
39/88
Supondo que o operador T1 :H0,N(div, )H0,N(div, )seja contnuo, temos
|a(u, v)| T1 uH0,N(div,) vH0,N(div,) (3.17)e
|b(u, p)| uH0,N(div,) pL2() . (3.18)Lema 2. Se g ImB e a forma bilineara(, ) coerciva em KerB, i.e., existe0 talque
a(u0, u0)
0
u0
2V
u0
KerB, (3.19)
ento existe um nico uV soluo de
a(u, v0) =f(v0)v0KerB (3.20)
que satisfaa
Bu = g emV. (3.21)
Demonstrao. Tomemos ug V tal que Bug = g. Pela coercividade dea em Ker B,podemos encontraru0Ker B satisfazendo
a(u0, v0) =f(v0) a(ug, v0)v0Ker B. (3.22).
Escrevendo u = u0+ ug, conclumos a existncia. Tomemos, ento, u1 e u2 solues de
(3.20) e (3.21). Temos
a(u1 u2, v0)v0Ker B. (3.23)
Como u1 u2Ker B, a coercividade de aem Ker B implica u1=u2.
claro que, se (3.13) tem uma soluo (u, p), ento u ser soluo de (3.20) e (3.21).
O lema acima implica, ento, que o primeiro componente da soluo (u, p)de (3.13) (se
39
7/26/2019 Dissertao Porosos 229
40/88
existir) nico. Notemos tambm que, pela prova do lema:
u ug + 10
(fV+ a ug) . (3.24)
Assim, se encontrarmos uma cota superior paraug, teremos feito o mesmo parau.Notemos tambm que a coercividade deapode ocorrer apenas em Ker B , sem que ocorra
em todo V.
Temos agora o problema de encontrarp. Para tal, mais hipteses tm de ser feitas sobre
o operador B. Mais precisamente, precisaremos que Im B seja fechada em Q. A razo
para tal est no seguinte resultado da anlise funcional, cuja prova pode ser vista em [56].
Lema 3. As afirmaes abaixo so equivalentes:
ImB fechada emQ.
ImBt fechada emV.
(KerB) = ImBt.
(KerBt) = ImB. Existek0>0 tal que, para todo gImB, existevgV comBvg =g e
vgV1/k0 gQ.
Existek0>0 tal que, para todo f ImBt, existeqf Q comBtqf=f eqfQ1/k0 fV.
Se qualquer uma das condies acima acontecer, dizemos queB admite um levantamentocontnuo deV aQ.
Podemos, finalmente, provar a existncia da soluo de (3.13).
Teorema 4. Sejag ImB e sejau a soluo do problema (3.20) e (3.21). Se ImB fechada emQ, existepQ tal que(u, p) soluo de (3.13). Alm disso, as seguintescotas so vlidas:
uV 1
0fV+
1
k0
1 +
a0
gQ, (3.25)
40
7/26/2019 Dissertao Porosos 229
41/88
pQ\KerBt 1
k0
1 +
a0
fV+
ak20
1 +
a0
gQ. (3.26)
Demonstrao. Para obter (3.25), podemos tomar ug satisfazendo (3.20), (3.21) e
ug 1k0
gQ,
portanto (3.24) nos d (3.25).
Considere o funcional em V:
L(v) =f(v) a(u, v).
Temos ento que L(v0) = 0 v0Ker B . Logo, pelo lema (3),LIm B t, mostrando queexistepQ que satisfaz
L(v) =b(v, p)vV,ou seja:
a(u, v) + b(v, p) =f(v)vV.
Por meio de projees ortogonais, podemos tomarp(KerB t). Sep0Ker B t, o lema(3) ainda nos diz
supvVv=0
b(v, p+p0)v = supvV
v=0
b(v, p)v = supvV
v=0
L(v)v infp0 Ker B t k0 p +p0Q= k0 pQ/ Ker B t.
Das desigualdades acima, e do fato que b(v, p) =L(v) =f(v) a(u, v)vV, obtemosa desigualdade (3.26).
Note que no conclumos a unicidade dep. De fato, pela primeira equao de (3.13), vemos
quepest definida a menos de um elemento de Ker B t. A unicidade estar garantida se,e somente se, Ker Bt = 0, ou seja, se B for sobrejetiva.
41
7/26/2019 Dissertao Porosos 229
42/88
Resta-nos verificar se as hipteses do teorema so satisfeitas no caso da Lei de Darcy.
Basta verificar que o operador B igual a -div restrito a H0,N(div, )e que
div:H0,N(div, )L2()
sobrejetivo se D=. Para tal, dado f L2(), basta utilizar a existncia da soluo([29], [24]) de
= f|D = 0n
|N= 0(3.27)
e tomar =, obtendo H0,N(div, )com div =f. No caso em que D =, asoluo estar definida a menos de uma constante arbitrria.
Vejamos como utilizar o Mtodo dos Elementos Finitos para construir uma aproxima-
o para o problema. Trabalharemos num domnio quadrado , dividido em quadrados
abertos disjuntos ij tais que i,j
ij = . (3.28)
Denotemos, num elemento ij da partio,
P0(ij) : o espao de funes constantes em ij
e
RT0(ij) :{(ax + b,cy+ d) :a, b, c, dP0} .
Denotemos por ho supremo dos dimetros (i.e, dos dimetros dos crculos circunscritos)
desses quadrados e
Vh:=
vhH0,N(div, ) :vh|ijRT0(ij)i, j
,
42
7/26/2019 Dissertao Porosos 229
43/88
Qh:=phL2(div, ) :ph|ijP0(ij)i, j
.
Note que div Vh=Qh. Procuramos ento (vh, ph)satisfazendo
a(u, vh) + b(u, ph) =f(u) uVhb(vh, q) =g(q) qQh
(3.29)
ou mesmo
Ahu + Bthp= f em V
h
Bhu= g em Qh.(3.30)
Aqui, o operadorBh a restrio de Ba Vh. ComoB =div, pela igualdade divVh=Qhvemos que Bh sobrejetiva e Ker Bh Ker B. Alm disso, como a(, ) coerciva emKer B, tambm o ser em Ker Bh. Assim, o problema discreto 3.29tem soluo nica.
Pode ser visto, tambm, em [10]:
Teorema 5. Sejam(v, p) soluo de (3.13) e(vh, ph) soluo de (3.29). Ento, existeK
independende deh tal que as seguintes estimativas so vlidas:
v vhH0,N(div,)+ p phL2()K h (3.31)
Isso significa que temos uma aproximao de primeira ordemao utilizar o mtodo misto.
Vimos que o mtodo misto uma alternativa interessante por permitir a computao
direta do campo de velocidades. Porm, pode-se mostrar ([10]) que a formulao mista
envolve a determinao de umponto de sela. Computacionalmente, problemas de ponto de
sela so mais difceis que problemas de minimizao. Alm disso, o espaoVh composto
de funes contnuas, o que gera um certo acoplamento no sistema associado, de forma
que a matriz de rigidez tenha menos zeros. Procuramos (se possvel) um mtodo que nos
leve a um problema de minimizao, preferencialmente que envolva funes descontnuas.Esse mtodo o chamado mtodo misto-hbrido, tema da prxima seo.
43
7/26/2019 Dissertao Porosos 229
44/88
3.3 A formulao mista-hbrida
A ideia desta seo trabalhar em um espao menos restritivo que H0,N(div, ). Come-
cemos definindo
X() ={v: v ijH(div, ij)} .
claro que H0,N(div, ) X(). Para trabalhar nesse espao, precisamos saber sobquais condies uma funo em X()pertence a H0,N(div, ).
Seja Eo conjunto das arestas dos ij e 0,D as funes em L2(E)que se anulam em D
(a parte da fronteira com condies de Dirichlet). Definimos a seguinte forma bilinear:
c: X() 0,D R(v, )
i,j
ij
v n ds (3.32)
Como usual, denotaremos por C :X()0,D o operador linear definido por c(, ). imediato ver que
Cu=i,j
v|ij nij
Esse operador sobrejetivo. De fato, dado 0,Dpodemos resolver em cada ij:
ij = 0ij = 0 em ij Dijn
=|ij em ij N
e escrevervcomo uma soma (direta) dos ij .
Em[10], temos as seguintes caracterizaes:
44
7/26/2019 Dissertao Porosos 229
45/88
Lema 4. Seja vX(). Ento
(c(v, ) = 00,D)(vH0,N(div, )). (3.33)
Lema 5. Se0,D ec(v, ) = 0vX() (3.34)
ento 0.
Podemos ver, ento, a condio v H0,N(div, ) como uma restrio em X() dadapor c(v, ) = 0
0,D.
Vejamos como uma nova formulao do problema pode ser apresentada. Escreveremos,
emX0,N(),v|ij =vij e, em L2(),p|ij =pij. Seja(v, p)soluo de (3.13) e considereo mapa
: uf(u) a(v, u) ij
b(uij, pij). (3.35)
claro que (u) = 0para todo uH0,N(div, ) =Ker C. Pelo lema (3), isso implica naexistncia de um 0tal que
c(u, 0) =(u)uX(). (3.36)
Ou seja, temos
a(v, u) +ij
b(uij, pij) + c(u, 0) =f(u). (3.37)
Para ver que 0 nico, basta usar o lema (5).
Podemos resumir os resultados obtidos no seguinte teorema:
Teorema 6. [Formulao mista-hbrida] Sejam (v, p) soluo de (3.13) e 0 definido
por (3.36). Ento o terno ordenado (v, p , 0) a unica soluo do seguinte problema:
45
7/26/2019 Dissertao Porosos 229
46/88
encontrar(v, p , 0)X() L2() 0,D tais que
a(v, u) +
i,jb(uij, pij) + c(u, 0) =f(u) uH()
i,jb(vij, qij) =g(q) qL2()
c(v, ) = 0 0,D.(3.38)
Denotaremos por Vh =
vX() :v|ij RT0(ij) i, j
,Ph =pL2() :p|ij P0()i, j
eh ={0,N; constante em cada aresta de E}. Esses espaos so normados (e deHilbert): Vhe Ph herdam a norma natural de X()e L2(). Alm disso:
2h =eE
|e| |e2L2(E) .
Podemos, ento, estender nosso resultado para o problema discreto:
Teorema 7. [Formulao mista-hbrida discreta] Existe um nico terno (vh, ph, h)Vh Ph h soluo de
a(v, u) +ijb(uij, pij) + c(u, ) =f(u) uVhijb(vij, qij) =g(q) qPh
c(v, ) = 0 h.(3.39)
Finalmente, em [10], encontramos o seguinte teorema para a convergncia do Mtodo dos
Elementos Finitos:
Teorema 8. [Convergncia dos multiplicadores de Lagrange] Seja (v, p , ) soluo de
(3.38), (vh, ph, h) soluo de (3.39) e a projeo deL2
() emh. Ento existe umaconstanteK independente deh tal que
p hh K h2 (3.40)
Assim, vemos que h converge para a projeo de p nas arestas das clulas. E mais:
h uma aproximao de segunda ordem de p, enquantop phPh somente O(h)(vide Equao (3.31)). Assim, os multiplicadores de Lagrange sero identificados com as
presses nas arestas.
As prximas sees so devotadas a estabelecer a implementao numrica da formulao
46
7/26/2019 Dissertao Porosos 229
47/88
mista-hbrida no escoamento de fluidos em meios porosos.
3.4 A formulao mista-hbrida e o escoamento em
meios porosos
Utilizaremos a abordagem descrita em [52] para a soluo da equao de Darcy. Seja
o domnio do problema uma coleo de retngulos =i,jij e seja Eo conjunto doslados desses retngulos interiores a . De acordo com a seo anterior, a formulao
mista-hbridaconsiste em encontrar v
ijH(div, ij), p
L2() e
0,D (funes
em L2(E)que se anulam em D) tais que:
u T1v i,j
ij
( u)p ij
(u n)
=
(ug)
(u n)p
u ijH(div, ij) (3.41)
i,j
ij
q v=
q wqL2() (3.42)
i,j
ij
v n= 00,N (3.43)
Os multiplicadores de lagrange sero associados s presses nos lados das clulas. De
fato, para solues fortes, o teorema da divergncia nos diz que o somatrio na equao
(3.41) deve ser igual a
u p, implicando em p|E = . Essa hiptese vai ao encontrodo teorema de convergncia (8), e assim passaremos a identificar os multiplicadores de
Lagrange com as presses nas arestas das clulas.
Para encontrar uma soluo aproximada, trabalharemos em subespaos Vh ijH(div, ij),PhL2()e hL2(E). Assim, procuramos vhVh, phPhe hhtais que
47
7/26/2019 Dissertao Porosos 229
48/88
u T1v
i,jij
( u)p ij
(u n)
=
(u g)
(n u)p
uVh (3.44)
i,j
ij
q.v=
q wqPh (3.45)
i,jij v n= 0h. (3.46)
Retomando a seo anterior, escolhemos Vh=
v ijH(div, ij) :v|ij RT0(ij) i, j
,
Ph=pL2() :p|ij P0()i, j
eh={0,N; constante em cada aresta de E}.
Veremos agora o formato do sistema linear associado s equaes do mtodo misto-hbrido
3.5 Construo do sistema linear
Seja{ui}uma base para o espao Vh,{qi}uma base para Ph e{i}uma base para h.Definimosaij =
uiT1uj,bij =
.ui qj,cij =
Ei(uin)+
p (uin),gi =
uig
e wi= qi w. Definamos, ento, as matrizes A= (aij), B= (bij), C= (cij), g= (gi)e
w= (wi). As equaes (3.44)-(3.46) se tornam, ento:
Av + Bp+ C= g (3.47)
Btv= w (3.48)
Ctv= 0 (3.49)
Notemos que so permitidas funes descontnuas nos espaosVhe Ph. A vantagem disso
o desacoplamento das funes em cada clula.
Mostraremos agora que o sistema de equaes para os multiplicadores de lagrange
48
7/26/2019 Dissertao Porosos 229
49/88
positivo-definido.
3.6 Positividade definida do sistema linear
Examinando as equaes (3.47) - (3.49), vemos que o sistema linear associado da seguinte
forma:
A B C
Bt 0 0
Ct 0 0
v
p
=
g
w
0
(3.50)
Sabemos que a matrizA diagonal em blocos e positivo-definida. Sendo assim, resolvendo
a primeira linha de blocos da equao (3.50) para v, obtemos:
v= A1(g Bp C)
Substituindo nas outras linhas, obtemos o seguinte sistema em pe :
BtA1B BtA1C
CtA1B CtA1C
p
=
w BtA1g
CtA1g
.
Isolamos pna primeira linha:
p= (BtA1B)1BtA1C + BtA1gw
Substituimospna segunda linha, obtendo um sistema na forma H= R, no qual:
H=CtA1C CtA1B(BtA1B)1BtA1C.
Ou seja:
49
7/26/2019 Dissertao Porosos 229
50/88
H=CtA1I B(BtA1B)1BtA1
C.
A quantidade entre colchetes uma projeo ortogonal (pois simtrica e idempotente),o que nos mostra que H no-negativa. Invocando a unicidade da soluo do sistema
de equaes do mtodo hbrido misto (Teorema6), conclumos que o sistema positivo-
definido.
3.7 Implementao numrica do mtodo hbrido-misto
No caso em que temos as presses especificadas em toda a fronteira (ou seja, N =),obtemos uma grande simplificao do problema. Vejamos mais a fundo o formato do
sistema nesse caso. Para fixar as ideias, consideremos a clula[0, x] [0, y]. Podemostomar uma base de Vh na qual cada elemento nulo exceto em uma clula ij. Em
[0, x] [0, y]h quatro elementos que no so identicamente nulos. So eles:
uL(t, s) =
1 t
x , 0
uR(t, s) = tx
, 0
uU(t, s) =
0, 1 sy
uD(t, s) =
0, s
y
Dessa forma, podemos escrever, nessa clula, a velocidade na forma v= vLuL+vRuR+
vUuU+ vDuD, em que vL, vR, vU e vD indicam, respectivamente, a velocidade na aresta
esquerda, direita, superior e inferior da clula. Analogamente, denotaremos os multipli-
cadores de lagrange por pL,pR,pU e pD.
Tomaremos como base de Ph os elementos que valem 1 em apenas uma clula, e 0 nas
demais clulas. A base de h obtida de forma anloga. Definimos o tensor transmissi-
bilidadeT= K/e a matriz auxiliar S:
S=
x 0
0 y
T11 T12
T21 T22
1
x 0
0 y
1
xy
Nessas condies, a as equaes (3.47) e (3.48) tm o seguinte formato:
50
7/26/2019 Dissertao Porosos 229
51/88
13
S1116
S1114
S1214
S12 1
16
S1113
S1114
S1214
S12 114
S21 14S2113
S22 16S22 1
14
S2114
S2116
S2213
S22 1
1 1 1 1 0
vLy
vRy
vUx
vDx
p
=
pL
pRpU
pD0
Simularemos o caso em que K uma matriz diagonal. Nessa situao, temos um sistema
ainda mais simples, pois S diagonal:
13
S1116
S11 0 0 1
16
S1113
S11 0 0 10 0 1
3S22
16
S22 1
0 0 16
S2213
S22 1
1 1 1 1 0
vLy
vRy
vUx
vDx
p
=
pL
pRpU
pD0
(3.51)
O sistema (3.51) pode ser invertido, de forma a obter:
vLy
vRy
vUx
vDx
p
=
M1 M2 M4 M7 m1
M2 M3 M5 M8 m2
M4 M5 M6 M9 m3
M7 M8 M9 M10 m4
m1 m2 m3 m4
pL
pRpU
pD
0
(3.52)
Agora, (3.49) nos d:
(vL)i+1,j = (vR)i,j
(vU)i,j+1 = (vD)i,j(3.53)
em que os ndices i, j indicam, respectivamente, a posio horizontal e vertical da clula.As equaes (3.53) nos dizem que o fluxo volumtrico contnuo entre clulas, i.e., todo
o fluido que vem de uma clula transportado para a clula vizinha.
51
7/26/2019 Dissertao Porosos 229
52/88
Combinando (3.52), (3.53) e (3.53), podemos eliminar as velocidades do sistema linear,
resultando num sistema que envolve apenas as presses nas arestas das clulas. Alm disso,
como vimos na seo anterior, esse sistema positivo-definido, permitindo a utilizao do
mtodo do gradiente conjugado para sua resoluo. O mtodo ser ainda mais eficientepor se acoplar a matriz de rigidez malha, pois aquela esparsa e requer alguns cuidados
no armazenamento de forma a garantir a otimizao dos recursos computacionais.
Uma vez obtidas as presses, podemos usar as equaes (3.52) para obter as velocidades
em cada clula, e proceder para o clculo do transporte de massa convectivo, que ser
visto na prxima seo.
No caso em que as condies de fronteira determinam a presso (condio do tipo de
Dirichlet) em uma parte, e da velocidade (que pode ser vista como condio do tipo
de Neumann, pois v proporcional ao gradiente de p) na outra parte da fronteira, o
sistema no envolve apenas as presses, e o mtodo acima no mais utilizvel. Tambm
no podemos garantir que o sistema seja simtrico ou positivo-definido. Dessa forma,
um mtodo diferente para a resoluo do problema deve ser utilizado. Como a matriz de
rigidez ainda esparsa, escolhemos aqui o mtodo do Gradiente Biconjugado Estabilizado.
As experincias feitas indicaram a necessidade de um precondicionamento do sistema e,
para tal, utilizamos a decomposio LU incompleta. Assim, em qualquer tipo de condio
de fronteira, j sabemos como obter as velocidades, e assim resta-nos apenas resolver asequaes convectivas.
52
7/26/2019 Dissertao Porosos 229
53/88
Captulo 4
Discretizao das equaes
convectivo-difusivas
Aps integrar a Eq.(1.4) no tempo, podemos escrever a conservao da massa do tracer
em uma clula da malha por
c(x, tn+1) dV =
c(x, tn) dV +
tn+1
tnF n dtds. (4.1)
Definimos
A=
dV (4.2)
e
cn=
c(x, tn) dV
1
A. (4.3)
Nas faces Sda clula, temos as integrais no tempo para o fluxo volumtrico
Vn+ 12S =
S
tn+1
tnv n dtds. (4.4)
53
7/26/2019 Dissertao Porosos 229
54/88
A equao da diferena conservativa para o tracer, Eq.(4.1) toma a forma
cn+1 =c
n+
1
A
Sfn+ 1
2
S (4.5)
em quefn+ 1
2
S a integral no tempo do fluxo dotracerna faceSda clula. Vejamos como
se d o clculo desse fluxo.
Dividiremos o fluxo em uma parcela hiperblica (correspondente contribuio convec-
tiva) e em uma parcela difusiva:
fn+ 1
2
S =fH,n+ 1
2
S + fD,n+ 1
2
S (4.6)
A parcela convectiva definida por:
fH,n+ 1
2
S =cn+ 1
2
S Vn+ 1
2
S =S
tn+1tn
cv.ndtds. (4.7)
Aqui, a concentrao dotracerna faceSda clula,cn+ 1
2
S , calculada de diferentes manei-
ras, dependendo do nmero de dimenses espaciais. J o fluxo difusivo ser calculado por
diferenas finitas explcitas. Trataremos do fluxo convectivo na seo seguinte. A parcela
difusiva ser tratada numa seo posterior.
4.1 Fluxo convectivo em 1D
Por razes ilustrativas, exibimos o clculo dos fluxos convectivos em 1D. Esse clculo ser
feito utilizando uma variao do conhecidoMUSCL Scheme, que consiste em uma recons-
truo linear por partes da concentrao do tracer, seguido do uso das retas caractersticas
para obter as concentraes nas faces da clula. A reconstruo detalhada a seguir.
Definimos
54
7/26/2019 Dissertao Porosos 229
55/88
muscl(a, b) =
sign (a)min {2 |a| , 2 |b| , |a + b| /2} se ab00 se ab
7/26/2019 Dissertao Porosos 229
56/88
fH,n+ 1
2
i+ 12
=cn+ 1
2
i+ 12
Vn+ 1
2
i+ 12
. (4.13)
possvel determinar as circunstncias nas quais esse mtodo ser estvel e convergente.
Lembramos que, na ausncia de poos em 1D, a velocidade independente de x (j que
divv = w = 0), e assim Vn+ 1
2
i+ 12
=vntn+1
2 . Utilizando o Lema de Harten [32], vemos que
esse mtodo ser TVD (Total Variation Diminishing) se:
i tn+12 tHi xi|vn| . (4.14)
Fisicamente, isso significa que o fluxo de fluido (i.e., vt) se restrinja a, no mximo, o
volume de vazio da clula (i.e., xi). fato conhecido que os mtodos TVD so estveis
e convergentes[32].
4.2 Fluxo convectivo em 2D
Em 2D, usaremos uma verso dodonor cell upwind scheme([53]). Nesse mtodo, a massa
do tracerna clula aps um intervalo de tempo igual ao valor no tempo anterior mais
ou menos correes nas regies trapezoidais associadas s arestas das clulas.
Como no caso unidimensional, utilizamos uma reconstruo linear por partes. Em cada
clula, calculamos os incrementos:
cn1,ij =muscl(cnij cni1,j, cni+1,j cnij)
e
cn2,ij =muscl(cnij cni,j1, cni,j+1 cnij)
alm do volume dos poros na clula:
56
7/26/2019 Dissertao Porosos 229
57/88
Aij =x1ix2j. (4.15)
Definindo a aproximao
Vn+ 1
2
i+ 12,j(vR)i,jx2jtn+ 12 (4.16)
e equaes anlogas para Vn+ 1
2
i,j+ 12
(lembrando que vL, vR, vU e vD foram definidos no
Captulo 3), computamos as concentraes nas arestas:
cni+ 12,j =
cnij+ 12
1 V
n+12
i+12,j
Aij
cn1,ij se Vn+ 12i+ 1
2,j0
cni+1,j 12
1 + V
n+12
i+12,j
Ai+1,j
cn1,i+1,j se Vn+ 12i+ 1
2,j
7/26/2019 Dissertao Porosos 229
58/88
cn+1ij =cnij
1
Aij
fn+ 1
2
i+ 12,j fn+ 12
i 12,j
+ fn+ 1
2
i,j 12
fn+ 12i,j 1
2
. (4.21)
As equaes (4.15) - (4.21) constituem o esquema 2D para o clculo das concentraes
em cada instante de tempo. De acordo com[52], o esquema ser estvel e de segunda
ordem caso
max
Vn+ 1
2+
i+ 12,j Vn+
1
2
i 12,j
Aij,
Vn+ 1
2+
i,j+ 12
Vn+1
2
i,j 12
Aij
1 (4.22)
em que V+ e V so as partes positiva e negativa de V:
V+ = max {V, 0} e V = min {V, 0} . (4.23)
No caso de uma malha uniforme, podemos equivalentemente definir os incrementos tem-
porais direcionais
1
t1,ij=vn+ 1
2+
i+ 12,j vn+
1
2
i+ 12,j
ijxij(4.24)
e
1t2,ij
= v
n+ 12+
i,j+ 12 vn+ 1
2
i,j+ 12
ijxij(4.25)
restringindo o incremento de tempotn+1
2 de forma que
tn+1
2 max
1
t1,ij,
1
t2,ij
1. (4.26)
Note que no exigimos que o mtodo bidimensional seja TVD, pois fato conhecido que
os esquemas TVD so, no mximo, de primeira ordem [31].
58
7/26/2019 Dissertao Porosos 229
59/88
7/26/2019 Dissertao Porosos 229
60/88
Captulo 5
Refinamento Adaptativo da Malha
5.1 Malhas Adaptativas
A soluo numrica de EDPs pode muitas vezes exigir uma malha no-uniforme, que se
adapte ao problema, sendo mais refinada onde a soluo ou as suas derivadas variem mais
rapidamente e menos refinada onde essa variao ocorra menos abruptamente. Em siste-mas dependentes do tempo, essa necessidade mais expressiva, uma vez que a natureza
dinmica desses problemas pode produzir regies de rpida variao na soluo.
H essencialmente dois mtodos de se modificar malhas. A primeira mov-las juntamente
com as frentes de onda (vide [4,14,15,38]). O segundo mtodo no movimenta a malha
no tempo, mas refina as regies de grande variao da soluo. Nessa situao, a cada
instante de tempo, a malha deve ser atualizada de forma a refletir a dinmica do problema
(veja[3,5,8,19,39,41, 43,44,45]). Alguns autores, inclusive, preferem combinar os dois
mtodos, como pode ser visto em[2].
Uma malha formada por clulas que cobrem todo o domnio estudado, cada uma delas
associada a um n da estrutura de dados. Dada uma malha inicial, que representada
como a raiz de uma rvore, ns-filho com graus crescentes de refinamento podem ser
criados. Cada n pode agora ser um n-paipara os refinamentos seguintes. Esse procedi-
mento pode ser executado em regies que atendam a um certo critrio de refinamento de
acordo com o problema em questo. Dessa forma, diferentes graus de refinamento podem
ser encontrados em regies adjacentes do domnio estudado [12].
60
7/26/2019 Dissertao Porosos 229
61/88
De acordo com [47], esse procedimento conhecido como quadtree, uma classe de es-
truturas de dados hierrquica cuja propriedade em comum a todas serem baseadas no
princpio da decomposio recursiva do espao. Em [36], uma nova estrutura de dados,
chamada de fully threaded tree(FTT), apresentada e aplicada a simulao do escoa-mento de fluidos. FFT representa um avano comparado s abordagens baseadas em
rvores mencionadas. Cada n (seja ele uma folhaou no) possui fcil acesso a seus
filhos, vizinhos e pais, por meio de ponteiros especiais. Mais ainda, no FFT, o acesso a
clulas vizinhas e mais rpido, pois no necessrio percorrer toda a rvore para tal.
Implementaes do tipoquadtreespermitem que regies adjacentes do domnio possuam
diferentes nveis de refinamento. Como consequncia, para alcanar clulas vizinhas, o
algoritmo realiza uma busca na rvore. FFT consegue evitar esses custos permitindo a
diferena de apenas um nvel de refinamento entre clulas vizinhas. No entanto, essa
restrio impe um efeito de mdio alcance na malha: se uma clula necessita de mais
refinamento, seus vizinhos tero de ser refinados tambm. H, ento, uma perda de
localidade no refinamento, impondo um custo computacional adicional.
Nesse trabalho, ser utilizada a estrutura de dados Autonomous Leaves Graph (ALG),
descrita em [12]. Trata-se de uma estrutura simples, embora flexvel, que nos propor-
ciona um refinamento local da malha, a um baixo custo computacional. Nele, em vez
da estrutura quadtreeconvencional, criado um grafo no qual os ns filho (os quais se-ro chamados folhas) tornam-se autnomos quando os respectivos ns paiso deletados.
Alm disso, essa estrutura permite a coexistncia de vizinhos com graus de refinamento
(arbitrariamente) diferentes por meio dos chamados ns de transio bem como a
comunicao rpida entre esses vizinhos ([12]).
O ordenamento das clulas feito por meio de uma space-filling curve. A vantagem desse
tipo de curva a utilizao de um mtodo simples e sistemtico para visitar todas as clulas
da malha. Alm disso, space-filling curvespossuem um alto grau de localidade (vide [46,47]). Nesse trabalho, ser utilizada a chamada Curva de Hilbert Modificada, uma verso da
curva de Hilbert tradicional, mas adaptada a malhas de diferentes nveis de refinamento.
Essa curva foi proposta em[13] e descrita em detalhes em [15]. Apresentaremos agora o
funcionamento do ALG.
61
7/26/2019 Dissertao Porosos 229
62/88
5.2 A estrutura de dadosAutonomous Leaves Graph
(ALG)
A explicao e as ilustraes dessa seo at o fim do captulo foram retirados de [ 12]. A
estrutura de grafo utilizada no ALG constituda por dois tipos de ns: oscell nodes, que
correspondem s clulas da malha, e os transition nodes, que so usados para conectar
cell nodesde diferentes nveis de refinamento.
No caso bidimensional, cadacell nodepossui quatro ponteiros, chamados norte, sul, leste
e oeste, orientados respectivamente na direo do eixo y positivo, y negativo, x positivo e
x negativo. Esses ponteiros podem apontar tanto para cell nodesquanto para transitionnodes. Cadatransition nodepossui apenas trs ponteiros, para conectarcelloutransition
nodes.
Os cell nodespossuem variveis adicionais correspondentes s coordenadas espaciais dos
seus centros e as propriedades fsicas pertinentes. Ambos os ns possuem informao sobre
o seu tipo (cell ou transition node) e sobre o seu nvel de refinamento. Por conveno,
representaremos oscell nodespor pontos negros, e ostransition nodespor pontos brancos.
5.2.1 A estrutura da malha
Para descrever a malha, mostraremos o caso do quadrado unitrio (embora a generaliza-
o possa ser feita para malhas com polgonos arbitrrios, em duas ou mais dimenses
espaciais). Comecemos com uma malha constituda por quatro clulas, cada uma delas
identificadas pelo seu centro: (0.25, 0.25), (0.75, 0.25), (0.75, 0.25) e (0.75, 0.75). Cada
um dos quatro cell nodespossuem elos orientados nas quatro direes: norte, sul, leste eoeste. Como resultado, temos o diagrama representado na Figura5.1.
Figura 5.1: Malha de 4 clulas
Os elos restantes, que no apontam para um dos quatro ns no quadrado so ento
62
7/26/2019 Dissertao Porosos 229
63/88
conectados aos quatro transition nodesmostrados na Figura5.2.
Figura 5.2: Definindo os ns de transio
Finalmente, como mostrado na Figura5.3, essestransition nodesso aterrados (ou seja,apontam para um ponteiro nulo).
Figura 5.3: Esquema completo
De forma a simplificar a representao grfica, cada par de setas conectando dois ns ser
substitudos por uma nica linha, como mostrado na Figura5.4.
Quatro clulas com o mesmo nvel de refinamento e originadas do mesmo vertic v so
chamadas decacho(bunch), com pai v. O n pai deletado, mas cada clula munida deum identificador de forma a determinar que as quatro clulas pertencem ao mesmo cacho.
5.3 Refinamento da malha
Para mostrar como se realiza o refinamento da malha, suponha que a clula centrada em
(0.25, 0.75)seja escolhida para ser refinada. Nesse caso, a configurao mostrada em naFigura5.5(a) ser criada:
63
7/26/2019 Dissertao Porosos 229
64/88
Figura 5.4: Esquema no-direcional
0101
01
001101
0
0
1
1
00
00
11
11
01 01
01
00
00
11
11
0
0
1
1
01
0
0
1
1
0
0
1
1
00
00
11
11
0
0
1
1
01001101
0011 001100
00
11
11
00
00
11
11
00
00
11
11
00
00
11
11
00
00
11
11
00
00
11
11
00
00
11
11
00
00
11
11
(a) (c)(b)
Figura 5.5: Sucesso de refinamentos
Esse refinamento implementado substituindo a estrutura bsica representando os elos do
n em questo (parte esquerda da Fig5.6) pela estrutura desenvolvida para o quadrado
unitrio (parte direita da Fig 5.7). Cada n (a, b) no centro da clula de lado c a ser
refinada substituido pelos ns (a c/4, b c/4), (a c/4, b+c/4), (a+c/4, b+c/4)e(a + c/4, b c/4). Os quatro elos partindo dos transition nodesso ento conectados aosvizinhos do n que est sendo substituido. O grafo resultante exibido na figura 5.7; o
nmero prximo a cada n indica o seu nvel de refinamento.
Figura 5.6: Elos de um n e estrutura bsica de refinamento
Consideremos agora um nvel a mais de refinamento, como mostrado na Figura5.5(b). A
aplicao do mesmo princpio para o refinamento resulta no grafo da Figura5.8:
64
7/26/2019 Dissertao Porosos 229
65/88
1
1
1
1
1
1
1
2
2
2
2
2
2
2
2
Figura 5.7: Refinando a clula noroeste
1
1
1
1
1
1
1
2
2
2
2
2
2
2
3
3
3
3
3
33
Figura 5.8: Mais um nvel de refinamento
5.4 Desrefino da malha
Antes de entrar em detalhes sobre o processo de desrefino, devemos ter em mente que,
para manter a estrutura bsica de blocos do grafo, necessrio que todas as quatro
clulas que sero desrefinadas pertenam ao mesmo bunch, permitindo a recuperao da
configurao anterior. Essa condio garante que configuraes antigas da malha possam
ser recuperadas executando-se mudanas locais mutuamente independentes, mesmo aps
um nmero arbitrariamente grade de refinos e desrefinos. O procedimento de derrefinopode ser esquematizado da seguinte forma:
1. Transformao do cacho a ser desrefinado em um nico n
2. Alimentao desse n com as variveis pertinentes (coordenadas espaciais, nvel e
parmetros fsicos relativos ao problema)
3. Conexo do novo n com seus vizinhos
4. Simplificao da malha, eliminando ns desnecessrios
65
7/26/2019 Dissertao Porosos 229
66/88
A conexo do novo n com seus vizinhos feita com a criao de quatro transition no-
des, com o mesmo nvel de refinamento do bunch original. Ilustremos a sequncia de
desrefinamento com as figuras abaixo:
0011
0011
0011
0
0
1
1
00
00
11
11
010011
0
0
1
1
00
00
11
11
0
0
1
1
01
0101
(a) Seleo do cacho para desrefinamento
0011
0011
010011
0
0
1
1
00
00
11
11
0
0
1
1
0101
0011
(b) Colapso do cacho
0011
0011
010011
0
0
1
1
00
00
11
11
0
0
1
1
0101
0011
01
(c) Criao dos ns de transio
0011
0011
01
0
0
1
1
00
00
11
11
0
0
1
1
0101
0011
01
0011
(d) Configurao final aps simplificao
Figura 5.9: Sequncia de desrefinamento
Adotando as estratgias acima, a configurao mostrada na Figura 5.9atorna-se aquela
mostrada na Figura5.9d. Note as trs simplificaes que ocorreram nesse exemplo.
5.5 Ordenao Total da Malha
A ordenao das clulas da malha feita utilizando uma Curva de Hilbert Modificada
(MHC). Trata-se de um algoritmo baseado na curva de Hilbert original, mas adaptado
diviso do domnio em uma malha no-uniforme. Essa ordenao implementada pormeio de uma lista, e a cada vez que um refinamento executado, as novas clulas so
inseridas nessa lista. O procedimento detalhado para a construo da MHC pode ser
encontrado em [12]. Podemos ver um exemplo da MHC na Figura 5.10abaixo.
66
7/26/2019 Dissertao Porosos 229
67/88
Figura 5.10: Exemplo de ordenao pela Curva de Hilbert Modificada
5.6 Refinando e Desrefinando
Descreveremos aqui como se d o refinamento e o desrefinamento da malha em termos
das propriedades fsicas armazenadas nas clulas.
Temos, no problema, um campo de velocidades esttico, em regime estacionrio. Esse
campo causa o transporte convectivo do tracerde clula para clula. Utilizamos como
critrio de refinamento a norma (da soma) do gradiente de concentrao; caso o seu valor
fosse superior a um limite predeterminado, a clula seria refinada.
Resta-nos saber como calcular os valores de presso e velocidade nas clulas refinadas.
Nossa primeira abordagem foi realizar um clculo local, resolvendo a Equao de Darcy
com condies de contorno iguais s propriedades da clula a ser refinada. Esse pro-
cedimento mostrou-se ineficaz, pois a soluo obtida se afastava muito daquela obtida
com a malha uniforme (como nosso problema no tem soluo analtica, a alternativa a
comparao com os resultados da aplicao do mesmo mtodo a uma malha uniforme).
Dessa forma, apelamos para um mtodo diferente. No incio da simulao, resolve-se aEquao de Darcy no maior nvel de refinamento permitido. Armazenamos esse resultado
parte na memria (numa malha auxiliar), e buscamos informao ali sempre que ne-
cessrio. Ao refinar uma clula da malha, percorremos essa malha auxiliar (que contm
a melhor informao fsica disponvel), e calculamos as mdias necessrias para definir as
propriedades das clulas refinadas. Denotaremos por Ie Ja posio espacial de uma c-
lula da malha do problema e por ie j a posio de uma clula na malha auxiliar. Quando
a aresta da malha auxiliar Si+ 12,j est contida na aresta SI+ 1
2,J, calculamos as mdias da
seguinte maneira:
67
7/26/2019 Dissertao Porosos 229
68/88
VI+ 12,J=
Si+1
2,jS
I+12,J
Vi+ 12,j (5.1)
pI+ 12,J=
1
yIJ
Si+1
2,jS
I+12,J
pi+ 12,jyij (5.2)
e se a clula ij est contida emIJ:
pIJ= 1
xIJyIJ
ijIJ
pijxijyij. (5.3)
Essas equaes garantem a conservao da massa do fluido. A Eq. (5.1) diz que o fluxo
volumtrico na aresta de uma clula menos refinada a soma dos fluxos volumtricos nas
arestas das clulas mais refinadas. Equaes (5.2) e (5.3) nos dizem que as presses, tanto
nas faces quanto nos centros das clulas, so mdias das presses nas respectivas clulas
mais refinadas. Quanto concentrao do tracer, as clulas refinadas recebem o mesmo
valor da clula menos refinada original, pois no h como se determinar qual a proporo
de tracer que deve ir para cada nova clula.
O procedimento de resolver a Equao de Darcy inicialmente para toda a malha no nvel
de refinamento mximo requer um certo esforo computacional. Porm, lembramos que
o sistema linear ser resolvido apenas uma vez: os clculos para a atualizao da concen-
trao do tracernas clulas usam um mtodo explcito. Assim, esse custo computacional
inicial compensado por no termos que resolver um novo sistema para determinar o
campo de velocidades a cada iterao na malha no-uniforme.
O desrefinamento inteiramente anlogo: caso um grupo de clulas com o mesmo nvel de
refinamento e pertencentes ao mesmo cacho tenham um gradiente de concentrao menor
que um valor limite, efetuamos o desrefinamento, como visto na seo 5.2. Nesse caso,
para esse grupo de clulas, calculamos as mdias pelas equaes (5.1) - (5.3). Alm disso,
calculamos a concentrao mdia do tracer:
cIJ= 1
xIJyIJ
ijIJ
cijxijyij. (5.4)
Esses valores mdios sero as propriedades fsicas da clula desrefinada. Note que a
68
7/26/2019 Dissertao Porosos 229
69/88
Equao (5.4) conserva a massa em cada clula, como era de se esperar.
Podemos resumir, ento, o procedimento adotado na simulao no pseudocodigo:
1: Definir o nvel mximo de refinamento
2: Definir as condies de contorno
3: Montar a matriz de rigidez da malha
4: ifPresso especificada em toda a fronteira then
5: Resoluo do sistema pelo CG
6: else
7: Fatorao LU incompleta
8: Resoluo do sistema pelo BiCGSTAB9: end if
10: Definir condio inicial das concentraes
11: whileMassa total for positiva do
12: Refinar\Desrefinar a malha13: Calcular tpara esta iterao
14: Atualizar as concentraes nas clulas
15: end while
69
7/26/2019 Dissertao Porosos 229
70/88
7/26/2019 Dissertao Porosos 229
71/88
6.1.1 Exemplo Numrico 1
Para ilustrar a aplicao do mtodo ao caso de presses especificadas, utilizamos uma
presso igual a 2 atm nas fronteiras norte e oeste, e igual a 1 atm nas fronteiras sul eleste de um domnio quadrado, de lado 1m. Foi utilizado um campo de permeabilidades
uniforme e isotrpico igual a 1000 milidarcy, e uma porosidade igual a 0, 27. O fluido
simulado foi a gua, com viscosidade igual a 1cP.
Para a soluo do sistema, aplicamos o Gradiente Conjugado (CG) ao sistema linear
obtido. Por envolver apenas multiplicaes matriz-vetor, uma forma bastante eficiente de
se aplicar o CG utilizar a prpria malha para o armazenamento da matriz e para realizar
as operaes necessrias. Dessa forma, nenhum pacote extra de lgebra linear precisa serutilizado, otimizando o uso de memria e de processamento.
Exibimos na Figura6.1o campo de presses obtido em diferentes nveis de refinamento.
6.1.2 Exemplo Numrico 2
Gostaramos de verificar a exatido do mtodo empregado. No caso geral, no existe
soluo analtica para a Lei de Darcy. Porm, condies de contorno bem-escolhidas nos
levam a um resultado conhecido. Sendo assim, realizamos a simulao com condio de
contorno mista: v= 0(no-flow) nas fronteiras norte e sul,p = 2atm na fronteira oeste e
p= 1atm na fronteira leste.
Os testes iniciais indicaram a necessidade de um precondicionamento: o solverno con-
vergia utilizando simplesmente o BiCGSTAB aplicado a esse sistema. Buscamos, ento,
alternativas que fossem satisfatrias em termos de eficincia e praticidade. H vrias bibli-otecas paraC++especializadas em lgebra linear, como por exemplo LAPACK,Amadillo
eEigen. Existem tambm aplicaes como o SciLab e MATLAB(R) que possuem mdulos
para a resoluo de sistemas.
O programa que demonstrou melhor eficincia na resoluo do sistema, alm de oferecer
as ferramentas necessrias para seu uso, foi o MATLAB(R). Sendo assim, utilizamos
uma aplicao desenvolvida em MATLAB(R) e compilada (via MATLAB(R) Compiler)
num arquivo executvel standalone, que pode ser executado mesmo em mquinas queno possuem MATLAB(R) instalado. Esse arquivo executvel l a matriz de rigidez e o
vetor de cargas, realiza a fatorao LU incompleta (iLU) e resolve o sistema utilizando
71
7/26/2019 Dissertao Porosos 229
72/88
(a) Nvel 4 (b) Nvel 6
(c) Nvel 7 (d) Nvel 8
Figura 6.1: Campo de presses em diferentes nveis de refinamento
72
7/26/2019 Dissertao Porosos 229
73/88
o BiCGSTAB. Os arquivos gerados pelo MATLAB(R) Compilergozam de livre uso e
distribuio, no implicando em problemas de licena na utilizao do simulador.
Utilizando esse programa misto, podemos resolver o sistema desse exemplo (e todos osque utilizam uma condio de fronteira mista). Mostramos na Figura6.2abaixo o campo
de presses obtido em vrios nveis de refinamento.
(a) Nvel 4 (b) Nvel 6
(c) Nvel 7 (d) Nvel 8
Figura 6.2: Campo de presses em diferentes nveis de refinamento
Para essas condies de contorno, fcil ver que a soluo analtica :
73
7/26/2019 Dissertao Porosos 229
74/88
p(x, y) = 2 xatm
v(x, y) =
K, 0
= (1000, 0)mD-atm
cP-m (8.64, 0)m/dia.
Dessa forma, podemos comparar a soluo obtida na simulao com a soluo analtica.
Realizamos a simula
Recommended