Análise de Estabilidade e de Consistência da Equação de Difusão Anômala com Fluxo Bimodal Empregando o Método de Volumes Finitos

RESUMO.

A análise de consistência e as condições de estabilidade da equação de difusão com fluxo bimodal são tratadas neste trabalho. Foram, também, apresentados diversos detalhes sobre a maneira como a equação diferencial foi discretizada no âmbito do método de Volumes Finitos incluindo a forma que as condições de contorno foram definidas nas equações discretizadas.

Palavras-chave:
difusão anômala; difusão de fluxo bimodal; método de Volumes Finitos

ABSTRACT

Consistency analysis and stability conditions of the bi-flux diffusion equation are discussed in this paper. Several details are presented on how the differential equation has been discretized based on the Finite Volume Method.

Key words
annomalous diffusion; bi-flux diffusion; Finite Volume method

1 INTRODUÇÃO

A equação da difusão descreve o movimento de matéria, momento ou energia em um meio sujeito a gradientes de matéria, momento ou energia, respectivamente^[11]. Para um problema unidimensional tal equação pode ser escrita como:

$\frac{\mathrm{\partial }\varphi }{\mathrm{\partial }t}=\frac{\mathrm{\partial }}{\mathrm{\partial }x}\left(\mathrm{\Gamma }\frac{\mathrm{\partial }\varphi }{\mathrm{\partial }x}\right)$ (1.1)

em que $\Gamma$ é o coeficiente de difusão.

O modelo matemático classicamente utilizado na modelagem de processos de difusão tem como característica básica a escala linear do deslocamento médio das partículas, $r$, com o tempo, isto é, $<r^2>\approx t^{\gamma }$, com $\gamma =1$ ^[23]. Para aquelas circunstâncias em que $\gamma \ne 1$ têm-se o que se denomina de difusão anômala. Existe uma infinidade de modelos de difusão anômala na literatura^[5],^[9],^[10],^[14],^[24],^[25], cada qual com suas características, equações e métodos apropriados para estudá-las. Neste trabalho consideraremos apenas a equação de difusão anômala com distribuição de fluxo bimodal, ou equação de difusão bi-fluxo, proposta por Bevilacqua e colaboradores^[3],^[4], definida pela equação:

$\frac{\partial \varphi }{\partial t}=\xi {\Gamma }_2\frac{\partial [}{\partial 2}]\varphi x-\xi (1-\xi ){\Gamma }_4\frac{\partial [}{\partial 4}]\varphi x$(1.2)

na qual $0<\xi \le 1$, ${\Gamma }_2>0$ e ${\Gamma }_4>0$ são, respectivamente, os coeficientes de difusão primária e secundária. O modelo de fluxo bimodal considera que se uma fração $(1-\xi )$ das partículas em difusão se atrasar em seu movimento por causa de alguma iteração mecânica, biológica, física, química ou físico-química com o meio, este atraso é definido pelo termo envolvendo a derivada de quarta ordem^[8]. Desta forma, $\xi$ fornece a fração da difusão primária neste modelo. Para $\xi =1$ esta equação passa a ser a equação da difusão clássica, Eq. (1.1), e para $\xi =0$ tem-se a equação estacionária.

Apesar de ser pouco conhecida, não faz parte do escopo deste trabalho apresentar como a Eq. (1.2) foi obtida ou mesmo questões sobre a sua física. Aspectos importantes sobre a sua obtenção foram tratados em Bevilacqua et al.^[3],^[4] e em Jiang^[8].

Enquanto a Eq. (1.1), ou a sua versão em duas ou três dimensões, já foram estudadas utilizando os mais diversos métodos de soluções de equações diferenciais, o mesmo não se pode dizer da Eq. (1.2). Existem relativamente poucas publicações^[13],^[19],^[20] que lidaram com esta equação, ou a sua versão em duas dimensões e, até onde os autores deste presente trabalho puderam determinar, não foram encontradas publicações em que se discutiu questões de estabilidade e convergência dos métodos de solução da equação de difusão de bi-fluxo, Eq. (1.2).

O propósito deste trabalho é apresentar uma formulação consistente da Eq. (1.2) utilizando o método de Volumes Finitos^[12],^[15] para discretizá-la. Também será apresentada uma análise de estabilidade para verificar em que condições esta equação é estável.

Por ser uma equação relativamente nova e por possuir um termo com uma derivada de 4ª ordem, derivada esta que não é usualmente encontrada em equações de fenômenos de transporte, é importante realizar análises de estabilidade e consistência para se garantir a acurácia dos resultados obtidos. Por exemplo, até o momento os trabalhos apresentados utilizam soluções totalmente implícitas^[21] ou semi-implícitas^[17] na modelagem do problema transiente sem garantia alguma que esta aproximação seja incondicionalmente estável como é a solução quando se considera somente o clássico modelo de difusão.

Para tanto este artigo está organizado da seguinte forma: na Seção 2 é apresentada a obtenção dos coeficientes da equação unidimensional discretizada, incluindo as condições de contorno. Na Seção 3 é realizada a análise de consistência e na Seção 4 a análise de estabilidade das discretizações propostas. Por fim, na Seção 5, são resumidas as principais conclusões obtidas neste trabalho.

2 EQUAÇÃO UNIDIMENSIONAL

Nesta seção será apresentada a discretização da Eq. (1.2) através do Método de Volumes Finitos. Para tanto, tal equação é inicialmente reescrita como:

$\frac{\partial \varphi }{\partial t}={\lambda }_2\frac{\partial [}{\partial 2}]\varphi x-{\lambda }_4\frac{\partial [}{\partial 4}]\varphi x$(2.1)

em que ${\lambda }_2=\xi {\Gamma }_2>0$ e ${\lambda }_4=\xi (1-\xi ){\Gamma }_4\ge 0$. O domínio, de comprimento $L$, é definido por $\{x\in ℝ|0<x<L\}$ com $L>0$. Esta equação está sujeita a uma condição inicial, $\varphi (x,t=0)=\Phi (x)$, a duas condições de contorno válidas em $x=0$ e a duas válidas em $x=L$.

Equações envolvendo derivadas de até 3ª ordem poderiam ser empregadas para definir as condições de contorno da Eq. (1.2). Contudo, neste trabalho nos restringiremos as condições de contorno que usualmente são adotadas em problemas de difusão, ou seja, as condições de Dirichlet, Neumann e Robin.

2.1 Discretização da Equação

A seguir é apresentada a integração de cada um dos termos da Eq. (2.1) no espaço e no tempo em uma malha unidimensional com espaçamento uniforme de tamanho $h$. A integração do termo envolvendo a derivada no tempo nesta malha assume a seguinte forma:

${\int }_{x_w}^{x_e}{\int }_t^{t+\Delta t}\frac{\partial \varphi }{\partial t}dtdx=({\varphi }_P-{\varphi }_P^0)h$(2.2)

em que $x_P$ é a coordenada do centro do $P$-ésimo volume sob análise, $x_w=x_P-h/2$ e $x_e=x_P+h/2$ são, respectivamente, as coordenadas da face esquerda e direita do volume finito típico, $h=x_e-x_w$, ${\varphi }_P=\varphi (x_P,t+\Delta t)$ e ${\varphi }_P^0=\varphi (x_P,t)$.

A integração do termo envolvendo a derivada de 2ª ordem foi realizada segundo a técnica usual do Método de Volumes Finitos^[12],^[15] da seguinte forma:

${\int }_t^{t+\Delta t}{\int }_{x_w}^{x_e}\frac{\partial [}{\partial 2}]\varphi xdxdt=\frac{\partial {\varphi }^{\theta }}{\partial x}|{}_{x_w}^{x_e}\Delta t$(2.3)

em que ${\varphi }^{\theta }=\theta {\varphi }^{t+\Delta t}+(1-\theta ){\varphi }^t$, com $0\le \theta \le 1$. Para $\theta =1$ tem-se a formulação totalmente implícita, para $\theta =0$, a formulação explícita^[12].

Por fim, o termo de 4ª ordem é integrado da seguinte forma:

${\int }_t^{t+\Delta t}{\int }_{x_w}^{x_e}\frac{\partial [}{\partial 4}]\varphi xdxdt=\frac{{\partial }^3{\varphi }^{\theta }}{\partial x^3}|{}_{x_w}^{x_e}\Delta t$ (2.4)

Substituindo os resultados das integrações definidas anteriormente, chega-se a equação básica do método de volumes finitos para a Eq. (2.1):

$({\varphi }_P-{\varphi }_P^0)\frac{h}{\Delta t}={\lambda }_2[\frac{\partial {\varphi }^{\theta }}{\partial x}|{}_{x_e}-\frac{\partial {\varphi }^{\theta }}{\partial x}\left|{}_{x_w}\right]-{\lambda }_4[\frac{{\partial }^3{\varphi }^{\theta }}{\partial x^3}|{}_{x_e}-\frac{{\partial }^3{\varphi }^{\theta }}{\partial x^3}\left|{}_{x_w}\right]$(2.5)

Para o caso unidimensional, ${\lambda }_2\frac{\partial *}{\partial \varphi }x$ representa o fluxo primário, e ${\lambda }_4\frac{\partial *}{\partial [}3]\varphi x$ o fluxo secundário ou como denominado nos primeiros trabalhos sobre esta equação^[4], fluxo de retenção.

2.2 Volumes Internos

Para os volumes internos da malha, as derivadas de primeira ordem obtidas da integração podem ser aproximadas por:

$\begin{array}{cc}\hfill \frac{\partial {\varphi }^{\theta }}{\partial x}|{}_{x_w}& =\frac{9}{8}\frac{{\varphi }_P^{\theta }-{\varphi }_W^{\theta }}{h}+\frac{1}{24}\frac{{\varphi }_{WW}^{\theta }-{\varphi }_E^{\theta }}{h}+𝒪\hfill \\ \hfill \frac{\partial {\varphi }^{\theta }}{\partial x}|{}_{x_e}& =\frac{9}{8}\frac{{\varphi }_E^{\theta }-{\varphi }_P^{\theta }}{h}+\frac{1}{24}\frac{{\varphi }_W^{\theta }-{\varphi }_{EE}^{\theta }}{h}+𝒪\hfill \end{array}$(2.6)

em que ${\varphi }_W=\varphi (x_P-h)$, ${\varphi }_E=\varphi (x_P+h)$, ${\varphi }_{WW}=\varphi (x_P-2h)$ e ${\varphi }_{EE}=\varphi (x_P+2h)$. A aproximação aqui adotada para as derivadas parciais de 1ª ordem não são as encontradas usualmente neste tipo de discretização. Vasconcellos et al.^[22] demonstraram que a aproximação definida pela Eq. (2.6) é mais adequada que a aproximação de 2ª ordem tradicional, considerando que as aproximações das derivadas de 3ª ordem usarão os mesmos volumes empregados na Eq. (2.6), como exposto a seguir.

As derivadas parciais de 3ª ordem são aproximadas por:

$\begin{array}{cc}\hfill \frac{{\partial }^3{\varphi }^{\theta }}{\partial x^3}|{}_{x_w}& =3\frac{{\varphi }_W^{\theta }-{\varphi }_P^{\theta }}{h^3}+\frac{{\varphi }_E^{\theta }-{\varphi }_{WW}^{\theta }}{h^3}+𝒪\hfill \\ \hfill \frac{{\partial }^3{\varphi }^{\theta }}{\partial x^3}|{}_{x_e}& =3\frac{{\varphi }_P^{\theta }-{\varphi }_E^{\theta }}{h^3}+\frac{{\varphi }_{EE}^{\theta }-{\varphi }_W^{\theta }}{h^3}+𝒪\hfill \end{array}$ (2.7)

Após substituir as aproximações das derivadas na Eq. (2.5) chega-se a equação discretizada válida para todos os volumes internos escrita na forma usual do Método de Volumes Finitos^[12],^[15]. Para $\theta =1$, método totalmente implícito, tem-se:

$A_P{\varphi }_P=A_W{\varphi }_W+A_E{\varphi }_E+A_{WW}{\varphi }_{WW}+A_{EE}{\varphi }_{EE}+S_P$(2.8)

e os seguintes coeficientes:

$A_P=\frac{h}{\Delta t}+\frac{9}{4}\frac{{\lambda }_2}{h}+6\frac{{\lambda }_4}{h^3},S_P=\frac{h}{\Delta t}{\varphi }_P^0$ (2.9)

e

$A_W=A_W=\frac{7}{6}\frac{{\lambda }_2}{h}+4\frac{{\lambda }_4}{h^3},A_{WW}=A_{EE}=\frac{1}{24}\frac{{\lambda }_2}{h}-\frac{{\lambda }_4}{h^3}$ (2.10)

Pode-se observar nos coeficientes da Eq. (2.9) e (2.10) que nem sempre se tem que $\left|A_P\right|>\left|A_W\right|+\left|A_E\right|+\left|A_{WW}\right|+\left|A_{EE}\right|$ o que torna este sistema de equações lineares inadequado para ser resolvido por muitos algoritmos que requerem que a matriz seja do tipo diagonal dominante^[6],^[18].

Para o caso em que $\theta =0$, método totalmente explícito, tem-se que:

$A_P{\varphi }_P=A_W{\varphi }_W^0+A_E{\varphi }_E^0+A_{WW}{\varphi }_{WW}^0+A_{EE}{\varphi }_{EE}^0+S_P$(2.11)

e

$A_P=\frac{h}{\Delta t},S_P=[\frac{h}{\Delta t}-\frac{9}{4}\frac{{\lambda }_2}{h}-6\frac{{\lambda }_4}{h^3}]{\varphi }_P^0$(2.12)

e os demais coeficientes são iguais aos da Eq. (2.10).

Esta formulação, como veremos neste trabalho, está sujeita a condições de estabilidade mais rígidas que a conhecida restrição de estabilidade para a clássica equação de difusão de clássica.

Para o problema unidimensional somente os dois primeiros e dois últimos volumes não são considerados volumes internos, uma vez que para estes quatro volumes, não estão definidos alguns dos pontos discretos empregados na discretização dos volumes internos. As equações particulares para as condições de contorno serão definidas na próxima seção.

2.3 Volumes do Contorno

Os dois primeiros e os dois últimos volumes não podem ser escritos segundo a Eq. (2.8), para cada um desses volumes serão escritas equações especiais utilizando as informações sobre as condições de contorno do problema.

Apesar de ser possível utilizar condições de contorno envolvendo derivadas de 2ª ou 3ª para a equação em análise, neste trabalho considerou-se como condições de contorno possíveis as condições de Dirichlet, Neumann e Robin, as quais podem ser resumidas pela equação a seguir:

$\alpha \varphi (x_f)+\beta \frac{\partial \varphi }{\partial x}|{}_{x=x_f}=\gamma$(2.13)

em que $x_f$ é a coordenada do contorno. Desta forma, por exemplo, a primeira das condições de contorno em $x=0$, face esquerda ($x_w$) do primeiro volume, é escrita como:

${\alpha }_w^1\varphi (x=0)+{\beta }_w^1\frac{\partial \varphi }{\partial x}|{}_{x=0}={\gamma }_w^1$(2.14)

e

${\alpha }_w^2\varphi (x=0)+{\beta }_w^2\frac{\partial \varphi }{\partial x}|{}_{x=0}={\gamma }_w^2$(2.15)

em que ${\alpha }_w^1$, ${\beta }_w^1$, ${\gamma }_w^1$ são os coeficientes que definem a primeira condição de contorno em $x=0$ e ${\alpha }_w^2$, ${\beta }_w^2$, ${\gamma }_w^2$ são os coeficientes que definem a segunda condição de contorno em $x=0$. Por exemplo, se a primeira condição de contorno em $x=0$ for $\frac{\partial *}{\partial \varphi }x=0$, condição de Neumann, tem-se ${\alpha }_w^1=0$, ${\beta }_w^1=1$ e ${\gamma }_w^1=0$.

A Eq. (2.14) será utilizada na obtenção dos fluxos da face localizada em $x_w$ do primeiro volume e a Eq. (2.15) será utilizada na obtenção dos fluxos da face localizada em $x_e$ do primeiro volume e na face localizada em $x_w$ do segundo volume.

Existem outras possibilidades de se informar as condições de contorno, como a adoção de volumes fictícios. A adoção desta técnica implicaria em aumentar em mais 4 o número de equações do sistema linear a ser resolvido. Porém, as equações específicas para os volumes fictícios seriam semelhantes as que apresentaremos a seguir.

É importante que se garanta a conservação dos fluxos, isto é, que o fluxo da face $x_e$ de um volume seja o mesmo fluxo da face $x_w$ do volume seguinte. Para todas as faces dos volumes internos isto acontece naturalmente, porém, dependendo de como for formulada as equações daqueles volumes que são influenciados diretamente pelas condições de contorno, isto pode não acontecer.

2.3.1 Equação para o Primeiro Volume

Após realizar uma expansão em série de Taylor chega-se à seguinte fórmula para as derivadas parciais de 1ª ordem:

$\begin{array}{cc}\hfill \frac{\partial {\varphi }^{\theta }}{\partial x}|{}_{x_w}& =\left[\frac{3675}{8}\frac{{\alpha }_w^1}{D_w}\right]{\varphi }_P^{\theta }+[-\frac{1225}{8}\frac{{\alpha }_w^1}{D_w}]{\varphi }_E^{\theta }+\hfill \\ \hfill & \left[\frac{441}{8}\frac{{\alpha }_w^1}{D_w}\right]{\varphi }_{EE}^{\theta }+[-\frac{75}{8}\frac{{\alpha }_w^1}{D_w}]{\varphi }_{EEE}^{\theta }-\frac{352{\gamma }_w^1}{D_w}\hfill \\ \hfill \frac{\partial {\varphi }^{\theta }}{\partial x}|{}_{x_e}& =[-\frac{1085}{8}\frac{{\alpha }_w^2}{D_e}+\frac{1154}{3h}\frac{{\beta }_w^2}{D_e}]{\varphi }_P^{\theta }+[\frac{1015}{8}\frac{{\alpha }_w^2}{D_e}-\frac{402}{h}\frac{{\beta }_w^2}{D_e}]{\varphi }_E^{\theta }+\hfill \\ \hfill & [-\frac{63}{8}\frac{{\alpha }_w^2}{D_e}+\frac{18}{h}\frac{{\beta }_w^2}{D_e}]{\varphi }_{EE}^{\theta }+[\frac{5}{8}\frac{{\alpha }_w^2}{D_e}-\frac{2}{3h}\frac{{\beta }_w^2}{D_e}]{\varphi }_{EEE}^{\theta }+\frac{16{\gamma }_w^2}{D_e}\hfill \end{array}$(2.16)

ao passo que para as derivadas 3ª ordem, obtemos:

$\begin{array}{cc}\hfill & \frac{{\partial }^3{\varphi }^{\theta }}{\partial x^3}|{}_{x_w}=[\frac{1575}{h^2}\frac{{\alpha }_w^1}{D_w}-\frac{1920}{h^3}\frac{{\beta }_w^1}{D_w}]{\varphi }_P^{\theta }+[-\frac{1365}{h^2}\frac{{\alpha }_w^1}{D_w}+\frac{3456}{h^3}\frac{{\beta }_w^1}{D_w}]{\varphi }_E^{\theta }+\hfill \\ \hfill & [\frac{693}{h^2}\frac{{\alpha }_w^1}{D_w}-\frac{1920}{h^3}\frac{{\beta }_w^1}{D_w}]{\varphi }_{EE}^{\theta }+[-\frac{135}{h^2}\frac{{\alpha }_w^1}{D_w}+\frac{384}{h^3}\frac{{\beta }_w^1}{D_w}]{\varphi }_{EEE}^{\theta }-\frac{768}{h^2}\frac{{\gamma }_w^1}{D_w}\hfill \\ \hfill & \frac{{\partial }^3{\varphi }^{\theta }}{\partial x^3}|{}_{x_e}=[\frac{735}{h^2}\frac{{\alpha }_w^2}{D_e}-\frac{784}{h^3}\frac{{\beta }_w^2}{D_e}]{\varphi }_P^{\theta }+[-\frac{525}{h^2}\frac{{\alpha }_w^2}{D_e}+\frac{1200}{h^3}\frac{{\beta }_w^2}{D_e}]{\varphi }_E^{\theta }+\hfill \\ \hfill & [\frac{189}{h^2}\frac{{\alpha }_w^2}{D_e}-\frac{432}{h^3}\frac{{\beta }_w^2}{D_e}]{\varphi }_{EE}^{\theta }+[-\frac{15}{h^2}\frac{{\alpha }_w^2}{D_e}+\frac{16}{h^3}\frac{{\beta }_w^2}{D_e}]{\varphi }_{EEE}^{\theta }-\frac{384}{h^2}\frac{{\gamma }_w^2}{D_e}\hfill \end{array}$(2.17)

com:

$\begin{array}{cc}\hfill D_w& =105{\alpha }_w^{1}h-352{\beta }_w^1\hfill \\ \hfill D_e& =105{\alpha }_w^{2}h-352{\beta }_w^2\hfill \end{array}$(2.18)

Substituindo os fluxos avaliados para as duas faces do primeiro volume na Eq. (2.5) chega-se a equação particular para este volume que é apresentada em sua forma compacta:

$A_P{\varphi }_P=A_E{\varphi }_E^{\theta }+A_{EE}{\varphi }_{EE}^{\theta }+A_{EEE}{\varphi }_{EEE}^{\theta }+S_P$(2.19)

em que ${\varphi }_{EEE}=\varphi (x_P+3h,t+\Delta t)$ e os coeficientes para uma formulação totalmente implícita, $\theta =1$ ficam:

$\begin{array}{cc}\hfill & A_P=\frac{h}{\Delta t}+\hfill \\ \hfill & \frac{{\lambda }_2}{D_P}[62475{\alpha }_w^1{\alpha }_w^2{h}^2-(47740{\alpha }_w^2{\beta }_w^1+202090{\alpha }_w^1{\beta }_w^2)h+\frac{406208}{3}{\beta }_w^1{\beta }_w^2]+\hfill \\ \hfill & \frac{{\lambda }_4}{D_P{h}^2}[-88200{\alpha }_w^1{\alpha }_w^2{h}^2+(-57120{\alpha }_w^2{\beta }_w^1+472080{\alpha }_w^1{\beta }_w^2)h-399872{\beta }_w^1{\beta }_w^2]\hfill \end{array}$(2.20)

$\begin{array}{cc}\hfill A_E& =\frac{{\lambda }_2}{D_P}[29400{\alpha }_w^1{\alpha }_w^2{h}^2-(44660{\alpha }_w^2{\beta }_w^1+96110{\alpha }_w^1{\beta }_w^2)h+\hfill \\ \hfill & 141504{\beta }_w^1{\beta }_w^2]+\frac{{\lambda }_4}{D_P{h}^2}[-88200{\alpha }_w^1{\alpha }_w^2{h}^2+\hfill \\ \hfill & (178080{\alpha }_w^2{\beta }_w^1+354480{\alpha }_w^1{\beta }_w^2)h-794112{\beta }_w^1{\beta }_w^2]\hfill \end{array}$ (2.21)

$\begin{array}{cc}\hfill & A_{EE}=\frac{{\lambda }_2}{D_P}[-6615{\alpha }_w^1{\alpha }_w^2{h}^2+(2772{\alpha }_w^2{\beta }_w^1+21294{\alpha }_w^1{\beta }_w^2)h-6336{\beta }_w^1{\beta }_w^2]+\hfill \\ \hfill & \frac{{\lambda }_4}{D_P{h}^2}[52920{\alpha }_w^1{\alpha }_w^2{h}^2-(135072{\alpha }_w^2{\beta }_w^1+198576{\alpha }_w^1{\beta }_w^2)h+523776{\beta }_w^1{\beta }_w^2]\hfill \end{array}$ (2.22)

$\begin{array}{cc}\hfill & A_{EEE}=\frac{{\lambda }_2}{D_P}[1050{\alpha }_w^1{\alpha }_w^2{h}^2-(220{\alpha }_w^2{\beta }_w^1+3370{\alpha }_w^1{\beta }_w^2)h+\frac{704}{3}{\beta }_w^1{\beta }_w^2]+\hfill \\ \hfill & \frac{{\lambda }_4}{D_P{h}^2}[-12600{\alpha }_w^1{\alpha }_w^2{h}^2+(35040{\alpha }_w^2{\beta }_w^1+45840{\alpha }_w^1{\beta }_w^2)h-129536{\beta }_w^1{\beta }_w^2]\hfill \end{array}$ (2.23)

$\begin{array}{cc}\hfill S_P& =\frac{h}{\Delta t}{\varphi }_P^0+\hfill \\ \hfill & \frac{{\lambda }_2}{D_P}[(36960{\alpha }_w^{2}h-123904{\beta }_w^2)h{\gamma }_w^1+(1680{\alpha }_w^{1}h-5632{\beta }_w^1)h{\gamma }_w^2]+\hfill \\ \hfill & \frac{{\lambda }_4}{D_P{h}^2}[(-80640{\alpha }_w^{2}h+27033{\beta }_w^2)h{\gamma }_w^1+(40320{\alpha }_w^{1}h-135168{\beta }_w^1)h{\gamma }_w^2]\hfill \end{array}$ (2.24)

e, por fim,

$D_P=D_w{D}_{e}h$ (2.25)

2.3.2 Equação para o Segundo Volume

Para escrever uma equação para o segundo volume finito devem-se avaliar os fluxos nas faces $x_w$ e $x_e$ deste volume em função das condições de contorno. O fluxo primário foi avaliado da seguinte maneira:

$\begin{array}{cc}\hfill \frac{\partial {\varphi }^{\theta }}{\partial x}|{}_{x_w}& =[\frac{1015}{8}\frac{{\alpha }_w^2}{D_e}-\frac{402}{h}\frac{{\beta }_w^2}{D_e}]{\varphi }_P^{\theta }+[-\frac{1085}{8}\frac{{\alpha }_w^2}{D_e}+\frac{1154}{3h}\frac{{\beta }_w^2}{D_e}]{\varphi }_W^{\theta }+\hfill \\ \hfill & [-\frac{63}{8}\frac{{\alpha }_w^2}{D_e}+\frac{18}{h}\frac{{\beta }_w^2}{D_e}]{\varphi }_E^{\theta }+[\frac{5}{8}\frac{{\alpha }_w^2}{D_e}-\frac{2}{3h}\frac{{\beta }_w^2}{D_e}]{\varphi }_{EE}^{\theta }+\frac{16{\gamma }_w^2}{D_e}\hfill \\ \hfill \frac{\partial {\varphi }^{\theta }}{\partial x}|{}_{x_e}& =\frac{9}{8}\frac{{\varphi }_E^{\theta }-{\varphi }_P^{\theta }}{h}+\frac{1}{24}\frac{{\varphi }_W^{\theta }-{\varphi }_{EE}^{\theta }}{h}\hfill \end{array}$(2.26)

e o fluxo secundário da seguinte forma:

$\begin{array}{cc}\hfill & \frac{{\partial }^3{\varphi }^{\theta }}{\partial x^3}|{}_{x_w}=[-\frac{525}{h^2}\frac{{\alpha }_w^2}{D_e}+\frac{1200}{h^3}\frac{{\beta }_w^2}{D_e}]{\varphi }_P^{\theta }+[\frac{735}{h^2}\frac{{\alpha }_w^2}{D_e}-\frac{784}{h^3}\frac{{\beta }_w^2}{D_e}]{\varphi }_W^{\theta }+\hfill \\ \hfill & [\frac{189}{h^2}\frac{{\alpha }_w^2}{D_e}-\frac{432}{h^3}\frac{{\beta }_w^2}{D_e}]{\varphi }_E^{\theta }+[-\frac{15}{h^2}\frac{{\alpha }_w^2}{D_e}+\frac{16}{h^3}\frac{{\beta }_w^2}{D_e}]{\varphi }_{EE}^{\theta }-\frac{384}{h^2}\frac{{\gamma }_w^2}{D_e}\hfill \\ \hfill & \frac{{\partial }^3{\varphi }^{\theta }}{\partial x^3}|{}_{x_e}=3\frac{{\varphi }_P^{\theta }-{\varphi }_E^{\theta }}{h^3}+\frac{{\varphi }_{EE}^{\theta }-{\varphi }_W^{\theta }}{h^3}\hfill \end{array}$(2.27)

Em ambas as equações anteriores têm-se que $D_e=105{\alpha }_w^{2}h-352{\beta }_w^2$.

De forma semelhante calculou-se a equação para o segundo volume finito, substituindo os fluxos na Eq. (2.5). Esta equação segue um padrão diferente da Eq. (2.8) e, por este motivo, deve ser reescrita como:

$A_P{\varphi }_P=A_W{\varphi }_W^{\theta }+A_E{\varphi }_E^{\theta }+A_{EE}{\varphi }_{EE}^{\theta }+S_P$(2.28)

em que os coeficientes para uma formulação totalmente implícita, $\theta =1$, são definidos por:

$A_P=\frac{h}{\Delta t}+\frac{{\lambda }_2}{D_{e}h}[245{\alpha }_w^{2}h-798{\beta }_w^2]+\frac{{\lambda }_4}{D_e{h}^2}[840{\alpha }_w^{2}h-2256{\beta }_w^2]$ (2.29)

$A_W=\frac{{\lambda }_2}{D_{e}h}[140{\alpha }_w^{2}h-\frac{1198}{3}{\beta }_w^2]+\frac{{\lambda }_4}{D_e{h}^2}[840{\alpha }_w^{2}h-1136{\beta }_w^2]$ (2.30)

$A_E=\frac{{\lambda }_2}{D_{e}h}[126{\alpha }_w^{2}h-414{\beta }_w^2]+\frac{{\lambda }_4}{D_e{h}^2}[540{\alpha }_w^{2}h-1488{\beta }_w^2]$ (2.31)

$A_{EE}=\frac{{\lambda }_2}{D_{e}h}[-5{\alpha }_w^{2}h+\frac{46}{3}{\beta }_w^2]+\frac{{\lambda }_4}{D_e{h}^2}[-120{\alpha }_w^{2}h+368{\beta }_w^2]$ (2.32)

$S_P=\frac{h}{\Delta t}{\varphi }_P^0+\frac{{\lambda }_2}{D_e}[-16{\gamma }_w^{2}h]+\frac{{\lambda }_4}{D_e}[-384{\gamma }_w^{2}h]$ (2.33)

2.3.3 Equação para o Penúltimo Volume

De modo análogo ao procedimento adotado para obtenção das equações das condições de contorno em $x=0$, Eqs. (2.14) e (2.15), as condições de contorno em $x=L$ seguirão as duas seguintes equações:

${\alpha }_e^1\varphi (x=L)+{\beta }_e^1\frac{\partial \varphi }{\partial x}|{}_{x=L}={\gamma }_e^1$(2.34)

e

${\alpha }_e^2\varphi (x=L)+{\beta }_e^2\frac{\partial \varphi }{\partial x}|{}_{x=L}={\gamma }_e^2$ (2.35)

Em função destas condições de contorno, os fluxos primários nas faces do penúltimo volume foram definidos segundo a seguinte expressão:

$\begin{array}{cc}\hfill \frac{\partial {\varphi }^{\theta }}{\partial x}|{}_{x_w}& =\frac{9}{8}\frac{{\varphi }_P^{\theta }-{\varphi }_W^{\theta }}{h}+\frac{1}{24}\frac{{\varphi }_{WW}^{\theta }-{\varphi }_E^{\theta }}{h}\hfill \\ \hfill \frac{\partial {\varphi }^{\theta }}{\partial x}|{}_{x_e}& =[-\frac{1015}{8}\frac{{\alpha }_w^2}{D_e}-\frac{402}{h}\frac{{\beta }_w^2}{D_e}]{\varphi }_P^{\theta }+[\frac{1085}{8}\frac{{\alpha }_w^2}{D_e}+\frac{1154}{3h}\frac{{\beta }_w^2}{D_e}]{\varphi }_E^{\theta }+\hfill \\ \hfill & [\frac{63}{8}\frac{{\alpha }_w^2}{D_e}+\frac{18}{h}\frac{{\beta }_w^2}{D_e}]{\varphi }_W^{\theta }+[\frac{5}{8}\frac{{\alpha }_w^2}{D_e}-\frac{2}{3h}\frac{{\beta }_w^2}{D_e}]{\varphi }_{WW}^{\theta }+\frac{16{\gamma }_w^2}{D_e}\hfill \end{array}$(2.36)

Já para o fluxo secundário tem-se as seguintes equações:

$\begin{array}{cc}\hfill & \frac{{\partial }^3{\varphi }^{\theta }}{\partial x^3}|{}_{x_w}=[\frac{525}{h^2}\frac{{\alpha }_w^2}{D_e}+\frac{1200}{h^3}\frac{{\beta }_w^2}{D_e}]{\varphi }_P^{\theta }+[-\frac{735}{h^2}\frac{{\alpha }_w^2}{D_e}-\frac{784}{h^3}\frac{{\beta }_w^2}{D_e}]{\varphi }_E^{\theta }+\hfill \\ \hfill & [-\frac{189}{h^2}\frac{{\alpha }_w^2}{D_e}-\frac{432}{h^3}\frac{{\beta }_w^2}{D_e}]{\varphi }_W^{\theta }+[\frac{15}{h^2}\frac{{\alpha }_w^2}{D_e}+\frac{16}{h^3}\frac{{\beta }_w^2}{D_e}]{\varphi }_{EE}^{\theta }+\frac{384}{h^2}\frac{{\gamma }_w^2}{D_e}\hfill \\ \hfill & \frac{{\partial }^3{\varphi }^{\theta }}{\partial x^3}|{}_{x_e}=3\frac{{\varphi }_P^{\theta }-{\varphi }_E^{\theta }}{h^3}+\frac{{\varphi }_{EE}^{\theta }-{\varphi }_W^{\theta }}{h^3}\hfill \end{array}$(2.37)

Nas equações acima tem-se:

$\begin{array}{cc}\hfill D_w& =105{\alpha }_e^{1}h-352{\beta }_e^1\hfill \\ \hfill D_e& =105{\alpha }_e^{2}h-352{\beta }_e^2\hfill \end{array}$ (2.38)

Reunindo as informações acima, pode-se escrever a equação para o penúltimo volume como:

$A_P{\varphi }_P=A_E{\varphi }_E+A_W{\varphi }_W+A_{WW}{\varphi }_{WW}+S_P$(2.39)

em que os coeficientes para uma formulação totalmente implícita, $\theta =1$ são definidos por:

$A_P=\frac{h}{\Delta t}+\frac{{\lambda }_2}{D_{e}h}[245{\alpha }_e^{2}h+798{\beta }_e^2]+\frac{{\lambda }_4}{D_e{h}^2}[840{\alpha }_e^{2}h+2256{\beta }_e^2]$ (2.40)

$A_E=\frac{{\lambda }_2}{D_{e}h}[140{\alpha }_w^{2}h+\frac{1198}{3}{\beta }_w^2]+\frac{{\lambda }_4}{D_e{h}^2}[840{\alpha }_w^{2}h+1136{\beta }_w^2]$ (2.41)

$A_W=\frac{{\lambda }_2}{D_{e}h}[126{\alpha }_w^{2}h+414{\beta }_w^2]+\frac{{\lambda }_4}{D_e{h}^2}[540{\alpha }_w^{2}h+1488{\beta }_w^2]$ (2.42)

$A_{WW}=\frac{{\lambda }_2}{D_{e}h}[-5{\alpha }_w^{2}h-\frac{46}{3}{\beta }_w^2]+\frac{{\lambda }_4}{D_e{h}^2}[-120{\alpha }_w^{2}h-368{\beta }_w^2]$ (2.43)

$S_P=\frac{h}{\Delta t}{\varphi }_P^0+\frac{{\lambda }_2}{D_e}[-16{\gamma }_w^{2}h]+\frac{{\lambda }_4}{D_e}[-384{\gamma }_w^{2}h]$ (2.44)

2.3.4 Equação para o Último Volume

Por fim, os fluxos primários nas faces do último volume foram modelados como:

$\begin{array}{cc}\hfill \frac{\partial {\varphi }^{\theta }}{\partial x}|{}_{x_w}& =[\frac{1085}{8}\frac{{\alpha }_e^2}{D_e}+\frac{1154}{3h}\frac{{\beta }_e^2}{D_e}]{\varphi }_P^{\theta }+[-\frac{1015}{8}\frac{{\alpha }_e^2}{D_e}-\frac{402}{h}\frac{{\beta }_e^2}{D_e}]{\varphi }_W^{\theta }+\hfill \\ \hfill & [\frac{63}{8}\frac{{\alpha }_e^2}{D_e}+\frac{18}{h}\frac{{\beta }_e^2}{D_e}]{\varphi }_{WW}^{\theta }+[-\frac{5}{8}\frac{{\alpha }_e^2}{D_e}-\frac{2}{3h}\frac{{\beta }_e^2}{D_e}]{\varphi }_{WWW}^{\theta }-\frac{16{\gamma }_e^2}{D_e}\hfill \\ \hfill \frac{\partial {\varphi }^{\theta }}{\partial x}|{}_{x_e}& =[-\frac{3675}{8}\frac{{\alpha }_e^1}{D_w}]{\varphi }_P^{\theta }+\left[\frac{1225}{8}\frac{{\alpha }_e^1}{D_w}\right]{\varphi }_W^{\theta }+\hfill \\ \hfill & [-\frac{441}{8}\frac{{\alpha }_e^1}{D_w}]{\varphi }_{WW}^{\theta }+\left[\frac{75}{8}\frac{{\alpha }_e^1}{D_w}\right]{\varphi }_{WWW}^{\theta }+\frac{352{\gamma }_e^1}{D_w}\hfill \end{array}$(2.45)

E os fluxos secundários foram modelados como:

$\begin{array}{cc}\hfill & \frac{{\partial }^3{\varphi }^{\theta }}{\partial x^3}|{}_{x_e}=[-\frac{735}{h^2}\frac{{\alpha }_w^2}{D_e}-\frac{784}{h^3}\frac{{\beta }_w^2}{D_e}]{\varphi }_P^{\theta }+[\frac{525}{h^2}\frac{{\alpha }_w^2}{D_e}+\frac{1200}{h^3}\frac{{\beta }_w^2}{D_e}]{\varphi }_E^{\theta }+\hfill \\ \hfill & [-\frac{189}{h^2}\frac{{\alpha }_w^2}{D_e}-\frac{432}{h^3}\frac{{\beta }_w^2}{D_e}]{\varphi }_{EE}^{\theta }+[\frac{15}{h^2}\frac{{\alpha }_w^2}{D_e}+\frac{16}{h^3}\frac{{\beta }_w^2}{D_e}]{\varphi }_{EEE}^{\theta }+\frac{384}{h^2}\frac{{\gamma }_w^2}{D_e}\hfill \\ \hfill & \frac{{\partial }^3{\varphi }^{\theta }}{\partial x^3}|{}_{x_w}=[-\frac{1575}{h^2}\frac{{\alpha }_w^1}{D_w}-\frac{1920}{h^3}\frac{{\beta }_w^1}{D_w}]{\varphi }_P^{\theta }+[\frac{1365}{h^2}\frac{{\alpha }_w^1}{D_w}+\frac{3456}{h^3}\frac{{\beta }_w^1}{D_w}]{\varphi }_E^{\theta }+\hfill \\ \hfill & [-\frac{693}{h^2}\frac{{\alpha }_w^1}{D_w}-\frac{1920}{h^3}\frac{{\beta }_w^1}{D_w}]{\varphi }_{EE}^{\theta }+[\frac{135}{h^2}\frac{{\alpha }_w^1}{D_w}+\frac{384}{h^3}\frac{{\beta }_w^1}{D_w}]{\varphi }_{EEE}^{\theta }+\frac{768}{h^2}\frac{{\gamma }_w^1}{D_w}\hfill \end{array}$(2.46)

A equação para o último volume também segue um padrão diferente da Eq. (2.8) e é escrita como:

$A_P{\varphi }_P=A_W{\varphi }_W+A_{WW}{\varphi }_{WW}+A_{WWW}{\varphi }_{WWW}+S_P$(2.47)

em que os coeficientes para uma formulação totalmente implícita, $\theta =1$ são:

$\begin{array}{cc}\hfill & A_P=\frac{h}{\Delta t}+\hfill \\ \hfill & \frac{{\lambda }_2}{D_P}[62475{\alpha }_e^1{\alpha }_e^2{h}^2+(47740{\alpha }_e^2{\beta }_e^1+202090{\alpha }_e^1{\beta }_e^2)h+\frac{406208}{3}{\beta }_e^1{\beta }_e^2]+\hfill \\ \hfill & \frac{{\lambda }_4}{D_P{h}^2}[-88200{\alpha }_e^1{\alpha }_e^2{h}^2+(57120{\alpha }_e^2{\beta }_e^1-472080{\alpha }_e^1{\beta }_e^2)h-399872{\beta }_e^1{\beta }_e^2]\hfill \end{array}$(2.48)

$\begin{array}{cc}\hfill & A_W=\frac{{\lambda }_2}{D_P}[29400{\alpha }_e^1{\alpha }_e^2{h}^2+(44660{\alpha }_e^2{\beta }_e^1+96110{\alpha }_e^1{\beta }_e^2)h+141504{\beta }_e^1{\beta }_e^2]+\hfill \\ \hfill & \frac{{\lambda }_4}{D_P{h}^2}[-88200{\alpha }_e^1{\alpha }_e^2{h}^2-(178080{\alpha }_e^2{\beta }_e^1+354480{\alpha }_e^1{\beta }_e^2)h-794112{\beta }_e^1{\beta }_e^2]\hfill \end{array}$ (2.49)

$\begin{array}{cc}\hfill & A_{WW}=\frac{{\lambda }_2}{D_P}[-6615{\alpha }_e^1{\alpha }_e^2{h}^2-(2772{\alpha }_e^2{\beta }_e^1+21294{\alpha }_e^1{\beta }_e^2)h-6336{\beta }_e^1{\beta }_e^2]+\hfill \\ \hfill & \frac{{\lambda }_4}{D_P{h}^2}[52920{\alpha }_e^1{\alpha }_e^2{h}^2+(135072{\alpha }_e^2{\beta }_e^1+198576{\alpha }_e^1{\beta }_e^2)h+523776{\beta }_e^1{\beta }_e^2]\hfill \end{array}$ (2.50)

$\begin{array}{cc}\hfill & A_{WWW}=\frac{{\lambda }_2}{D_P}[1050{\alpha }_e^1{\alpha }_e^2{h}^2+(220{\alpha }_e^2{\beta }_e^1+3370{\alpha }_e^1{\beta }_e^2)h+\frac{704}{3}{\beta }_e^1{\beta }_e^2]+\hfill \\ \hfill & \frac{{\lambda }_4}{D_P{h}^2}[-12600{\alpha }_e^1{\alpha }_e^2{h}^2-(35040{\alpha }_e^2{\beta }_e^1+45840{\alpha }_e^1{\beta }_e^2)h-129536{\beta }_e^1{\beta }_e^2]\hfill \end{array}$ (2.51)

$\begin{array}{cc}\hfill S_P& =\frac{h}{\Delta t}{\varphi }_P^0+\hfill \\ \hfill & \frac{{\lambda }_2}{D_P}[(36960{\alpha }_e^{2}h+123904{\beta }_e^2)h{\gamma }_e^1+(1680{\alpha }_e^{1}h+5632{\beta }_e^1)h{\gamma }_e^2]+\hfill \\ \hfill & \frac{{\lambda }_4}{D_P{h}^2}[(-80640{\alpha }_e^{2}h-27033{\beta }_e^2)h{\gamma }_e^1+(40320{\alpha }_e^{1}h+135168{\beta }_e^1)h{\gamma }_e^2]\hfill \end{array}$ (2.52)

Observa-se que as equações dos dois primeiros volumes são bastante semelhantes, mas não iguais, às equações dos dois últimos volumes. Este procedimento foi intencional e tem o objetivo de facilitar a elaboração do código computacional.

As equações discretizadas da Eq. (2.1), segundo o Método de Volumes Finitos, e que definem um sistema de equações lineares, para o caso de $\theta =1$, foram detalhadas nesta seção. Para resolver este sistema de equações lineares, estes autores indicam o uso de algoritmos mais robustos, como o GMRES^[16] ou outros simulares e a biblioteca PETSC^[1],^[2], a qual possui um grande acervo de metodologias amplamente testadas, como uma boa e eficiente ferramenta computacional.

3 ANÁLISE DE CONSISTÊNCIA

Na Seção 2 foram detalhadas todas as equações necessárias para se resolver numericamente a Eq. (2.1). Obviamente, mesmo no âmbito do Método de Volumes Finitos, a formulação apresentada neste trabalho não é a única forma de discretizar tal equação. Outras equações podem ser propostas, em especial para os dois primeiros e dois últimos volumes. Contudo, nem toda a aproximação adotada para calcular os fluxos, mesmo que correta, levará a uma equação discretizada consistente.

Seja uma equação diferencial qualquer $f=0$ e seja uma discretização desta equação empregando um método numérico qualquer, $\hat{f}=0$. Uma aproximação será dita consistente se $f-\hat{f}\to 0$ quando $limh\to 0$ ^[7]. No caso da equação objeto deste trabalho, as aproximações empregadas serão consistentes se quando $lim\Delta t\to 0$ e $limh\to 0$ as Eqs. (2.8), (2.19), (2.28), (2.39) e (2.47) tenderem a Eq. (2.1)^[27].

A seguir é apresentada a metodologia para a avaliação de consistência da discretização do último volume finito.

Primeiramente, são feitas expansões em série de Taylor em torno de $x_P$ para cada uma das incógnitas da equação. Sendo assim:

$\begin{array}{cc}\hfill {\varphi }_W& ={\varphi }_P-\frac{\partial \varphi }{\partial x}|{}_{{}_{\begin{array}{c}\hfill x=x_P,\hfill \\ \hfill t=t+\Delta t\hfill \end{array}}}h+\frac{1}{2}\frac{\partial [}{\partial 2}]\varphi x|{}_{\begin{array}{c}\hfill x=x_P,\hfill \\ \hfill t=t+\Delta t\hfill \end{array}}{h}^2-\frac{1}{6}\frac{{\partial }^3\varphi }{\partial x^3}|{}_{\begin{array}{c}\hfill x=x_P,\hfill \\ \hfill t=t+\Delta t\hfill \end{array}}{h}^3+\cdots \hfill \\ \hfill {\varphi }_{WW}& ={\varphi }_P-2\frac{\partial \varphi }{\partial x}|{}_{{}_{\begin{array}{c}\hfill x=x_P,\hfill \\ \hfill t=t+\Delta t\hfill \end{array}}}h+2\frac{\partial [}{\partial 2}]\varphi x|{}_{\begin{array}{c}\hfill x=x_P,\hfill \\ \hfill t=t+\Delta t\hfill \end{array}}{h}^2-\frac{4}{3}\frac{{\partial }^3\varphi }{\partial x^3}|{}_{\begin{array}{c}\hfill x=x_P,\hfill \\ \hfill t=t+\Delta t\hfill \end{array}}{h}^3+\cdots \hfill \\ \hfill {\varphi }_{WWW}& ={\varphi }_P-3\frac{\partial \varphi }{\partial x}|{}_{{}_{\begin{array}{c}\hfill x=x_P,\hfill \\ \hfill t=t+\Delta t\hfill \end{array}}}h+\frac{9}{2}\frac{\partial [}{\partial 2}]\varphi x|{}_{\begin{array}{c}\hfill x=x_P,\hfill \\ \hfill t=t+\Delta t\hfill \end{array}}{h}^2-\frac{9}{2}\frac{{\partial }^3\varphi }{\partial x^3}|{}_{\begin{array}{c}\hfill x=x_P,\hfill \\ \hfill t=t+\Delta t\hfill \end{array}}{h}^3+\cdots \hfill \end{array}$(3.1)

Os valores de ${\gamma }_e^1$ e ${\gamma }_e^2$ são também recalculados em função de expansões em série de Taylor em torno de $x_P$. Estas equações ficam da seguinte maneira:

$\begin{array}{cc}\hfill {\gamma }_e^1=& {\alpha }_e^1({\varphi }_P+\frac{1}{2}\frac{\partial \varphi }{\partial x}|{}_{{}_{\begin{array}{c}\hfill x=x_P,\hfill \\ \hfill t=t+\Delta t\hfill \end{array}}}h+\frac{1}{8}\frac{\partial [}{\partial 2}]\varphi x|{}_{\begin{array}{c}\hfill x=x_P,\hfill \\ \hfill t=t+\Delta t\hfill \end{array}}{h}^2+\cdots )\hfill \\ \hfill & {\beta }_e^1(\frac{\partial \varphi }{\partial x}|{}_{{}_{\begin{array}{c}\hfill x=x_P,\hfill \\ \hfill t=t+\Delta t\hfill \end{array}}}+\frac{1}{2}\frac{\partial [}{\partial 2}]\varphi x|{}_{\begin{array}{c}\hfill x=x_P,\hfill \\ \hfill t=t+\Delta t\hfill \end{array}}h+\frac{1}{8}\frac{{\partial }^3\varphi }{\partial x^3}|{}_{\begin{array}{c}\hfill x=x_P,\hfill \\ \hfill t=t+\Delta t\hfill \end{array}}{h}^2+\cdots )\hfill \end{array}$(3.2)

e

$\begin{array}{cc}\hfill {\gamma }_e^2=& {\alpha }_e^1({\varphi }_P+\frac{1}{2}\frac{\partial \varphi }{\partial x}|{}_{{}_{\begin{array}{c}\hfill x=x_P,\hfill \\ \hfill t=t+\Delta t\hfill \end{array}}}h+\frac{1}{8}\frac{\partial [}{\partial 2}]\varphi x|{}_{\begin{array}{c}\hfill x=x_P,\hfill \\ \hfill t=t+\Delta t\hfill \end{array}}{h}^2+\cdots )\hfill \\ \hfill & {\beta }_e^2(\frac{\partial \varphi }{\partial x}|{}_{{}_{\begin{array}{c}\hfill x=x_P,\hfill \\ \hfill t=t+\Delta t\hfill \end{array}}}+\frac{1}{2}\frac{\partial [}{\partial 2}]\varphi x|{}_{\begin{array}{c}\hfill x=x_P,\hfill \\ \hfill t=t+\Delta t\hfill \end{array}}h+\frac{1}{8}\frac{{\partial }^3\varphi }{\partial x^3}|{}_{\begin{array}{c}\hfill x=x_P,\hfill \\ \hfill t=t+\Delta t\hfill \end{array}}{h}^2+\cdots )\hfill \end{array}$(3.3)

Substituindo as Eqs. (3.1) a (2.47) na Eq. (2.47) e realizando $h\to 0$ chega-se a equação original do problema, Eq. (2.1).

Procedimento semelhante foi feito para cada um dos outros três volumes do contorno e também para os volumes internos. Para todos estes casos, conclui-se que a formulação é consistente.

É possível obter outras aproximações para o primeiro volume como, por exemplo, uma que não inclua o termo ${\varphi }_{EEE}$ e que utilize as duas condições de contorno para determinar os fluxos em $x_e$ e em $x_w$. Neste caso, a consistência da discretização dependeria da ordem de aproximação empregada nas expansões em Séries de Taylor. Testes, que não apresentaremos aqui, demonstraram isto.

Desta forma, pode-se afirmar que dependendo da metologia empregada para modelar os fluxos nas faces do volume finito pode-se obter aproximações inconsistentes. Por esta razão, esta análise é obrigatória para se obter soluções numéricas acuradas.

4 ANÁLISE DE ESTABILIDADE

Uma vez verificada a consistência, a próxima etapa consiste em realizar a análise de estabilidade do comportamento do esquema numérico. Nos diversos trabalhos que estes autores tiveram acesso até o momento e que empregam a Eq. (1.2) foram empregadas formulações totalmente implícitas ou semi-implícitas a modelagem da variação temporal de $\varphi$. No entanto, em nenhum desses trabalhos pode-se encontrar uma análise da estabilidade dos esquemas numéricos apresentados, garantindo que tais formulações são incondicionalmente estáveis, como acontece na discretização da Eq. (1.1).

A estabilidade é uma condição da solução numérica que estabelece que, para valores finitos de $\Delta t$ e $h$, todo o erro definido como a diferença entre a solução numérica e a solução exata do esquema numérico, erro de arredondamento, deve permanecer limitado quando o processo de iteração avança^[2].

Definindo-se o erro de arredondamento como:

$$ϵ=\Phi -\varphi$$ (4.1)

em que $\Phi$ é o valor calculado em computador com precisão infinita e $\varphi$ é o valor em um computador real com precisão finita e substituindo $\varphi$ na equação discretizada para os volumes internos, chega-se a:

$$A_P{ϵ}_P^{\theta }=A_W{ϵ}_W^{\theta }+A_E{ϵ}_E^{\theta }+A_{WW}{ϵ}_{WW}^{\theta }+A_{EE}{ϵ}_{EE}^{\theta }+S_P$$(4.2)

com ${ϵ}^{\theta }=\theta ϵ+(1-\theta ){ϵ}^0$, e $0\le \theta \le 1$, e os valores de $A_P$, $S_P$ dependem do valor de $\theta$.

Neste trabalho o método empregado para a análise de estabilidade foi o método de von Neumann^[3]. Inicialmente, o erro $ϵ$ é expresso em cada um dos volumes finitos utilizando-se as seguintes equações:

$$\begin{array}{cccc}\hfill & {ϵ}_P=e^{\alpha (t+\Delta t)}{e}^{iK{x}_P};\hfill & \hfill & {ϵ}_P^0=e^{\alpha t}{e}^{iK{x}_P}\hfill \\ \hfill & {ϵ}_W=e^{\alpha (t+\Delta t)}{e}^{iK(x_P-h)};\hfill & \hfill & {ϵ}_W^0=e^{\alpha t}{e}^{iK(x_P-h)}\hfill \\ \hfill & {ϵ}_E=e^{\alpha (t+\Delta t)}{e}^{iK(x_P+h)};\hfill & \hfill & {ϵ}_E^0=e^{\alpha t}{e}^{iK(x_P+h)}\hfill \\ \hfill & {ϵ}_{WW}=e^{\alpha (t+\Delta t)}{e}^{iK(x_P-2h)};\hfill & \hfill & {ϵ}_{WW}^0=e^{\alpha t}{e}^{iK(x_P-2h)}\hfill \\ \hfill & {ϵ}_{EE}=e^{\alpha (t+\Delta t)}{e}^{iK(x_P+2h)};\hfill & \hfill & {ϵ}_{EE}^0=e^{\alpha t}{e}^{iK(x_P+2h)}\hfill \end{array}$$ (4.3)

A seguir são apresentados os resultados da análise de estabilidade para as três aproximações mais populares empregadas na modelagem da variação de $\varphi$ com o tempo.

4.1 Formulação Totalmente Implícita

Para esta formulação tem-se que $\theta =1$. A equação discretizada Eq. (2.8) é reescrita como:

$$A_P{ϵ}_P=A_W{ϵ}_W+A_E{ϵ}_E+A_{WW}{ϵ}_{WW}+A_{EE}{ϵ}_{EE}+S_P$$(4.4)

em que os coeficientes $A_P$ e $S_P$ são definidos pela Eq. (2.9) e os demais coeficientes pela Eq. (2.10).

Substituindo os valores de $ϵ$ adequados da Eq. (4.2) na equação anterior, e fazendo as simplificações usuais^[3], e, por fim, definido-se o fator de amplitude como:

$$G=\frac{{ϵ}_P}{{ϵ}_P^0}$$(4.5)

chega-se a seguinte expressão:

$G=\frac{3}{3+14\frac{{\lambda }_2\mathrm{\Delta }t}{h^2}+48\frac{{\lambda }_4\mathrm{\Delta }t}{h^4}}$ (4.6)

Considerando que a condição de estabilidade é $|G|\le 1$ tem-se que tal condição sempre será satisfeita, visto que ${\lambda }_2\ge 0$ e ${\lambda }_4\ge 0$.

Com o resultado acima pode-se confirmar que a formulação totalmente implícita é incondicionalmente estável para a discretização da Eq. (1.2).

4.2 Formulação Explícita

Para esta formulação tem-se que $\theta =0$. A equação discretizada pode ser escrita como:

$$A_P{ϵ}_P=A_W{ϵ}_W^0+A_E{ϵ}_E^0+A_{WW}{ϵ}_{WW}^0+A_{EE}{ϵ}_{EE}^0+S_P$$(4.7)

em que os coeficientes $A_P$ e $S_P$ são definidos pela Eq. (1.2) e os demais coeficientes pela Eq. (2.10).

Novamente, substituindo os valores de $ϵ$ adequados da Eq. (4.2) na equação anterior e realizando-se uma série de simplificações usuais desta metodologia chega-se a:

$$G=1-\frac{7}{3}\frac{{\lambda }_2\Delta t}{h^2}-8\frac{{\lambda }_4\Delta t}{h^4}$$ (4.8)

Para satisfazer a condição de estabilidade para este caso tem-se que:

$$\Delta t\le \frac{3h^4}{7{\lambda }_2{h}^2+24{\lambda }_4}$$(4.9)

ou seja, esta formulação é condicionalmente estável.

Este resultado é ligeiramente diferente do que é habitualmente encontrado na literatura para a equação da difusão clássica, isto é,

$$\Delta t\le \frac{h^2}{2{\lambda }_2}$$ (4.10)

e, como se pode observar, fazendo-se ${\lambda }_4=0$ na Eq. (4.9) não se obtém a Eq. (4.10). Isto se dá porque a Eq. (4.10) foi definida empregando-se uma aproximação de $2^{\underset{\_}{a}}$ ordem para o termo $\frac{\partial *}{\partial [}2]\varphi x$ e a utilizada neste trabalho foi uma aproximação de $4^{\underset{\_}{a}}$ ordem, como explicitado na Eq. (2.6).

A Figura1 apresenta o comportamento de $\Delta t$ calculado segundo a Eq. (4.9) em função de $h$ para três valores da relação ${\lambda }_2/{\lambda }_4$ distintos. Os valores obtidos para o caso em que ${\lambda }_2/{\lambda }_4={10}^{-4}$ diferem minimamente dos obtidos para ${\lambda }_2/{\lambda }_4=1$ e, por esta razão, suas curvas estão sobrepostas.

Figure 1
$\Delta t$ para se obter uma solução estável em função de $h$ para três valores de ${\lambda }_2/{\lambda }_4$.

Nestes três casos, como seria de se esperar, quanto menor for o valor de $h$, menor é o $\Delta t$ necessário para se obter uma solução estável.

Observando na Figura1 conclui-se que para $h=0.01$ tem-se $\Delta t<{10}^{-10}$. Valores de $\Delta t$ desta magnitude certamente inviabilizam o uso desta formulação, pois seria necessário um tempo computacional exorbitantemente alto.

Por fim, tem-se uma discreta diferença para o caso em que ${\lambda }_2/{\lambda }_4={10}^4$, isto é, para problemas que se aproximam da equação de difusão clássica. Esta diferença aparece nos casos em que $h$ é grande. Mas, mesmo assim, a conclusão que se pode obter da Figura1 é que valores de $\Delta t$ que permite uma solução estável são tão pequenos que soluções totalmente implícitas, que não possuem limites dos valores de $\Delta t$, levarão bem menos tempo computacional para resolver o mesmo problema.

4.3 Formulação de Crank–Nicolson

No texto não foram apresentados os coeficientes para quando $\theta =1/2$, isto é, para a formulação de Crank–Nicolson^[1]. Contudo, por ser relativamente comum nos livros de solução numérica de equações diferenciais, também será apresentada a parte final da demonstração de estabilidade para essa formulação.

Seguindo o mesmo procedimento apresentado anteriormente chega-se a expressão do fator de amplitude para esta formulação em especial:

$G=\frac{7\frac{{\lambda }_2\mathrm{\Delta }t}{h^2}+24\frac{{\lambda }_4\mathrm{\Delta }t}{h^4}-3}{7\frac{{\lambda }_2\mathrm{\Delta }t}{h^2}+24\frac{{\lambda }_4\mathrm{\Delta }t}{h^4}+3}$ (4.11)

Para ser estável tem-se que $|G|\le 1$ e, como pode ser facilmente determinado, tem-se que esta condição é sempre satisfeita, independentemente do valor de $\Delta t$, logo a formulação de Crank-Nicolson é, também, incondicionalmente estável.

5 COMENTÁRIOS E CONCLUSÕES

Neste trabalho apresentou-se uma formulação para Volumes Finitos da equação de difusão de fluxo bimodal proposta por Bevilacqua e colaboradores^[3],^[4].

A EDP que governa tal fenômeno, Eq. (1.2), possui um termo com uma diferencial de 4ª ordem. Diferenciais desta ordem não são usualmente encontradas em equações de Fenômenos de Transporte. Como consequência deste fato, não há muito na literatura sobre como proceder para discretizar tais diferenciais no âmbito do Método de Volumes ou Diferenças Finitas, principalmente como proceder com os volumes que terão seus fluxos calculados em função das condições de contorno.

Diferentemente de outros trabalhos congêneres, em que se apresenta somente a equação discretizada para os volumes internos, neste trabalho pode se encontrar a dedução de equações para todos os volumes do domínio.

Outro aspecto importante do trabalho está relacionado com a análise de consistência da equação discretizada. Ao demonstrar que a formulação apresentada é consistente, se está se assegurando que a solução numérica gerada por esta formulação possuirá uma acurácia que não seria possível obter caso a equação discretizada não fosse consistente.

Por fim, foi apresentada uma análise de estabilidade utilizando a metodologia de von Neumann. Entendem estes autores que este é o resultado mais relevante deste trabalho. As condições de estabilidade para a equação de difusão clássica são amplamente conhecidas, contudo, para a equação de fluxo bimodal, Eq. (1.2), esta é a primeira vez que é apresentada, até onde estes autores puderam determinar.

AGRADECIMENTOS

Os autores agradecem as sugestões apresentadas pelos revisores deste texto e pelo apoio financeiro fornecido da FAPERJ, Fundação Carlos Chagas Filho de Amparo à Pesquisa do Estado do Rio de Janeiro, do CNPq, Conselho Nacional de Desenvolvimento Científico e Tecnológico, e da CAPES, Fundação Coordenação de Aperfeiçoamento de Pessoal de Nível Superior.

Consistency analysis and stability conditions of the bi-flux diffusion equation are discussed in this paper. Several details are presented on how the differential equation has been discretized based on the Finite Volume Method.

annomalous diffusion, bi-flux diffusion, Finite Volume method.

REFERENCES

• ^[1]
  S. Balay, J. Brown, K. Buschelman, W.D. Gropp, D. Kaushik, M.G. Knepley, L.C. McInnes, B.F. Smith & H. Zhang. PETSc web page (2011).
• ^[2]
  S. Balay, W.D. Gropp, L.C. McInnes & B.F. Smith. Efficient management of parallelism in object oriented numerical software libraries. In E. Arge, A.M. Bruaset & H.P. Langtangen (editors), “Modern software tools in scientific computing”. Birkhäuser Press (1997), pp. 163–202.
• ^[3]
  L. Bevilacqua, A.C.N.R. Gale˜ao & F.P. Costa. A new analytical formulation of retention effects on particle diffusion processes. Anais da Academia Brasileira de Ciências, 83 (2011), 1443–1464. URL http://www.scielo.br/scielo.php?script=sci_arttext&pid=S0001-37652011000400031&nrm=iso
  » http://www.scielo.br/scielo.php?script=sci_arttext&pid=S0001-37652011000400031&nrm=iso
• ^[4]
  L. Bevilacqua, A.C.N.R. Gale˜ao & F.P. Costa. On the significance of higher order differential terms in diffusion processes. Journal of the Brazilian Society of Mechanical Sciences and Engineering, 33(2) (2011), 166–175. doi:10.1590/S1678-58782011000200007. URL http://www.scielo.br/scielo.php?script=sci_arttext&pid=S1678-58782011000200007&lng=en&nrm=iso&tlng=en.
• ^[5]
  J. Crank & P. Nicolson. A practical method for numerical evaluation of solutions of partial differential equations of the heat-conduction type. Advances in Computational Mathematics, 6(1) (1996), 207–226. doi:10.1007/BF02127704. URL http://link.springer.com/10.1007/BF02127704.
• ^[6]
  JV. Ganti, M.M. Meerschaert, E. Foufoula-Georgiou, E. Viparelli & G. Parker. Normal and anomalous diffusion of gravel tracer particles in rivers. Journal of Geophysical Research: Earth Surface, 115(F2) (2010). doi:10.1029/2008JF001222. URL http://doi.wiley.com/10.1029/2008JF001222.
• ^[7]
  G.H. Golub & C.F. Van Loan. “Matrix computations”. Johns Hopkins University Press, Baltimore, 4th ed. ed. (2013), 784 pp.
• ^[8]
  C. Hirsch. Consistency, Stability and Error Analysis of Numerical Schemes. In “Numerical Computation of Internal and External Flows”. Elsevier (2007), pp. 283–335. doi:10.1016/B978-075066594-0/ 50049-7. URL http://linkinghub.elsevier.com/retrieve/pii/B9780750665940500497.
• ^[9]
  M. Jiang. “The fourth order diffusion model for a bi-flux mass transfer”. Ph.D. thesis, Universidade Federal do Rio de Janeiro, (2017). URL www.coc.ufrj.br/pt/documents2/doutorado/2017/2922-jiang-m-td-17/file.
• ^[10]
  J. Klafter & I.M. Sokolov. Anomalous diffusion spreads its wings. Physics World, 18(8) (2005), 29–32. doi:10.1088/2058-7058/18/8/33. URL http://stacks.iop.org/2058-7058/18/i=8/a=33?key=crossref.49d9670701f966532fc6c0805ce1fde2.
• ^[11]
  J.Y. Lin. The non-Fourier effect on the fin performance under periodic thermal conditions. Applied Mathematical Modelling, 22(8) (1998), 629–640. doi:10.1016/S0307-904X(98)10061-6. URL https://linkinghub.elsevier.com/retrieve/pii/S0307904X98100616.
• ^[12]
  E.H. Macdonald. Sedimentation and detrital gold. In “Handbook of Gold Exploration and Evaluation”. Elsevier (2007), pp. 195–266. doi:10.1533/9781845692544.195. URL https://linkinghub.elsevier.com/retrieve/pii/B9781845691752500040.
• ^[13]
  C.R. Maliska. “Transferência de calor e mecânica dos fluidos computacional”. LTC Editora (1994), 472 pp.
• ^[14]
  G.M. Marinho. “Análise numérica do problema de difus˜ao anômala unidimensional”. Ph.D. thesis, Universidade do Estado do Rio de Janeiro (2014).
• ^[15]
  P.F. Nealey, R.E. Cohen & A.S. Argon. Limited-supply non-Fickian diffusion in glassy polymers. Polymer, 36(19) (1995), 3687–3695. doi:10.1016/0032-3861(95)93771-D. URL http://www.sciencedirect.com/science/article/pii/003238619593771D.
• ^[16]
  S.V. Patankar. “Numerical heat transfer”. CRC Press, Inc., New York (1980).
• ^[17]
  Y. Saad & M. Schultz. GMRES: A generalized minimal residual algorithm for solving nonsymmetric linear systems. SIAM Journal on Scientific and Statistical Computing, 7(3) (1986), 856–869. doi:10.1137/0907058.
• ^[18]
  M. Santos, F.P. Costa, A.C.N.R. Gale˜ao & L. Bevilacqua. O Método de Crank-Nicolson aplicado ao modelo de difusão do conhecimento: Uma simulação para o processo de transmissão do conhecimento. International Journal of Knowledge Engineering and Management, 4(9) (2015), 129–146. URL http://incubadora.periodicos.ufsc.br/index.php/IJKEM/article/viewFile/3310/4164.
• ^[19]
  J.B. Scarborough. “Numerical mathematical analysis”. The Johns Hopkins Press, Baltimore (1966).
• ^[20]
  L.G. Silva. “Problemas inversos em processos difusivos com retenção”. Ph.D. thesis, Universidade do Estado do Rio de Janeiro (2013).
• ^[21]
  J.G. Simas. “Modelagem computacional do problema de difusão com retenção”. Ph.D. thesis, Laboratório Nacional de Computação Científica (2012).
• ^[22]
  J.F. Vasconcellos, D.C. Knupp & G.M. Marinho. Uma comparação entre o Método de Volumes Finitos e a Técnica de Transformada Integral Generalizada par a solução de uma equação de difusão bidimensional. Revista Mundi Engenharia, Tecnologia e Gest˜ao, 4(3) (2019). doi:10.21575/25254782rmetg2019vol4n3850. URL http://periodicos.ifpr.edu.br/index.php?journal=MundiETG&page=article&op=view&path[]=850.
• ^[23]
  J.F. Vasconcellos, G.M. Marinho & J.H. Zani. Análise numérica da equação da difusão anômala com fluxo bimodal. Revista Internacional de Métodos Numéricos para Cálculo y Diseño en Ingeniería, 33(3-4) (2017), 242–249. doi:10.1016/j.rimni.2016.05.001. URL https://www.scipedia.com/public/Vasconcellos_et_al_2016a.
• ^[24]
  L. Vlahos, H. Isliker, Y. Kominis & K. Hizanidis. Normal and Anomalous Diffusion: A Tutorial. (2008). URL http://arxiv.org/abs/0805.0419.
• ^[25]
  P.A. Wankhade, B. Kundu & R. Das. Establishment of non-Fourier heat conduction model for an accurate transient thermal response in wet fins. International Journal of Heat and Mass Transfer, 126 (2018), 911–923. doi:10.1016/j.ijheatmasstransfer.2018.05.094. URL https://linkinghub.elsevier.com/retrieve/pii/S0017931018315977.
• ^[26]
  B. Xu, G. Gong, Y. Fan, B. Wu & J.H. Gao. Directional sensitivity of anomalous diffusion in human brain assessed by tensorial fractional motion model. Magnetic Resonance Imaging, 42 (2017), 74–81. doi:10.1016/j.mri.2017.05.006. URL https://linkinghub.elsevier.com/retrieve/pii/S0730725X17301030.
• ^[27]
  O.C. Zienkiewicz. Numerical computation of internal and external flow vol. 1. Fundamentals of numerical discretisation, C. Hirsch, Wiley, Chichester, 1988, ISBN 0471 917621. International Journal for Numerical Methods in Engineering, 28(10) (1989), 2465–2465. doi:10.1002/nme.1620281016. URL http://doi.wiley.com/10.1002/nme.1620281016.

Datas de Publicação

• Publicação nesta coleção
  20 Dez 2019
• Data do Fascículo
  Sep-Dec 2019

Histórico

• Recebido
  06 Abr 2019
• Aceito
  01 Jul 2019