UM MÉTODO CONSISTENTE DE GALERKIN DESCONTÍNUOESTABILIZADO PARA PROBLEMAS DE QUARTA ORDEM COM

RESTRIÇÃO DE INCOMPRESSIBILIDADE

Antônio Guilherme Barbosa da Cruz

Tese de Doutorado apresentada ao Programade Pós-graduação em Engenharia Mecânica,COPPE, da Universidade Federal do Rio deJaneiro, como parte dos requisitos necessáriosà obtenção do título de Doutor em EngenhariaMecânica.

Orientadores: Fernando Pereira DudaEduardo Gomes Dutra doCarmo

Rio de JaneiroMarço de 2011

UM MÉTODO CONSISTENTE DE GALERKIN DESCONTÍNUOESTABILIZADO PARA PROBLEMAS DE QUARTA ORDEM COM

RESTRIÇÃO DE INCOMPRESSIBILIDADE

Antônio Guilherme Barbosa da Cruz

TESE SUBMETIDA AO CORPO DOCENTE DO INSTITUTO ALBERTO LUIZCOIMBRA DE PÓS-GRADUAÇÃO E PESQUISA DE ENGENHARIA (COPPE)DA UNIVERSIDADE FEDERAL DO RIO DE JANEIRO COMO PARTE DOSREQUISITOS NECESSÁRIOS PARA A OBTENÇÃO DO GRAU DE DOUTOREM CIÊNCIAS EM ENGENHARIA MECÂNICA.

Examinada por:

Prof. Fernando Pereira Duda, D.Sc.

Prof. Eduardo Gomes Dutra do Carmo, D.Sc.

Prof. José Herskovits Norman, D.Sc.

Prof. Webe João Mansur, Ph.D.

Prof. Antônio José da Silva Neto, Ph.D.

RIO DE JANEIRO, RJ – BRASILMARÇO DE 2011

Cruz, Antônio Guilherme Barbosa daUm método consistente de Galerkin descontínuo

estabilizado para problemas de quarta ordem com restriçãode incompressibilidade/Antônio Guilherme Barbosa daCruz. – Rio de Janeiro: UFRJ/COPPE, 2011.

XI, 53 p.: il.; 29, 7cm.Orientadores: Fernando Pereira Duda

Eduardo Gomes Dutra do CarmoTese (doutorado) – UFRJ/COPPE/Programa de

Engenharia Mecânica, 2011.Referências Bibliográficas: p. 47 – 53.1. Método de elementos finitos. 2. Segundo

gradiente. 3. Problema de escoamento incompressívelde quarta ordem. I. Duda, Fernando Pereira et al.II. Universidade Federal do Rio de Janeiro, COPPE,Programa de Engenharia Mecânica. III. Título.

iii

PARAminha mãe Nazaré, Lourdinha,

Sheila e Maiara.

iv

Agradecimentos

A minha mãe, por acreditar em meus sonhos. Por enxergar que, quanto melhor éadquirir a sabedoria do que o ouro, e quanto mais excelente é escolher o entendi-mento do que a prata.

A minha companheira Lourdinha, pela caminhada de apoio a cada palavraescrita nesta tese.

Aos meus orientadores, Eduardo Gomes Dutra do Carmo e Fernando PereiraDuda, pela amizade e encorajamento nos caminhos da pesquisa científica.

Aos meus colegas do Laboratório de Mecânica do Sólidos, entre eles, Gabriel,Fábio e Zio.

A CAPES, pela concessão da bolsa de doutorado, e a Fundação COPPE-TEC por intermédio do Projeto Galileu da Petrobrás, pela ajuda financeira quepossibilitou o término deste doutoramento.

v

Resumo da Tese apresentada à COPPE/UFRJ como parte dos requisitos necessáriospara a obtenção do grau de Doutor em Ciências (D.Sc.)

UM MÉTODO CONSISTENTE DE GALERKIN DESCONTÍNUOESTABILIZADO PARA PROBLEMAS DE QUARTA ORDEM COM

RESTRIÇÃO DE INCOMPRESSIBILIDADE

Antônio Guilherme Barbosa da Cruz

Março/2011

Orientadores: Fernando Pereira DudaEduardo Gomes Dutra do Carmo

Programa: Engenharia Mecânica

Este trabalho apresenta uma formulação consistente e estabilizada de elementosfinitos para problema de escoamento incompressível de quarta ordem. A formulaçãoé baseada no método de penalização C0-interior e no esquema de Galerkin mínimosquadrados (GLS), que assegura que a formulação seja fracamente coerciva para es-paços que não satisfazem a condição inf-sup, e considera interpolações descontínuasda pressão. Experimentos são realizados para ilustrar a capacidade da formulação delidar com interpolações arbitrárias para a velocidade e pressão, e de estabilizar gra-dientes de pressão elevados. A análise da estabilidade é apresentada através de umlema que engloba a condição inf-sup, e confirma a robustez da formulação proposta.Este lema indica que existe, a nível de elemento, um parâmetro de estabilidade GLSótimo ou quase ótimo que depende do grau do polinômio usado para interpolar oscampos de velocidade e pressão, da geometria do elemento da malha, e do termo deviscosidade do fluido.

vi

Abstract of Thesis presented to COPPE/UFRJ as a partial fulfillment of therequirements for the degree of Doctor of Science (D.Sc.)

A CONSISTENT STABILIZED DISCONTINUOUS GALERKIN METHOD FORFOURTH-ORDER PROBLEMS WITH INCOMPRESSIBILITY CONSTRAINT

Antônio Guilherme Barbosa da Cruz

March/2011

Advisors: Fernando Pereira DudaEduardo Gomes Dutra do Carmo

Department: Mechanical Engineering

This paper presents a consistent and stabilized finite-element formulation forfourth-order incompressible flow problems. The formulation is based on the C0-interior penalty method, the Galerkin least-square (GLS) scheme, which assuresthat the formulation is weakly coercive for spaces that fail to satisfy the inf-sup con-dition, and considers discontinuous pressure interpolations. Numerical experimentsare carried out to illustrate the ability of the formulation to deal with arbitraryinterpolations for velocity and pressure, and to stabilize large pressure gradients. Astability analysis is presented via a lemma that encompasses the inf-sup conditionand confirms the robustness of the proposed formulation. This lemma indicatesthat, at the element level, there exists an optimal or quasi–optimal GLS stabilityparameter that depends on the polynomial degree used to interpolate the velocityand pressure fields, the geometry of the finite element, and the fluid viscosity term.

vii

Sumário

Lista de Figuras ix

1 Introdução 1

2 Problema de quarta ordem com restrição de incompressibilidade 52.1 Problema modelo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52.2 Formulação variacional associada . . . . . . . . . . . . . . . . . . . . 6

3 Aproximação por elementos finitos 83.1 Métodos de elementos finitos para problemas de quarta ordem . . . . 83.2 Método de penalização C0-interior para o problema incompressível de

quarta ordem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12

4 Estabilização da formulação de Galerkin contínuo/descontínuo paraescoamento incompressível 174.1 Oscilações espúrias e a necessidade de estabilização em problemas

incompressíveis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 174.2 Estabilização pelo método de Galerkin mínimos quadrados . . . . . . 19

5 Análise de estabilidade do problema de quarta ordem 21

6 Resultados numéricos 296.1 Caso 1: Escoamento plano de Poiseuille . . . . . . . . . . . . . . . . . 296.2 Caso 2: Escoamento em cavidade . . . . . . . . . . . . . . . . . . . . 35

7 Conclusões e sugestões de trabalhos futuros 45

Referências Bibliográficas 47

viii

Lista de Figuras

3.1 Elementos cúbicos de Hermite P3 (4 graus de liberdade). . . . . . . . 93.2 Elemento triangular P5 Argyris (21 graus de liberdade); adaptado de

BRENNER e SUNG [1]. . . . . . . . . . . . . . . . . . . . . . . . . . 103.3 Elemento retangular Q3 Bogner-Fox-Schmit (16 graus de liberdade);

adaptado de BRENNER e SUNG [1]. . . . . . . . . . . . . . . . . . . 10

6.1 Domínio do escoamento em canal plano e condições de contorno. . . . 306.2 Escoamento em canal: campo de pressão estabilizado obtido pela

formulação CDGLS, α=1.5, elementos finitos de Lagrange Q2/Q1 compressão descontínua. . . . . . . . . . . . . . . . . . . . . . . . . . . . 31

6.3 Escoamento em canal: campo de pressão estabilizado obtido pelaformulação CDGLS, α=0,05, elementos finitos com Lagrange Q2/Q1

com pressão descontínua. . . . . . . . . . . . . . . . . . . . . . . . . . 316.4 Escoamento em canal: campo de pressão estabilizado obtido pela

formulação CDGLS, α=0,005, elementos finitos de Lagrange Q2/Q1

com pressão descontínua. . . . . . . . . . . . . . . . . . . . . . . . . . 326.5 Escoamento em canal: perfis de velocidade exatos e numéricos através

da linha de centro x/10d = 0, 5 para diferentes valores do parâmetrode estabilização α, comparados à solução do problema de Poiseuilleclássico. Todos os resultados numéricos foram obtidos com uma ma-lha construída com 100 elementos Lagrange Q2/Q1 com pressão des-contínua. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32

6.6 Escoamento em canal: perfis de velocidade exatos e numéricos atravésde x/10d = 0, 5 para a razão l/ℓ iguais a 0, 1, 1, 0 e 3, 0. Todos osresultados numéricos foram obtidos com uma malha construída com100 elementos finitos de Lagrange Q2/Q1 com pressão descontínua. . 33

6.7 Escoamento em canal: campo de pressão oscilante obtidos pela for-mulação CDGLS, α=0,000015, e elementos finitos de Lagrange Q2/Q2

com pressão contínua. . . . . . . . . . . . . . . . . . . . . . . . . . . 34

ix

6.8 Escoamento em canal: campo de pressão oscilante obtidos pela for-mulação CDGLS, α=0,000015, e elementos finitos de Lagrange Q2/Q2

com pressão descontínua. . . . . . . . . . . . . . . . . . . . . . . . . . 346.9 Escoamento em canal: campo de pressão estabilizado obtido pela

formulação CDGLS, α=0,05, e elementos finitos de Lagrange Q2/Q2

com pressão descontínua. . . . . . . . . . . . . . . . . . . . . . . . . . 356.10 Domínio do escoamento em cavidade e condições de contorno. . . . . 366.11 Escoamento em cavidade: campo de pressão oscilante obtido com

elementos de Lagrange Q2/Q1 com pressão contínua, e parâmetroα=0,0005 para o método CDGLS. . . . . . . . . . . . . . . . . . . . . 36

6.12 Escoamento em cavidade: campo de pressão oscilante obtido comelementos de Lagrange Q2/Q1 com pressão descontínua, e parâmetroα=0,0005 para o método CDGLS. . . . . . . . . . . . . . . . . . . . . 37

6.13 Escoamento em cavidade: campo de pressão estabilizado obtido comelementos de Lagrange Q2/Q1 com pressão descontínua, e parâmetroα=0,01 para o método CDGLS. . . . . . . . . . . . . . . . . . . . . . 37

6.14 Escoamento em cavidade: campo de pressão estabilizado obtido comelementos de Lagrange Q2/Q1 com pressão descontínua, e parâmetroα=0,1 para o método CDGLS. . . . . . . . . . . . . . . . . . . . . . . 38

6.15 Escoamento em cavidade: distribuição de pressão (descontínua) aolongo da linha y = 0, 5 mm, malha contendo 20 × 20 elementos. . . . 38

6.16 Escoamento em cavidade: perfis de velocidade normalizados ao longoda linha x/d=0,5, malha contendo 20 × 20 elementos. . . . . . . . . . 39

6.17 Escoamento em cavidade: perfis de velocidade normalizados ao longoda linha y/d=0,5; malha contendo 20 × 20 elementos. . . . . . . . . . 39

6.18 Escoamento em cavidade: campo de pressão oscilante obtido comelementos de Lagrange Q2/Q2 com pressão contínua, e parâmetroα=0,00025 para o método CDGLS; malha contendo 10 × 10 elementos. 41

6.19 Escoamento em cavidade: campo de pressão oscilante obtido comelementos de Lagrange Q2/Q2 com pressão descontínua, e parâmetroα=0,00025 para o método CDGLS; malha contendo 10 × 10 elementos. 41

6.20 Escoamento em cavidade: campo de pressão estabilizado obtido comelementos de Lagrange Q2/Q2 com pressão descontínua, e parâmetroα=0,005 para o método CDGLS; malha contendo 10 × 10 elementos. 42

6.21 Escoamento em cavidade: campo de pressão estabilizado obtido comelementos de Lagrange Q2/Q2 com pressão descontínua, e parâmetroα=0,035 para o método CDGLS; malha contendo 10 × 10 elementos. 42

x

6.22 Escoamento em cavidade: campo de pressão estabilizado obtido comelementos de Lagrange Q2/Q2, e parâmetro de estabilização α=0,1para o método CDGLS; malha contendo 10 × 10 elementos. . . . . . . 43

6.23 Escoamento em cavidade: distribuição da pressão (descontínua) aolongo da linha central y = 0, 5 mm; malha contento 10 × 10 elementos. 43

6.24 Escoamento em cavidade: perfil de velocidade normalizados ao longoda linha x/d = 0, 5; malha contendo 10 × 10 elementos. . . . . . . . . 44

6.25 Escoamento em cavidade: perfil de velocidade normalizados ao longoda linha y/d = 0, 5; malha contendo 10 × 10 elementos. . . . . . . . . 44

xi

Capítulo 1

Introdução

O desenvolvimento de formulações de elementos finitos para aproximar operadoresdiferenciais de quarta ordem requer a introdução de espaços de elementos finitoscom continuidade C1 afim de satisfazer a conformidade, uma característica quecausa dificuldades em termos de construção –por envolver um grande número degraus de liberdade– e implementação computacional, cf. [1–3]. Isto motivou odesenvolvimento de métodos de penalização C0-interior (também conhecido comométodos de Galerkin contínuos/decontínuos) por ENGEL et al. [4] para problemaselípticos de quarta ordem, e que foram posteriormente investigados em diferentestrabalhos [1, 5–8]. Esses métodos usam somente elementos finitos de Lagrange usuaispara problemas de segunda ordem e, portanto, são mais simples que os métodosC1 conformais. Além disso, são uma alternativa atrativa aos métodos mistos queapresentam um elevado custo computacional e requerem a aproximação separada demúltiplos campos [3].

Por outro lado, formulações usuais para problemas de escoamentos incompressí-veis requerem, devido sua natureza mista (cf. [9]), a escolha de espaços de interpola-ções para o par velocidade e pressão que satisfaça a condição de inf-sup ou condiçãode Banach-Necas-Babuska [10] para garantir soluções estáveis e convergentes. Istoproíbe, por exemplo, o uso de aproximações de igual ordem para a velocidade epressão. Dificuldades com a inf-sup também ocorrem quando espaços de elemen-tos finitos de baixa ordem são usados [11]. Não obstante, como é conhecido (veja,e.g., GRESHO e SANI [12]), a violação da condição de inf-sup introduz, dentreoutras patologias, oscilações espúrias na solução da pressão e, possivelmente, umadeterioração ou travamento (locking) do campo de velocidade. Portanto, a condi-ção de inf-sup deve ser satisfeita ou contornada para evitar soluções indesejáveis doproblema discreto (cf. [10, 12–14]).

Uma abordagem eficiente para contornar as dificuldades fixadas pela condição deinf-sup consiste em usar o método de Galerkin mínimos quadrados ou GLS (cf. [15–18]), em que termos residuais das equações governantes são contabilizados, tornando

1

a equação variacional fracamente coerciva para espaços que falham ao satisfazer acondição de inf-sup [10, 19, 20]. Em particular, o uso do método GLS torna possívelescolhas arbitrárias de espaços de interpolações para a velocidade e pressão, que égeralmente desejável do ponto de vista computacional, cf. [21].

Neste trabalho, nós propomos uma formulação consistente e estabilizada paraproblemas de escoamentos incompressíveis de quarta ordem. Motivado pelas discus-sões mencionadas anteriormente, a formulação é baseada nos métodos de Galerkincontínuo/descontínuo e estabilização por mínimos quadrados. Além disso, interpo-lações por pressão descontínua são adotadas, que são fisicamente mais apropriadaspara forçar a condição de incompressibilidade de uma maneira forte a nível de ele-mento [12, 13] e, ademais, permitir uma condensação parcial dos graus de liberdade(veja, e.g., [22]). Uma análise de estabilidade usando um lema englobando a con-dição de inf-sup e pressão descontínua é também apresentado, mostrando que aformulação proposta é fracamente coerciva e, portanto, robusta. Esse lema sugereque existe um parâmetro de estabilização de mínimos quadrados ótimo ou quaseótimo, que não é necessariamente o mesmo para todos os elementos da malha.

Experimentos numéricos são fornecidos para ilustrar a eficiência da formulaçãoem lidar com qualquer combinação de elementos de velocidade-pressão, bem comosua capacidade de estabilizar altos gradientes de pressão. Além disso, os resultadosnuméricos também dão suporte à análise de estabilidade e indicam que o parâmetrode estabilização –até o momento obtido globalmente por experimentação numérica–depende do grau do polinômio usado para interpolar os campos de velocidade, dageometria dos elementos da malha e do termo de viscosidade do fluido.

Uma formulação de elementos finitos para o problema considerado aqui foi desen-volvido inicialmente em KIM et al. [23], porém envolvendo o gradiente do rotacionaldo campo de velocidade. Similar ao presente trabalho, KIM et al. [23] baseiam a for-mulação deles no método de penalização C0-interior proposto em ENGEL et al. [4].Entretanto, o tratamento deles para a restrição de incompressibilidade não leva emconta a condição de inf-sup. Na formulação deles, a pressão é estabilizada usandouma técnica de estabilização baseada em HUGHES e FRANCA [15], impondo acontinuidade do campo de pressão de maneira fraca entre os elementos e introduzum parâmetro de estabilidade global adicional. Todavia, o método tem apelo li-mitado, ao ganhar estabilidade a incompressibilidade é sacrificada, o que torna ométodo inconsistente em sua forma variacional e pode comprometer a convergên-cia do mesmo. Por essa razão, não fica evidente que eles possam evitar as bemconhecidas patologias associadas a violação da condição de inf-sup.

O presente trabalho é motivado pela necessidade de resolver problemas de es-coamento incompressível em micro e nanoescalas. Por exemplo, um entendimentodo escoamento de materiais fluidos em pequenas escalas de comprimento é de fun-

2

damental importância em uma ampla variedade de investigações, incluindo biologia[24, 25], química [26, 27], e em particular nos campos recentes da micro e nanotec-nologia [28, 29]. Experimentos recentes indicam que as equações de Navier-Stokese suas condições de contorno usuais descrevem corretamente o escoamento de líqui-dos com comportamento linear através de canais com dimensões transversais tãopequenas quanto 50µm [30, 31]. Simulações atomísticas sustentam tais conclusões,e indicam que a teoria convencional tem falhado para diâmetros menores [32–34].

Uma teoria de escoamento de fluido incompressível envolvendo um segundo gra-diente no campo de velocidade foi desenvolvida em FRIED e GURTIN [35]. Essateoria foi sugerida como um meio para descrever comportamentos em escalas decomprimento suficientemente pequenas para que o desvio da teoria convencionalNavier-Stokes seja concebível [23]. Além da possibilidade de descrever escoamentosem dispositivos em microescalas e problemas de nanolubrificação, a teoria tambémpermite descrever escoamentos em superfícies rugosas. Uma variação dessa mesmateoria envolve o gradiente do rotacional do campo de velocidade, e foi introduzidoem [36] ao construir uma formulação generalizada que fornece um alicerce contínuo-mecânico alternativo para o modelo de turbulência baseado na média LagrangeanaNavier-Stokes-α (LANS-α), dando origem ao modelo de turbulência Navier-Stokes-αβ (cf. [37, 38]). O modelo LANS-α confere algumas vantagens sobre os mode-los de turbulência LES e RANS, incluindo a preservação de estruturas coerentesnão-lineares [23] (veja HOLM et al. [39] para uma sinopse). Ambos os trabalhospropostos em [35] e [36] dão origem à equações de escoamentos idênticas sob certashipóteses constitutivas, e envolvem derivadas de quarta ordem no campo de veloci-dade. Além disso, a teoria também fornece condições de contorno associadas comformulações de alta ordem.

Por outro lado, dada a complexidade dessas equações, apenas soluções analíticaspara problemas de escoamento simples unidimensionais –análogos aos problemasde escoamento de Poiseuille–, estão disponíveis e, portanto, métodos alternativasde solução são claramente desejáveis. O desenvolvimento de um método numéricorobusto possibilita explorar as potencialidades e capacidade de predição de escoa-mentos complexos em geometrias microfluídicas que essas equações de escoamentode alta ordem ordem podem fornecer.

Este trabalho é composto dos seguintes capítulos.O Capítulo 2 introduz o problema modelo de quarta ordem com restrição interna

de incompressibilidade, que tem sua origem em uma teoria de segundo gradientepara escoamento de fluido incompressível. Em seguida é introduzido o problemavariacional associado.

No Capítulo 3 revisaremos algumas técnicas numéricas usuais dentro da estru-tura dos métodos dos elementos finitos convencionais para solução de equações de

3

quarta ordem, considerando as principais dificuldades em se discretizar essas equa-ções diferenciais. O problema de valor de contorno é então aproximado pelo métodode Galerkin contínuo/descontínuo.

Dentro do Capítulo 4 é considerada a estabilização do problema aproximadopor meio do método de Galerkin mínimos quadrados combinado com interpolaçãodescontínua da pressão. Uma discussão inicial no que concerne as oscilações numé-ricas que podem ocorrer sempre que processos incompressíveis são considerados eda necessidade de algum tipo de estabilidade.

O Capítulo 5 apresenta uma análise de estabilidade da nova formulação a partirde um lema estendido a problemas de quarta ordem com restrição interna, queengloba a condição de inf-sup e indica a existência de uma parâmetro de mínimosquadrados ótimo para cada elemento da malha.

No Capítulo 6, alguns exemplos numéricos dão suporte a análise de estabilidadee ilustram o bom desempenho e capacidade de aproximação da formulação proposta.

Finalmente no Capítulo 7, considerações finais e sugestões para trabalhos futurossão apresentadas.

4

Capítulo 2

Problema de quarta ordem comrestrição de incompressibilidade

2.1 Problema modelo

Seja um domínio Ω ⊂ Rn (n = 2 ou 3) com contorno Γ Lipschitz contínuo. Seja ΓD

e ΓN tais queΓ = ΓD ∪ ΓN e meas(ΓD ∩ ΓN) = 0, (2.1)

onde meas( · ) denota medida positiva de Lebesgue.Seja η e µ constantes reais e considere a função f ∈ L2(Ω)n (L2(Ω)n =

s = (s1, . . . , sn); si ∈ L2(Ω)), g0 ∈ H32 (ΓD)n ∩ C0(ΓD)n (H v

2 (ΓD)n = g =(g1, . . . , gn); gi ∈ H

v2 (ΓD) (v = 1, 2, 3...)), g1 ∈ H

12 (ΓD)n, h0 ∈ L2(ΓN)n e h1 ∈

L2(ΓN)n, onde os espaços são espaços de Sobolev como definidos em [40].Nós consideramos ainda os seguintes espaços de SobolevHv(Ω) eHv(Ω)n = w =

(w1, . . . , wn);wi ∈ Hv(Ω) (v = 1, 2, 3 . . .) todos definidos conforme a referência [40].Nosso problema modelo de valor de contorno consiste em encontrar o par (u, p) ∈

H2(Ω)n × L2(Ω) se η > 0 ou o par (u, p) ∈ H1(Ω)n × L2(Ω) se η = 0 e satisfazendo

η∆(∆u) − µ∆u + ∇p = f em Ω, (2.2)

∇ · u = 0 em Ω, (2.3)

u = g0 em ΓD, (2.4)∂u∂n

= (∇u)n = g1 em ΓD se η > 0, (2.5)

(µ∂u∂n

− pn)

− η∂(∆u)∂n

= h0 em ΓN, (2.6)

η∆u = h1 em ΓN se η > 0, (2.7)

5

onde ∆u denota o Laplaciano do vetor u, ∇·u denota o divergente de u, ∇u denotao gradiente de u, n denota o vetor unitário normal ao contorno Γ orientado saindode Ω e definido quase em toda parte de Γ, e “v · w” denota o usual produto escalarentre os vetores v e w do Rn:

v · w =n∑

i=1viwi (2.8)

Note que aqui e em todo o texto nós usamos notação direta para simplificar nossoscálculos. Contudo, para maior clareza, também apresentamos definições básicas eresultados na forma de componente.

Observe que para η = 0 o problema de valor de contorno dado em (2.2) a (3.15)é reduzido ao problema de escoamento de Stokes convencional. Do contrário, paraη = µℓ2 > 0, ele é análogo as equações governantes para segundo-gradiente e fluidosincompressível, cf. [35, 36]. Nesse caso, é introduzido um comprimento de escalaℓ =

√η/µ característico do material, o que abre caminho para contabilizar os efeitos

em pequena escala [35].Problemas de valor de contorno formulados em termos de teorias generalizadas

de alta ordem requerem condições de contorno adicionais, e determiná-las é umatarefa delicada e não trivial. Embora essas condições sejam de simples tratamentona definição variacional, o significado físico, aplicabilidade e valores no contornopara um dado problema nem sempre são evidentes [41]. Sem perda de generalidade,nós consideramos as condições de contorno essenciais (3.4) e (3.13) e naturais (3.14)e (3.15) prescritas no contorno, onde h0 e h1 denotam trações no contorno.

Veja [35, 36] para uma justificativa compreensiva e discussão concernente àsequações governantes (2.2) e condições de contorno formuladas com base em umaestrutura para teorias contínuo-mecânicas envolvendo dependências de gradientes.Em [12] é apresentada uma exaustiva discussão de condições de contorno e seusefeitos em simulações numéricas de escoamentos incompressíveis.

2.2 Formulação variacional associada

Para definir o problema variacional associado ao problema de valor de contorno dadona seção 2.1, nós precisamos definir os conjuntos solução e os espaços das variaçõesadmissíveis.

Como as condições associadas com g0 e g1 são essenciais, nós definimos o conjuntosolução para o campo de velocidade u como segue:

Su =

u ∈ H2(Ω)n; u = g0 e ∂u

∂n= g1 em ΓD se η > 0

u ∈ H1(Ω)n; u = g0 em ΓD se η = 0(2.9)

6

e definimos o espaço correspondente das variações admissíveis como segue:

Vu =

v ∈ H2(Ω)n; v = 0 e ∂v

∂n= 0 em ΓD se η > 0

v ∈ H1(Ω)n; v = 0 em ΓD se η = 0(2.10)

O conjunto solução para a pressão p e o espaço correspondente das variaçõesadmissíveis são dados como segue:

Sp =

L2(Ω) e meas(ΓN) > 0

q ∈ L2(Ω);∫

Ω q dΩ = 0 se meas(ΓN) = 0(2.11)

Vp = Sp. (2.12)

Para se obter a formulação variacional associada nós definimos a forma multili-near

A∗(u, p,v, q) =∫

Ω(η∆u · ∆v + µ∇u : ∇v) dΩ

+∫

Ωq (∇ · u) dΩ −

∫Ωp (∇ · v) dΩ,

(2.13)

e o funcional linear

l(v) =

∫

Ω f · v dΓ +∫

ΓNh0 · v dΓ +

∫ΓN

h1 · ∂v∂ndΓ se η > 0

∫Ω f · v dΓ +

∫ΓN

h0 · v dΓ se η = 0(2.14)

O problema variacional associado consiste em encontrar o par (u, p) ∈ Su × Sp

satisfazendo a seguinte equação variacional

A∗(u, p,v, q) = l(v) ∀(v, q) ∈ Vu × Vp. (2.15)

O desafio em construir aproximações de elementos finitos para o problema devalor de contorno dado na Seção 2.1 é encontrar uma formulação consistente nos es-paços de dimensões finitas com uma solução aproximada estável tão próxima quantopossível da solução exata correspondente em dimensão infinita. No presente traba-lho buscamos essa meta para o problema de quarta ordem com restrição interna deincompressibilidade, e consiste em incorporar adequadamente as derivadas de altaordem do campo de velocidade e ao mesmo tempo estabilizar, de modo consistente,o campo de pressão.

7

Capítulo 3

Aproximação por elementos finitos

O método de elementos finitos (MEF) é um método numérico para resolver proble-mas complexos em engenharia e ciências aplicadas. Esse método tem proporcionadovantagens quando aplicado a geometrias complexas, modelos mais complicados, eo comportamento da solução é incerto [11] (cf., por exemplo, [10]). A robustez econfiabilidade do MEF é um resultado direto do seu sólido fundamento matemáticoque conecta (via sua arquitetura matemática) os princípios físicos fundamentais doproblema contínuo às equações discretas que o MEF produz para aproximação doproblema.

3.1 Métodos de elementos finitos para problemasde quarta ordem

O problema de valor de contorno dado na Seção 2.1 pode ser aproximado por méto-dos usuais de elementos finitos. O material que apresentamos a seguir é fortementebaseado em ENGEL et al. [4], BRENNER e GUDI [8] e BRENNER e SUNG [1] ereferências citadas.

Dentro da estrutura dos métodos usuais de elementos finitos, esse problema podeser resolvido numericamente por métodos de elementos finitos C1 conformais. Es-ses métodos são chamados conformais porque o espaço de elementos finitos paraaproximar o problema discreto é um subespaço do espaço onde o problema contínuoestá posto. Dada uma partição de Ω em elementos quadriláteros ou triangulares,os espaços de elementos finitos são construídos juntando funções de base polinomi-ais contínuas por partes de modo que eles sejam globalmente C1. Tais construçõesrequerem, portanto, elementos C1. Em particular, métodos conformais sempre fun-cionam e produzem a melhor aproximação dos espaços de elementos finitos. Con-tudo, elementos finitos C1 são bastantes complicados de construir, em particularquando interfaces entre diferentes elementos são consideradas, o que envolve um

8

grande numero de graus de liberdade (veja, por exemplo, [1–3, 42, 43] e referênciascitadas).

Em problemas unidimensionais, o requerimento de continuidade C1 é trivial,pois cada elemento tem um único ponto nodal em comum. Se considerarmos, porexemplo, os elementos (ou intervalos) da figura 3.1, podemos notar que a continui-dade da função no ponto nodal “•” e a continuidade da derivada no ponto nodal“⃝” impõem quatro condições (ou graus de liberdade) em cada elemento. Estascondições podem ser combinadas usando polinômios cúbicos, que são elementos deHermite com quatro graus de liberdade. E estes polinômios são determinados porseus valores nodais e os valores de suas derivadas em cada ponto. Pela especificaçãodesses valores podemos juntar polinômios cúbicos contínuos por partes para formarfunções globalmente C1.

Figura 3.1: Elementos cúbicos de Hermite P3 (4 graus de liberdade).

Casos bidimensionais envolvem muitos elementos vizinhos e a continuidade C1

é mais complicada de ser satisfeita. Existem muitas condições fixadas nos vérticese arestas de um elemento, nos quais a continuidade C1 tem de ser garantida. Con-sequentemente são necessários muitos graus de liberdade para satisfazer todas essascondições. Elementos especiais desenvolvidos para aproximar problemas de placase cascas, por exemplo, foram adaptados com algum sucesso para problemas comdependências de gradientes. Esses elementos são construídos com arestas retas, compontos nodais nos vértices e usam como graus de liberdade a variável primal (i.e.,deslocamento) e suas derivadas de primeira e segunda ordem [3].

Por exemplo, nos elementos triangulares Argyris (TUBA) [43] e elementos re-tangulares Bogner-Fox-Schmit [42], os polinômios são determinados pelos valores desuas derivadas até ordem dois nos vértices e os valores de suas derivadas normaisnos pontos médios e derivadas segunda mista nos pontos nodais. Esses elementossão mostrados nas figuras 3.2 e 3.3, onde o ponto sólido “•” denota a função deinterpolação avaliada em cada ponto do elemento, o circulo pequeno “⃝” e o cir-culo grande “⃝” denotam as derivadas de primeira e segunda ordem das funçõesde forma avaliadas em cada ponto, a seta “↑” denota a derivada normal avaliada noponto médio, e “” a derivada de segunda ordem mista avaliada em cada ponto [1].Através da especificação desse valores, nós podemos então juntar esses polinômiosde alta ordem contínuos por partes para formar funções que são globalmente C1.

Por outro lado, elementos isoparamétricos C1 têm sua construção dificultadapela necessidade de mapear –do sistema paramétrico ao sistema de coordenadasglobais– graus de liberdade que são derivadas espaciais da variável primaria. Além

9

Figura 3.2: Elemento triangular P5 Argyris (21 graus de liberdade); adaptado deBRENNER e SUNG [1].

Figura 3.3: Elemento retangular Q3 Bogner-Fox-Schmit (16 graus de liberdade);adaptado de BRENNER e SUNG [1].

disso, estender os elementos unidimensionais cúbicos de Hermite para duas dimen-sões não é trivial. Elementos de Hermite retangulares isoparamétricos C1 para duasdimensões foram desenvolvidos em [44] e posteriormente estendidos para três dimen-sões em [41], obtendo um elemento isoparamétrico hexaédrico de Hermite. Devidoà formulação isoparamétrica, esses elementos embora sejam de fácil construção emmalhas estruturadas, eles apresentam dificuldades em malhas não-estruturadas ecertas partições não são possíveis.

Em resumo, métodos C1 conformais são as opções mais naturais de aproximaçãode problemas de quarta ordem e sempre funcionam, porém, a construção dessesespaços de elementos finitos é bastante complicada e mais ainda em três dimensões.

Alternativamente, o problema de valor de contorno dado na Seção 2.1 pode serresolvido numericamente por elementos finitos não-conformais, cf., por exemplo,[45–47]. Contudo, a construção de espaços de elementos finitos não-conformais con-vergentes é uma tarefa delicada, especialmente para problemas de quarta ordem comrestrição interna. Além disso, elementos não-conformais são de baixa ordem –nãoapresentam nenhuma hierarquia natural–, os quais não são muito eficientes parasoluções suaves. E pouco se sabe a respeito de elementos não-conformais tridimen-sionais para problemas de quarta ordem.

Uma outra alternativa é usar uma formulação de elementos finitos misto, cf. [9],em que o problema de quarta ordem é dividido e aproximações separadas para oscampos primários e secundários são requeridas. A relação entre os campos é esta-belecida através do uso de multiplicadores de Lagrange, que constitui um campo

10

adicional a ser discretizado. Nessas formulações são usados apenas espaços de ele-mentos finitos que são subespaços de H1(Ω) e, portanto, elementos C0 de Lagrange.Contudo, construir uma formulação de elementos finitos misto estável e convergenteé uma tarefa delicada e altamente não-trivial para cada problema de quarta ordem[8], e mais ainda quando apresenta restrição interna de incompressibilidade. Em for-mulações mistas são usados diferentes espaços de elementos finitos para aproximaros campos primários e secundários. O método é funcional somente se a combinaçãodesses espaços satisfaz a condição de inf-sup para escolhas particulares desses espa-ços discretos. E encontrar tais espaços não é uma tarefa simples. Ao final, resolvero problema de ponto de sela resultando é ainda mais mais complicado que resolverum problema simétrico positivo definido [8].

Um método de penalização C0-interior que preserva a forma simétrica positivadefinida de problemas elípticos de quarta ordem e que usa ao mesmo tempo apenaselementos finitos de Lagrange usuais para problemas de segunda ordem foi propostoem ENGEL et al. [4]. Esse método permite contornar os principais inconvenientesdos métodos convencionais de elementos finitos para problemas de quarta ordem.Métodos de penalização C0-interior pertencem à família dos métodos de Galerkindescontínuos. E os espaços de elementos finitos usados são subespaços de H1(Ω),ou seja, espaços das funções contínuas C0.

Ao contrário dos métodos de elementos finitos mistos e não-conformais, o desen-volvimento de métodos de penalização C0-interior quase ótimo é simples mesmo paraproblemas de quarta ordem complexos, e usa apenas integração por partes, simetri-zação e penalização [8]. Além disso, possui ainda as seguintes vantagens em relaçãoa outros métodos de elementos finitos (como investigado em [1, 4, 8, 48]): (i) é maissimples que os métodos de elementos finitos conformais C1. (ii) A complexidadedos métodos de penalização C0-interior de baixa ordem é comparável a dos méto-dos não-conformais usuais, porém vêm de uma hierarquia natural. (iii) Ao mesmotempo os métodos de penalização C0-interior de alta ordem, que são apropriadospara capturar soluções suaves, são mais simples que os métodos de elementos finitosC1. (iv) O fato de os espaços de elementos finitos usados em aproximações por pena-lização C0-interior serem subespaços de H1 pode ser explorado no desenvolvimentode algoritmos multigrid e métodos de decomposição de domínios [48].

O método de penalização C0-interior de ENGEL et al. [4] foi posteriormente in-vestigado em [1] e empregado com sucesso na solução de diversos problemas de valorde contorno de quarta ordem de natureza elíptica e parabólica, incluindo equação deCahn-Hilliard para separação de fases [6], problemas de dano com dependência degradientes [5, 49], e na discretização de teorias de segundo gradiente para escoamentode fluido incompressível [23, 37]. Mais recente, um estimador de erro a posterioribaseado em resíduo para um método de penalização C0-interior quadrático para um

11

problema biharmônico foi proposto em BRENNER e GUDI [8], fornecendo novasferramentas para análise de outros métodos de elementos finitos para problemas dequarta ordem.

O estudo de métodos de penalização interior não é atual e remonta aos anos de1970. Métodos de penalização interior para problemas biharmônicos foram estu-dados em BABUSKA e ZLÁMAL [50] e BAKER [51]. Métodos com penalizaçãointerior para equações biharmônicos usando funções polinomiais descontínuas foraminvestigados em SULI e MOZOLEVSKI [52]. Apesar da capacidade de refinamentohp que os métodos deles apresentam, possuem o inconveniente de dobrar os graus deliberdade no contorno entre elementos [7] e, portanto, tornam-se caros computacio-nalmente. Em comparação, o método de penalização C0-interior proposto em [4], eque é usado no presente trabalho para aproximar o problema de valor de contornodado na Seção 2.1, envolve menos graus de liberdade devido os valores nodais seremcompartilhados nos contornos dos elementos interiores ao domínio. Isto é possívelpela inexistência de hanging nodes nas formulações [1].

Para uma profunda discussão do método de penalização C0-interior ou Galerkincontínuo/descontínuo e seu desempenho numérico, o leitor é reportado ao trabalhode ENGEL et al. [4]. No trabalho deles é apresentada uma revisão das técnicasconvencionais de elementos finitos tais como métodos de Galerkin contínuo, Galerkindescontínuo e esquemas de estabilização, e que fornece as bases para explorar ométodo de Galerkin contínuo/descontínuo para resolver problemas de quarta ordem.WELLS et al. [5] e BRENNER e SUNG [1] proveram investigações e análises de errodo método de ENGEL et al. [4] para problemas elípticos de quarta ordem.

3.2 Método de penalização C0-interior para o pro-blema incompressível de quarta ordem

Nesta seção nós adaptamos o método de penalização C0-interior ou método de Galer-kin contínuo/descontínuo proposto em ENGEL et al. [4] para aproximar o problemade valor de contorno dado na Seção 2.1 e, assim, contornar a necessidade de usarfunções de interpolação com continuidade C1 e evitar formulações de elementos fini-tos mistos. A idéia básica por trás do método é o uso de funções de interpolação C0

de Lagrange, e impor de maneira fraca a continuidade das derivadas primeira e dealta ordem através dos contornos entre elementos usando uma variação do métodooriginalmente proposto em NITSCHE [53] (veja [4] para uma revisão).

Para construirmos as bases para a formulação, nós vamos considerar Mh =Ω1, · · · ,Ωne uma partição regular do domínio Ω em ne elementos finitos Ωe, tal queΩe pode ser mapeado em elementos convencionais por mapeamentos isoparamétricos

12

e Ωe ∩ Ωe′ = ∅ se e = e′. Ademais, Ω ∪ Γ = ∪nee=1

(Ωe − Γe

)e Γee′ = Γe ∩ Γe′ , onde

Γe denota o contorno de Ωe. Nós definimos ainda os espaços quebrados de Sobolevassociados à partição Mh como segue:

Hv,b(Mh) = φ : Ω 7−→ R ; ∀e ∈ 1, . . . , ne, φe ∈ Hv(Ωe),

Hv,b,1(Mh) = φ ∈ Hv,b ∩H10 (Ω),

onde φe é a restrição de φ ao elemento Ωe. Hv(Ωe) e Hv(Ωe)n são espaços de Sobolevde ordem v; H0(Ωe) ≡ L2(Ωe) são espaços de Hilbert. (Nós seguimos a definiçãousual dos espaços de Sobolev em [40].)

Seja k > 1 e k > l > 0 inteiros e considere Pr(Ωe) o espaço dos polinômios degrau menor que ou igual a r restrito ao elemento Ωe. Introduzindo os espaços deelementos finitos

Hh,k = φ ∈ H2,b,1(Mh);φe ∈ Pk(Ωe),

Hh,k,n = φ = (φi)1≤ i ≤ n; φi ∈ Hh,k,

Lh,l = φ ∈ H0,b(Mh);φe ∈ P l(Ωe),

nós podemos então definir os seguinte espaços discretos:

Sh,ku = uh ∈ Hh,k,n; uh = gh

0 em ΓD,

V h,ku = vh ∈ Hh,k,n; vh = 0 em ΓD,

Shp = V h

p = Sp ∩ Lh,l,

onde gh0 é o usual interpolante de g0.

Note que são contempladas funções que são contínuas em todo o domínio, masdescontínuas nas derivadas primeira e de segunda ordem através dos contornos doselementos interiores. Vale ressaltar que esses espaços de funções possuem menos re-gularidade do que seria necessário em abordagens convencionais para as equações dequarta ordem [6]. Em uma aproximação usual de Galerkin, o natural seria procurarsoluções em um subespaço de H2(Ω), ao invés de fazê-lo em subespaços de H1(Ω).Consequentemente, Sh,k

u , V h,ku são espaços de elementos finitos C0 de Lagrange. Em

particular, estamos assumindo que Shp e V h

p são espaços de elementos finitos C−1.Isto é, nós interpolamos o campo de pressão por funções descontínuas. Além disso,nós também estamos interessados em casos de interpolação de igual ordem k = l,que é atrativo do ponto de vista de construção e implementação computacional, cf.[21].

A aproximação de Galerkin contínuo/descontínuo para o problema variacional

13

dado em (2.15) é encontrar (uh, ph) ∈ Sh,ku × Sh

p tal que

A(uh, ph,vh, qh) = l(vh) ∀(vh, qh) ∈ V h,ku × V h

p (3.1)

onde

A(uh, ph,vh, qh) =ne∑

e=1

∫Ωe

µ∇uhe : ∇vh

e dΩ +∫

Ωe

µℓ2∆uhe · ∆vh

e dΩ

−∫

Ωe

phe (∇ · vh

e ) dΩ +∫

Ωe

qhe (∇ · uh

e ) dΩ

+∑e′>e

(−∫

Γee′

12(µℓ2∆uh

e + µℓ2∆uhe′

)·(∂vh

e

∂ne

+ ∂vhe′

∂ne′

)dΓ

+∫

Γee′

12(∂uh

e

∂ne

+ ∂uhe′

∂ne′

)·(µℓ2∆vh

e + µℓ2∆vhe′

)dΓ

+∫

Γee′τu

(∂uhe

∂ne

+ ∂uhe′

∂ne′

)·(∂vh

e

∂ne

+ ∂vhe′

∂ne′

)dΓ)

−∫

Γe∩ΓDµℓ2∆uh

e · ∂vhe

∂ne

dΓ +∫

Γe∩ΓD

∂uhe

∂ne

· µℓ2∆vhe dΓ

+∫

Γe∩ΓDτD∂uh

e

∂ne

· ∂vhe

∂ne

dΓ,

(3.2)

l(vh) =ne∑

e=1

∫Ωe

fe · vhe dΓ +

∫Γe∩ΓN

he,0 · vhe dΓ +

∫Γe∩ΓN

he,1 · ∂vhe

∂ne

dΓ

+∫

Γe∩ΓDge,1 ·

(µℓ2∆vh

e

)dΓ +

∫Γe∩ΓD

τD ge,1 · ∂vhe

∂ne

dΓ,

(3.3)

e τu > 0, τD > 0 são parâmetros de penalização a serem fixados.De análise dimensional, nós supomos que os parâmetros de penalização τu e τD

devam ser tais que τu = τD ∼ Cµℓ2/hee′ , onde hee′ é um comprimento característicoda aresta do elemento e C é uma constante positiva independente. Essa constantetem de ser cuidadosamente escolhida para assegurar acurácia suficiente e boa con-vergência do método, cf. [4, 7, 23, 37]. Em [4] foi sugerida uma abordagem analíticapara determinar ou limitar a constante C, porém técnicas alternativas são desejáveis[37]. Para os propósitos deste trabalho, a constante C é até o momento determinadopor experimentações numéricas. Vale ressaltar que um parâmetro de penalização devalor elevado afeta negativamente a precisão do método de penalização C0-interior[1, 37].

14

Consistência

As equações de Euler-Lagrange associadas à partição Mh do domínio Ω são dadaspor:

− ∇pe + µ∆ue − µℓ2∆∆ue + fe = 0 em Ωe, (3.4)

∇ · ue = 0 em Ωe, (3.5)

ue = g0 em Γe ∩ ΓD se meas(Γe ∩ ΓD) > 0, (3.6)

∂ue

∂ne

= (∇ue)n = g1 em Γe ∩ ΓD se η > 0 e meas(Γe ∩ ΓD) > 0, (3.7)

(µ∂ue

∂ne

− pene

)− ∂(µℓ2∆ue)

∂ne

= he,0 se meas(Γe ∩ ΓN) > 0, (3.8)

µℓ2∆ue = he,1 em ΓN, se meas(Γe ∩ ΓN) > 0, (3.9)

∂ue

∂ne

+ ∂ue′

∂ne′= 0 em Γee′ , (3.10)

[(µ∂ue

∂ne

− pene

)− ∂(µℓ2∆ue)

∂ne

]+[(µ∂ue′

∂ne′− pe′ne′

)− ∂(µℓ2∆ue′)

∂ne′

]= 0 em Γee′ , (3.11)

µℓ2∆ue − µℓ2∆ue′ = 0 em Γee′ , (3.12)

onde ue e pe são restrições ao elemento Ωe. As equações governantes são assegura-das por (3.4) e (3.5) nos elementos da partição Mh, (3.6) a (3.9) são condições decontorno essenciais e naturais, (3.10) garante a continuidade dos gradientes de velo-cidade através dos contornos dos elementos, e (3.11) e (3.12) garante a continuidadedos fluxos através dos elementos.

Para (ue, pe) ∈ H4,b(Mh)n × H1,b(Mh) e (vhe , q

he ) ∈ V h,k

u × V hp , integração por

partes conduz as seguintes expressões:

∫Ωe

µℓ2∆ue · ∆vhe dΩ =

∫Ωe

(µℓ2∆∆ue) · vhe dΩ −

∫Γe

∂(µℓ2∆ue)∂ne

· vhe dΓ

+∫

Γe

µℓ2∆ue · vhe dΓ,

(3.13)

∫Ωe

µ∇ue · ∇vhe dΩ = −

∫Ωe

(µ∆ue) · vhe dΩ +

∫Γe

µ∂ue

∂ne

· vhe dΓ, (3.14)

−∫

Ωe

pe∇ · vhe dΩ =

∫Ωe

∇pe · vhe dΩ −

∫Γe

pene · vhe dΓ. (3.15)

15

Somando (3.13), (3.14) e (3.15) em todos os elementos Ωe em Mh e usando o fatode ve ser contínuo através dos contornos dos elementos interiores ao dominínio, i.e.,ve − ve′ = 0, e se anula em ΓD, temos que

0 = A(u, p,vh, qh) − l(vh)

=ne∑

e=1

∫Ωe

(µ∆ue − µℓ2∆∆ue − ∇pe + fe

)· vh

e dΩ +∫

Ωe

qhe (∇ · ue) dΩ

+∑e′>e

(−∫

Γee′

12(∂ue

∂ne

+ ∂ue′

∂ne′

)·(µℓ2∆vh

e + µℓ2∆vhe′

)dΓ

+∫

Γee′

12(µℓ2 ∆ue − µℓ2∆ue′

)·(∂vh

e

∂ne

+ ∂vhe′

∂ne′

)dΓ

+∫

Γee′τu(∂ue

∂ne

+ ∂ue′

∂ne′

)·(∂vh

e

∂ne

+ ∂vhe′

∂ne′

)dΓ)

−∫

Γe∩ΓN

12

([(µ∂ue

∂ne

− pene

)− ∂(µℓ2∆ue)

∂ne

]+[(µ∂ue′

∂ne′− pe′ne′

)− ∂(µℓ2∆ue′)

∂ne′

])· (vh

e + vhe′) dΓ

−∫

Γe∩ΓN

([(µ∂ue

∂ne

− pene

)− ∂(µℓ2∆ue)

∂ne

]− he,0

)· vh

e dΓ

+∫

Γe∩ΓN

(µℓ2∆ue − he,1

)· ∂vh

e

∂ne

dΓ

−∫

Γe∩ΓD

(∂ue

∂ne

− ge,1)

· µℓ2∆vhe dΓ

+∫

Γe∩ΓDτD(∂ue

∂ne

− ge,1)

· ∂vhe

∂ne

dΓ.

(3.16)

Uma vez que essa equação é assegurada para todo (vh, qh) ∈ V h,ku × V h

p , a aproxi-mação (3.1) é portanto consistente no sentido que

A(u, p,vh, qh) = l(vh) ∀(vh, qh) ∈ V h,ku × V h

p . (3.17)

16

Capítulo 4

Estabilização da formulação deGalerkin contínuo/descontínuopara escoamento incompressível

Neste capítulo será discutido brevemente a conhecida necessidade de estabilização,considerando as principais razões da ocorrência de oscilações não físicas da soluçãoaproximada sempre que um processo incompressível é considerado, e as vantagens dese contornar tais patologias. A forma variacional discreta é modificada pela adiçãode termos de mínimos quadrados para obter estabilidade.

O foco é a formulação de um método de elementos finitos estabilizado baseadonos métodos de penalização C0-interior e GLS, tal que nenhuma restrição adicionalà estabilidade ocorra.

4.1 Oscilações espúrias e a necessidade de estabi-lização em problemas incompressíveis

É bem conhecido que formulações de elementos finitos com base nos métodos deGalerkin usuais para escoamentos incompressíveis podem falhar drasticamente. Istoé em parte devido à condição de incompressibilidade e o papel desempenhado pelapressão (que não é uma variável termodinâmica).

Um olhar mais atento sobre a aproximação de Galerkin contínuo/descontínuo(3.1) revela que a variável pressão atua como um multiplicador de Lagrange parafazer cumprir a condição de incompressibilidade. Isto evita a necessidade de cons-truir funções de interpolação de elementos finitos que satisfaçam de maneira exataa condição de divergência livre sobre o campo de velocidade. Uma vez que as fun-ções de interpolação para a velocidade não são divergente livre, a formulação tem,portanto, um caráter misto –ambas as variáveis u e p são aproximadas–, e requer a

17

solução de problemas de ponto de sela, cf. [9].Por outro lado, a escolha dos espaços de aproximação para velocidade e pressão

têm de satisfazer a condição de inf-sup (ou condição de Banach-Necas-Babuska)para que as soluções aproximadas uh e ph sejam estáveis e convergentes, com acu-rácia ao menos quase ótima. Satisfazer a condição de inf-sup implica em geral naestabilidade da pressão. A violação dessa condição conduz, dentre outras coisas,à não-unicidade da solução da pressão, manifestando patologias caracterizadas poroscilações espúrias (não-físicas) no campo de pressão e/ou degradação do campo develocidade (cf. GRESHO e SANI [12], DONEA e HUERTA [54]). A exigência decompatibilidade entre os espaços decorre do fato que, ao contrário dos problemasfortemente elípticos, a restrição interna de incompressibilidade causa dificuldades deaproximação porque a estabilidade dos espaços onde é definido o problema contínuonão é automaticamente transferida para o modelo discreto [13]. Desta forma, con-dições de compatibilidade têm de ser consideradas para cada escolha particular dosespaços discretos.

Da discussão acima exposto, o problema discreto (3.1) é estável se os espaçosde elementos finitos Sh,k

u , Shp , V h,k

u , V hp e, portanto, os elementos são escolhidos

tais que a versão discreta da condição de inf-sup seja satisfeita. Contudo, essacondição é violada para algumas combinações desses espaços que são desejáveis paraaproximar as variáveis discretas velocidade uh e pressão ph. Por exemplo, elementosC0 com interpolações de igual ordem definidos com relação à mesma malha ou paresde elementos com interpolações de baixa ordem. Estes são elementos amplamenteusados para problemas de escoamento complexo [11]. Poucas são as combinaçõesdos espaços de elementos finitos Sh,k

u , Shp , V h,k

u e V hp capazes de satisfazer a condição

de inf-sup dentro da formulação de Galerkin contínuo/descontínuo.Existem excelentes combinações de pares de elementos para velocidade e pres-

são que satisfazem a condição de inf-sup a priori [14, 55], tais como o elementoTaylor-Hood de velocidade quadrática e pressão linear. Esse elemento usa espaçode interpolação contínua para pressão, que implica em uma discretização global daincompressibilidade. Contudo, para um valor do comprimento característico ℓ ele-vado, o elemento de Taylor-Hood pode não garantir a mesma estabilidade da pressãocomo acontece em problemas de segunda ordem (cf. [23]). Por exemplo, os altosgradientes do campo de pressão não são estabilizados, como ilustram os resultadosnuméricos presentes no Capítulo 5. Portanto, é necessário o uso de uma estabilizaçãoadicional para a pressão.

Uma maneira eficiente de estabilização consiste na adição de termos residuaisque tornam a equação variacional coerciva, ou fracamente coerciva, para espaçosque falham ao satisfazer a condição de inf-sup [19] (cf., por exemplo, [10, 20]). Istoé, permite contornar as restrições impostas pela condição de inf-sup e possibilita

18

ampliar o conjunto de elementos possíveis para simulações de escoamentos incom-pressíveis, tal que a escolha dos espaços de interpolação para velocidade e pressãopode ser feito livremente [56]. Além disso, métodos de elementos finitos estabiliza-dos têm mostrado ser mais eficientes que o uso de interpolações mistas estáveis, cf.[21].

4.2 Estabilização pelo método de Galerkin míni-mos quadrados

Como discutido na seção anterior, é possível trabalhar sem a necessidade de usarinterpolações de elementos finitos mistos [10] desde que a aproximação de Galerkincontínuo/descontínuo (3.1) seja ligeiramente modificada.

Para aliviar a necessidade de satisfazer a condição de inf-sup, e ao mesmo tempotirar vantagem do uso de interpolações descontínuas para a pressão (cf. [21]), nósadotamos o método GLS. A técnica consiste na inclusão de termos adicionais naforma de Galerkin de maneira a reforçar sua estabilidade. Esses termos são obtidospela minimização do quadrado da norma-L2 do residual discreto dentro de cadaelemento [57], multiplicado por parâmetros de regularização (ou estabilização) ele-mental adequados que irá controlar a contribuição da parcela de mínimos quadradosna sentença variacional. O método GLS para problemas de Stokes foi originalmenteintroduzido em [15], a análise de convergência foi investigada em [18]; cf. também[10, 17, 20, 58, 59].

A aproximação de Galerkin contínuo/descontínuo (3.1) escrita no sentido dosmínimos quadrados (de agora em diante referenciado como formulação CDGLS)consiste em encontrar (uh, ph) ∈ Sh,k

u × Shp tal que

Ah(uh, ph,vh, qh) = ltot(vh, qh) ∀(vh, qh) ∈ V h,ku × V h

p , (4.1)

comAh(uh, ph,vh, qh) = A(uh, ph,vh, qh) +

ne∑e=1

AeLS(uh, ph,vh, qh), (4.2)

ltot(vh, qh) = l(vh, qh) +ne∑

e=1l eLS(vh, qh), (4.3)

onde os termos a frente são os termos de mínimos quadrados

AeLS(w, q, δw, δq) =

∫Ωe

δGLS(he)( n∑

i=1Ri(w, q)Ri(δw, δq)

)dΩ, (4.4)

Ri(w, q) = ∂p

∂xi

− µ∆wi + µℓ2∆∆wi, (4.5)

19

l eLS(vh, qh) =

∫Ωe

δGLS(he)( n∑

i=1fi Ri(vh, qh)

)dΩ, (4.6)

com (w, q) ∈ H 4,b(Mh)n ×H 1,b(Mh) e (δw, δq) ∈ H 4,b(Mh)n ×H 1,b(Mh), e δGLS(he)é o parâmetro de estabilização, que será definido aqui por [17]:

δGLS(he) = αh 2e

µ, (4.7)

onde α é um parâmetro adimensional a ser fixado e he é o comprimentos caracterís-tico do elemento usualmente dado por he =

( ∫ΩedΩ)1/n

, cf. [15].A adição desses termos fornece a estabilização necessária e ao mesmo tempo se

anulam para a solução exata, preservando a acurácia natural do método de Galerkin[20], mantendo a consistência da aproximação (4.1). Ressaltamos também que paraos espaços de elementos finitos triangulares e tetraédricos até ordem três, elementosquadriláteros quadráticos e hexaédricos quadráticos, o operador diferencial de quartaordem se anula (∆∆ui = 0) e, portanto, nenhum cálculo de derivadas de ordemmaior que dois é efetuado para esses elementos, e que são amplamente usados naprática.

A aproximação CDGLS, Eq. (4.1), adapta o método de penalização C0-interiorcom a técnica de estabilização GLS, o que faz da formulação uma alternativa atrativaem relação aos métodos convencionais de elementos finitos para problemas de esco-amento incompressível de quarta ordem. Além de usar apenas espaços de elementosfinitos usuais em problemas de segunda ordem, ela permite satisfazer a condição deinf-sup para uma ampla faixa de combinações de espaços de elementos finitos para opar velocidade e pressão, incluindo o emprego de interpolações de igual ordem paraambas as variáveis. Ademais, é capaz de estabilizar altos gradientes de pressão,como claramente é mostrado pelos resultados numéricos discutidos no Capítulo 5.

O parâmetro adimensional α em (4.7) precisa ser cuidadosamente determinadopara assegurar ao método acurácia suficiente e boa taxa de convergência [60]. Esteparâmetro possui um valor ótimo ou quase ótimo para cada tipo de problema, erequer testes numéricos para sua calibração, sendo em geral definido globalmentepara todos os elementos da malha. Contudo, experimentos numéricos revelam umefeito significativo da magnitude do parâmetro de estabilização na acurácia dosresultados discretos da pressão, cf. [61].

Um procedimento que permite determinar, para cada elemento da malha, umparâmetro α ótimo ou quase ótimo de formulações estabilizadas GLS para problemasincompressíveis de segunda ordem foi investigado em CARMO et al. [62] e [63].Contudo, estender esses resultados aos problemas de quarta ordem não é uma tarefatrivial e requer uma análise de estabilidade em uma norma mais adequada para adeterminação desses parâmetros. Isso é investigado no capítulo seguinte.

20

Capítulo 5

Análise de estabilidade doproblema de quarta ordem

Neste capítulo nós apresentamos uma análise de estabilidade do método CDGLSestabilizado quando aplicado ao problema de valor de contorno dado na Seção 2.1.A análise sugere a existência de um parâmetro de estabilização α que é dependentedo grau do polinômio usado para interpolar os campos de velocidade u e pressão p,da geometria dos elementos da malha e da viscosidade dinâmica do fluido µ. E ésustentado por experimentos numéricos. O resultado que permite obter a metodolo-gia apropriada para determinar esse parâmetro é apresentado no Lema 1, que é umaextensão do proposto em CARMO et al. [62, 63] para problemas incompressíveis desegunda ordem (problemas de Stokes e elasticidade incompressível).

Nós consideramos que os espaços de elementos finitos são construídos usandomalhas que satisfazem a condição de regularidade de malha CM1, em que a distorçãodos elementos é controlada: existem constantes reais C1,distor > C0,distor > 0 tais que

C0,distor <ρe

he

< C1,distor ∀Ωe, (5.1)

onde he = diam(Ωe) e ρe = supdiam(S), S é uma hiperesfera contida em Ωe ∪ Γe;cf. [10]. Devido este controle de distorção da malha nós encontramos c′

0 > 0 ec′

1 > 0 para inferir sobre c′0he < hee′ < c′

1he, c′0he′ < hee′ < c′

1he, c′0he < he′ < c′

1he ec′

0he′ < he < c′1he′ ; onde hee′ = min he, he′.

Note que ao longo desse texto nós usamos C (com ou sem subscrito) para denotaruma constante genérica, que pode tomar diferentes valores em diferentes aspectos.

Segue do bem conhecido teorema que diz que todas as normas em um espaço dedimensão finita são equivalentes, e da condição de regularidade de malha CM1 acima,as seguintes estimativas inversas são validas: existem constantes reais CLap,1 > 0 eCLap,2 > 0 e Ccont,1 > 0 independentes do parâmetro da malha he, dependendosomente do grau do polinômio e da distorção da malha que é controlada pelas

21

constantes C0,distor, C1,distor tais que∫

Γe

heψ2 dΩ 6 Ccont,1

∫Ωe

ψ2 dΩ ∀ψ ∈ Pk(Ωe) (k > 1), (5.2)

∫Ωe

ℓ2|∆φ|2 dΩ 6 ℓ2

(he)2CLap,1

∫Ωe

|∇φ|2 dΩ ∀φ ∈ Pk(Ωe) (k > 1), (5.3)

∫Ωe

µ2ℓ4|∆∆φ|2 dΩ 6 ℓ2

(he)2CLap,2 µ

(he)2

∫Ωe

µℓ2|∆φ|2 dΩ ∀φ ∈ Pk(Ωe) (k > 1). (5.4)

Nós assumimos para todos we ∈ H1(Ωe)n e qe ∈ L2(Ωe) a seguinte decomposiçãoa nível de elemento

we = we + we′ e qe = qe + qe, (5.5)

com as partes constantes

we =∫

Ωewe dΩ∫

ΩedΩ

e qe =∫

Ωeqe dΩ∫

ΩedΩ

. (5.6)

Nós consideramos ainda as seguintes desigualdades de Poincaré

∀we ∈ H1(Ωe)n, ∥we∥L2(Ωe)n 6 CeP oinc he|we|H1(Ωe)n , (5.7)

∀qe ∈ H1(Ωe)n, ∥qe∥L2(Ωe)n 6 CeP oinc he|qe|H1(Ωe). (5.8)

onde | · |H1(Ωe)n e | · |H1(Ωe) denotam as seminormas de H1(Ωe)n e H1(Ωe), respecti-vamente, e a constante real positiva Ce

P oinc independe do parâmetro de malha he edepende somente da razão de aspecto usual da malha de elementos finitos e que écontrolada (condição CM1).

O resultado que segue é fortemente inspirado em CARMO et al. [62, 63] e repre-senta uma extensão do Lemma (4.38) da referência [10]. Além de considerar as duasdiferenças fundamentais presentes em [62, 63], isto é: (i) o parâmetro de estabilizaçãodepende do elemento Ωe e não é necessariamente o mesmo para todos os elementosda malha, (ii) a análise é derivada com uma norma diferente e mais apropriada paraa determinação do parâmetro de estabilização, nós estendemos estas duas diferençasaos problemas de quarta ordem com restrição interna de incompressibilidade.

Para formalizar essa idéia, para 0 < β < 1, α > αmin > 0 ∀Ωe e k > 1 (sendo

22

k um inteiro), nós introduzimos a seguinte norma em V h,ku × V h

p :

∥|(vh, qh)|∥2β,h,Ω =

ne∑e=1

β

(1 + β)|qh|2αe,he,µ,H1(Ωe) + (1 − β)

(|vh|2µ,H1(Ωe)n

+ |vh|2µ,ℓ,Lap +∑e′>e

∫Γee′

τv

∣∣∣(∂vhe

∂ne

+ ∂vhe′

∂ne′

)∣∣∣2 dΓ+

n∑i=1

∫Ωe

αe (he)2

µ|Ri(vh, qh)|2 dΩ

)+∥∥∥ qh

e

µ1/2

∥∥∥2

L2(Ωe)

,

(5.9)

|vh|2µ,H1(Ωe)n =∫

Ωe

( n∑i=1

µ|∇vhi |2)dΩ, (5.10)

|vh|2µ,ℓ,Lap =∫

Ωe

( n∑i=1

µℓ2|∆vhi |2)dΩ, (5.11)

|qh|2αe,he,µ,H1(Ωe) =∫

Ωe

αe (he)2

µ|∇qh|2 dΩ, (5.12)

que foi estendida aqui para os problemas de quarta ordem.Desde modo nós resumimos a discussão acima no seguinte lema.

Lema 1. Sob pressupostos que (MC1) e (5.2) a (5.4), se o parâmetro αe é tal que

αe > αmin > 0 ∀Ωe, (5.13)

2∫

Ωe

αe (he)2

µ(µ∆ui)2 dΩ 6 β

∫Ωe

µ |∇ui|2 dΩ ∀Ωe e ∀i, (5.14)

2∫

Ωe

αe (he)2

µ(µℓ2∆∆ui)2 dΩ 6 β

∫Ωe

µℓ2 |∆ui|2 dΩ ∀Ωe e ∀i, (5.15)

então existe uma constante CInfSup > 0 dependente de (ℓ/he) tal que

inf(uh,ph)∈V h,k

u ×V hp

sup(vh,qh)∈V h,k

u ×V hp

Ah(uh, ph,vh, qh)∥|(uh, ph)|∥h,Ω∥|(vh, qh)|∥β,h,Ω

> CInfSup. (5.16)

Pode ser observado da condição (5.16) e da demonstração do Lema 1 que aformulação CDGLS (4.1) é fracamente coerciva com relação à norma (5.9) para oproblema discreto. Ademais, os resultados no Lema 1 indicam um modo de obterum procedimento para determinar ou limitar parâmetros de GLS ótimos ou quaseótimos.

Nota 1. Para a demonstração do Lema 1 nós assumimos interpolação por pressãodescontínua. Formulações com pressão descontínua pode impor a condição de in-compressibilidade de maneira forte a nível de elemento. Contudo, pode-se provar,sem dificuldades adicionais, a estabilidade para pressão contínua, cf. [63].

Demonstração. Seja (uh, ph) ∈ V h,ku × V h

p . A demonstração procede em três etapas.

23

A seguinte desigualdade elementar será usada na demonstração:

± ab 6 12

(γa2 + 1γb2) ∀a, b ∈ R and ∀γ ∈ R (γ > 0). (5.17)

Seja λ definido como segue

λ = max(

ℓ

he

)2; e = 1, · · · , ne

. (5.18)

1) Um cálculo direto conduz a

Ah(uh, ph,uh, ph) =ne∑

e=1

(1 − β)

(|uh|2µ,H1(Ωe)n

+ |uh|2µ,ℓ,Lap +∑e′>e

∫Γee′

τu

∣∣∣(∂uhe

∂ne

+ ∂uhe′

∂ne′

)∣∣∣2 dΓ

+n∑

i=1

∫Ωe

αe (he)2

µ|Ri(uh, ph)|2 dΩ

)+ β

(|uh|2µ,H1(Ωe)n

+ |uh|2µ,ℓ,Lap +∑e′>e

∫Γee′

τu

∣∣∣(∂uhe

∂ne

+ ∂uhe′

∂ne′

)∣∣∣2 dΓ+

n∑i=1

∫Ωe

αe (he)2

µ|Ri(uh, ph)|2 dΩ

).

(5.19)

Usando a desigualdade (5.17) com γ = β na identidade seguinte

β

(1 + β)

∫Ωe

αe(he)2

µ

(∂ph

∂xi

)2

= β

(1 + β)

∫Ωe

αe(he)2

µ|Ri(uh, ph) + µ∆uh

i − µℓ2∆∆uhi |2 dΩ,

(5.20)

obtemos

β

(1 + β)

∫Ωe

αe(he)2

µ

(∂ph

∂xi

)26 β

(1 + β)

( ∫Ωe

αe(he)2

µ(1 + β) |Ri(uh, ph)|2 dΩ

+∫

Ωe

2 αe(he)2

µ

(1 + 1

β

)(|µ∆uh

i |2 + |µℓ2 ∆∆uhi |2) dΩ

).

(5.21)

Resulta de (5.10), (5.12), (5.14), (5.15) e (5.21) que

β

(1 + β)|ph|2αe,he,µ,H1(Ωe)

6 β(

|uh|2µ,H1(Ωe)n + |uh|2µ,ℓ,Lap

+∑e′>e

∫Γee′

τu

∣∣∣(∂uhe

∂ne

+ ∂uhe′

∂ne′

)∣∣∣2 dΓ +∫

Ωe

αe(he)2

µ|Ri(uh, ph)|2

).

(5.22)

24

Combinando (5.9) e (5.19) com (5.22), encontramos

Ah(uh, ph,uh, ph) > ∥|(uh, ph)|∥2β,h,Ω −

ne∑e=1

(∥∥∥ phe

µ1/2

∥∥∥2

L2(Ωe)

). (5.23)

2) k > 1 (Shp ⊂ L2(Ω) e V h

p ⊂ L2(Ω))Do resultado apresentando em FORTIN [64] (cf. também [65]) e Lema (4.19)

dado na referência [10], nós temos que que para k > 1 existe uma constante realpositiva C∗,F ortin tal que, para todo ph ∈ V h

p , existe vh ∈ V h,ku satisfazendo

ne∑e=1

−(∇ · vhe , p

he )L2(Ωe) =

ne∑e=1

∥∥∥ phe

µ1/2

∥∥∥2

L2(Ωe), (5.24)

∥vh∥2H1(Ωe)n 6 C∗,F ortin

( ne∑e=1

∥∥∥ phe

µ1/2

∥∥∥2

L2(Ωe)

) 12. (5.25)

Ademais, usando as estimativas inversas (5.2) e (5.3) juntamente com condição deregularidade de malha CM1, a desigualdade (5.17) com um adequado valor de γ eo resultado de Fortin (5.25), após algumas manipulações algébricas, obtemos

ne∑e=1

( ∫Ωe

µℓ2(∆v)2 dΩ +∫

Ωe

µ|∇v|2 dΩ +∑e′>e

∫Γee′

τu

(∂ve

∂ne

+ ∂ve′

∂ne′

)2dΓ)

6 C∗∗ne∑

e=1

(∥ve∥2

H1(Ωe)n +∑e′>e

∥ve′∥2H1(Ωe)n

)

6 C∗∗ (Nface + 1)( ne∑

e=1∥ve∥2

H1(Ωe)n

)

6 C∗∗ (Nface + 1)C∗,F ortin( ne∑

e=1

∥∥∥ pe

µ1/2

∥∥∥2

L2(Ωe)

)

= C2( ne∑

e=1

∥∥∥ pe

µ1/2

∥∥∥2

L2(Ωe)

)

(5.26)

onde Nface(Ωe) denota o número de faces do elemento Ωe.Um cálculo simples conduz à relação

Ah(uh, ph,vh, 0) = A′(uh,vh) +ne∑

e=1

(− (∇ · vh

e , phe + ph

e )L2(Ωe)

+n∑

i=1

∫Ωe

αe (he)2

µRi(uh, ph)(µ∆vh

i ) dΩe

),

(5.27)

25

A′(uh,vh) =ne∑

e=1

n∑1=1

( ∫Ωe

µ∇uhe,i · ∇vh

e,i dΩe +∫

Ωe

µℓ2∆uhe,i∆vh

e,i dΩe

)

+∑e′>e

(−∫

Γee′

12

(µℓ2∆uhe + µℓ2∆uh

e′) ·(∂vh

e

∂ne

+ ∂vhe′

∂ne′

)dΓ

+∫

Γee′

12(∂uh

e

∂ne

+ ∂uhe′

∂ne′

)· (µℓ2∆vh

e + µℓ2∆vhe′) dΓ

+∫

Γee′τu

(∂uhe

∂ne

+ ∂uhe′

∂ne′

)·(∂vh

e

∂ne

+ ∂vhe′

∂ne′

)dΓ).

(5.28)

Uma vez quene∑

e=1∥∇ · vh

e ∥2L2(Ωe) 6 n

ne∑e=1

∥v∥2H1(Ωe)n , (5.29)

pela desigualdade de Cauchy-Schwartz em L2(Ωe), (5.8), (5.13), (5.24) e (5.25) epela desigualdade de Cauchy-Schwartz Rne, temos

ne∑e=1

−(∇ · vhe , p

he + ph

e )L2(Ωe) >( ne∑

e=1

∥∥∥ phe

µ1/2

∥∥∥2

L2(Ωe)

)

− C2

( ne∑e=1

∥∥∥ phe

µ1/2

∥∥∥2

L2(Ωe)

) 12( ne∑

e=1

β

1 + β|ph|2αe,he,µ,H1(Ωe)

) 12,

(5.30)

C2 = C∗,F ortin sup(

nµ (1 + β)β αe

) 12Ce

P oinc; e = 1, · · · , ne. (5.31)

Similarmente, pela desigualdade de Cauchy-Schwartz em L2(Ωe), (5.14), (5.25) epela desigualdade de Cauchy-Schwartz inequality em Rne, temos

ne∑e=1

n∑i=1

∫Ωe

αe(he)2

µRi(uh, ph)(µ∆vh

i )

> −(µ

β

1 − β

) 12C∗,F ortin

( ne∑e=1

∥∥∥ phe

µ1/2

∥∥∥2

L2(Ωe)

) 12

×( ne∑

e=1(1 − β)

n∑i=1

∫Ωe

αe(he)2

µ|Ri(uh, ph)|2

) 12.

(5.32)

Novamente, usando a desigualdade de Cauchy-Schwartz em L2(Ωe), as relações dedesigualdades definidas em (MC1), (5.2) e pela desigualdade de Cauchy-Schwartzem Rne, nós podemos determinar a constante Cc

1 > 0 tal que

A′(uh,vh) > −Cc1

ne∑e=1

( n∑1=1

∫Ωe

(µ(∇uhe,i)2 dΩe + µℓ2(∆uh

e,i)2) dΩe

+∑e′>e

∫Γee′

τu

(∂uhe

∂ne

+ ∂uhe′

∂ne′

)2dΓ) 1

2 ne∑

e=1

( n∑1=1

∫Ωe

(µ(∇vhe,i)2

+ µℓ2(∆vhe,i)2) dΩe +

∑e′>e

∫Γee′

τu

(∂vhe

∂ne

+ ∂vhe′

∂ne′

)2dΓ) 1

2

(5.33)

26

Segue do resultado (5.26) que

A′(uh,vh) > −Cc1 C

((1 − β)

ne∑e=1

( n∑1=1

∫Ωe

(µ (∇uhe,i)2 + µℓ2(∆uh

e,i)2) dΩe

+∑e′>e

∫Γee′

τu

(∂uhe

∂ne

+ ∂uhe′

∂ne′

)2dΓ)) 1

2( 1

(1 − β)

ne∑e=1

∥∥∥ phe

µ1/2

∥∥∥2

L2(Ωe)

) 12.

(5.34)

De (5.27), (5.30)-(5.34) e pela desigualdade Cauchy-Schwartz em R3, nós inferimosque

Ah(uh, ph,vh, 0) >( ne∑

e=1

∥∥∥ phe

µ1/2

∥∥∥2

L2(Ωe)

)− C3

( ne∑e=1

∥∥∥ phe

µ1/2

∥∥∥2

L2(Ωe)

) 12

×( ne∑

e=1(1 − β)

(|uh|2µ,H1(Ωe) + |uh|2µ,ℓ,Lap +

n∑i=1

∫Ωe

αe(he)2

µ|Ri(uh, ph)|2 dΩ

+∑e′>e

∫Γee′

τu

(∂uhe

∂ne

+ ∂uhe′

∂ne′

)2dΓ)

+(

β

1 − β

)|ph|2αe,he,µ,H1(Ωe)

) 12,

(5.35)

C3 = (3)12 sup

(µ

β

1 − β

) 12C∗,F ortin, Cc

1 C( 1

1 − β

) 12, C2

. (5.36)

Combinando a desigualdade elementar (5.17) com γ = 1 e (5.35) encontramos

Ah(uh, ph,vh, 0) > 12

( ne∑e=1

∥∥∥ phe

µ1/2

∥∥∥2

L2(Ωe)

)− C4

( ne∑e=1

(1 − β)(

|uh|2µ,H1(Ωe)

+ |uh|2µ,ℓ,Lap +n∑

i=1

∫Ωe

αe(he)2

µ|Ri(uh, ph)|2 dΩ

+∑e′>e

∫Γee′

τu

(∂uhe

∂ne

+ ∂uhe′

∂ne′

)2dΓ)

+(

β

1 − β

)|ph|2αe,he,µ,H1(Ωe)

),

(5.37)

onde C4 = (C3)12/2. Segue da norma (5.9) que

Ah(uh, ph,uh, ph) > 12

ne∑e=1

(∥∥∥ phe

µ1/2

∥∥∥2

L2(Ωe)

)

− C4

(∥|(uh, ph)|∥2

β,h,Ω −ne∑

e=1

(∥∥∥ phe

µ1/2

∥∥∥2

L2(Ωe)

)).

(5.38)

3) Escolhendoθ = 1

1 + C4 + 14, (5.39)

e fixandowh = (1 − θ)uh + θvh e qh = (1 − θ)ph, (5.40)

27

e combinando (5.23) com (5.38) (k > 1), encontramos

Ah(uh, ph,uh, ph) > −((1 − θ) − θC4)(

∥|(uh, ph)|∥2β,h,Ω −

( ne∑e=1

∥∥∥ phe

µ1/2

∥∥∥2

L2(Ωe)

))

+ θ

2

( ne∑e=1

∥∥∥ phe

µ1/2

∥∥∥2

L2(Ωe)

),

(5.41)

e portanto segue que

Ah(uh, ph,uh, ph) > θ

4∥|(uh, ph)|∥2

β,h,Ω. (5.42)

Uma vez que(wh, ph) = (1 − θ)(uh, ph) + θ(vh, 0), (5.43)

temos pela desigualdade triangular, (5.9), (5.14) e (5.25) que

∥|(wh, qh)|∥β,h,Ω 6 (1 − θ)∥|(uh, ph)|∥β,h,Ω + θ∥|(vh, 0)|∥β,h,Ω

6 ((1 − β2)µ)12 C5 ∥|(uh, ph)|∥β,h,Ω,

(5.44)

C5 = C∗,F ortin (1 + λ(ℓ,Mh)). (5.45)

FixandoCInfSup = θ

4((1 − β2)µ) 14 C5

, (5.46)

finalmente conduz a

Ah(uh, ph,wh, qh) > CInfSup ∥|(uh, ph)|∥β,h,Ω∥|(wh, qh)|∥β,h,Ω, (5.47)

e o resultado segue imediatamente, completando a demonstração.

28

Capítulo 6

Resultados numéricos

Neste capítulo apresentamos simulações numéricas que ilustram o desempenho daformulação CDGLS.

6.1 Caso 1: Escoamento plano de Poiseuille

Neste caso teste, é considerado de um fluido Newtoniano com viscosidade µ, e escoa-mento laminar permanente entre duas placas paralelas separadas por uma distânciad. Este é um dos tipos mais comuns de escoamento observado em canais longos eestreitos (e.g., dispositivos de microfuídica), cf. [28]. Uma solução analítica paraeste problema foi desenvolvida em [35] para o escoamento de um fluido com segundogradiente e condições de contorno generalizadas.

Na solução, o campo de pressão é conhecido a menos de uma constante aditivaarbitrária com gradiente

∇p = −Pex (6.1)

e P = constant (P > 0); sem perda de generalidade, nós assumimos que a pressãodecresce com o crescimento de x. Ao passo que a solução da velocidade u = u(y)ex

é dada por

u(y) = Pd 2

2µ yd

(1 − y

d

)− blℓ

d sinh(d/ℓ)(

sinh dℓ

− sinh yℓ

− sinh d− y

ℓ

), (6.2)

onde a constante positiva

bl = (2ℓ/d) + (l/ℓ)1 + (l/ℓ) tanh(d/2ℓ)

(6.3)

é um número adimensional não-negativa, e l é um comprimento de scala de aderência;cf. [35].

Como na referência KIM et al. [23], usamos esta solução exata que é essen-

29

cialmente unidimensional, para construir um problema bidimensional para verifi-car os resultados obtidos com a formulação CDGLS estabilizada. Esse problemade escoamento foi resolvido considerando um domínio que é o canal retangularΩ = [0, 10d] × [0, d] com d = 0, 5 mm, prescrevendo as condições de contornou = u(y)ex e (∇u)n = 0 sobre o contorno de entrada x = 0, 0<y<d, satisfazendoa velocidade exata (6.2), e condições de contorno u livre, v = 0, (∇u)n = 0 na regiãode saída do escoamento x = 10d, 0<y<d. Condições de contorno de parede fixau = 0 e (∇u)n = g1 são especificadas em ambos os contornos y = 0, 0<x< 10de y = d, 0<x<10d. A figura 6.1 descreve o domínio e as condições de contorno.

Figura 6.1: Domínio do escoamento em canal plano e condições de contorno.

Inicialmente, examinamos a solução discreta e sua sensibilidade ao parâmetrode estabilização δGLS experimentando valores diferentes para constante α. Para isto,nós fixamos o comprimento de escala ℓ igual a d/4 e a condição de aderência fraca(l = 0) é considerada para a teoria de gradiente. A malha espacial é construída com100 elementos.

As Figuras 6.2 a 6.4 mostram soluções do campo de pressão para elementos finitosde Lagrange Q2/Q1 com pressão descontínua. Note que o elemento finito Q2/Q1 cominterpolação descontínua da pressão não satisfaz a condição de inf-sup e, portanto,instável para o método de Galerkin (cf., por exemplo, [10]). A estabilização pelométodo GLS alivia as restrições impostas pela condição de inf-sup e permite o usoadequadamente desse elemento para aproximar a velocidade e pressão. A figura 6.5mostra os perfis de velocidade através da linha central x/10d = 0.5, juntamente comas soluções exatas e solução de referência para o problema de segunda ordem (cf.[35]), e são similares aqueles obtidos por KIM et al. [23]. Note que a estabilidadeé alcançada para diferentes valores de α. Contudo, podemos observar que para umvalor de α pequeno, a pressão exibe pequenas pertubações nos cantos do domínio.

30

Figura 6.2: Escoamento em canal: campo de pressão estabilizado obtido pela for-mulação CDGLS, α=1.5, elementos finitos de Lagrange Q2/Q1 com pressão descon-tínua.

Figura 6.3: Escoamento em canal: campo de pressão estabilizado obtido pela for-mulação CDGLS, α=0,05, elementos finitos com Lagrange Q2/Q1 com pressão des-contínua.

31

Figura 6.4: Escoamento em canal: campo de pressão estabilizado obtido pela for-mulação CDGLS, α=0,005, elementos finitos de Lagrange Q2/Q1 com pressão des-contínua.

0 0.05 0.1 0.15 0.2 0.250

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

velocidade adimensional 2µu/Pd2

dist

anci

a ad

imen

sion

al y

/d

CDGLS, α=1.5

CDGLS, α=0.5

CDGLS, α=0.05

CDGLS, α=0.005sol. exatasol. classica

Figura 6.5: Escoamento em canal: perfis de velocidade exatos e numéricos através dalinha de centro x/10d = 0, 5 para diferentes valores do parâmetro de estabilizaçãoα, comparados à solução do problema de Poiseuille clássico. Todos os resultadosnuméricos foram obtidos com uma malha construída com 100 elementos LagrangeQ2/Q1 com pressão descontínua.

A figura 6.6 compara as soluções exatas e numéricas ao longo da linha x/10d =0, 5, para o caso de condições de aderência generalizada usando três razões diferentesde comprimento de aderência pelo comprimento de escala do material, l/ℓ. Os

32

resultados indicam uma correspondência satisfatória entre os campos de velocidadenumérico e o exato. Como discutido no trabalho [35], é esperado que os perfis develocidade em escala de comprimento suficientemente pequenas sejam menores queaqueles da teoria Navier-Stokes convencional para este tipo de escoamento.

0 0.05 0.1 0.15 0.2 0.250

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

velocidade adimensional 2µu/Pd2

dist

anci

a ad

imen

sion

al y

/d

0.11.03.0

classico

exata

numerica

Figura 6.6: Escoamento em canal: perfis de velocidade exatos e numéricos através dex/10d = 0, 5 para a razão l/ℓ iguais a 0, 1, 1, 0 e 3, 0. Todos os resultados numéricosforam obtidos com uma malha construída com 100 elementos finitos de LagrangeQ2/Q1 com pressão descontínua.

As figuras 6.7 e 6.8 mostram o campo de pressão obtido com elementos biqua-dráticos de Lagrange Q2/Q2 com pressão contínua e descontínua, respectivamente.E claramente descrevem os efeitos adversos da violação da condição de inf-sup. Afigura 6.9 mostra que a pressão é estabilizada, contudo, apresenta pequenas pertu-bações na solução discreta. Isto pode sugerir que existe um valor de α apropriadoque não é necessariamente o mesmo para cada elemento da malha. No momento estevalor é escolhido globalmente para todos os elementos. Estes resultados confirmama importância da estabilização da pressão pelo método GLS em nossa formulação.

33

Figura 6.7: Escoamento em canal: campo de pressão oscilante obtidos pela for-mulação CDGLS, α=0,000015, e elementos finitos de Lagrange Q2/Q2 com pressãocontínua.

Figura 6.8: Escoamento em canal: campo de pressão oscilante obtidos pela for-mulação CDGLS, α=0,000015, e elementos finitos de Lagrange Q2/Q2 com pressãodescontínua.

34

Figura 6.9: Escoamento em canal: campo de pressão estabilizado obtido pela for-mulação CDGLS, α=0,05, e elementos finitos de Lagrange Q2/Q2 com pressão des-contínua.

6.2 Caso 2: Escoamento em cavidade

Neste segundo caso, nós resolvemos o problema do escoamento em cavidade em umdomínio quadrado Ω = [0, d]2 com d = 1mm em que as condições de parede deparede u = 0 e (∇u)n = g1 são fixadas em todo o domínio exceto no contornode topo y = d, 0<x< d, que move com velocidade u = Udex com Ud constante;(∇u)n = 0. A figura 6.10 mostra o domínio e as condições de contorno. Nósfixamos o comprimento de escala material ℓ igual a 1, 5d, mesma ordem de grandezado comprimento de escala geométrico d. Um valor de pressão nodal com média nulaé fixado para evitar algum modo constante e, uma vez que escoamentos tipo Stokesgovernam fenômenos viscosos, o número de Reynolds foi igual 2 × 10−3.

Inicialmente, consideramos uma malha espacial feita de 20 × 20 elementos. Mos-trados nas Figuras 6.11 e 6.12 estão os campos de pressão com oscilações obtidaspelo método CDGLS usando elementos de Lagrange Q2/Q1 com pressão contínuae descontínua, respectivamente. As figuras 6.13 e 6.14 mostram que os altos gradi-entes de pressão podem ser gradualmente estabilizados usando valores apropriadosde α. É importante observar que para um α suficientemente pequeno a pressão époluída por oscilações. Para valores de α altos existe o risco de uma estabilizaçãoexcessiva pelo amortecimento dos resultados do campo de pressão, mudando a físicado problema e pode não capturar corretamente a pressão nos cantos do domínio.(Lembramos que a pressão é singular nos cantos superiores do domínio para esteproblema.) Estas observações também são evidenciadas em [61]. A distribuição

35

Figura 6.10: Domínio do escoamento em cavidade e condições de contorno.

de pressão ao longo da linha horizontal y = 0, 5 mm é mostrada na figura 6.15 eacompanhada por aproximações razoáveis do campo de velocidade como mostradopelos perfis de velocidade através das linhas x/d = 0, 5 e y/d = 0, 5 das figuras 6.16e 6.17, comparadas com a solução de referência clássica para esse problema.

Esses resultados indicam claramente que a acurácia da solução do campo depressão é sensível ao parâmetro de estabilização.

0

0.5

1x 10

−3

0

0.5

1x 10−3

−1.5

−1

−0.5

0

0.5

1

1.5

x−coordenaday−coordenada

Pre

ssao

Figura 6.11: Escoamento em cavidade: campo de pressão oscilante obtido comelementos de Lagrange Q2/Q1 com pressão contínua, e parâmetro α=0,0005 para ométodo CDGLS.

36

0

0.5

1x 10

−3

0

0.5

1x 10−3

−1.5

−1

−0.5

0

0.5

1

1.5

x−coordenaday−coordenada

Pre

ssao

Figura 6.12: Escoamento em cavidade: campo de pressão oscilante obtido comelementos de Lagrange Q2/Q1 com pressão descontínua, e parâmetro α=0,0005 parao método CDGLS.

0

0.5

1x 10

−3

0

0.5

1x 10−3

−1.5

−1

−0.5

0

0.5

1

1.5

x−coordenaday−coordenada

Pre

ssao

Figura 6.13: Escoamento em cavidade: campo de pressão estabilizado obtido comelementos de Lagrange Q2/Q1 com pressão descontínua, e parâmetro α=0,01 parao método CDGLS.

37

0

0.5

1x 10

−3

0

0.5

1x 10−3

−1.5

−1

−0.5

0

0.5

1

1.5

x−coordenaday−coordenada

Pre

ssao

Figura 6.14: Escoamento em cavidade: campo de pressão estabilizado obtido comelementos de Lagrange Q2/Q1 com pressão descontínua, e parâmetro α=0,1 para ométodo CDGLS.

0 0.2 0.4 0.6 0.8 1

x 10−3

−0.08

−0.06

−0.04

−0.02

0

0.02

0.04

x−coordenada

Pre

ssao

CDGLS, α=0.0005CDGLS, α=0.01CDGLS, α=0.1

Figura 6.15: Escoamento em cavidade: distribuição de pressão (descontínua) aolongo da linha y = 0, 5 mm, malha contendo 20 × 20 elementos.

38

−0.4 −0.2 0 0.2 0.4 0.6 0.8 10

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

velocidade adimensional u/Ud

altu

ra a

dim

ensi

onal

y/d

CDGLS, α=0.0005CDGLS, α=0.01CDGLS, α=0.1sol. classica

Figura 6.16: Escoamento em cavidade: perfis de velocidade normalizados ao longoda linha x/d=0,5, malha contendo 20 × 20 elementos.

0 0.2 0.4 0.6 0.8 1−0.4

−0.3

−0.2

−0.1

0

0.1

0.2

0.3

velo

cida

de a

dim

ensi

onal

v/U

d

largura adimensional x/d

CDGLS, α=0.0005CDGLS, α=0.01CDGLS, α=0.1sol. classica

Figura 6.17: Escoamento em cavidade: perfis de velocidade normalizados ao longoda linha y/d=0,5; malha contendo 20 × 20 elementos.

Nós concluímos as simulações considerando aproximações de igual ordem para oscampos de velocidade e pressão. Para tal, consideramos uma malha contendo 10×10elementos. Os campos de pressão são mostrados nas Figuras 6.18 e 6.19 para os ele-mentos biquadráticos de Lagrange Q2/Q2 com pressões contínuas e descontínuasrespectivamente. O campo de pressão obtido com elementos Q2/Q2 é poluído porfortes oscilações espúrias. Esses efeitos adversos são causados pela violação da con-

39

dição de inf-sup, cf. [10]. A estabilidade do campo de pressão discreto pode seralcançada gradualmente usando uma parâmetro α apropriado como mostrado nasfiguras de 6.20 a 6.22. Novamente podemos observar que para valores elevados de αas soluções discretas da pressão são bastante difusivas devido à estabilização exces-siva dos altos gradientes de pressão. Isso indica que um parâmetro de estabilizaçãoótimo ou quase ótimo tem de ser identificado quando usando o método de estabiliza-ção GLS. A figura 6.23 mostram a distribuição da pressão ao longo da linha y = 0, 5mm, sem qualquer perda significativa na aproximação do campo de velocidade comomostrado pelos perfis de velocidade através da linha central da cavidade x/d = 0, 5e y/d = 0, 5 das figuras 6.24 e 6.25, respectivamente.

Estes resultados numéricos reforçam a importância da estabilização pelo métodoGLS na formulação proposta, e confirmam ainda a necessidade de se determinar umparâmetro de estabilização ótimo ou quase ótimo. O lema 1 apresenta, portanto,um caminho para a determinação do parâmetro de mínimos quadrados ótimo a nívelde elemento.

40

0

0.5

1

0

0.5

1

−0.2

−0.1

0

0.1

0.2

0.3

x−coordenaday−coordenada

Pre

ssao

Figura 6.18: Escoamento em cavidade: campo de pressão oscilante obtido comelementos de Lagrange Q2/Q2 com pressão contínua, e parâmetro α=0,00025 parao método CDGLS; malha contendo 10 × 10 elementos.

0

0.5

1x 10

−3

0

0.5

1x 10−3

−0.2

−0.1

0

0.1

0.2

0.3

x−coordenaday−coordenada

Pre

ssao

Figura 6.19: Escoamento em cavidade: campo de pressão oscilante obtido comelementos de Lagrange Q2/Q2 com pressão descontínua, e parâmetro α=0,00025para o método CDGLS; malha contendo 10 × 10 elementos.

41

0

0.5

1x 10

−3

0

0.5

1x 10−3

−0.2

−0.1

0

0.1

0.2

0.3

x−coordenaday−coordenada

Pre

ssao

Figura 6.20: Escoamento em cavidade: campo de pressão estabilizado obtido comelementos de Lagrange Q2/Q2 com pressão descontínua, e parâmetro α=0,005 parao método CDGLS; malha contendo 10 × 10 elementos.

0

0.5

1x 10

−3

0

0.5

1x 10−3

−0.2

−0.1

0

0.1

0.2

0.3

x−coordenaday−coordenada

Pre

ssao

Figura 6.21: Escoamento em cavidade: campo de pressão estabilizado obtido comelementos de Lagrange Q2/Q2 com pressão descontínua, e parâmetro α=0,035 parao método CDGLS; malha contendo 10 × 10 elementos.

42

0

0.5

1x 10

−3

0

0.5

1x 10−3

−0.2

−0.1

0

0.1

0.2

0.3

x−coordenaday−coordenada

Pre

ssao

Figura 6.22: Escoamento em cavidade: campo de pressão estabilizado obtido comelementos de Lagrange Q2/Q2, e parâmetro de estabilização α=0,1 para o métodoCDGLS; malha contendo 10 × 10 elementos.

0 0.2 0.4 0.6 0.8 1

x 10−3

−0.02

−0.015

−0.01

−0.005

0

0.005

0.01

x−coordenada

Pre

ssao

CDGLS, α=0.00025CDGLS, α=0.005CDGLS, α=0.035CDGLS, α=0.1

Figura 6.23: Escoamento em cavidade: distribuição da pressão (descontínua) aolongo da linha central y = 0, 5 mm; malha contento 10 × 10 elementos.

43

−0.4 −0.2 0 0.2 0.4 0.6 0.8 10

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

velocidade adimensional u/Ud

altu

ra a

dim

ensi

onal

y/d

CDGLS, α=0.00025CDGLS, α=0.005CDGLS, α=0.035CDGLS, α=0.1sol. classica

Figura 6.24: Escoamento em cavidade: perfil de velocidade normalizados ao longoda linha x/d = 0, 5; malha contendo 10 × 10 elementos.

0 0.2 0.4 0.6 0.8 1−0.4

−0.3

−0.2

−0.1

0

0.1

0.2

0.3

largura adimensional x/d

velo

cida

de a

dim

ensi

onal

v/U

d

CDGLS, α=0.00025CDGLS, α=0.005CDGLS, α=0.035CDGLS, α=0.1sol. classica

Figura 6.25: Escoamento em cavidade: perfil de velocidade normalizados ao longoda linha y/d = 0, 5; malha contendo 10 × 10 elementos.

44

Capítulo 7

Conclusões e sugestões detrabalhos futuros

Neste trabalho nós obtemos uma formulação de elementos finitos estabilizada e va-riacionalmente consistente para solução de problemas de escoamento incompressívelde quarta ordem. O método proposto é baseado no método de Galerkin contí-nuo/descontínuo e estabilização por mínimos quadrados. Em particular, a formu-lação é fracamente coerciva para os espaços de elementos finitos que falham aosatisfazer as restrições fixadas pela condição de inf-sup, e usa espaços de elementosfinitos C0 de Lagrange, em vez de elementos com continuidade C1 para aproximaro campo de velocidade. A continuidade das derivadas de alta ordem da velocidadeé imposta de maneira fraca entre os contornos dos elementos interiores, preservandoao mesmo tempo a consistência e mantendo a estabilidade. Além disso, nós adota-mos interpolações por pressão descontínua através dos elementos, de maneira que aincompressibilidade é localizada a nível de elemento.

Usando o método proposto, nós consideramos como caso teste o escoamento dePoiseuille plano com a teoria de segundo gradiente. Foi obtida uma concordânciasatisfatória entre os resultados numéricos e a solução analítica. Um segundo casoteste consistiu da simulação do escoamento em cavidade em escala de comprimentosuficientemente pequena. Os resultados numéricos mostram que a formulação pro-posta é capaz de estabilizar altos gradientes de pressão, e indicam que os espaçosde elementos finitos para pressão e velocidade podem ser escolhidos livremente. Poroutro lado, existe uma sensibilidade dos resultados à escolha do parâmetro de mí-nimos quadrados, e a acurácia do cálculo da pressão é estabelecida de acordo comessa escolha.

Uma metodologia que torna possível determinar parâmetros de estabilidade óti-mos ou quase ótimos do método GLS para problemas de escoamentos de Stokes eelasticidade incompressível foi investigado em CARMO et al. [62, 63]. Entretanto,como ficou evidenciado neste trabalho, estender esses resultados para problemas de

45

quarta ordem com restrição interna é uma tarefa delicada. Os resultados da aná-lise de inf-sup presentes no Lema 1 sugere que existe um parâmetro de mínimosquadrados ótimo ou quase ótimo, que não é necessariamente o mesmo para cadaelemento da malha. O lema estabelece ainda que esse parâmetro depende do graudo polinômio usado para interpolar os campos de velocidade e pressão, da geometriados elementos da malha e do termo de viscosidade do fluido.

Com base nesses resultados, sugerimos como trabalho futuro, estender o procedi-mento proposto em [62, 63] para identificar ou limitar parâmetros ótimos ou quaseótimos, a nível de elemento, do método GLS aplicado a problemas de quarta or-dem com restrição interna de incompressibilidade. Além disso, é desejável tambémobter uma análise de estabilidade cuja constante de inf-sup independa do parâme-tro de malha (ℓ/he). Uma possível abrangência do resultado de FORTIN [64, 65]para o problema de quarta ordem poderia ser uma condição de eliminação dessadependência. Mais testes numéricos do método CDGLS devem explorar suas taxasde convergências e estimativas de erros para problemas modelo mais sofisticados.Como a acurácia da solução aproximada e taxas de convergência dependem forte-mente do parâmetro de penalização global, um esquema alternativo para identificarou limitar esse parâmetro é claramente desejável.

46

Referências Bibliográficas

[1] BRENNER, S. C., SUNG, L.-Y. “C0 Interior Penalty Methods for Fourth Or-der Elliptic Boundary Value Problems on Polygonal Domains”, J. Sci.Comput., v. 22-23, pp. 83–118, 2005.

[2] DOUGLAS, J. J., DUPONT, T., PERCELL, P., et al. “A family of C1 fi-nite elements with optimal approximations properties for various Galerkinmethod for 2nd and 4th order problems”, R.A.I.R.O. Modél. Math. Anal.Numér., v. 13, pp. 227–255, 1979.

[3] PAPANICOLOPULOS, S.-A., ZERVOS, A., VARDOULAKIS, I. “Discretizationof Gradient Elasticity Problems Using C1 Finite Elements”. In: Maugin,G. A., Metrikine, A. V. (Eds.), Mechanics of Generalized Continua, v. 21,Advances in Mechanics and Mathematics, Springer New York, pp. 269–277, 2010.

[4] ENGEL, G., GARIKIPATI, K., HUGHES, T. J. R., et al. “Continu-ous/discontinuous finite element approximations of fourth-order ellipticproblems in structural and continuum mechanics with applications to thinbeams and plates, and strain gradient elasticity”, Journal of Elasticity,v. 191, n. 34, pp. 3669–3750, 2002.

[5] WELLS, G. N., GARIKIPATI, K., MOLARI, L. “A discontinuous Galerkinformulation for a strain gradient-dependent damage model”, Compu-ter Methods in Applied Mechanics and Engineering, v. 193, n. 33–35,pp. 3633–3645, 2004.

[6] WELLS, G. N., KUHL, E., GARIKIPATI, K. “A discontinuous Galerkin methodfor the Cahn-Hilliard equation”, Journal of Computational Physics, v. 218,n. 2, pp. 860–877, 2006.

[7] WELLS, G. N., DUNG, N. T. “A C0 discontinuous Galerkin formulation for Kir-chhoff plates”, Computer Methods in Applied Mechanics and Engineering,v. 196, n. 35–36, pp. 3370–3380, 2007.

47

[8] BRENNER, S. C., GUDI, T. “An a posteriori error estimator for a quadra-tic C0-interior penalty method for the biharmonic problem”, Journal ofNumerical Analysis, v. 30, pp. 777–798, 2010.

[9] BREZZI, F., FORTIN, M. Mixed and hybrid finite element methods. 1 ed. NewYork, Springer, 1991.

[10] GUERMOND, J. L., ERN, A. A theory and pratictice of finite elements. 1 ed.New York, Springer-Verlag, 2004.

[11] LAYTON, W. Introduction to the Numerical Analysis of Incompressible Vis-cous Flows. Philadelphia, PA, USA, Society for Industrial and AppliedMathematics, 2008.

[12] GRESHO, P. M., SANI, R. L. Incompressible Flow and the Finite ElementMethod, Volume 2, Isothermal Laminar Flow. 1 ed. New York, Wiley,2000.

[13] KARAM, J., LOULA, A. F. D. “A non-standard application of Babuska-Brezzitheory to finite element analysis of Stokes problem”, Mat. Aplic. Comp.,v. 10, pp. 243–262, 1991.

[14] BREZZI, F., FALK, R. “Stability of Higher-Order Hood-Taylor Methods”,SIAM Journal on Numerical Analysis, v. 28, n. 3, pp. 581–590, 1991.

[15] HUGHES, T. J. R., FRANCA, L. P. “A new finite element formulation forcomputational fluid dynamics: VII. The Stokes problem with various well-posed boundary conditions: symmetric formulations that converge for allvelocity/pressure spaces”, Computer Methods in Applied Mechanics andEngineering, v. 65, pp. 85–96, 1987.

[16] HUGHES, T. J. R., FRANCA, L. P., BALESTRA, M. A. “A new finite elementformulation for computational fluid dynamics: VIII. The Galerkin/leastsquare methods for advective-diffusive equations”, Computer Methods inApplied Mechanics and Engineering, v. 73, pp. 173–189, 1989.

[17] FRANCA, L. P., FREY, S. L. “Stabilized finite element methods: II. Theincompressible Navier-Stokes equations”, Computer Methods in AppliedMechanics and Engineering, v. 99, n. 2-3, pp. 209–233, 1992.

[18] FRANCA, L. P., HUGHES, T. J. R. “Convergence analyses of Galerkin least-squares methods for symmetric advective-diffusive forms of the Stokes andincompressible Navier-Stokes”, Computer Methods in Applied Mechanicsand Engineering, v. 105, pp. 285–298, 1993.

48

[19] BOCHEV, P. B., GUNZBURGER, M. D., LEHOUCQ, R. B. “On stabilizedfinite element methods for the Stokes problem in the small time steplimit”, International Journal for Numerical Methods in Fluids, v. 53, n. 4,pp. 573–597, 2007.

[20] BARTH, T., BOCHEV, P., GUNZBURGER, M., et al. “A Taxonomy of Consis-tently Stabilized Finite Element Methods for the Stokes Problem”, SIAMJournal on Numerical Analysis, v. 25, pp. 1585–1607, May 2004.

[21] NORBURN, S., SILVESTER, D. “Stabilised vs. stable mixed methods forincompressible flow”, Computer Methods in Applied Mechanics and Engi-neering, v. 166, n. 1–2, pp. 131–141, 1998.

[22] ELMAN, H. C., SILVESTER, D. J., WATHEN, A. J. Finite Elements and FastIterative Solvers: with Applications in Incompressible Fluid Dynamics.Numerical Mathematics and Scientific Computation. Oxford UniversityPress, 2005.

[23] KIM, T. Y., DOLBOW, J., FRIED, E. “A numerical method for a second-gradient theory of incompressible fluid flow”, Journal of ComputationalPhysics, v. 223, n. 2, pp. 551–570, 2007.

[24] BEEBE, D. J., MENSING, G. A., WALKER, G. M. “Physics and applicationsof microfluidics in biology”, Annual Review of Biomedical Engineering,v. 4, n. 1, pp. 261–286, 2002.

[25] MINC, N., VIOVY, J.-L. “Microfluidique et applications biologiques: enjeuxet tendances”, Comptes Rendus Physique, v. 5, n. 5, pp. 565–575, 2004.

[26] GARDENIERS, H., DEN BERG, A. V. “Micro- and nanofluidic devices forenvironmental and biomedical applications”, International Journal of En-vironmental Analytical Chemistry, v. 84, pp. 809–819, 2004.

[27] JENSEN, K. “Chemical kinetics: Smaller, faster chemistry”, Nature, v. 393,pp. 735–737, 1998.

[28] KIRBY, B. J. Micro- and nanoscale fluid mechanics: transport in microfluidicdevices. 2 ed. New York, Cambridge University Press, 2010.

[29] STONE, H. A., STROOCK, A. D., AJDARI, A. “Engineering flows in smalldevices: Microfluidics towards a lab-on-a-chip”, Annual Review of FluidMechanics, v. 36, pp. 381–411, 2004.

49

[30] HETSRONI, G., MOSYAK, A., POGREBNYAK, E., et al. “Fluid flow inmicro-channels”, International Journal of Heat and Mass Transfer, v. 48,pp. 1982–1998, 2005.

[31] SHARP, K. V., ADRIAN, R. A. “Transition from turbulent to laminar flow inliquid filled mircotubes”, Experiments in Fluids, v. 36, pp. 741–747, 2004.

[32] BITSANIS, I., SOMERS, S. A., DAVIS, H. T., et al. “Microscopic dynamicsof flow in molecularly narrow pores”, Journal of Chemical Physics, v. 93,pp. 3427–3431, 1990.

[33] KOPLIK, J., BANAVAR, J. R. “Continuum deductions from molecular hy-drodynamics”, Annual Review of Fluid Mechanics, v. 27, pp. 257–292,1995.

[34] OKAMURA, H., HEYES, D. M. “Comparisons between molecular dynamicsand hydrodynamics treatment of nonstationary thermal processes in aliquid”, Physical Review E, v. 70, 2004.

[35] FRIED, E., GURTIN, M. “Tractions, Balances, and Boundary Conditions forNonsimple Materials with Application to Liquid Flow at Small-LengthScales”, Archive for Rational Mechanics and Analysis, v. 182, pp. 513–554, 2006.

[36] FRIED, E., GURTIN, M. “A continuum mechanical theory for turbulence: agereralized Navier-Stokes-α equation with boundary conditions”, Theore-tical and Computational Fluid Dynamics, v. 22, pp. 433–470, 2008.

[37] KIM, T.-Y., CASSIANI, M., ALBERTSON, J. D., et al. “Impact of the inherentseparation of scales in the Navier–Stokes-αβ equations”, Phys. Rev. E,v. 79, n. 4, pp. 045307, Apr 2009.

[38] FRIED, E., GURTIN, M. E. “Turbulent kinetic-energy and a possible hierarchyof length scales in a generalization of the Navier-Stokes-α theory”, Phys.Rev. E, v. 82, n. 2, 2010.

[39] HOLM, D. D., JEFFREY, C., KURIEN, S., et al. “The LANS-α model forcomputing turbulence: origins, results, and open problems”, Los AlamosSci., v. 29, pp. 152–171, 2005.

[40] ADAMS, R., FOURNIER, J. Sobolev spaces. New York, Academic Press, 2003.

[41] PAPANICOLOPULOS, S.-A., , ZERVOS, A., et al. “A three-dimensional C1finite element for gradient elasticity”, International Journal for NumericalMethods in Engineering, v. 77, n. 10, pp. 1396–1415, 2009.

50

[42] BOGNER, F. K., FOX, R. L., SCHMIT, L. A. “The generation of interele-ment compatible stiffness and mass matrices by the use of interpolationformula”. In: Proceedings Conference on Matrix Methods in StructuralMechanics, pp. 397–444, Wright Patterson A.F.B., Dayton, OH, 1965.

[43] ARGYRIS, J. H., FRIED, I., SCHARPF, D. W. “The TUBA family of plateelements for the matrix displacement method”, J. Royal Aeronaut. Soc.,v. 72, pp. 701–709, 1968.

[44] PETERA, J., PITTMAN, J. “Isoparametric hermite elements”, InternationalJournal for Numerical Methods in Engineering, v. 37, n. 20, pp. 3489–3519, 1994.

[45] MORLEY, L. S. D. “The triangular equilibrium problem in the solution ofplate bending problems”, Aerosp. Q., v. 19, pp. 149–169, 1968.

[46] SHI, Z. “On the convergence of the incomplete biquadratic nonconforming plateelement”, Math. Numer. Sinica, v. 8, pp. 53–62, 1979.

[47] MING, W., XU, J. “The Morley element for fourth order elliptic equations inany dimensions”, Numerische Mathematik, v. 103, pp. 155–169, 2006.

[48] BRENNER, S. C., SUNG, L.-Y. “Multigrid algorithms for C0 interor penaltymethods”, SIAM J. Numer. Anal., v. 44, pp. 199–223, 2006.

[49] MOLARI, L., WELLS, G. N., GARIKIPATI, K., et al. “A discontinuous Ga-lerkin method for strain gradient-dependent damage: Study of interpo-lations and convergence”, Computer Methods in Applied Mechanics andEngineering, v. 195, n. 13–16, pp. 1480–1498, 2006.

[50] BABUSKA, I., ZLÁMAL, M. “Nonconforming elements in the finite elementmethod with penalty”, SIAM J. Numer. Anal., v. 10, pp. 863–875, 1973.

[51] BAKER, G. “Finite element methods for elliptic equations using noncoformingelements”, Math. Comp., v. 31, pp. 45–89, 1977.

[52] SULI, E., MOZOLEVSKI, I. “hp-version interior penalty DGFEMs for bihar-monic equation”, Comp. Methods Appl. Mech. Eng., v. 196, pp. 1851–1863, 2007.

[53] NITSCHE, J. “Über ein Variationsprinzip zur Lösung von Dirichlet-Problemenbei Verwendung von Teilräumen, die keinen Randbedingungen unterwor-fen sind”. In: Abhandlungen aus dem Mathematischen Seminar der Uni-versität Hamburg, pp. 9–15, Hamburg, 1970.

51

[54] DONEA, J., HUERTA, A. Finite elements methods for flow problems. NewYork, John Wiley & Sons, 2003.

[55] ARNOLD, D. N., BREZZI, F., FORTIN, M. “A stable finite element for theStokes equations”, Calcolo, v. 21, n. 4, pp. 337–344, 1984.

[56] FÖRSTER, C. Robust methods for fluid-structure interaction with stabilisedfinite elements. Tese de Doutorado, Institut für Baustatik und Baudyna-mik, Unirsittat Stuttgart.

[57] BLASCO, J. “An anisotropic GLS-stabilized finite element method for incom-pressible flow problems”, Computer Methods in Applied Mechanics andEngineering, v. 197, n. 45–48, pp. 3712–3723, 2008.

[58] HUGHES, T. J. R., FRANCA, L. P., BALESTRA, M. A. “A new finite ele-ment formulation for computational fluid dynamics: V. Circumventingthe Babuska–Brezzi condition: a stable Petrov–Galerkin formulation ofthe stokes problem accommodating equal-order interpolations”, Compu-ter Methods in Applied Mechanics and Engineering, v. 59, pp. 85–99, 1986.

[59] FRANCA, L. P., STEMBERG, R. “Error analysis of some Galerkin least squa-res methods for the elasticity equations”, SIAM Journal on NumericalAnalysis, v. 28, pp. 1680–1697, 1991.

[60] ZIENKIEWICZ, O. C., TAYLOR, R. L., ZHU, J. Z. The Finite ElementMethod: Its Basis and Fundamentals. 6 ed. London, Elsevier Butterworth-Heinemann, 2005.

[61] XIA, K., YAO, H. “A Galerkin/least-square finite element formulation fornearly incompressible elasticity/stokes flow”, Applied Mathematical Mo-delling, v. 31, n. 3, pp. 513–529, 2007.

[62] CARMO, E. G. D., SANTOS, J. P. L., MANSUR, W. J. “Uma determinaçãosistemática do parâmetro de calibração do método Galerkin Least Square(GLS)”. In: Proceedings of 30th Congresso Ibero-Latino-Americano de Mé-todos Computacionais em Engenharia, pp. 237–253, Armação dos Búzios,Brazil, ago. 2009.

[63] CARMO, E. G. D., SANTOS, J. P. L., MANSUR, W. J. “A systematic deter-mination at element level of the parameter of stabilization of the methodGLS for Stokes flow and nearly incompressible linear elasticity”, Numeri-cal Methods for Partial Differential Equations, submitted.

52

[64] FORTIN, M. “An analysis of the convergence of mixed finite element methods”,RAIRO Analyse Numérique, v. 11, pp. 341–354, 1977.

[65] FORTIN, M. “Old and new finite elements for incompressible flows”, Interna-tional Journal for Numerical Methods in Fluids, v. 1, pp. 347–364, 1981.

53