Elementos finitos em formulação mista de mínimos quadrados para a simulação da convecção-difusão em regime transiente

Dedicado a Eugenio Oñate no seu 60°. aniversário

Resumo

Um esquema de elementos finitos mistos de tipo mínimos quadrados é estudado para resolver as equações da convecção-difusão transientes expressas em função tanto do campo incógnito primário como do seu fluxo, incorporando ou não um termo reativo. Uma vez efetuada uma discretização temporal à la Crank-Nicholson, o sistema de equações resultante permite uma aproximação estável desses 2 campos, com elementos finitos contínuos de Lagrange clássicos de grau arbitrário em geometria de tipo simplex ou não, em dimensão espacial qualquer. O esquema é convergente em espaço no sentido da média quadrática no que tange ao campo incógnito primário, a seu gradiente, à variável de fluxo e à divergência desta última, e no tempo num sentido apropriado para cada um desses campos. Os resultados numéricos atestam o bom desempenho do esquema, para quaisquer números de Péclet, confirmando as previsões teóricas, pelo menos no caso de camadas limite estreitas em que o método se mostra fracassado. A técnica é também comparada com outros métodos para resolver essas equações, incluindo 2 propostos recentemente pelo primeiro autor et al.

Abstract

A mixed least-squares finite element scheme designed for solving the transient convection-diffusion equations expressed in terms of both the primal unknown and its flux, incorporating or not a reaction term, is studied. Once a time discretization of the Crank-Nicholson type is performed, the resulting system of equations allows for a stable approximation of both fields, by means of classical Lagrange continuous piecewise polynomial functions of arbitrary degree, in both simplicial and non-simplicial geometry, in any space dimension. The scheme is also convergent in space in the mean-square sense as far as the primal unknown field, its gradient, the flux variable and its divergence are concerned, and in time in an appropriate sense for each one of these four fields. Numerical results certify that the scheme performs well for any Péclet number, thereby allowing to confirm theoretical predictions, at least in the case where there is no narrow boundary layer. In the latter case however the method fails to produce reliable results. The technique is also compared with three existing methods to solve the convection-diffusion equations in the transient case. These include two recent ones proposed by the first author and collaborators.

Palavras-chave

Convecção-difusão ; Elementos finitos ; Mínimos quadrados ; Mistos ; Transiente

Keywords

Convection-diffusion ; Finite elements ; Least-squares ; Lixed ; Transient

1. Introdução

A resolução numérica das equações da convecção-difusão pode ser um problema delicado por diferentes razões, mesmo quando elas são lineares. Em particular, este é o caso de simulações com números de Péclet elevados, devido à necessidade de capturar fortes gradientes da solução próximos de camadas limite. No âmbito de uma resolução via elementos finitos, foram adotadas diversas, abordagens para lidar com esse problema de maneira satisfatória já nos anos 70. De entre elas, gostaríamos de citar as contribuições pioneiras de Heinrich et al. (ver p. ex. [1] ) e de Baba e Tabata [2] . O procedimento introduzido no início da década de 80 por Hughes e Brooks [3] continua a ser amplamente utilizado até hoje. Baseia-se numa modificação da formulação padrão de Galerkin mediante a introdução de difusão numérica estabilizadora na direção das linhas de corrente. Essa técnica, que é comumente conhecida sob a sigla SUPG (de streamline upwind Petrov-Galerkin), deu lugar a diversas variantes desde então, tais como as formulações de Galerkin de mínimos quadrados. Outros métodos estáveis se adequam à resolução explícita de problemas transientes, tais como os esquemas explorados por Kawahara et al. (ver p. ex. [4] ) nos anos 80 e 90. Essa técnica pode ser vista como adaptações ao contexto de elementos finitos, dos esquemas de Lax ou de Lax-Wendroff para discretizações espaciais e temporais por diferenças finitas. Mais especificamente, devem ser usadas no âmbito da modelagem numérica com elementos finitos lineares contínuos clássicos em formulação de Galerkin padrão. Recentemente, o primeiro autor et al. exploraram essa ideia modificando o método de maneira a fazê-lo gerar aproximações convergentes na norma do máximo, mesmo tratando-se de malhas não uniformes [5] .

Mais tarde, diversas formulações de Petrov-Galerkin ou de mínimos quadrados foram estudadas para lidar com essas equações na formulação mista que resulta da introdução da variável de fluxo. Isso foi feito principalmente para o caso estacionário (ver p. ex. [6] ). A propósito, uma alternativa de um campo a esses métodos foi proposta mais recentemente por Carneiro de Araújo e o primeiro autor (ver p. ex. [8] ). A principal característica desse método é uma interpolação quadrática de Hermite do campo incógnito primário, incorporando os fluxos médios através das interfaces dos elementos da malha como graus de liberdade. Experiências numéricas com esses elementos e métodos mistos clássicos de ordem comparável revelaram que os primeiros são globalmente mais precisos.

Outra contribuição original relativamente recente que merece ser registada aqui, por permitir resoluções numéricas estáveis com qualquer tipo de interpolação, é a técnica baseada em cálculo finito, introduzida por Oñate et al. em [7] .

No que tange a problemas transientes, exceto por um pequeno número de artigos com análises parciais tal como [9] , o estudo da convergência de métodos mistos parece não ter sido atacado a fundo, especialmente no caso em que a convecção tem um papel significativo. Nesse contexto preciso, a contribuição do segundo autor et al. baseia-se numa formulação mista de mínimos quadrados em espaço e em tempo [10] mostrou-se bastante eficaz do ponto de vista computacional, inclusive por permitir o uso de interpolações de ordem arbitrária para os 2 campos incógnitos. Neste artigo, os autores apresentam precisamente os resultados de estabilidade e convergência que podem ser esperados para esta última técnica na sua versão ligeiramente modificada, de forma a que a dependência temporal seja tratada por um esquema de Crank-Nicholson clássico, mantendo a mesma formulação de mínimos quadrados no que diz respeito à dependência espacial. Outro ponto interessante da presente contribuição é que nem a estabilidade nem a convergência do método está subordinada ao facto de que um termo reativo figure nas equações. Na realidade, o caso das equações de convecção-difusão-reação é tratado aqui como mero sub produto do estudo para problemas de convecção-difusão.

Uma das contribuições mais relevantes deste artigo é a apresentação de quantidade apreciável de resultados numéricos, destinados a testar e/ou validar a metodologia de mínimos quadrados aqui tratada, obtidos para uma série de casos-teste representativos de diversas situações correntemente encontradas na prática. De forma sucinta, podemos afirmar que esses resultados confirmam as nossas previsões teóricas quanto à estabilidade e a convergência, qualquer que seja o número de Péclet considerado, salvo na presença de camadas limite demasiado estreitas, em que o método necessitaria de alguma adaptação para as poder captar corretamente.

O artigo está organizado da forma seguinte: na secção 2 descrevemos o problema a ser resolvido e introduzimos as notações empregues no trabalho. Na secção 3 apresentamos as discretizações espacial e temporal relativas ao nosso método. Na secção 4, são fornecidos os principais resultados de estabilidade e convergência que se aplicam a este último. Em seguida, na secção 5, exploramos o esquema numérico computacionalmente para o caso de elementos finitos lineares, comparando-o com outros baseados num tratamento espacial pelo método de Galerkin. Mais especificamente, trata-se do esquema de Kawahara citado acima [4] com a modificação de [5] , do método misto de Raviart-Thomas de mais baixa ordem [11] e do método quadrático de Hermite com um só campo introduzido em [8] . A propósito, isso faz deste artigo o primeiro a incluir resultados computacionais com esse último método, pelo menos de publicações sobre o mesmo que não sejam do conhecimento dos autores. Por fim concluimos com alguns comentários na secção 6.

2. Descrição das equações e notações

Considera-se as equações transientes da convecção-difusão, incluindo ou não um termo reativo, definidas num domínio Ω × (0, T   ), onde Ω é um sub conjunto limitado de ${\textstyle R^{N}}$ , N  = 1, 2 ou 3 de fronteira ∂Ω, e T é um tempo finito, descritas como se segue:

Pretende-se determinar uma função escalar u (x , t ) satisfazendo:

${\begin{array}{l}{\partial }_{t}u-\nabla \cdot {\boldsymbol {\mbox{K}}}\nabla u+{\boldsymbol {\mbox{w}}}\cdot \nabla u+\sigma u=f\quad em\quad \Omega \times (0,T)\\u=g\quad sobre\quad {\Gamma }_{0}\times (0,T)\\{\boldsymbol {\mbox{K}}}\nabla u\cdot {\boldsymbol {\nu }}=0\quad sobre\quad {\Gamma }_{1}\times (0,T)\\u=u^{0}\quad em\quad \Omega \quad parat=0\end{array}}$

( 1)

onde ν é o vetor normal externo unitário sobre ∂Ω, Γ₀ e Γ₁ são 2 porções disjuntas de ∂Ω, ${\textstyle {\partial }_{t}G}$ representa derivada parcial de primeira ordem de uma função escalar ou vetorial ${\textstyle G}$ , e ∇ designa o operador gradiente. A medida de Γ₁ pode ser nula mas a de Γ₀ é supostamente estritamente positiva.

Ao longo de todo o texto, a notação correspondente a um ponto   intercalado entre 2 campos ou vetores é usada para designar o produto interno padrão de ${\textstyle R^{N}}$ . Num termo do tipo ∇ · a serve para representar a divergência de uma função vetorial a , e num termo do tipo a · ∇ significa o operator ${\textstyle {\sum }_{i=1}^{N}a_{i}{\frac {\partial }{\partial x_{i}}}}$ .

Em (1)K é um tensor de difusividade simétrico, que se supõe ser constante e positivo-definido, e w é uma velocidade de transporte de divergência nula pertencente a ${\textstyle [C^{m}({\overline {\Omega }})]^{N}}$ para certo ${\textstyle m\in {\mathbb {N} }^{\ast }}$ . Supomos ainda que w satisfaz a condição w  · ν  ≥ 0 sobre Γ₁  ∀ t . Nesse sentido, Γ₁ pode ser vista como uma parte de ∂Ω que contém somente porções, quer da saída, quer de paredes com deslizamento da região Ω, na qual um fluido incompressível escoa com velocidade w . σ , por sua vez, é um coeficiente constante não negativo relativo a um eventual fenómeno reativo associado ao processo convectivo-difusivo em estudo.

O dados f e g são, respetivamente, uma função de forças pertencente a L^∞ [(0, T ) ; L² (Ω)] (cf. [12] ) e um valor imposto sobre Γ₀  × (0, T ). Em benefício da simplicidade, mas sem nenhum prejuízo essencial para os resultados que apresentaremos a seguir, tomaremos neste trabalho g  ≡ 0. Suporemos ainda que u⁰  ∈ H¹ (Ω) (cf. [13] ). Vamos ainda requerer alguma regularidade adicional tanto para u⁰ como para f , a ser especificada mais adiante.

Em seguida, definimos 2 espaços de Hilbert para as normas naturais (cf. [14] e [15] ), os quais vão representar um papel fundamental neste estudo, a saber, ${\textstyle V:=\lbrace v/v\in H^{1}(\Omega ),v_{/{\Gamma }_{0}}=}$$0\rbrace$ e ${\textstyle {\boldsymbol {\mbox{Q}}}:=\lbrace {\boldsymbol {\mbox{q}}}/{\boldsymbol {\mbox{q}}}\in {\boldsymbol {\mbox{H}}}(\Omega ,div),{\boldsymbol {\mbox{q}}}_{/{\Gamma }_{1}}\cdot {\boldsymbol {\nu }}=}$$0\rbrace$ (cf. [15] ). Introduzimos ainda a variável de fluxo p  : = − K  ∇ u , a qual pertence a Q por hipótese. Assim sendo, dada uma constante α estritamente positiva, colocamos (1) na forma variacional mista equivalente de tipo mínimos quadrados (2) abaixo, onde (A , B ) designa o produto interno padrão de 2 funções escalares ou vetoriais A e B de L² (Ω), e (A , B )_K representa (KA , B ), sendo A e B 2 funções vetoriais. Assim, sendo I o operador identidade, temos:

${\begin{array}{l}Encontrar\quad u(\cdot ,t)\in V\quad bem\quad como\quad {\boldsymbol {\mbox{p}}}\in {\boldsymbol {\mbox{Q}}},\\com\quad {\partial }_{t}u(\cdot ,t)\in L^{2}(\Omega )\quad e\quad u=u^{0}\quad em\quad \Omega \quad para\quad t=0,\\e\quad {\boldsymbol {\mbox{p}}}=-{\boldsymbol {\mbox{K}}}\nabla u^{0}\quad em\quad \Omega \quad para\quad t=0,\\tais\quad que\quad \forall v\in V\quad e\quad \forall {\boldsymbol {\mbox{q}}}\in {\boldsymbol {\mbox{Q}}}\quad tenhamos\quad \forall t\in (0,T):\\({\partial }_{t}u+\nabla \cdot {\boldsymbol {\mbox{p}}}+[{\boldsymbol {\mbox{w}}}\cdot \nabla +\sigma I]u,\\v+\alpha \lbrace \nabla \cdot {\boldsymbol {\mbox{q}}}+[{\boldsymbol {\mbox{w}}}\cdot \nabla +\sigma I]v\rbrace )\\+(\nabla u+{\boldsymbol {\mbox{K}}}^{-1}{\boldsymbol {\mbox{p}}},\nabla v+{\boldsymbol {\mbox{K}}}^{-1}{\boldsymbol {\mbox{q}}})_{\boldsymbol {\mbox{K}}}\\=(f,v+\alpha \lbrace \nabla \cdot {\boldsymbol {\mbox{q}}}+[{\boldsymbol {\mbox{w}}}\cdot \nabla +\sigma I]v\rbrace ).\end{array}}$

( 2)

Observe-se que (· , ·)_K é um produto interno sobre L² (Ω)^N cuja norma associada denotada por ∥ ·  ∥ _K é equivalente à norma padrão desse espaço, a qual será representada daqui em diante por ∥·  ∥, inclusive no caso escalar. Mais especificamente, vale

$\lambda \parallel \quad {\boldsymbol {\mbox{A}}}\quad {\parallel }^{2}\leq \parallel \quad {\boldsymbol {\mbox{A}}}\quad {\parallel }_{\boldsymbol {\mbox{K}}}^{2}\leq \mu \parallel \quad {\boldsymbol {\mbox{A}}}\quad {\parallel }^{2},\quad \forall {\boldsymbol {\mbox{A}}}\in L^{2}(\Omega )^{N}$

( 3)

onde λ e μ são, respetivamente, o menor autovalor e o maior autovalor de K .

Observação 1: Caso a equação (1) esteja escrita na forma adimensional, λ⁻¹ representa o número de Péclet.   □

No restante deste artigo, denotamos por ∥ ·  ∥ _{m ,p ,Ω} a norma padrão do espaço de Sobolev W^{m ,p} (Ω) (cf. [13] ) com ${\textstyle m\in \mathbb {N} }$ e ${\textstyle p\in R}$ , p  ≥ 1. W^{m ,2} (Ω) é representado por H^m (Ω) para m  ≠ 0. A norma padrão de H^m (Ω) é simplesmente denotada por ∥ ·  ∥ _{m ,Ω} exceto para m  = 0, sendo que todas essas notações se aplicam a versões tanto escalares como vetoriais dos espaços correspondentes.

3. Discretização no espaço e no tempo

Para a discretização temporal da equação (1) , empregamos um esquema do tipo Crank-Nicolson. Mais especificamente, dado um número inteiro M , M  > 1, definimos um passo de tempo Δt  = T /M . Em seguida, definindo f^r  = f (r Δt ) para todo r  ∈ [0, M ] e partindo de u⁰ e p⁰ , determinamos sucessivamente para n  = 1, 2, …, M uma aproximação de (u [n Δt ] ; p [n Δt ]) designada por (uⁿ  ; pⁿ ), como solução do problema de mínimos quadrados seguinte, em tudo análogo a (2) :

${\begin{array}{l}Dado\quad u^{0}\in V\quad e\quad tomando\quad {\boldsymbol {\mbox{p}}}^{0}:=-K\nabla u^{0}\in {\boldsymbol {\mbox{Q}}},\\para\quad n=1,2,\ldots ,M,\quad encontrar\quad u^{n}\in V\quad e\quad {\boldsymbol {\mbox{p}}}^{n}\in {\boldsymbol {\mbox{Q}}}\\tais\quad que\quad \forall v\in V\quad e\quad \forall {\boldsymbol {\mbox{q}}}\in {\boldsymbol {\mbox{Q}}}:\\(u^{n}+{\frac {\Delta t}{2}}\lbrace \nabla \cdot {\boldsymbol {\mbox{p}}}^{n}+[{\boldsymbol {\mbox{w}}}^{n-1/2}\cdot \nabla +\sigma I]u^{n}\rbrace ,\\v+\alpha \lbrace \nabla \cdot {\boldsymbol {\mbox{q}}}+[{\boldsymbol {\mbox{w}}}^{n-1/2}\cdot \nabla +\sigma I]v\rbrace )\\+\Delta t(\nabla u^{n}+K^{-1}{\boldsymbol {\mbox{p}}}^{n},\nabla v+K^{-1}{\boldsymbol {\mbox{q}}})_{K}/2\\=\quad \Delta t(f^{n-1/2},v+\alpha \lbrace \nabla \cdot {\boldsymbol {\mbox{q}}}+[{\boldsymbol {\mbox{w}}}^{n-1/2}\cdot \nabla +\sigma I]v\rbrace )\\+(u^{n-1}-{\frac {\Delta t}{2}}\lbrace \nabla \cdot {\boldsymbol {\mbox{p}}}^{n-1}+[{\boldsymbol {\mbox{w}}}^{n-1/2}\cdot \nabla +\sigma I]u^{n-1}\rbrace ,\\v+\alpha \lbrace \nabla \cdot {\boldsymbol {\mbox{q}}}+[{\boldsymbol {\mbox{w}}}^{n-1/2}\cdot \nabla +\sigma I]v\rbrace )\\-\Delta t(\nabla u^{n-1}+K^{-1}{\boldsymbol {\mbox{p}}}^{n-1},\nabla v+K^{-1}{\boldsymbol {\mbox{q}}})_{K}/2\\\end{array}}$

( 4)

Daqui em diante, em benefício da simplicidade, suporemos que Ω é um intervalo se N  = 1, um polígono se N  = 2 ou um poliedro se N  = 3. Assim fazendo, vamos considerar um sistema discretizado no espaço análogo a (4) definido como se segue:

Seja ${\textstyle T_{h}}$ uma partição de Ω em intervalos para N  = 1, em triângulos ou quadrângulos convexos para N  = 2, e em tetraedros ou hexaedros convexos com faces quadrangulares para N  = 3, de tamanho de lado máximo igual a h   . Faremos a hipótese de que ${\textstyle T_{h}}$ satisfaz as condições de compatibilidade habituais que se aplicam a malhas de elementos finitos, e que ela pertence a uma família quasiuniforme de partições de Ω. Suporemos ainda que tanto Γ₀ como Γ₁ são tais que podem ser completamente cobertas pela união de lados para N  = 2 ou faces para N  = 3, de elementos pertencentes a ${\textstyle T_{h}}$ . Se ${\textstyle T_{h}}$ é constituído de N   -simplex, para cada ${\textstyle E\in T_{h}}$R_k (E ) é definido como o espaço de polinómios de grau menor ou igual a k e, senão R_k (E ) é o espaço de funções definidas como as transformadas de polinómios de grau menor ou igual a k em cada uma das N   variáveis espaciais, definidos num quadrado ou cubo unitário ${\textstyle {\overset {\mbox{ˆ}}{E}}}$ conforme o caso, pela aplicação N   -linear bijetiva que mapeia ${\textstyle {\overset {\mbox{ˆ}}{E}}}$ em E   . Dessa forma, para todo ${\textstyle k\in {\mathbb {N} }^{\ast }}$ , introduzimos os seguintes espaços vetoriais associados a ${\textstyle T_{h}}$ :

$S_{h,k}:=\left\{v\quad \vert \quad v\in C^{0}({\overline {\Omega }})\quad e\quad v_{/E}\in R_{k}(E),\quad \forall E\in T_{h}\right\},$

${\boldsymbol {\mbox{Q}}}_{h}:=\left\{{\boldsymbol {\mbox{q}}}\in {\boldsymbol {\mbox{Q}}}\quad \vert \quad \forall i\quad q_{i}\in S_{h,l}\right\}\quad para\quad l\in {\mathbb {N} }^{\ast },$

$V_{h}:=S_{h,j}\cap V\quad para\quad j\in {\mathbb {N} }^{\ast }.$

Sejam agora ${\textstyle u_{h}^{0}}$ e ${\textstyle {\boldsymbol {\mbox{p}}}_{h}^{0}}$ os interpolantes padrão de u⁰ em V_h e de p⁰ em Q_h , respectivamante. Além disso, para efeitos práticos, em princípio temos de utilizar integração numérica nos termos que involvem w e f . Neste trabalho, isto será levado a cabo substituindo w ou f^r pelos seus interpolantes padrão w_h e ${\textstyle f_{h}^{r}}$ em [S_{h ,l +1} ]^N e S_{h ,j} , respetivamente. Com isso, a nossa aproximação espacial e temporal do problema (1) na forma mista formula-se como se segue:

${\begin{array}{l}Partindo\quad de\quad u_{h}^{0}\in V_{h}\quad e\quad {\boldsymbol {\mbox{p}}}_{h}^{0}\in {\boldsymbol {\mbox{Q}}}_{h},\\para\quad n=1,2,\ldots ,M,encontrar\quad u_{h}^{n}\in V_{h}\quad e\quad {\boldsymbol {\mbox{p}}}_{h}^{n}\in {\boldsymbol {\mbox{Q}}}_{h}\\tais\quad que\quad \forall v\in V_{h}\quad e\quad \forall {\boldsymbol {\mbox{q}}}\in {\boldsymbol {\mbox{Q}}}_{h}:\\(u_{h}^{n}+{\frac {\Delta t}{2}}\lbrace \nabla \cdot {\boldsymbol {\mbox{p}}}_{h}^{n}+[{\boldsymbol {\mbox{w}}}_{h}^{n-1/2}\cdot \nabla +\sigma I]u_{h}^{n}\rbrace ,\\v+\alpha \lbrace \nabla \cdot {\boldsymbol {\mbox{q}}}+[{\boldsymbol {\mbox{w}}}_{h}^{n-1/2}\cdot \nabla +\sigma I]v\rbrace )\\+\Delta t(\nabla u_{h}^{n}+K^{-1}{\boldsymbol {\mbox{p}}}_{h}^{n},\nabla v+K^{-1}{\boldsymbol {\mbox{q}}})_{K}/2\\=\quad \Delta t(f^{n-1/2},v+\alpha \lbrace \nabla \cdot {\boldsymbol {\mbox{q}}}+[{\boldsymbol {\mbox{w}}}_{h}^{n-1/2}\cdot \nabla +\sigma I]v\rbrace )\\+(u_{h}^{n-1}-{\frac {\Delta t}{2}}\lbrace \nabla \cdot {\boldsymbol {\mbox{p}}}_{h}^{n-1}+[{\boldsymbol {\mbox{w}}}_{h}^{n-1/2}\cdot \nabla +\sigma I]u_{h}^{n-1}\rbrace ,\\v+\alpha \lbrace \nabla \cdot {\boldsymbol {\mbox{q}}}+[{\boldsymbol {\mbox{w}}}_{h}^{n-1/2}\cdot \nabla +\sigma I]v\rbrace )\\-\Delta t(\nabla u_{h}^{n-1}+K^{-1}{\boldsymbol {\mbox{p}}}_{h}^{n-1},\nabla v+K^{-1}{\boldsymbol {\mbox{q}}})_{K}/2\\\end{array}}$

( 5)

Uma simples aplicação do Teorema de Lax-Milgram, permite afirmar que o problema (5) admite solução única.

4. Estabilidade e convergência

Observação 2: No restante deste artigo, a letra C combinada ou não com outros símbolos representará diversas constantes estritamente positivas independentes de Δt e de h .  □

Primeiramente, resumimos no Teorema 1 abaixo os resultados de estabilidade que se aplicam ao esquema (5) . Com esse fim, lembramos primeiramente que ∇ · w  = 0. Assim, supondo que vale o resultado de regularidade w  ∈ [H^{l +2} (Ω)]^N , onde m é o maior número inteiro tal que 1 ≤ m  < l  + 2 − N   /2, temos ${\textstyle {\boldsymbol {\mbox{w}}}\in [C^{m}({\overline {\Omega }})]^{N}}$ (cf. [13] ). Então, de acordo com resultados de aproximação clássicos (cf. [14] ), conclui-se facilmente que existe uma constante C_w tal que,

$D:=\vert \vert \nabla \cdot {\boldsymbol {\mbox{w}}}_{h}\vert \vert _{0,\infty ,\Omega }\leq C_{\boldsymbol {\mbox{w}}}h^{m-1}\vert \vert {\boldsymbol {\mbox{w}}}\vert \vert _{m,\infty ,\Omega }.$

( 6)

Então, sendo C_l uma constante para a qual vale ||w  − w_h ||_0,∞,Ω  ≤ (C_l  − 1)||w ||_0,∞,Ω (cf. [14] ), temos:

Teorema 1.

Tomando α  = Δt   /2, e definindo ${\textstyle \gamma :=max\lbrace max[2\sigma ,4]+}$${\frac {D}{4}},{\frac {19W^{2}}{2\delta \lambda }},4+{\frac {2W^{2}}{\lambda }},{\frac {3W^{2}}{\lambda }}\rbrace$ com ${\textstyle \delta =(3-{\sqrt {5}})/4}$ , onde ${\textstyle W:=C_{l}{\underset {{\boldsymbol {\mbox{x}}}\in \Omega }{sup}}\vert {\boldsymbol {\mbox{w}}}({\boldsymbol {\mbox{x}}})\vert }$ , supondo que ${\textstyle \gamma \Delta t\leq {\frac {1}{2}}}$ , o resultado de estabilidade seguinte aplica-se ao esquema (5) , para qualquer n  ≤ M :

${\begin{array}{c}\parallel \quad u_{h}^{n}\quad {\parallel }^{2}+{\frac {\lambda \Delta t}{2}}\parallel \quad \nabla u_{h}^{n}\quad {\parallel }^{2}\\+\Delta t\sum _{i=1}^{n}{\frac {\parallel \quad u_{h}^{i}-u_{h}^{i-1}\quad {\parallel }^{2}}{8\Delta t}}+{\frac {\delta }{\mu }}\parallel \quad {\boldsymbol {\mbox{p}}}_{h}^{i}+{\boldsymbol {\mbox{p}}}_{h}^{i-1}\quad {\parallel }^{2}\\+{\frac {\Delta t}{24}}\parallel \quad \nabla \cdot ({\boldsymbol {\mbox{p}}}_{h}^{i}+{\boldsymbol {\mbox{p}}}_{h}^{i-1})\quad {\parallel }^{2}\leq (3e^{2})^{\gamma T}\\\left[\parallel \quad u_{h}^{0}\quad {\parallel }^{2}+{\frac {\mu \Delta t}{2}}\parallel \quad \nabla u_{h}^{0}\quad {\parallel }^{2}+\Delta t\sum _{i=1}^{n}\parallel \quad f_{h}^{i-1/2}\quad {\parallel }^{2}\right].\end{array}}$

( 7)

□

Em seguida, fechamos esta secção dando um resultado de convergência para o esquema (5) nas normas naturais delineadas em (7) . Esses resultados são garantidos pelo menos para o caso em que Ω é um domínio convexo. Além disso, sendo τ um número real estritamente positivo, fazemos a hipótese de que M é escolhido de tal forma que Δt  = C_Δh^τ .

Além da norma de L^∞ [(0, T ), L² (Ω)] (cf. [12] ) empregamos a norma discreta ∥ ·  ∥ _M de L² [(0, T ) ; L² (Ω)] definida por ${\textstyle \parallel \quad G\quad {\parallel }_{M}:=}$$[{\sum }_{n=1}^{M}\Delta t\parallel \quad G^{n}\quad {\parallel }^{2}]^{1/2}$ , onde ${\textstyle G^{n}=G(n\Delta t),n=1,2,\ldots ,M}$ , sendo ${\textstyle G}$ um campo de C⁰ {[0, T ] ; L² (Ω)}^L com L  = 1 ou L  = N . Além disso, é conveniente definir uma função p_h (x , t ) constante em cada intervalo I_n  : = ([n  − 1]Δt , n Δt ) para cada x , de valor igual a ${\textstyle [{\boldsymbol {\mbox{p}}}_{h}^{n}({\boldsymbol {\mbox{x}}})+}$${\boldsymbol {\mbox{p}}}_{h}^{n-1}({\boldsymbol {\mbox{x}}})]/2$ , n  = 1, 2, …, M . Uma estimativa de erro na norma ∥ ·  ∥ _M é dada também para a derivada temporal de u aproximada por uma função u_h obtida com o nosso método. Mais especificamente, u_h (x , t ) é definida como a função que varia linearmente com t em cada intervalo I_n , e cujo valor para t  = n Δt   é ${\textstyle u_{h}^{n}({\boldsymbol {\mbox{x}}})}$ ∀x  ∈ Ω, n  = 1, 2, …, M   . Desta forma, fazendo ${\textstyle ({\tilde {u}}^{n};{\tilde {\boldsymbol {\mbox{p}}}}^{n}):=}$$(u[n\Delta t];{\boldsymbol {\mbox{p}}}[n\Delta t])$ , temos:

Teorema 2.

Suponhamos que f  ∈ L^∞ [(0, T ) ; H^{j +1} (Ω)], w  ∈ H^{l +2} (Ω)}^N , ∂_tttu  ∈ L^∞ [(0, T ) ; H¹ (Ω)], u , ∂_tu ∈L^∞ [(0, T ) ; H^k (Ω)], u⁰  ∈ H^k (Ω) para k  = max {l , j } + 1, que Δt satisfaz a relação Δt  = C_Δh^τ e que h é suficientemente pequeno para que valha a desigualdade γ Δt  ≤ 1/2. Nessas condições, existe uma constante ${\textstyle {\overline {C}}}$ que depende somente de τ , l , j , Ω, T , K , w e σ tal que ∀n  ≤ M a seguinte estimativa de erro a priori é verdadeira:

${\begin{array}{l}\parallel \quad {\tilde {u}}^{n}-u_{h}^{n}\quad \parallel +h^{\tau /2}\parallel \quad \nabla ({\tilde {u}}^{n}-u_{h}^{n})\quad \parallel \\\leq {\overline {C}}h^{\eta }[G(u)+\parallel \quad u^{0}\quad {\parallel }_{k,\Omega }],\\\parallel \quad {\boldsymbol {\mbox{p}}}-{\boldsymbol {\mbox{p}}}_{h}\quad {\parallel }_{M}+h^{\tau /2}\parallel \quad \nabla \cdot ({\boldsymbol {\mbox{p}}}-{\boldsymbol {\mbox{p}}}_{h})\quad {\parallel }_{M}\\\leq {\overline {C}}h^{\eta -\tau /2}[G(u)+\parallel \quad u^{0}\quad {\parallel }_{k,\Omega }]\\\parallel \quad {\partial }_{t}[u-u_{h}]\quad {\parallel }_{M}\leq {\overline {C}}h^{\rho }\quad [G(u)+\parallel \quad u^{0}\quad {\parallel }_{k,\Omega }],\end{array}}$

( 8)

onde η  : = min {2τ , l  + min [2, τ ], j  + min [1, τ /2]}, ρ  : = min {τ , l  + min [2 − τ /2, τ /2], j  + min [1 − τ   /2, 0]}, e ${\textstyle G(u):={\underset {s\in (0,T)}{sup}}ess\lbrace \parallel \quad f(s)\quad {\parallel }_{j+1,\Omega }+}$$\parallel \quad {\partial }_{t}u(s)\quad {\parallel }_{k,\Omega }+$$\parallel \quad u(s)\quad {\parallel }_{k,\Omega }+$$\parallel \quad \nabla {\partial }_{ttt}u(s)\quad \parallel \rbrace$ .

□

A demonstração dos Teoremas 1 e 2 pode ser encontrada em [16] .

5. Experiências numéricas

Passamos a fornecer alguns resultados numéricos representativos do desempenho do método estudado nas secções 3 e 4.

5.1. Testes para verificação da convergência

No intuito de observar a validade das estimativas de erro apresentadas na secção anterior, levamos a cabo algumas experiências com o método de mínimos quadrados em questão, nas quais resolvemos problemas de solução analítica conhecida, tomando representações lineares por triângulo para ambos os campos incógnitos.

Teste n. 1: Neste primeiro teste, o domínio de definição da equação (1) é Ω × (0, T ), onde Ω é o quadrado unitário (0, 1)² e T  = 1. Apresentamos abaixo alguns resultados relevantes para o caso em que Γ₀ é a porção de ∂Ω de equação xy  = 0, tomando K  = I e w (x , y ) = (y  ; x   )/2, o que corresponde a um número de Péclet igual a ${\textstyle {\sqrt {2}}/2}$ . Consideramos uma solução analítica exata dada por:

${\begin{array}{ccc}u(x,y,t)&=&(x-x^{2}/2)(y-y^{2}/2)e^{-t};\\{\boldsymbol {\mbox{p}}}(x,y,t)&=&([x-1][y-y^{2}/2];[y-1][x-x^{2}/2])e^{-t}.\end{array}}$

Assim fazendo, a função de força f é definida por

$f:=u_{t}+\nabla \cdot {\boldsymbol {\mbox{p}}}-{\boldsymbol {\mbox{w}}}\cdot {\boldsymbol {\mbox{p}}}.$

( 9)

Resolvemos este problema com malhas uniformes triangulares obtidas mediante uma primeira subdivisão de Ω em L² quadrados iguais com lado de comprimento h  = 1/L , e, em seguida, subdividindo cada um deles em 2 triângulos pela sua diagonal paralela à reta x  = y . Mostramos nas tabelas 1 , 2 e 3 abaixo os erros absolutos da aproximação para valores crescentes de L , tomando M  = 2 × L , i.e. Δt  = h /2. Na tabela 1 fornecemos erros na norma de L² (Ω) dos valores approximados de u , ∇u , p , ∇ · p e u_t para t  = 1. Na tabela 2 provemos os erros de p , ∇ · p e u_t na norma de L² [(0, T ) ; L² (Ω)], em versão tanto discreta como contínua, dependendo do campo (cf. secção 4). Em seguida, na tabela 3 , fornecemos os erros absolutos de u na norma de L^∞ (Ω) para t  =0,5 e t =1.

Tabela 1. Erros absolutos de u , ∇u , p , ∇ · p e u_t na norma de L² (Ω) para t =1, com Pé=0,7071.

L	u	∇u	p	∇ · p	ut
8	0,51623 × 10−3	0,10008 × 10−1	0,14413 × 10−2	0,12313 × 10−1	0,59391 × 10−3
16	0,12767 × 10−3	0,49162 × 10−2	0,36909 × 10−3	0,61624 × 10−2	0,14355 × 10−3
32	0,31657 × 10−4	0,24350 × 10−2	0,95238 × 10−4	0,30742 × 10−2	0,35225 × 10−4
64	0,78814 × 10−5	0,12116 × 10−2	0,24893 × 10−4	0,15321 × 10−2	0,87263 × 10−5

Tabela 2. Erros absolutos de p , ∇ · p e u_t na norma de L² [(0, 1) ; L² (Ω)], com Pé=0,7071.

L	p	∇ · p	ut
8	0,23292 × 10−2	0,21268 × 10−1	0,14317 × 10−2
16	0,60978 × 10−3	0,10835 × 10−1	0,38081 × 10−3
32	0,15944 × 10−3	0,54507 × 10−2	0,10111 × 10−3
64	0,42068 × 10−4	0,27276 × 10−2	0,26818 × 10−4

Tabela 3. Erro absoluto máximo de u (x , t ) em Ω para t =0,5 e t =1, com Pé=0,7071.

L	∥ [u  − uh ](0, 5) ∥ 0,∞,Ω	∥ [u  − uh ](1) ∥ 0,∞,Ω
8	0,51960 × 10−3	0,34118 × 10−3
16	0,16323 × 10−3	0,10422 × 10−3
32	0,62399 × 10−4	0,39235 × 10−4
64	0,21037 × 10−4	0,13119 × 10−4

Como se pode inferir da tabela 1 , a taxa de convergência observada a partir dos resultados numéricos concordam perfeitamente com as previsões teóricas, no que tange o gradiente de u . No entanto, de acordo com a mesma tabela, uma convergência quadrática é observada para a norma L² dos erros de u , p , ∇ · p e u_t a cada tempo t . No que diz respeito a u , isto é melhor do que a taxa de convergência de 3/2 que estimamos, mas todos esses resultados poderiam explicar-se pelo fenómeno de superconvergência próprio de malhas uniformes.

Já a tabela 2 indica taxas de convergência na norma de L² [(0 ; 1) ; L² (Ω)] de cerca de 2, 1 e 2 para p , ∇ · p e u_t , respectivamente. Aqui também isto representa uma melhoria em relação às, nossas estimadas taxas de 1, 1/2 e 1, mas uma vez mais, ela pode ter origem na superconvergência típica de malhas uniformes. Enfim, a tabela 3 que não corresponde a nenhuma de nossas estimativas a priori de erro, mostra uma taxa de convergência para u na norma de L^∞ [(0 ; 1) ; L^∞ (Ω)] menos bem definida, mas oscilando em torno de 3/2.

Teste n. 2: Com o objetivo de verificar o comportamento do método no caso de fronteiras curvas, tomamos agora w (x , y ) = 2(x /r²  ; y /r² ), com ${\textstyle r={\sqrt {x^{2}+y^{2}}}}$ e novamente K  = I . A equação (1) é considerada no domínio Ω × (0, 1), onde Ω é o quarto de coroa circular dado por {(x  ; y )/0, 2 < r  < 1, x  > 0, y  > 0}, o que dá um número de Péclet igual a 10. Apresentamos abaixo alguns resultados relevantes para o caso em que Γ₀ é a união das porções de ∂Ω formadas pelas circunferências interna e externa da coroa (i.e. para r  = 0, 2 e r  = 1).

Consideramos uma solução analítica exata dada por:

${\begin{array}{ccc}u(x,y,t)&=&(1-x^{2}-y^{2})(2t^{2}-t^{3});\\{\boldsymbol {\mbox{p}}}(x,y,t)&=&(x;y)[4t^{2}-2t^{3}].\end{array}}$

e, para isso, tomamos de novo a função de força f expressa por (9) , o que, neste caso particular, se reduz a fazer f  = u_t .

Resolvemos este problema com malhas quasi-uniformes triangulares obtidas pelo mapeamento da malha uniforme de um quadrado unitário descrita no Teste n. 1, em Ω de tal forma que as coordenadas polares dos vértices dos triângulos da malha deste último sejam dadas por θ  = πy ′/2 e ρ  = 0, 2(1 − x ′) + x ′, sendo (x ′ ; y ′) as coordenadas cartesianas dos vértices da malha do quadrado unitário.

De forma idêntica às tabelas 1 , 2 e 3 , apresentamos nas tabelas 4 , 5 e 6 abaixo os erros absolutos das aproximações determinadas para valores crescentes de L , tomando novamente M  = 2 × L .

Tabela 4. Erros absolutos de u , ∇u , p , ∇ · p e u_t na norma de L² (Ω) para t = 1, com Pé = 10.

L	u	∇u	p	∇ · p	ut
8	0,58714 × 10−2	0,96409 × 10−1	0,10442 × 10−1	0,31405 × 10−1	0,71557 × 10−2
16	0,14758 × 10−2	0,47999 × 10−1	0,27637 × 10−2	0,12091 × 10−1	0,17624 × 10−2
32	0,36766 × 10−3	0,23924 × 10−1	0,72462 × 10−3	0,46634 × 10−2	0,43382 × 10−3
64	0,91531 × 10−4	0,11939 × 10−1	0,19051 × 10−3	0,18316 × 10−2	0,10731 × 10−3

Tabela 5. Erros absolutos de p , ∇ · p e u_t na norma de L² [(0, 1) ; L² (Ω)], com Pé=10.

L	p	∇ · p	ut
8	0,56371 × 10−2	0,16351 × 10−1	0,66105 × 10−2
16	0,14533 × 10−2	0,61521 × 10−2	0,16610 × 10−2
32	0,37594 × 10−3	0,23792 × 10−2	0,41373 × 10−3
64	0,98229 × 10−4	0,94291 × 10−3	0,10299 × 10−3

Tabela 6. Erro absoluto máximo de u (x , t ) em Ω para t =0,5 e t =1, com Pé=10.

L	∥ [u  − uh ](0, 5) ∥ 0,∞,Ω	∥ [u  − uh ](1) ∥ 0,∞,Ω
8	0,16996 × 10−2	0,54675 × 10−2
16	0,40829 × 10−3	0,13321 × 10−2
32	0,98540 × 10−4	0,32295 × 10−3
64	0,24078 × 10−4	0,79155 × 10−4

As tabelas 4 e 5 indicam taxas de convergência semelhantes às constatadas no Teste n. 1. Isto mostra não haver nenhuma perda na qualidade da solução numérica devido à aproximação da fronteira curva por uma poligonal. Ao contrário, a tabela 6 permite observar uma taxa de convergência quadrática para u na norma do máximo. No entanto, isto poderia ser decorrente da forma particular da solução exata.

5.2. Testes de comparação com outros métodos

Nesta subsecção, resolvemos com o mesmo tipo de aproximação usada na anterior, problemas que admitem solução analítica exata, no intuito de comparar essa técnica de resolução numérica com outros métodos de elementos finitos de ordem equivalente, propostos para resolver a equação (1) .

Numa primeira série de experiências, consideramos apenas a difusão pura tomando w  = 0 , ou seja Pé=0. Comparamos o desempenho dos elementos finitos lineares clássicos para aproximar os 2 campos incógnitos, mediante o emprego da nossa formulação de mínimos quadrados, com o elemento misto de Raviart-Thomas de mais baixa ordem [11] , e também com o elemento de Hermite quadrático - ou equivalentemente, P2 de Hermite -, com um só campo incógnito, introduzido em [8] , ambos em formulação de Galerkin padrão, mista ou clássica. A integração no tempo é efetuada também com o esquema de Crank-Nicholson para os 2 métodos.

Nos conjuntos de tabelas Tabela 7 , Tabela 8 , Tabela 9  and Tabela 10 e Tabela 12 , Tabela 13 , Tabela 14  and Tabela 15 , apresentamos os erros absolutos no tempo final T escolhido, medidos na norma de L² (Ω) e calculados para os 3 métodos, sucessivamente de u , ∇u , u_t e Δu para cada um dos 2 testes desta subsecção. No que tange o gradiente de u , é o valor do campo p_h que é utilizado para o calcular para os métodos de mínimos quadrados e de Raviart-Thomas. Quanto ao laplaciano de u , trata-se, na realidade, da divergência de p_h para esses 2 métodos mistos. Como no caso do elemento de Hermite, o laplaciano de u_h não é computa´vel, para o tipo de interpolação particular utilizado, adotamos o seguinte procedimento para a aproximação de Δu : os 3 fluxos através dos lados são interpolados com o mesmo tipo de representação que o usado para o campo p no caso do elemento de Raviart-Thomas de mais baixa ordem (vale dizer com campos da forma (cx  + a  ; cy  + b ) por triângulo). Em seguida, aproximamos o laplaciano de u pela divergência do campo interpolante elemento a elemento, sendo, esta sim, computável, por ser o elemento de Raviart-Thomas conforme em H (Ω, div ).

Tabela 7. Erros absolutos de u na norma de L² (Ω) para t =2, com Pé=0.

L	Mín. quadrados	Raviart-Thomas	P2 de Hermite
8	0,74186 × 10−2	0,70430 × 10−1	0,24029 × 10−2
16	0,18750 × 10−2	0,35196 × 10−1	0,60755 × 10−3
32	0,46958 × 10−3	0,17596 × 10−1	0,15220 × 10−3
64	0,11735 × 10−3	0,87976 × 10−2	0,38065 × 10−4

Tabela 8. Erros absolutos de ∇u (ou p ) na norma de L² (Ω) para t =2, com Pé=0.

L	Mín. quadrados	Raviart-Thomas	P2 de Hermite
8	0,21285 × 10−1	0,13731 × 10+0	0,12452 × 10+0
16	0,59689 × 10−2	0,68856 × 10−1	0,62450 × 10−1
32	0,16144 × 10−2	0,34438 × 10−1	0,31243 × 10−1
64	0,43448 × 10−3	0,17217 × 10−1	0,15622 × 10−1

Tabela 9. Erros absolutos de u_t na norma de L² (Ω) para t =2, com Pé=0.

L	Mín. quadrados	Raviart-Thomas	P2 de Hermite
8	0,65171 × 10−2	0,92591 × 10−1	0,58182 × 10−1
16	0,16981 × 10−2	0,46534 × 10−1	0,29920 × 10−1
32	0,43178 × 10−3	0,23324 × 10−1	0,15176 × 10−1
64	0,10872 × 10−3	0,11676 × 10−1	0,76430 × 10−2

Tabela 10. Erros absolutos de Δu (ou ∇ · p ) na norma de L² (Ω) para t =2, com Pé=0.

L	Mín. quadrados	Raviart-Thomas	P2 de Hermite
8	0,21831 × 10+0	0,21218 × 10+0	0,21107 × 10+0
16	0,11584 × 10+0	0,10828 × 10+0	0,10697 × 10+0
32	0,59636 × 10−1	0,56156 × 10−1	0,54893 × 10−1
64	0,30209 × 10−1	0,29963 × 10−1	0,28812 × 10−1

Tabela 12. Erros absolutos de u na norma de L² (Ω) para t =1, com Pé=0.

L	Mín. quadrados	Raviart-Thomas	P2 de Hermite
8	0,34024 × 10−3	0,33929 × 10−2	0,17325 × 10−3
16	0,84995 × 10−4	0,16977 × 10−2	0,43393 × 10−4
32	0,21198 × 10−4	0,84897 × 10−3	0,10854 × 10−4
64	0,52896 × 10−5	0,42451 × 10−3	0,27130 × 10−5

Tabela 13. Erros absolutos de ∇u (ou p ) na norma de L² (Ω) para t =1, com Pé=0.

L	Mín. quadrados	Raviart-Thomas	P2 de Hermite
8	0,51724 × 10−3	0,42261 × 10−2	0,13665 × 10−3
16	0,15157 × 10−3	0,14797 × 10−2	0,34674 × 10−4
32	0,44395 × 10−4	0,52048 × 10−3	0,87319 × 10−5
64	0,12986 × 10−4	0,18354 × 10−3	0,21922 × 10−5

Tabela 14. Erros absolutos de u_t na norma de L² (Ω) para t =1, com Pé=0.

L	Mín. quadrados	Raviart-Thomas	P2 de Hermite
8	0,36462 × 10−3	0,38883 × 10−2	0,13397 × 10−2
16	0,89528 × 10−4	0,19013 × 10−2	0,70337 × 10−3
32	0,22142 × 10−4	0,93845 × 10−3	0,35977 × 10−3
64	0,55037 × 10−5	0,46598 × 10−3	0,18187 × 10−3

Tabela 15. Erros absolutos de Δu na norma de L² (Ω) para t =1, com Pé=0.

L	Mín. quadrados	Raviart-Thomas	P2 de Hermite
8	0,30261 × 10−2	0,79777 × 10−1	0,18590 × 10−2
16	0,11917 × 10−2	0,56347 × 10−1	0,88593 × 10−3
32	0,48081 × 10−3	0,39858 × 10−1	0,42961 × 10−3
64	0,19739 × 10−3	0,28194 × 10−1	0,21100 × 10−3

Nas tabelas 11 e 16 são mostrados os erros absolutos máximos de u no domínio Ω no instante t  = T , para cada um dos 2 casos-teste.

Tabela 11. Erro absoluto máximo de u nos baricentros dos triângulos para t =2, com Pé=0.

L	Mín. quadrados	Raviart-Thomas	P2 de Hermite
8	0,97375 × 10−2	0,73056 × 10−2	0,80272 × 10−2
16	0,24267 × 10−2	0,18436 × 10−2	0,21745 × 10−2
32	0,60173 × 10−3	0,46199 × 10−3	0,56464 × 10−3
64	0,14921 × 10−3	0,11557 × 10−3	0,14379 × 10−3

Teste n. 3: Este teste visa observar o comportamento dos 3 métodos no caso em que a solução evolui exponencialmente no tempo de forma crescente. Tomando ainda K  = I , aqui a equação da difusão transiente é definida no domínio Ω × (0, T ) onde T  = 2 e Ω é o quadrado unitário do Teste n. 1. Γ₀ é a mesma porção de fronteira do que naquele teste, sendo agora a solução exata dada por:

${\begin{array}{ccc}u(x,y,t)&=&(x-x^{2}/2)(y-y^{2}/2)e^{t};\\{\boldsymbol {\mbox{p}}}(x,y,t)&=&([x-1][y-y^{2}/2];[y-1][x-x^{2}/2])e^{t}.\end{array}}$

Calculamos com o mesmo tipo de malhas que no Teste n. 1, e tomamos Δt  = 1/L .

Os resultados da tabela 8 acima indicam que o esquema de mínimos quadrados é bem superior aos 2 outros no que tange o gradiente de u como esperado. Por outro lado o método P2 de Hermite é bem mais preciso do que o primeiro, e mais ainda do que o de Raviart-Thomas, para efeitos de aproximação da função u propriamente dita no sentido da média (quadrática), como se pode depreender da tabela 7 . Por outro lado, no caso de valores em pontos da malha, os 3 métodos revelam-se sensivelmente equivalentes, com ligeira vantagem para o de Raviart-Thomas, conforme indica a tabela 11 . Por outro lado, como mostra a tabela 9 , a superioridade do método de mínimos quadrados aparece muito claramente na aproximação da derivada temporal de u , posto que apresenta uma taxa de convergência quadrática, em contraste com uma taxa linear para o método P2 de Hermite e o de Raviart-Thomas. Enfim, no que diz respeito ao laplaciano de u , os 3 métodos testados parecem ser igualmente competitivos, de acordo com a tabela 10 .

Teste n. 4: Trata-se agora de comparar os 3 métodos no caso em que o domínio tem uma fronteira curva. Um dos objetivos do teste é avaliar o efeito da aproximação dessa fronteira por uma poligonal, o que implica transladar condições de contorno para a função u , que, na realidade, incidem sobre Γ₀ , para os lados de elementos que têm suas 2 extremidades sobre essa porção da fronteira. Observa-se que tal translação só é passível de erodir a qualidade das aproximações dos 2 métodos de Galerkin, dado que somente eles têm graus de liberdade associados aos lados dos elementos.

Tomamos novamente K  = I , sendo que agora a equação da difusão transiente é considerada no domínio Ω × (0, T ), onde T  = 1 e Ω é o quarto de disco definido por {(x  ; y )/r  < 1, x  > 0, y  > 0}. Apresentamos abaixo resultados para o caso em que Γ₀ é a borda do disco, e a solução exata é dada por:

${\begin{array}{ccc}u(x,y,t)&=&(1-x^{2}-y^{2})e^{-t}/4;\\{\boldsymbol {\mbox{p}}}(x,y,t)&=&(x;y)e^{-t}/2.\end{array}}$

Fixamos Δt em 1/2L para todos os 3 métodos, onde L é um parâmetro da malha utilizada constituida de 2L² triângulos obtidos pelo mapeamento baseado em coordenadas polares, da malha uniforme usada para um quadrado unitário no Teste n. 1 da forma descrita em [17] .

Na aproximação do gradiente de u , a tabela 13 não aponta mais para uma superioridade do esquema de mínimos quadrados em relação ao método P2 de Hermite, e nem tanto para o elemento de Raviart-Thomas. No entanto, isto deve-se possivelmente à forma particular da solução exata em termos das variáveis espaciais. Aqui, o método P2 de Hermite continua a ser mais preciso do que o primeiro, e ainda mais do que o de Raviart-Thomas, no que diz respeito à aproximação da função u propriamente dita na norma de L² (Ω), como mostra a tabela 12 . No entanto, no caso de valores em pontos da malha, os 3 métodos são comparáveis, embora com vantagem para o método P2 de Hermite, e desvantagem para o de mínimos quadrados, conforme se pode observar da tabela 16 . De novo, o método de mínimos quadrados revela a mesma excelente capacidade de aproximação da derivada temporal de u , quando comparado com os 2 outros bem menos precisos, segundo a tabela 14 . Mas agora, no que diz respeito ao laplaciano de u , os 3 métodos testados não podem mais ser considerados equivalentes. De facto, observando-se a tabela 15 , embora o método de mínimos quadrados apareça como um pouco mais competitivo do que o método P2 de Hermite, com uma taxa de convergência estimada em cerca de 5/4, contra uma taxa linear para este último, uma erosão significativa é constatada para o elemento de Raviart-Thomas, o qual apresenta uma taxa de convergência sublinear - de cerca de 1/2, para ser mais preciso. Se lembrarmos os resultados para o Teste n. 3 dados na tabela 10 , isto poderia significar que a aproximação da fronteira curva por uma poligonal afeta o elemento de Raviart-Thomas de alguma forma, diferentemente dos 2 outros métodos testados.

Tabela 16. Erro absoluto máximo de u nos baricentros dos triângulos para t =1, com Pé=0.

L	Mín. quadrados	Raviart-Thomas	P2 de Hermite
8	0,28506 × 10−3	0,36076 × 10−3	0,24068 × 10−3
16	0,94938 × 10−4	0,91001 × 10−4	0,60377 × 10−4
32	0,29559 × 10−4	0,22799 × 10−4	0,15111 × 10−4
64	0,88337 × 10−5	0,57152 × 10−5	0,37801 × 10−5

O objetivo das 2 experiências seguintes é observar o comportamento do esquema de mínimos quadrados com elementos finitos lineares, em simulações com números de Péclet de moderados a elevados. Voltamos então a considerar a equação (1) completa, e comparamos essa técnica com o esquema explícito estudado em [5] , de um tipo ostensivamente explorado por Kawahara et al. (ver p. ex. [4] ), em que unicamente a função incógnita u é representada por elementos finitos lineares por triângulos.

Teste n. 5: Escolhemos o mesmo quadrado unitário como domínio Ω e a mesma discretização espacial do Teste n. 3, mas desta vez com T  = 0, 1 e K  = 10^−kI , para k  = 2 e k  = 5. Tomamos novamente w  = (y  ; x )/2, sendo a solução exata:

$u(x,y,t)=xye^{-t},$

o que corresponde a um número de Péclet igual a ${\textstyle {10}^{k}/{\sqrt {2}}}$ .

Consideramos que Γ₀ é toda a fronteira de Ω, impondo assim condições de Dirichlet não homogéneas sobre uma parte de ∂Ω.

Tomamos valores diferentes de Δt para cada método, como seja, Δt  = 0, 1/M com M  = L /4 para o esquema de mínimos quadrados e M  = 5L para o esquema explícito. Observa-se que, neste último caso, tal valor de Δt é requerido para que o esquema seja estável (cf. [5] ). Nas tabelas 17 e 18 são apresentados os erros relativos da solução numérica medidos nas normas de L² (Ω) e de L^∞ (Ω), respectivamente, para o caso k  = 2 e os valores crescentes de L indicados. As tabelas 19 e 20 correspondem aos mesmos tipos de resultados para k  = 5.

Tabela 17. Erro relativo de u na norma de L² (Ω) para t =0,1 com Pé = 0,7071 × 10² .

L	Mínimos quadrados	Método explícito
8	0,69701 × 10−2	0,16150 × 10−1
16	0,17454 × 10−2	0,46988 × 10−2
32	0,43656 × 10−3	0,24474 × 10−2
64	0,10915 × 10−3	0,12479 × 10−2

Tabela 18. Erro relativo máximo de u em Ω para t =0,1 com Pé = 0,7071 × 10² .

L	Mínimos quadrados	Método explícito
8	0,20766 × 10−3	0,38972 × 10−1
16	0,47550 × 10−4	0,11084 × 10−1
32	0,93816 × 10−5	0,57410 × 10−2
64	0,21402 × 10−5	0,27483 × 10−2

Tabela 19. Erro relativo de u na norma de L² (Ω) para t =0,1 com Pé = 0,7071 × 10⁵ .

L	Mínimos quadrados	Método explícito
8	0.69766 × 10−2	0,92090 × 10−2
16	0,17482 × 10−2	0,48018 × 10−2
32	0,43754 × 10−3	0,24667 × 10−2
64	0,10943 × 10−3	0,12561 × 10−2

Tabela 20. Erro relativo máximo de u em Ω para t =0,1 com Pé = 0,7071 × 10⁵ .

L	Mínimos quadrados	Método explícito
8	0,22926 × 10−3	0,54364 × 10−2
16	0,66230 × 10−4	0,23567 × 10−2
32	0,20337 × 10−4	0,10837 × 10−2
64	0,58799 × 10−5	0,50561 × 10−3

Como se pode depreender dos resultados acima, na simulação da convecção-difusão, o método de mínimos quadrados é capaz de produzir resultados convergentes com a mesma ordem do que no caso da difusão pura ou dominante, e isto mesmo a altos números de Péclet. No entanto, observa-se que, neste exemplo, está ausente qualquer camada limite, o que não ocorre no exemplo seguinte.

Teste n. 6: Escolhemos novamente o quadrado unitário como domínio Ω e tomamos ainda K  = νI com ν  = 10^−k para k  = 5, T =0,1 e a mesma discretização espacial do que no Teste n. 5. No entanto, agora w  = (1/π  ; 1/π ), o que implica que o número de Péclet seja igual a 10^k /π . Na realidade, vamos considerar o mesmo problema testado em [5] com o esquema explícito mencionado no Teste n. 5, ou seja, com uma solução exata possuindo uma dupla camada limite na vizinhança das abcissas x  = 1 e y  = 1. Mais especificamente, tomamos

$u(x,y,t)=e^{-t}[s(x)\quad sen(\pi y)+s(y)\quad sen(\pi x)]$

onde para z  ∈ [0, 1],

$s(z):={\frac {[e^{\frac {{\pi }^{2}\nu -1}{\pi \nu }}-1]e^{\pi z}+(e^{\pi }-1)e^{\frac {\left({\pi }^{2}\nu -1)(1-z\right)}{\pi \nu }}}{e^{\frac {2{\pi }^{2}\nu -1}{\pi \nu }}-1}}.$

A equação (1) possui efetivamente a solução acima, desde que tenhamos f (x , y , t ) = e^−t [s (x ) cos(πy ) + s (y ) cos(πx )], e condições iniciais e de contorno de Dirichlet compatíveis com u .

Como se pode observar nas tabelas 21 e 22 , o método de mínimos quadrados é capaz de simular corretamente camadas limite mais espessas, isto é, para um número de Péclet moderado da ordem de 10² . No entanto, falhou completamente no caso de uma camada limite mais fina com Pé da ordem de 10⁵ , como se deduz claramente das tabelas 23 e 24 . A propósito, vale a pena comentar que, com base nesses resultados extraídos de [5] , ainda assim, o método explícito comportou-se satisfatoriamente neste exemplo. Esse comentário aplica-se em particular ao caso de Pé elevado, levando em conta que os valores máximos do erro ocorrem justamente dentro da camada limite, totalmente fora do alcance de qualquer método, pelo menos para o grau de refinamento de malhas utilizado nos testes acima. Em resumo, em virtude dos grandes erros em normas do máximo e L 2, parece ser problemática a captação de camadas limite muito estreitas com o método de mínimos quadrados estudado neste trabalho. A notar que essa dificuldade foi constatada apesar do facto bem conhecido de que o método de mínimos quadrados introduz por si só uma estabilização da convecção dominante, à maneira do método SUPG (cf. [3] ). No entanto, como se sabe, a estabilidade está longe de ser o único critério para se avaliar a qualidade de um método de resolução numérica de equações diferenciais.

Tabela 21. Erro relativo de u na norma de L² (Ω) para t =0,1 com Pé = 0,3183 × 10² .

L	Mínimos quadrados	Método explícito
16	0,74477 × 10−1	0,27550 × 10+0
32	0,20806 × 10−1	0,13929 × 10+0
64	0,53658 × 10−2	0,68628 × 10−1
128	0,13526 × 10−2	0,50245 × 10−1

Tabela 22.

Erro absoluto máximo de u em Ω para t =0,1 com Pé = 0,3183 × 10² .

L	Mínimos quadrados	Método explícito
0,10739 × 10+0	0,84379 × 10+0	16
32	0,33655 × 10−1	0,40427 × 10+0
64	0,78437 × 10−2	0,21973 × 10+0
128	0,19133 × 10−2	0,16351 × 10+0

Tabela 23. Erro relativo de u na norma de L² (Ω) para t  =0,1 com Pé = 0,3183 × 10⁵ .

L	Mínimos quadrados	Método explícito
16	0,18919 × 10+0	0,16682 × 10+0
32	0,20459 × 10+0	0,14754 × 10+0
64	0,27201 × 10+0	0,12560 × 10+0
128	0,29258 × 10+0	0,99616 × 10−1

Tabela 24. Erro absoluto máximo de u em Ω para t =0,1 com Pé = 0,3183 × 10⁵ .

L	Mínimos quadrados	Método explícito
16	0,88493 × 10+0	0,33093 × 10+0
32	0,16663 × 10+1	0,32787 × 10+0
64	0,23959 × 10+1	0,32539 × 10+0
128	0,23286 × 10+1	0,32244 × 10+0

6. Conclusões

Pode-se concluir este trabalho com as seguintes considerações.

Uma primeira observação importante é a de que, na resolução de problemas de difusão transientes, o método misto de mínimos quadrados com elementos finitos lineares é tão ou mais preciso do que métodos de um só campo com capacidade de representar o vetor fluxo, tal o método P2 de Hermite testado neste trabalho, e do que métodos mistos de Galerkin consagrados como o de Raviart-Thomas de mais baixa ordem. Este é um ponto que merece ser enfatizado, ainda mais porque o primeiro método é bem menos custoso em termos computacionais do que os últimos, o que se deve ao facto de os seus graus de liberdade serem definidos unicamente por valores nos vértices dos elementos.

De modo geral, o esquema de mínimos quadrados aqui proposto para tratar a convecção-difusão transiente parece ser amplamente competitivo com relação a outros métodos existentes, pelo memos em condições não muito extremas. Essa expressão encontra todo seu sentido no caso de simulações com alto número de Péclet, as quais parecem exigir modificações no método, de forma a aumentar a sua precisão. Uma possibilidade nesse sentido poderia ser a de combinar esta formulação de mínimos quadrados para a discretização espacial, com o esquema explícito estudado em [5] para a integração no tempo. Um investimento nesse tipo de metodologia deverá ser realizado pelos autores no futuro próximo.

Agradecimentos

O autor Vitoriano Ruas gostaria de agradecer o apoio financeiro concedido pelo CNPq, Brasil, vinculado ao processo de n. 307996/2008-5.

References

1. [1] J. Heinrich, P. Huyakorn, O.C. Zienkiewicz, A. Mitchell; An Upwind Finite Element Scheme for Two-dimensional Convective Transport Equation; Int. J. Num. Meth. Engng., 12 (1978), pp. 187–190
2. [2] M. Baba, M. Tabata; On a Conservative Upwind Finite Element Scheme for Convective-diffusion Equations; RAIRO Analyse Numérique, 15 (1) (1981), pp. 3–25
3. [3] A.N. Brooks, T.J.R. Hughes; The Streamline Upwind/Petrov-Galerkin Formulation for Convection Dominated Flows with Particular Emphasis on the Navier-Stokes Equations; Comp. Meth. Appl. Mechanics & Engng., 32 (1982), pp. 199–259
4. [4] M. Kawahara, H. Hirano; A finite element method for high Reynolds number viscous fluid flow using two step explicit scheme; Int. J. Num. Meth. Fluids, 3 (1983), pp. 137–163
5. [5] V. Ruas, A.P. Brasil Jr., P.R. Trales; An explicit method for convection-diffusion equations; Japan J. Industrial & Appl. Maths., 26 (1–4) (2009), pp. 65–91
6. [6] A.I. Pehlivanov, G.F. Carey, R.D. Lazarov; Leat-squares Mixed Finite Elements for Second-order Elliptic Problems; SIAM J. Numer. Anal., 31 (5) (1994), pp. 1368–1377
7. [7] E. Oñate, R.L. Taylor, O.C. Zienkiewicz, J. Rojek; A residual correction method based on finite calculus; Engineering Computations, 20 (5-6) (2003), pp. 629–658
8. [8] V. Ruas, J.H. Carneiro de Araujo; A quadratic triangle of the Hermite type for second order elliptic problems; ZAMM, 89 (6) (2009), pp. 445–453
9. [9] D.P. Yang; Least-squares Finite Element Methods for Nonlinear Parabolic Problems; J. Comp. Maths., 20 (2) (2002), pp. 153–164
10. [10] C. Novo, R.C. Leal-Toledo, E.M. Toledo e L.S. Martins, Discontinuous Mixed Space-time Least-squares Formulation for Transient Advection-diffusion-reaction Equations. In: Mecánica Computacional, A.Cardona, N.Nigro, V.Sonzogni and M.Storti eds., XXV: 1113-1125, 2006.
11. [11] P.A. Raviart, J.M. Thomas; Mixed Finite Element Methods for Second Order Elliptic Problems; Lecture Notes in Mathematics, Springer Verlag, 606 (1977), pp. 292–315
12. [12] H. Fujita, N. Saito, T. Suzuki; Operator Theory and Numerical Methods; North Holland, N.Y (2001)
13. [13] R.A. Adams; Sobolev Spaces; Academic Press, N.Y (1975)
14. [14] P.G. Ciarlet; The Finite Element Method for Elliplic Problems; North-Holland, Amsterdam (1977)
15. [15] V. Girault, P.A. Raviart; Finite Element Methods for Navier-Stokes Equations; Springer-Verlag, Berlin (1986)
16. [16] R.C.P. Leal-Toledo, V. Ruas; Numerical analysis of a least-squares finite element method for the time-dependent advection-diffusion equation; Journal of Computational and Applied Mathematics, 235 (2011), pp. 3615–3631
17. [17] V. Ruas; Automatic generation of triangular finite element meshes; Comp. & Maths. with Applications, 5 (1979), pp. 125–140