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. [[#bib0005|[1]]] ) e de Baba e Tabata [[#bib0010|[2]]] . O procedimento introduzido no início da década de 80 por Hughes e Brooks [[#bib0015|[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. [[#bib0020|[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 [[#bib0025|[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. [[#bib0030|[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. [[#bib0040|[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 [[#bib0035|[7]]] .      

No que tange a problemas transientes, exceto por um pequeno número de artigos com análises parciais tal como [[#bib0045|[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 [[#bib0050|[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 [[#bib0020|[4]]]  com a modificação de [[#bib0025|[5]]] , do método misto de Raviart-Thomas de mais baixa ordem [[#bib0055|[11]]]  e do método quadrático de Hermite com um só campo introduzido em [[#bib0040|[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 <math display="inline">R^N</math> , ''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:

<span id='eq0005'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\begin{array}{l}
{\partial }_tu-\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}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 1)
|}

onde '''''ν'''''  é o vetor normal externo unitário sobre ∂Ω, Γ<sub>0</sub>  e Γ<sub>1</sub>  são 2 porções disjuntas de ∂Ω, <math display="inline">{\partial }_tG</math>  representa derivada parcial de primeira ordem de uma função escalar ou vetorial <math display="inline">G</math> , e ∇ designa o operador gradiente. A medida de Γ<sub>1</sub>  pode ser nula mas a de Γ<sub>0</sub>  é 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 <math display="inline">R^N</math> . 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 <math display="inline">{\sum }_{i=1}^Na_i\frac{\partial }{\partial x_i}</math> .      

Em [[#eq0005|(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 <math display="inline">[C^m(\overline{\Omega })]^N</math>  para certo <math display="inline">m\in {\mathbb{N}}^\ast </math> . Supomos ainda que '''w'''  satisfaz a condição '''w'''  · '''''ν'''''  ≥ 0 sobre Γ<sub>1</sub>  ∀ ''t'' . Nesse sentido, Γ<sub>1</sub>  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''<sup>∞</sup> [(0, ''T'' ) ; ''L''<sup>2</sup> (Ω)] (cf. [[#bib0060|[12]]] ) e um valor imposto sobre Γ<sub>0</sub>  × (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''<sup>0</sup>  ∈ ''H''<sup>1</sup> (Ω) (cf. [[#bib0065|[13]]] ). Vamos ainda requerer alguma regularidade adicional tanto para ''u''<sup>0</sup>  como para ''f'' , a ser especificada mais adiante.      

Em seguida, definimos 2 espaços de Hilbert para as normas naturais (cf. [[#bib0070|[14]]]  e [[#bib0075|[15]]] ), os quais vão representar um papel fundamental neste estudo, a saber, <math display="inline">V:=\lbrace v/v\in H^1(\Omega ),v_{/{\Gamma }_0}=</math><math>0\rbrace </math>  e <math display="inline">\boldsymbol{\mbox{Q}}:=\lbrace \boldsymbol{\mbox{q}}/\boldsymbol{\mbox{q}}\in \boldsymbol{\mbox{H}}(\Omega ,div),\boldsymbol{\mbox{q}}_{/{\Gamma }_1}\cdot \boldsymbol{\nu }=</math><math>0\rbrace </math>  (cf. [[#bib0075|[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  [[#eq0005|(1)]]  na forma variacional mista equivalente de tipo mínimos quadrados [[#eq0010|(2)]]  abaixo, onde (''A'' , ''B'' ) designa o produto interno padrão de 2 funções escalares ou vetoriais ''A''  e ''B''  de ''L''<sup>2</sup> (Ω), e ('''A''' , '''B''' )<sub>'''K'''</sub>  representa ('''KA''' , '''B''' ), sendo '''A'''  e '''B'''  2 funções vetoriais. Assim, sendo ''I''  o operador identidade, temos:

<span id='eq0010'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\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 }_tu(\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 }_tu+\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}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 2)
|}

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

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\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</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 3)
|}

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

'''Observação 1:'''  Caso a equação [[#eq0005|(1)]]  esteja escrita na forma adimensional, ''λ''<sup>−1</sup>  representa o número de Péclet.   □      

No restante deste artigo, denotamos por ∥ ·  ∥ <sub>''m'' ,''p'' ,Ω                            </sub>  a norma padrão do espaço de Sobolev ''W''<sup>''m'' ,''p''</sup> (Ω) (cf. [[#bib0065|[13]]] ) com <math display="inline">m\in \mathbb{N}</math>  e <math display="inline">p\in R</math> , ''p''  ≥ 1. ''W''<sup>''m'' ,2                            </sup> (Ω) é representado por ''H''<sup>''m''</sup> (Ω) para ''m''  ≠ 0. A norma padrão de ''H''<sup>''m''</sup> (Ω) é simplesmente denotada por ∥ ·  ∥ <sub>''m'' ,Ω                            </sub>  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 [[#eq0005|(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''<sup>''r''</sup>  = ''f'' (''r'' Δ''t'' ) para todo ''r''  ∈ [0, ''M'' ] e partindo de ''u''<sup>0</sup>  e '''p'''<sup>0</sup> , determinamos sucessivamente para ''n''  = 1, 2, …, ''M''  uma aproximação de (''u'' [''n'' Δ''t'' ] ; '''p''' [''n'' Δ''t'' ]) designada por (''u''<sup>''n''</sup>  ; '''p'''<sup>''n''</sup> ), como solução do problema de mínimos quadrados seguinte, em tudo análogo a [[#eq0010|(2)]] :

<span id='eq0020'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\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}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 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 [[#eq0020|(4)]]  definido como se segue:      

Seja <math display="inline">T_h</math>  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 <math display="inline">T_h</math>  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 Γ<sub>0</sub>  como Γ<sub>1</sub>  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 <math display="inline">T_h</math> . Se <math display="inline">T_h</math>  é constituído de ''N  '' -simplex, para cada <math display="inline">E\in T_h</math>''R''<sub>''k''</sub> (''E'' ) é definido como o espaço de polinómios de grau menor ou igual a ''k''  e, senão ''R''<sub>''k''</sub> (''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 <math display="inline">\overset{\mbox{ˆ}}{E}</math>  conforme o caso, pela aplicação ''N  '' -linear bijetiva que mapeia <math display="inline">\overset{\mbox{ˆ}}{E}</math>  em ''E  '' . Dessa forma, para todo <math display="inline">k\in {\mathbb{N}}^\ast </math> , introduzimos os seguintes espaços vetoriais associados a <math display="inline">T_h</math> :

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>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\},</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | 
|}

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\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 ,</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | 
|}

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>V_h:=S_{h,j}\cap V\quad para\quad j\in {\mathbb{N}}^\ast .</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | 
|}

Sejam agora <math display="inline">u_h^0</math>  e <math display="inline">\boldsymbol{\mbox{p}}_h^0</math>  os interpolantes padrão de ''u''<sup>0</sup>  em ''V''<sub>''h''</sub>  e de '''p'''<sup>0</sup>  em '''Q'''<sub>''h''</sub> , 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''<sup>''r''</sup>  pelos seus interpolantes padrão '''w'''<sub>''h''</sub>  e <math display="inline">f_h^r</math>  em [''S''<sub>''h'' ,''l'' +1                            </sub> ]<sup>''N''</sup>  e ''S''<sub>''h'' ,''j''</sub> , respetivamente. Com isso, a nossa aproximação espacial e temporal do problema [[#eq0005|(1)]]  na forma mista formula-se como se segue:

<span id='eq0040'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\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}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 5)
|}

Uma simples aplicação do Teorema de Lax-Milgram, permite afirmar que o problema [[#eq0040|(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 [[#enun0005|Teorema 1]]  abaixo os resultados de estabilidade que se aplicam ao esquema [[#eq0040|(5)]] . Com esse fim, lembramos primeiramente que ∇ · '''w'''  = 0. Assim, supondo que vale o resultado de regularidade '''w'''  ∈ [''H''<sup>''l'' +2                            </sup> (Ω)]<sup>''N''</sup> , onde ''m''  é o maior número inteiro tal que 1 ≤ ''m''  < ''l''  + 2 − ''N  '' /2, temos <math display="inline">\boldsymbol{\mbox{w}}\in [C^m(\overline{\Omega })]^N</math>  (cf. [[#bib0065|[13]]] ). Então, de acordo com resultados de aproximação clássicos (cf. [[#bib0070|[14]]] ), conclui-se facilmente que existe uma constante ''C''<sub>'''w'''</sub>  tal que,

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>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 }.</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 6)
|}

Então, sendo ''C''<sub>''l''</sub>  uma constante para a qual vale ||'''w'''  − '''w'''<sub>''h''</sub> ||<sub>0,∞,Ω</sub>  ≤ (''C''<sub>''l''</sub>  − 1)||'''w''' ||<sub>0,∞,Ω</sub>  (cf. [[#bib0070|[14]]] ), temos:

<span id='enun0005'></span>

==Teorema 1.                     ==

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

<span id='eq0050'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\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}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 7)
|}

□      

Em seguida, fechamos esta secção dando um resultado de convergência para o esquema [[#eq0040|(5)]]  nas normas naturais delineadas em [[#eq0050|(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''<sub>Δ</sub>''h''<sup>''τ''</sup> .      

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

<span id='enun0010'></span>

==Teorema 2.                     ==

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

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\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}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 8)
|}

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

□      

A demonstração dos Teoremas [[#enun0005|1]]  e [[#enun0010|2]]  pode ser encontrada em [[#bib0080|[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 [[#eq0005|(1)]]  é Ω × (0, ''T'' ), onde Ω é o quadrado unitário (0, 1)<sup>2</sup>  e ''T''  = 1. Apresentamos abaixo alguns resultados relevantes para o caso em que Γ<sub>0</sub>  é 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 <math display="inline">\sqrt{2}/2</math> . Consideramos uma solução analítica exata dada por:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\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}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | 
|}

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

<span id='eq0065'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>f:=u_t+\nabla \cdot \boldsymbol{\mbox{p}}-\boldsymbol{\mbox{w}}\cdot \boldsymbol{\mbox{p}}.</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 9)
|}

Resolvemos este problema com malhas uniformes triangulares obtidas mediante uma primeira subdivisão de Ω em ''L''<sup>2</sup>  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  [[#tbl0005|1]] , [[#tbl0010|2]]  e [[#tbl0015|3]]  abaixo os erros absolutos da aproximação para valores crescentes de ''L'' , tomando ''M''  = 2 × ''L'' , i.e. Δ''t''  = ''h'' /2. Na tabela  [[#tbl0005|1]]  fornecemos erros na norma de ''L''<sup>2</sup> (Ω) dos valores approximados de ''u'' , ∇''u'' , '''p''' , ∇ · '''p'''  e ''u''<sub>''t''</sub>  para ''t''  = 1. Na tabela [[#tbl0010|2]]  provemos os erros de '''p''' , ∇ · '''p'''  e ''u''<sub>''t''</sub>  na norma de ''L''<sup>2</sup> [(0, ''T'' ) ; ''L''<sup>2</sup> (Ω)], em versão tanto discreta como contínua, dependendo do campo (cf. secção 4). Em seguida, na tabela [[#tbl0015|3]] , fornecemos os erros absolutos de ''u''  na norma de ''L''<sup>∞</sup> (Ω) para ''t''  =0,5 e ''t'' =1.

<span id='tbl0005'></span>

{| class="wikitable" style="min-width: 60%;margin-left: auto; margin-right: auto;"
|+

Tabela 1.

Erros absolutos de ''u'' , ∇''u'' , '''p''' , ∇ · '''p'''  e ''u''<sub>''t''</sub>  na norma de ''L''<sup>2</sup> (Ω) para ''t'' =1, com Pé=0,7071.                  

|-

! ''L''
! ''u''
! ∇''u''
! '''p'''
! ∇ · '''p'''
! ''u''<sub>''t''</sub>
|-

| 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
|}

<span id='tbl0010'></span>

{| class="wikitable" style="min-width: 60%;margin-left: auto; margin-right: auto;"
|+

Tabela 2.

Erros absolutos de '''p''' , ∇ · '''p'''  e ''u''<sub>''t''</sub>  na norma de ''L''<sup>2</sup> [(0, 1) ; ''L''<sup>2</sup> (Ω)], com Pé=0,7071.                  

|-

! ''L''
! '''p'''
! ∇ · '''p'''
! ''u''<sub>''t''</sub>
|-

| 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
|}

<span id='tbl0015'></span>

{| class="wikitable" style="min-width: 60%;margin-left: auto; margin-right: auto;"
|+

Tabela 3.

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

|-

! ''L''
! ∥ [''u''  − ''u''<sub>''h''</sub> ](0, 5) ∥ <sub>0,∞,Ω</sub>
! ∥ [''u''  − ''u''<sub>''h''</sub> ](1) ∥ <sub>0,∞,Ω</sub>
|-

| 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 [[#tbl0005|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''<sup>2</sup>  dos erros de ''u'' , '''p''' , ∇ · '''p'''  e ''u''<sub>''t''</sub>  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 [[#tbl0010|2]]  indica taxas de convergência na norma de ''L''<sup>2</sup> [(0 ; 1) ; ''L''<sup>2</sup> (Ω)] de cerca de 2, 1 e 2 para '''p''' , ∇ · '''p'''  e ''u''<sub>''t''</sub> , 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 [[#tbl0015|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''<sup>∞</sup> [(0 ; 1) ; ''L''<sup>∞</sup> (Ω)] 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''<sup>2</sup>  ; ''y'' /''r''<sup>2</sup> ), com <math display="inline">r=\sqrt{x^2+y^2}</math>  e novamente ''K''  = ''I'' . A equação  [[#eq0005|(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 Γ<sub>0</sub>  é 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:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\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}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | 
|}

e, para isso, tomamos de novo a função de força ''f''  expressa por  [[#eq0065|(9)]] , o que, neste caso particular, se reduz a fazer ''f''  = ''u''<sub>''t''</sub> .      

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 [[#tbl0005|1]] , [[#tbl0010|2]]  e [[#tbl0015|3]] , apresentamos nas tabelas [[#tbl0020|4]] , [[#tbl0025|5]]  e [[#tbl0030|6]]  abaixo os erros absolutos das aproximações determinadas para valores crescentes de ''L'' , tomando novamente ''M''  = 2 × ''L'' .

<span id='tbl0020'></span>

{| class="wikitable" style="min-width: 60%;margin-left: auto; margin-right: auto;"
|+

Tabela 4.

Erros absolutos de ''u'' , ∇''u'' , '''p''' , ∇ · '''p'''  e ''u''<sub>''t''</sub>  na norma de ''L''<sup>2</sup> (Ω) para ''t'' = 1, com Pé = 10.                  

|-

! ''L''
! ''u''
! ∇''u''
! '''p'''
! ∇ · '''p'''
! ''u''<sub>''t''</sub>
|-

| 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
|}

<span id='tbl0025'></span>

{| class="wikitable" style="min-width: 60%;margin-left: auto; margin-right: auto;"
|+

Tabela 5.

Erros absolutos de '''p''' , ∇ · '''p'''  e ''u''<sub>''t''</sub>  na norma de ''L''<sup>2</sup> [(0, 1) ; ''L''<sup>2</sup> (Ω)], com Pé=10.                  

|-

! ''L''
! '''p'''
! ∇ · '''p'''
! ''u''<sub>''t''</sub>
|-

| 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
|}

<span id='tbl0030'></span>

{| class="wikitable" style="min-width: 60%;margin-left: auto; margin-right: auto;"
|+

Tabela 6.

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

|-

! ''L''
! ∥ [''u''  − ''u''<sub>''h''</sub> ](0, 5) ∥ <sub>0,∞,Ω</sub>
! ∥ [''u''  − ''u''<sub>''h''</sub> ](1) ∥ <sub>0,∞,Ω</sub>
|-

| 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 [[#tbl0020|4]]  e [[#tbl0025|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 [[#tbl0030|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 [[#eq0005|(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 [[#bib0055|[11]]] , e também com o elemento de Hermite quadrático - ou equivalentemente, P2 de Hermite -, com um só campo incógnito, introduzido em [[#bib0040|[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 [[#tbl0035|Tabela 7]] , [[#tbl0040|Tabela 8]] , [[#tbl0045|Tabela 9]]  and [[#tbl0050|Tabela 10]]  e [[#tbl0060|Tabela 12]] , [[#tbl0065|Tabela 13]] , [[#tbl0070|Tabela 14]]  and [[#tbl0075|Tabela 15]] , apresentamos os erros absolutos no tempo final ''T''  escolhido, medidos na norma de ''L''<sup>2</sup> (Ω) e calculados para os 3 métodos, sucessivamente de ''u'' , ∇''u'' , ''u''<sub>''t''</sub>  e Δ''u''  para cada um dos 2 testes desta subsecção. No que tange o gradiente de ''u'' , é o valor do campo '''p'''<sub>''h''</sub>  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'''<sub>''h''</sub>  para esses 2 métodos mistos. Como no caso do elemento de Hermite, o laplaciano de ''u''<sub>''h''</sub>  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'' ).

<span id='tbl0035'></span>

{| class="wikitable" style="min-width: 60%;margin-left: auto; margin-right: auto;"
|+

Tabela 7.

Erros absolutos de ''u''  na norma de ''L''<sup>2</sup> (Ω) 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
|}

<span id='tbl0040'></span>

{| class="wikitable" style="min-width: 60%;margin-left: auto; margin-right: auto;"
|+

Tabela 8.

Erros absolutos de ∇''u''  (ou '''p''' ) na norma de ''L''<sup>2</sup> (Ω) 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
|}

<span id='tbl0045'></span>

{| class="wikitable" style="min-width: 60%;margin-left: auto; margin-right: auto;"
|+

Tabela 9.

Erros absolutos de ''u''<sub>''t''</sub>  na norma de ''L''<sup>2</sup> (Ω) 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
|}

<span id='tbl0050'></span>

{| class="wikitable" style="min-width: 60%;margin-left: auto; margin-right: auto;"
|+

Tabela 10.

Erros absolutos de Δ''u''  (ou ∇ · '''p''' ) na norma de ''L''<sup>2</sup> (Ω) 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
|}

<span id='tbl0060'></span>

{| class="wikitable" style="min-width: 60%;margin-left: auto; margin-right: auto;"
|+

Tabela 12.

Erros absolutos de ''u''  na norma de ''L''<sup>2</sup> (Ω) 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
|}

<span id='tbl0065'></span>

{| class="wikitable" style="min-width: 60%;margin-left: auto; margin-right: auto;"
|+

Tabela 13.

Erros absolutos de ∇''u''  (ou '''p''' ) na norma de ''L''<sup>2</sup> (Ω) 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
|}

<span id='tbl0070'></span>

{| class="wikitable" style="min-width: 60%;margin-left: auto; margin-right: auto;"
|+

Tabela 14.

Erros absolutos de ''u''<sub>''t''</sub>  na norma de ''L''<sup>2</sup> (Ω) 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
|}

<span id='tbl0075'></span>

{| class="wikitable" style="min-width: 60%;margin-left: auto; margin-right: auto;"
|+

Tabela 15.

Erros absolutos de Δ''u''  na norma de ''L''<sup>2</sup> (Ω) 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 [[#tbl0055|11]]  e [[#tbl0080|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.

<span id='tbl0055'></span>

{| class="wikitable" style="min-width: 60%;margin-left: auto; margin-right: auto;"
|+

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. Γ<sub>0</sub>  é a mesma porção de fronteira do que naquele teste, sendo agora a solução exata dada por:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\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}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | 
|}

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

Os resultados da tabela [[#tbl0040|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  [[#tbl0035|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 [[#tbl0055|11]] . Por outro lado, como mostra a tabela [[#tbl0045|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  [[#tbl0050|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 Γ<sub>0</sub> , 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 Γ<sub>0</sub>  é a borda do disco, e a solução exata é dada por:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\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}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | 
|}

Fixamos Δ''t''  em 1/2''L''  para todos os 3 métodos, onde ''L''  é um parâmetro da malha utilizada constituida de 2''L''<sup>2</sup>  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 [[#bib0085|[17]]] .      

Na aproximação do gradiente de ''u'' , a tabela  [[#tbl0065|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''<sup>2</sup> (Ω), como mostra a tabela [[#tbl0060|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 [[#tbl0080|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  [[#tbl0070|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  [[#tbl0075|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 [[#tbl0050|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.

<span id='tbl0080'></span>

{| class="wikitable" style="min-width: 60%;margin-left: auto; margin-right: auto;"
|+

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 [[#eq0005|(1)]]  completa, e comparamos essa técnica com o esquema explícito estudado em [[#bib0025|[5]]] , de um tipo ostensivamente explorado por Kawahara et al. (ver p. ex. [[#bib0020|[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<sup>−''k''</sup>''I'' , para ''k''  = 2 e ''k''  = 5. Tomamos novamente '''w'''  = (''y''  ; ''x'' )/2, sendo a solução exata:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>u(x,y,t)=xye^{-t},</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | 
|}

o que corresponde a um número de Péclet igual a <math display="inline">{10}^k/\sqrt{2}</math> .      

Consideramos que Γ<sub>0</sub>  é 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''  = 5''L''  para o esquema explícito. Observa-se que, neste último caso, tal valor de Δ''t''  é requerido para que o esquema seja estável (cf.  [[#bib0025|[5]]] ). Nas tabelas [[#tbl0085|17]]  e [[#tbl0090|18]]  são apresentados os erros relativos da solução numérica medidos nas normas de ''L''<sup>2</sup> (Ω) e de ''L''<sup>∞</sup> (Ω), respectivamente, para o caso ''k''  = 2 e os valores crescentes de ''L''  indicados. As tabelas  [[#tbl0095|19]]  e [[#tbl0100|20]]  correspondem aos mesmos tipos de resultados para ''k''  = 5.

<span id='tbl0085'></span>

{| class="wikitable" style="min-width: 60%;margin-left: auto; margin-right: auto;"
|+

Tabela 17.

Erro relativo de ''u''  na norma de ''L''<sup>2</sup> (Ω) para ''t'' =0,1 com Pé = 0,7071 × 10<sup>2</sup> .                  

|-

! ''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
|}

<span id='tbl0090'></span>

{| class="wikitable" style="min-width: 60%;margin-left: auto; margin-right: auto;"
|+

Tabela 18.

Erro relativo máximo de ''u''  em Ω para ''t'' =0,1 com Pé = 0,7071 × 10<sup>2</sup> .                  

|-

! ''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
|}

<span id='tbl0095'></span>

{| class="wikitable" style="min-width: 60%;margin-left: auto; margin-right: auto;"
|+

Tabela 19.

Erro relativo de ''u''  na norma de ''L''<sup>2</sup> (Ω) para ''t'' =0,1 com Pé = 0,7071 × 10<sup>5</sup> .                  

|-

! ''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
|}

<span id='tbl0100'></span>

{| class="wikitable" style="min-width: 60%;margin-left: auto; margin-right: auto;"
|+

Tabela 20.

Erro relativo máximo de ''u''  em Ω para ''t'' =0,1 com Pé = 0,7071 × 10<sup>5</sup> .                  

|-

! ''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<sup>−''k''</sup>  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<sup>''k''</sup> /''π'' . Na realidade, vamos considerar o mesmo problema testado em  [[#bib0025|[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

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>u(x,y,t)=e^{-t}[s(x)\quad sen(\pi y)+s(y)\quad sen(\pi x)]</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | 
|}

onde para ''z''  ∈ [0, 1],

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>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}.</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | 
|}

A equação [[#eq0005|(1)]]  possui efetivamente a solução acima, desde que tenhamos ''f'' (''x'' , ''y'' , ''t'' ) = ''e''<sup>−''t''</sup> [''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 [[#tbl0105|21]]  e [[#tbl0110|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<sup>2</sup> . No entanto, falhou completamente no caso de uma camada limite mais fina com Pé da ordem de 10<sup>5</sup> , como se deduz claramente das tabelas [[#tbl0115|23]]  e [[#tbl0120|24]] . A propósito, vale a pena comentar que, com base nesses resultados extraídos de [[#bib0025|[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.  [[#bib0015|[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.

<span id='tbl0105'></span>

{| class="wikitable" style="min-width: 60%;margin-left: auto; margin-right: auto;"
|+

Tabela 21.

Erro relativo de ''u''  na norma de ''L''<sup>2</sup> (Ω) para ''t'' =0,1 com Pé = 0,3183 × 10<sup>2</sup> .                  

|-

! ''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
|}

<span id='tbl0110'></span>

Tabela 22.

Erro absoluto máximo de ''u''  em Ω para ''t'' =0,1 com Pé = 0,3183 × 10<sup>2</sup> .  

{| class="wikitable" style="min-width: 60%;margin-left: auto; margin-right: auto;"

|-

! ''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
|}

<span id='tbl0115'></span>

{| class="wikitable" style="min-width: 60%;margin-left: auto; margin-right: auto;"
|+

Tabela 23.

Erro relativo de ''u''  na norma de ''L''<sup>2</sup> (Ω) para ''t''  =0,1 com Pé = 0,3183 × 10<sup>5</sup> .                  

|-

! ''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
|}

<span id='tbl0120'></span>

{| class="wikitable" style="min-width: 60%;margin-left: auto; margin-right: auto;"
|+

Tabela 24.

Erro absoluto máximo de ''u''  em Ω para ''t'' =0,1 com Pé = 0,3183 × 10<sup>5</sup> .                  

|-

! ''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 [[#bib0025|[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==

<ol style='list-style-type: none;margin-left: 0px;'><li><span id='bib0005'></span>
[[#bib0005|[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</li>
<li><span id='bib0010'></span>
[[#bib0010|[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</li>
<li><span id='bib0015'></span>
[[#bib0015|[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</li>
<li><span id='bib0020'></span>
[[#bib0020|[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</li>
<li><span id='bib0025'></span>
[[#bib0025|[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</li>
<li><span id='bib0030'></span>
[[#bib0030|[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</li>
<li><span id='bib0035'></span>
[[#bib0035|[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</li>
<li><span id='bib0040'></span>
[[#bib0040|[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</li>
<li><span id='bib0045'></span>
[[#bib0045|[9]]] D.P. Yang; Least-squares Finite Element Methods for Nonlinear Parabolic Problems; J. Comp. Maths., 20 (2) (2002), pp. 153–164</li>
<li><span id='bib0050'></span>
[[#bib0050|[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.</li>
<li><span id='bib0055'></span>
[[#bib0055|[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</li>
<li><span id='bib0060'></span>
[[#bib0060|[12]]] H. Fujita, N. Saito, T. Suzuki; Operator Theory and Numerical Methods; North Holland, N.Y (2001)</li>
<li><span id='bib0065'></span>
[[#bib0065|[13]]] R.A. Adams; Sobolev Spaces; Academic Press, N.Y (1975)</li>
<li><span id='bib0070'></span>
[[#bib0070|[14]]] P.G. Ciarlet; The Finite Element Method for Elliplic Problems; North-Holland, Amsterdam (1977)</li>
<li><span id='bib0075'></span>
[[#bib0075|[15]]] V. Girault, P.A. Raviart; Finite Element Methods for Navier-Stokes Equations; Springer-Verlag, Berlin (1986)</li>
<li><span id='bib0080'></span>
[[#bib0080|[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</li>
<li><span id='bib0085'></span>
[[#bib0085|[17]]] V. Ruas; Automatic generation of triangular finite element meshes; Comp. & Maths. with Applications, 5 (1979), pp. 125–140</li>
</ol>