Numerical analysis of an anomalous diffusion with a bimodal flux distribution

==Resumo==

Um modelo numérico da equação de difusão unidimensional com fluxo bimodal em regime permanente e malhas uniformes foi desenvolvido utilizando‐se o método de volumes finitos. Esta equação possui, além do fluxo de difusão, fluxo de Fick, um outro fluxo denominado fluxo secundário ou fluxo de retenção. No modelo desenvolvido, analisaram‐se as intensidades normalizadas, parceladas e distintas para ambos os fluxos. Neste trabalho, também foram estudados os comportamentos dos erros numéricos, relativos aos diferentes parâmetros, e.g., amplitude e frequência da solução como também para condições de contorno distintas. Os resultados numéricos mostram o bom desempenho do modelo proposto.

==Abstract==

A numerical study of unidimensional and steady state bimodal flux diffusion equation has been developed using the Finite Volume Method. In addition, with the Ficks flux, this equation has another flux named retention or secondary flux. The influence of each flux has been analysed in this work. The influence of several different parameters, such as boundary conditions or amplitude and frequency of the solution, in the numerical errors has been studied. Numerical results show the performance of the proposed model.

==Palavras‐chave==

Difusão anômala;Método de volumes finitos;Análise numérica;Difusão de quarta ordem

==Keywords==

Bimodal diffusion;Finite volume method;Numerical analysis;Fourth order diffusion

==1. Introdução==

O modelo universalmente adotado para descrever a difusão de um soluto em um meio contínuo baseia‐se na lei de Fick. Contudo, há na literatura relatos que descrevem processos difusivos que não seguem esta lei [[#bib0080|[1]]]; [[#bib0085|[2]]]; [[#bib0090|[3]]]; [[#bib0095|[4]]] ;  [[#bib0100|[5]]]. Bevilacqua et al. [[#bib0105|[6]]]; [[#bib0110|[7]]] ;  [[#bib0115|[8]]] apresentaram uma formulação analítica para descrever o processo difusivo em que ocorre uma forma de retenção do fluxo de Fick em um meio reativo sob condições termodinâmicas estáveis. Baseado nestas e outras considerações, estes autores propuseram a seguinte equação para descrever tal fenômeno para uma situação unidimensional:

<span id='eq0005'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\frac{\partial \phi }{\partial t}=\left(1-\beta \right)k_{22}\frac{{\partial }^2\phi }{\partial x^2}-</math><math>\beta \left(1-\beta \right)k_{44}\frac{{\partial }^4\phi }{\partial x^4}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (1)
|}

O termo de quarta ordem da Eq. [[#eq0005|(1)]] foi denominado originalmente de termo de difusão de retenção [[#bib0105|[6]]] ;  [[#bib0110|[7]]]. Em trabalhos posteriores, houve uma mudança na nomenclatura para termo de difusão secundária [[#bib0115|[8]]]. Nesta nova nomenclatura a difusão de Fick passou a ser denominada de difusão primária. Adotaremos a segunda opção para diferenciar os fluxos envolvidos no processo difusivo bimodal em análise neste trabalho. Desta forma, na Eq. [[#eq0005|(1)]], <math display="inline">\left(1-\beta \right)k_{22}</math> representa o coeficiente de difusão de Fick, e que denominaremos de difusão primária, e <math display="inline">\beta \left(1-\beta \right)k_{44}</math> o coeficiente de difusão de retenção, que denominaremos de difusão secundária. Tem‐se, ainda, que 0 ≤ ''β'' < 1 e neste trabalho os valores de ''k''<sub>22</sub> e ''k''<sub>44</sub> são constantes. Nas situações em que ''β'' = 0 ou ''k''<sub>44</sub> = 0, a Eq. [[#eq0005|(1)]] se reduzirá a forma clássica da lei de Fick, i. e., uma equação para um fluxo de difusão unimodal.

Os trabalhos de Bevilacqua et al. [[#bib0105|[6]]]; [[#bib0110|[7]]] ;  [[#bib0115|[8]]] abordaram quase que exclusivamente a apresentação da teoria que levou ao desenvolvimento da Eq. [[#eq0005|(1)]], enquanto a apresentação de uma metodologia para resolver esta equação foi tratada de forma secundária. Neste trabalho, trataremos exclusivamente da utilização do método de volumes finitos (MVF) para resolver numericamente o sistema composto pela Eq. [[#eq0005|(1)]] e suas respectivas condições de contorno. Neste ponto, para melhor compreensão das hipóteses que levaram à elaboração da Eq. [[#eq0005|(1)]], recomenda‐se a leitura dos artigos supracitados.

O grande desafio na modelagem da Eq. [[#eq0005|(1)]] é descrever de maneira adequada o termo envolvendo a derivada de quarta ordem. Isto porque enquanto o termo transiente e aquele que contém a derivada de segunda ordem são encontrados frequentemente na literatura, derivadas de quarta ordem são bem menos encontradas na modelagem de problemas que envolvem fenômenos de transporte.

O objetivo principal deste trabalho é o de apresentar a discretização da Eq. [[#eq0005|(1)]] através do MVF, de forma que tanto o termo com a derivada de segunda ordem quanto o de quarta ordem estejam corretamente modelados, independente dos valores de ''k''<sub>22</sub> ou ''k''<sub>44</sub>. Portanto, desconsiderou‐se o termo transiente da Eq. [[#eq0005|(1)]], com o propósito de destacar os problemas que possam vir a surgir nas soluções numéricas relativas aos sistemas de equações discretizadas de balanço espacial. Os problemas associados à discretização temporal serão objeto de um futuro trabalho. Desta forma, a equação objeto deste trabalho restringe‐se a Eq. [[#eq0010|(2)]]:

<span id='eq0010'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>{\Gamma }_{22}\frac{d{\phi }^2}{dx^2}-{\Gamma }_{44}\frac{d{\phi }^4}{dx^4}+</math><math>B\left(x\right)=0</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (2)
|}

Discretizando‐se a Eq. [[#eq0010|(2)]], c.f., utilizando técnicas do MVF [[#bib0120|[9]]], obtém‐se:

<span id='eq0015'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>{\Gamma }_{22}\left({\frac{d\phi }{dx}}_{x_i+\frac{1}{2}}-\right. </math><math>\left. {\frac{d\phi }{dx}}_{x_i-\frac{1}{2}}\right)-</math><math>{\Gamma }_{44}\left({\frac{d^3\phi }{dx^3}}_{x_i+\frac{1}{2}}-\right. </math><math>\left. {\frac{d^3\phi }{dx^3}}_{x_i-\frac{1}{2}}\right)+</math><math>B\left(x_i\right)h=0</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (3)
|}

Na Eq. [[#eq0015|(3)]], ''x''<sub>''i''</sub> representa a coordenada do centro do i‐ésimo volume de tamanho ''h'', detalhado na [[#fig0005|figura 1]].

<span id='fig0005'></span>

{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; max-width: 100%;" 
|-
|


[[Image:draft_Content_428379833-1-s2.0-S0213131516300232-gr1.jpg|center|301px|Esquema dos volumes internos.]]


|-
| <span style="text-align: center; font-size: 75%;">

Figura 1.

Esquema dos volumes internos.

</span>
|}

A seguir, são apresentadas e discutidas formas diferentes para avaliar os fluxos primário e secundário. Inicialmente, trata‐se da modelagem dos fluxos que serão utilizados na definição das equações do MVF, isto para os volumes internos. Em seguida, para avaliar os fluxos dos volumes adjacentes, trata‐se às fronteiras do contorno. Após detalhar a equação discretizada, apresentam‐se os resultados obtidos com esta metodologia.

==2. Modelo simplificado==

Podem‐se escrever os valores das variáveis nos centros dos volumes em função de expansões em série de Taylor em torno da interface. Desta forma, para a interface esquerda do volume genérico «i» tem‐se que:

<span id='eq0020'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>{\phi }_{i-1}={\phi }_{x_i-\frac{1}{2}}-\frac{h}{2}{\frac{d\phi }{dx}}_{x_i-\frac{1}{2}}+</math><math>\frac{h^2}{8}{\frac{d^2\phi }{dx^2}}_{x_i-\frac{1}{2}}-</math><math>\frac{h^3}{48}{\frac{d^2\phi }{dx^2}}_{x_i-\frac{1}{2}}+</math><math>\frac{h^3}{384}{\frac{d^3\phi }{dx^3}}_{x_i-\frac{1}{2}}-</math><math>\cdots </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (4)
|}

<span id='eq0025'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>{\phi }_i={\phi }_{x_i-\frac{1}{2}}+\frac{h}{2}{\frac{d\phi }{dx}}_{x_i-\frac{1}{2}}+</math><math>\frac{h^2}{8}{\frac{d^2\phi }{dx^2}}_{x_i-\frac{1}{2}}+</math><math>\frac{h^3}{48}{\frac{d^2\phi }{dx^2}}_{x_i-\frac{1}{2}}+</math><math>\frac{h^3}{384}{\frac{d^3\phi }{dx^3}}_{x_i-\frac{1}{2}}+</math><math>\cdots </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (5)
|}

Os índices presentes nas Eqs. [[#eq0020|(4)]] a [[#eq0035|(7)]] estão representados na [[#fig0005|figura 1]]. Resolvendo‐se o sistema linear de equações formado pelas Eqs. [[#eq0020|(4)]] e [[#eq0025|(5)]], obtém‐se uma aproximação para o fluxo primário na interface a esquerda do volume «i»:

<span id='eq0030'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>{\frac{d\phi }{dx}}_{x_i-\frac{1}{2}}=\frac{{\phi }_i-{\phi }_{i-1}}{h}-</math><math>\frac{h^2}{24}{\frac{d^3\phi }{dx^3}}_{x_i-\frac{1}{2}}+</math><math>\cdots </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (6)
|}

Utilizando‐se o mesmo procedimento empregado na obtenção da Eq. [[#eq0030|(6)]], pode‐se obter o fluxo primário para a interface direita do volume «i»:

<span id='eq0035'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>{\frac{d\phi }{dx}}_{x_i+\frac{1}{2}}=\frac{{\phi }_{i+1}-{\phi }_i}{h}-</math><math>\frac{h^2}{24}{\frac{d^3\phi }{dx^3}}_{x_i+\frac{1}{2}}+</math><math>\cdots </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (7)
|}

Para a avaliação do fluxo secundário foram realizadas expansões em série de Taylor em torno das faces dos volumes, da mesma maneira utilizada na obtenção das Eqs. [[#eq0030|(6)]] e [[#eq0035|(7)]]. Obviamente que para se obter uma aproximação de uma derivada de terceira ordem serão expansões em 4 pontos. Nesse caso, além dos pontos definidos pelas Eqs. [[#eq0020|(4)]] e [[#eq0025|(5)]], utilizou‐se as seguintes expansões na definição do fluxo secundário na face esquerda do volume:

<span id='eq0040'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>{\phi }_{i-2}={\phi }_{x_i-\frac{1}{2}}-\frac{3h}{2}{\frac{d\phi }{dx}}_{x_i-\frac{1}{2}}+</math><math>\frac{9h^2}{8}{\frac{d^2\phi }{dx^2}}_{x_i-\frac{1}{2}}-</math><math>\frac{9h^3}{16}{\frac{d^3\phi }{dx^3}}_{x_i-\frac{1}{2}}+</math><math>\frac{27h^4}{128}{\frac{d^4\phi }{dx^4}}_{x_i-\frac{1}{2}}-</math><math>\frac{81h^5}{1280}{\frac{d^5\phi }{dx^5}}_{x_i-\frac{1}{2}}+</math><math>\cdots </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (8)
|}

<span id='eq0045'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>{\phi }_{i+1}={\phi }_{x_i-\frac{1}{2}}+\frac{3h}{2}{\frac{d\phi }{dx}}_{x_i-\frac{1}{2}}+</math><math>\frac{9h^2}{8}{\frac{d^2\phi }{dx^2}}_{x_i-\frac{1}{2}}+</math><math>\frac{9h^3}{16}{\frac{d^3\phi }{dx^3}}_{x_i-\frac{1}{2}}+</math><math>\frac{27h^4}{128}{\frac{d^4\phi }{dx^4}}_{x_i-\frac{1}{2}}+</math><math>\frac{81h^5}{1280}{\frac{d^5\phi }{dx^5}}_{x_i-\frac{1}{2}}+</math><math>\cdots </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (9)
|}

A solução do sistema linear de equações envolvendo as Eqs. [[#eq0020|(4)]], [[#eq0025|(5)]], [[#eq0040|(8)]] e [[#eq0045|(9)]] leva a uma aproximação do fluxo secundário na face esquerda do volume:

<span id='eq0050'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>{\frac{d^3\phi }{dx^3}}_{x_i-\frac{1}{2}}=\frac{{\phi }_{i+1}-{\phi }_{i-2}}{h^3}+</math><math>3\frac{{\phi }_{i-1}-{\phi }_i}{h^3}-{\frac{1}{8}\frac{d^5\phi }{dx^5}}_{x_i-\frac{1}{2}}h^2+</math><math>\cdots </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (10)
|}

Utilizando‐se o mesmo procedimento empregado para a obtenção da Eq. [[#eq0050|(10)]] para se obter o fluxo secundário na face direita do volume, chega‐se a seguinte equação:

<span id='eq0055'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>{\frac{d^3\phi }{dx^3}}_{x_i+\frac{1}{2}}=\frac{{\phi }_{i+2}-{\phi }_{i-1}}{h^3}+</math><math>3\frac{{\phi }_i-{\phi }_{i+1}}{h^3}-{\frac{1}{8}\frac{d^5\phi }{dx^5}}_{x_i+\frac{1}{2}}h^2+</math><math>\cdots </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (11)
|}

A Eq. [[#eq0050|(10)]] não pode ser utilizada para calcular os fluxos secundários das faces esquerdas dos seguintes volumes: primeiro, segundo e último. Para estas faces devem‐se calcular os fluxos de outra maneira, como mostraremos mais adiante. O mesmo se pode dizer para o cálculo do fluxo secundário das faces direita do primeiro, do penúltimo e do último volume. Para estes volumes, a Eq. [[#eq0055|(11)]] não pode ser empregada.

===2.1. Equação dos volumes internos===

As Eqs. [[#eq0030|(6)]], [[#eq0035|(7)]], [[#eq0050|(10)]] e [[#eq0055|(11)]] definem os fluxos em cada uma das faces dos volumes. Substituindo‐as na Eq. [[#eq0015|(3)]], chega‐se a equação típica para os volumes internos do modelo denominado modelo simplificado. Esta equação tem a forma dada pela Eq. [[#eq0045|(9)]]:

<span id='eq0060'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\left(\frac{2{\Gamma }_{22}}{h}+\frac{6{\Gamma }_{44}}{h^3}\right){\phi }_i-</math><math>\left(\frac{{\Gamma }_{22}}{h}+\frac{4{\Gamma }_{44}}{h^3}\right){\phi }_{i-1}-</math><math>\left(\frac{{\Gamma }_{22}}{h}+\frac{4{\Gamma }_{44}}{h^3}\right){\phi }_{i+1}+</math><math>\left(\frac{{\Gamma }_{44}}{h^3}\right){\phi }_{i-2}+</math><math>\left(\frac{{\Gamma }_{44}}{h^3}\right){\phi }_{i+2}=</math><math>B_ih+\left[\frac{{\Gamma }_{22}}{24}{\frac{d^4\phi }{dx^4}}_{x_i}-\right. </math><math>\left. {\frac{{\Gamma }_{44}}{8}\frac{d^6\phi }{dx^6}}_{x_i}\right]h^3+</math><math>\cdots </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (12)
|}

e <math display="inline">B_i=B\left(x_i\right)</math>. Esta equação é válida para todos os volumes, com exceção dos 2 primeiros e dos 2 últimos volumes. Para estes volumes, em ao menos uma das faces, os fluxos secundários não poderão ser calculados através das Eqs. [[#eq0050|(10)]] ou [[#eq0055|(11)]]. Observe‐se que nas Eqs. [[#eq0030|(6)]], [[#eq0035|(7)]], [[#eq0050|(10)]] e [[#eq0055|(11)]] as derivadas foram avaliadas nas interfaces e na Eq. [[#eq0060|(12)]] as derivadas foram avaliadas no centro do volume «i». Para isto ser possível, uma expansão em série de Taylor somente das derivadas das Eqs. [[#eq0030|(6)]], [[#eq0035|(7)]], [[#eq0050|(10)]] e [[#eq0055|(11)]] foi realizada, recalculando‐as em torno do centro do volume «i».

===2.2. Condição de contorno===

Para se obter a solução numérica da Eq. [[#eq0010|(2)]] é necessário utilizar as 4 condições de contorno e desenvolver equações diferentes da Eq. [[#eq0060|(12)]], de modo que as condições de contorno estejam presentes nas equações discretizadas. No modelo desenvolvido empregaram‐se as condições de contorno do tipo: Dirichlet, Neumann ou Robin. Definiram‐se as 2 condições de contorno em ''x'' = 0, através dos valores dos coeficientes das equações apresentadas a seguir:

<span id='eq0065'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>w_{10}{\phi }_{x=0}+{w_{11}\frac{d\phi }{dx}}_{x=0}=</math><math>w_{12};{ }\mbox{e};{ }w_{20}{\phi }_{x=0}+</math><math>{w_{21}\frac{d\phi }{dx}}_{x=0}=w_2</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (13)
|}

Na Eq. [[#eq0065|(13)]], <math display="inline">w_{10}</math>, <math display="inline">w_{11}</math>, <math display="inline">w_{12}</math>, <math display="inline">w_{20}</math>, <math display="inline">w_{21}</math> e <math display="inline">w_{22}</math> são valores conhecidos ''a priori  '' e determinam a condição de contorno que será utilizada na solução da equação. Por exemplo, para empregar uma condição de contorno tipo Neumann deve‐se definir os seguintes valores: <math display="inline">w_{10}=0</math> e <math display="inline">w_{11}\not =0</math>.

A maneira selecionada para que as condições de contorno fossem incluídas no modelo se mostrou fundamental na obtenção de resultados confiáveis para a da Eq. [[#eq0010|(2)]]. A primeira forma idealizada para implementar as condições de contorno baseou‐se na utilização da Eq. [[#eq0065|(13)]], em conjunto com expansões em série de Taylor utilizando‐se os 2 primeiros volumes, Eqs. [[#eq0025|(5)]] e [[#eq0045|(9)]]. Nessa forma, os fluxos primário e secundário da face esquerda do primeiro volume, conforme mostra a [[#fig0010|figura 2]], assumem a forma:

<span id='eq0070'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>{\frac{d\phi }{dx}}_{x=0}=\frac{w_{10}w_{22}-w_{20}w_{12}}{w_{10}w_{21}-w_{20}w_{11}}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (14)
|}

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>{\frac{d^3\phi }{dx^3}}_{x=0}=-\frac{24{\phi }_0}{h^3}+</math><math>\frac{8{\phi }_1}{3h^3}+\frac{8}{h^2}\left(\frac{w_{10}w_{22}-w_{20}w_{12}}{w_{10}w_{21}-w_{20}w_{11}}\right)+</math><math>\frac{64}{3h^3}\left(\frac{w_{12}w_{21}-w_{22}w_{11}}{w_{10}w_{21}-w_{20}w_{11}}\right)+</math><math>{\frac{1}{2}\frac{d^4\phi }{dx^4}}_{x=x_0}h+\cdots </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (15)
|}

e para a face direita deste mesmo volume utilizaram‐se expansões em série de Taylor em torno desta face para os 4 primeiros volumes representados na [[#fig0010|figura 2]]. Desta forma, para esta face, os fluxos foram assim determinados:

<span id='eq0080'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>{\frac{d\phi }{dx}}_{x=h}=\frac{{\phi }_1-{\phi }_0}{h}-</math><math>\frac{1}{24}{\frac{d^3\phi }{dx^3}}_{x=h}h^2+\cdots </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (16)
|}

<span id='eq0085'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>{\frac{d^3\phi }{dx^3}}_{x=h}=\frac{{\phi }_3-{\phi }_0+3\left({\phi }_1-{\phi }_2\right)}{h^3}-</math><math>{\frac{d^4\phi }{dx^4}}_{x=x_0}h+\cdots </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (17)
|}

<span id='fig0010'></span>

{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; max-width: 100%;" 
|-
|


[[Image:draft_Content_428379833-1-s2.0-S0213131516300232-gr2.jpg|center|224px|Malha na região do contorno em x=0.]]


|-
| <span style="text-align: center; font-size: 75%;">

Figura 2.

Malha na região do contorno em ''x'' = 0.

</span>
|}

Substituindo as Eqs. [[#eq0070|(14)]] a [[#eq0085|(17)]] na Eq. [[#eq0015|(3)]], chega‐se a equação discretizada para o primeiro volume:

<span id='eq0090'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\left(\frac{{\Gamma }_{22}}{h}+23\frac{{\Gamma }_{44}}{h^3}\right){\phi }_0+</math><math>\left(-\frac{{\Gamma }_{22}}{h}+\frac{{\Gamma }_{44}}{3h^3}\right){\phi }_1-</math><math>3\frac{{\Gamma }_{44}}{h^3}{\phi }_2+\frac{{\Gamma }_{44}}{h^3}{\phi }_3=</math><math>B_ih-{\Gamma }_{22}\left(\frac{w_{10}w_{22}-w_{20}w_{12}}{w_{10}w_{21}-w_{20}w_{11}}\right)+</math><math>{\Gamma }_{44}\frac{8}{h^2}\left(\frac{w_{10}w_{22}-w_{20}w_{12}}{w_{10}w_{21}-w_{20}w_{11}}\right)+</math><math>\frac{64}{3h^3}\left(\frac{w_{12}w_{21}-w_{22}w_{11}}{w_{10}w_{21}-w_{20}w_{11}}\right)+</math><math>\frac{{\Gamma }_{22}}{24}{\frac{d^3\phi }{dx^3}}_{x=x_0}h^2+</math><math>{\frac{{\Gamma }_{44}}{2}\frac{d^4\phi }{dx^4}}_{x=x_0}h+</math><math>\cdots </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (18)
|}

Este modelo não chegou a ser utilizado, pois uma análise prévia demonstrou que ao se escrever os valores de ''ϕ''<sub>0</sub>, ''ϕ''<sub>1</sub>, ''ϕ''<sub>2</sub> e ''ϕ''<sub>3</sub> em termos de série de Taylor e substituindo‐se estes termos na Eq. [[#eq0090|(18)]] para o limite em que ''h'' → 0 chega‐se a equação:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>{\Gamma }_{22}\frac{d{\phi }^2}{dx^2}-\frac{3}{2}{\Gamma }_{44}\frac{d{\phi }^4}{dx^4}+</math><math>B\left(x\right)=0</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (19)
|}

Conclui‐se desta forma que a Eq. [[#eq0090|(18)]] é inconsistente, i.e., quando ''h'' → 0, a equação discretizada não tende para a equação diferencial original. A razão desta inconsistência está associada às ordens diferentes das aproximações empregadas nas Eqs. [[#eq0080|(16)]] e [[#eq0085|(17)]].

A solução encontrada para se resolver o problema de inconsistência descrita anteriormente foi a de incluir um ponto a mais nas expansões em série de Taylor. Nesta nova abordagem, os fluxos secundários foram calculados como:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>{\frac{d^3\phi }{dx^3}}_{x=0}=-\frac{48{\phi }_0}{h^3}+</math><math>\frac{8{\phi }_1}{h^3}-\frac{24{\phi }_2}{25h^3}+\frac{72}{5h^2}\left(\frac{w_{10}w_{22}-w_{20}w_{12}}{w_{10}w_{21}-w_{20}w_{11}}\right)+</math><math>\frac{1024}{25h^3}\left(\frac{w_{12}w_{21}-w_{22}w_{11}}{w_{10}w_{21}-w_{20}w_{11}}\right)+</math><math>{\frac{23}{80}\frac{d^5\phi }{dx^5}}_{x=x_0}h^2+\cdots </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (20)
|}

e

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>{\frac{d^3\phi }{dx^3}}_{x=h}=\frac{7{\phi }_1-2{\phi }_0-9{\phi }_2+5{\phi }_3-{\phi }_4}{h^3}-</math><math>\frac{7}{8}{\frac{d^5\phi }{dx^5}}_{x=x_0}h^2+\cdots </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (21)
|}

Os fluxos primários para este modelo ficam inalterados, Eqs. [[#eq0070|(14)]] e [[#eq0080|(16)]].

A equação para o primeiro volume fica reescrita como:

<span id='eq0110'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\left(\frac{{\Gamma }_{22}}{h}\mbox{+46}\frac{{\Gamma }_{44}}{h^3}\right){\phi }_0-</math><math>\left(\frac{{\Gamma }_{22}}{h}\mbox{+}\frac{{\Gamma }_{44}}{h^3}\right){\phi }_1-</math><math>\frac{201}{25}\frac{{\Gamma }_{44}}{h^3}{\phi }_2+</math><math>5\frac{{\Gamma }_{44}}{h^3}{\phi }_3-\frac{{\Gamma }_{44}}{h^3}{\phi }_4=</math><math>B_ih-{\Gamma }_{22}\left(\frac{w_{10}w_{22}-w_{20}w_{12}}{w_{10}w_{21}-w_{20}w_{11}}\right)+</math><math>{\Gamma }_{44}\frac{72}{5h^2}\left(\frac{w_{10}w_{22}-w_{20}w_{12}}{w_{10}w_{21}-w_{20}w_{11}}\right)+</math><math>\frac{1024}{25h^3}\left(\frac{w_{12}w_{21}-w_{22}w_{11}}{w_{10}w_{21}-w_{20}w_{11}}\right)-</math><math>\frac{{\Gamma }_{22}}{24}{\frac{d^3\phi }{dx^3}}_{x=x_0}h^2-</math><math>{{\Gamma }_{44}\frac{47}{80}\frac{d^5\phi }{dx^5}}_{x=x_0}h^2+</math><math>\cdots </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (22)
|}

A análise que determinou a inconsistência da Eq. [[#eq0090|(18)]] foi realizada também na Eq. [[#eq0110|(22)]] e demonstrou que a Eq. [[#eq0110|(22)]] é consistente. Equações similares foram obtidas para o último volume do domínio, baseadas na mesma metodologia que apresentamos para o primeiro volume.

Para a determinação de uma equação para o segundo volume, foram utilizadas as seguintes equações:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>{\frac{d\phi }{dx}}_{x=h}=\frac{{\phi }_1-{\phi }_0}{h}-</math><math>\frac{1}{24}{\frac{d^3\phi }{dx^3}}_{x=h}h^2+\cdots </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (23)
|}

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>{\frac{d^3\phi }{dx^3}}_{x=h}=\frac{7{\phi }_1-2{\phi }_0-9{\phi }_2+5{\phi }_3-{\phi }_4}{h^3}+</math><math>{\frac{7}{8}\frac{d^5\phi }{dx^5}}_{x=h}h^2+\cdots </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (24)
|}

e

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>{\frac{d\phi }{dx}}_{x=2h}=\frac{{\phi }_2-{\phi }_1}{h}-</math><math>\frac{1}{24}{\frac{d^3\phi }{dx^3}}_{x=2h}h^2+\cdots </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (25)
|}

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>{\frac{d^3\phi }{dx^3}}_{x=2h}=\frac{{\phi }_3-{\phi }_0+3\left({\phi }_1-{\phi }_2\right)}{h^3}-</math><math>{\frac{1}{8}\frac{d^5\phi }{dx^5}}_{x=2h}h^2+\cdots </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (26)
|}

Substituindo estas equações na Eq. [[#eq0015|(3)]], chega‐se a:

<span id='eq0135'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\left(-\frac{{\Gamma }_{22}}{h}+\frac{{\Gamma }_{44}}{h^3}\right){\phi }_0+</math><math>\left(\frac{{\Gamma }_{22}}{h}-\frac{4{\Gamma }_{44}}{h^3}\right){\phi }_1+</math><math>\left(-\frac{{\Gamma }_{22}}{h}+\frac{6{\Gamma }_{44}}{h^3}\right){\phi }_2-</math><math>4\frac{{\Gamma }_{44}}{h^3}{\phi }_3+\frac{{\Gamma }_{44}}{h^3}{\phi }_4=</math><math>B_ih+\frac{{\Gamma }_{22}}{24}{\frac{d^3\phi }{dx^3}}_{x=x_1}h^2-</math><math>{{\Gamma }_{44}\frac{d^5\phi }{dx^5}}_{x=x_1}h^2+\cdots </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (27)
|}

que é a equação discretizada para o segundo volume. De forma análoga, obteve‐se uma equação semelhante para o penúltimo volume.

Apesar de as Eqs. [[#eq0110|(22)]] e [[#eq0135|(27)]] serem consistentes, os resultados numéricos obtidos utilizando estas equações foram invariavelmente insatisfatórios para as situações em que ''Γ''<sub>22</sub>/''Γ''<sub>44</sub> < 1. Os erros, medidos pela Eq. [[#eq0190|(38)]], foram da ordem ou superior a 10<sup>−1</sup> e, como será mostrado, erros desta ordem são muito superiores ao obtido com o modelo que será descrito a seguir. Estudos estão sendo elaborados para justificar corretamente as razões que levaram os erros a assumir valores altos assim.

A alternativa adotada foi elaborar uma equação especial para os 2 primeiros e para os 2 últimos volumes. Para isso, foram incluídos 2 volumes fictícios antes do primeiro volume e mais 2 após o último volume, como apresentado na [[#fig0015|figura 3]] e na [[#fig0020|figura 4]], respectivamente. Para cada um destes volumes foram elaboradas equações especiais que satisfazem apenas as condições de contorno.

<span id='fig0015'></span>

{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; max-width: 100%;" 
|-
|


[[Image:draft_Content_428379833-1-s2.0-S0213131516300232-gr3.jpg|center|225px|Esquema dos volumes fictícios em x=0.]]


|-
| <span style="text-align: center; font-size: 75%;">

Figura 3.

Esquema dos volumes fictícios em ''x'' = 0.

</span>
|}

<span id='fig0020'></span>

{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; max-width: 100%;" 
|-
|


[[Image:draft_Content_428379833-1-s2.0-S0213131516300232-gr4.jpg|center|225px|Esquema dos volumes fictícios em x=L.]]


|-
| <span style="text-align: center; font-size: 75%;">

Figura 4.

Esquema dos volumes fictícios em ''x'' = ''L''.

</span>
|}

O grande diferencial desta técnica é que, ao se empregar estes volumes fictícios, a Eq. [[#eq0050|(10)]] passa a ser válida para todos os volumes que não os fictícios.

A partir da Eq. [[#eq0065|(13)]], e utilizando‐se séries de Taylor, pode‐se obter as equações para os volumes fictícios, de tal forma que estes volumes satisfaçam exclusivamente as condições de contorno. Para o primeiro volume teremos que:

<span id='eq0140'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>{\phi }_{-2}=90{\phi }_0-20{\phi }_1+\frac{27}{5}{\phi }_2-</math><math>\frac{36}{49}{\phi }_3+\frac{192}{245}\left(\frac{\left(35w_{20}h-94w_{21}\right)w_{12}}{w_{10}w_{21}-w_{20}w_{11}}\right)-</math><math>\frac{192}{245}\left(\frac{\left(35w_{10}h-94w_{11}\right)w_{22}}{w_{10}w_{21}-w_{20}w_{11}}\right)</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (28)
|}

e para o segundo volume teremos:

<span id='eq0145'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>{\phi }_{-1}=4{\phi }_0-\frac{2}{3}{\phi }_1+\frac{4}{25}{\phi }_2-</math><math>\frac{1}{49}{\phi }_3+\frac{64}{3675}\left(\frac{\left(105w_{20}h-142w_{21}\right)w_{12}}{w_{10}w_{21}-w_{20}w_{11}}\right)-</math><math>\frac{64}{3675}\left(\frac{\left(105w_{10}h-142w_{11}\right)w_{22}}{w_{11}w_{21}-w_{20}w_{11}}\right)</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (29)
|}

De forma análoga, podem‐se obter as equações para os 2 últimos volumes fictícios considerando que as 2 condições de contorno podem ser definidas utilizando a Eq. [[#eq0150|(30)]]:

<span id='eq0150'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
| <math>e_{10}{\phi }_{x=L}+{e_{11}\frac{d\phi }{dx}}_{x=L}=</math><math>|-e_{12};{ }\mbox{e};{ }e_{20}{\phi }_{x=L}+</math><math>{e_{21}\frac{d\phi }{dx}}_{x=L}=e_{22}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (30)
|}

Desta forma, a equação do último volume, como observado na [[#fig0020|figura 4]], pode ser escrita como:

<span id='eq0155'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>{\phi }_{N+1}=90{\phi }_{N-1}-20{\phi }_{N-2}+\frac{27}{5}{\phi }_{N-3}-</math><math>\frac{36}{49}{\phi }_{N-4}-\frac{192}{245}\left(\frac{\left(35e_{20}h+94e_{21}\right)e_{12}}{e_{10}e_{21}-e_{20}e_{11}}\right)+</math><math>\frac{192}{245}\left(\frac{\left(35e_{10}h+94e_{11}\right)e_{22}}{e_{10}e_{21}-e_{20}e_{11}}\right)</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (31)
|}

e para o penúltimo volume teremos:

<span id='eq0160'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>{\phi }_N=4{\phi }_{N-1}-\frac{2}{3}{\phi }_{N-2}+</math><math>\frac{4}{25}{\phi }_{N-3}-\frac{1}{49}{\phi }_{N-4}-</math><math>\frac{64}{3675}\left(\frac{\left(105e_{20}h+142e_{21}\right)e_{12}}{e_{10}e_{21}-e_{20}e_{11}}\right)+</math><math>\frac{64}{3675}\left(\frac{\left(105e_{10}h+142e_{11}\right)e_{22}}{e_{10}e_{21}-e_{20}e_{11}}\right)</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (32)
|}

Com os valores das variáveis dos volumes fictícios calculados segundo as expressões apresentadas, uma única expressão precisará ser adotada para todos os demais volumes na montagem do sistema de equações lineares para a determinação de ''ϕ''<sub>''i''</sub>, a Eq. [[#eq0060|(12)]]. Os resultados utilizando volumes fictícios para esta equação foram bastante satisfatórios, demonstrando o acerto ao se escolher esse procedimento para este método.

==3. Modelo completo==

Observando as Eqs. [[#eq0030|(6)]] e [[#eq0035|(7)]], percebe‐se que as derivadas que foram desconsideradas na aproximação em série de Taylor são de terceira ordem, a mesma ordem das derivadas envolvidas no fluxo secundário. Pode‐se imaginar que dependendo do valor da relação ''Γ''<sub>22</sub>/''Γ''<sub>44</sub> esta aproximação para o fluxo primário pode comprometer a avaliação do fluxo secundário.

Como a modelagem do fluxo secundário exige que mais pontos sejam utilizados, visto que este fluxo envolve derivadas de terceira ordem, empregar mais pontos também na avaliação do fluxo primário não altera a forma da matriz resultante do sistema linear de equações envolvido na solução deste problema. Desta forma, poder‐se‐ia modelar o fluxo primário na forma:

<span id='eq0165'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>{\frac{d\phi }{dx}}_{x_i-\frac{1}{2}}=\frac{9}{8}\frac{{\phi }_i-{\phi }_{i-1}}{h}-</math><math>\frac{1}{24}\frac{{\phi }_{i+1}-{\phi }_{i-2}}{h}+</math><math>\frac{3}{640}{\frac{d^5\phi }{dx^5}}_{x_i-\frac{1}{2}}h^4+</math><math>\cdots </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (33)
|}

e

<span id='eq0170'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>{\frac{d\phi }{dx}}_{x_i+\frac{1}{2}}=\frac{9}{8}\frac{{\phi }_{i+1}-{\phi }_i}{h}-</math><math>\frac{1}{24}\frac{{\phi }_{i+2}-{\phi }_{i-1}}{h}+</math><math>\frac{3}{640}{\frac{d^5\phi }{dx^5}}_{x_i+\frac{1}{2}}h^4+</math><math>\cdots </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (34)
|}

Ao se descrever o fluxo primário utilizando‐se as Eqs. [[#eq0165|(33)]] e [[#eq0170|(34)]], estão sendo desconsideradas derivadas de quinta ordem, o que, em princípio, eliminaria as possibilidades da avaliação deste fluxo comprometer na modelagem do fluxo secundário. Com os fluxos primário e secundário calculados para as 2 faces do volume, chega‐se a equação de volumes finitos para o modelo completo:

<span id='eq0175'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\left(\frac{9}{4}\frac{{\Gamma }_{22}}{h}+6\frac{{\Gamma }_{44}}{h^3}\right){\phi }_i-</math><math>\left(\frac{7}{6}\frac{{\Gamma }_{22}}{h}+4\frac{{\Gamma }_{44}}{h^3}\right){\phi }_{i-1}-</math><math>\left(\frac{7}{6}\frac{{\Gamma }_{22}}{h}+4\frac{{\Gamma }_{44}}{h^3}\right){\phi }_{i+1}+</math><math>\left(\frac{1}{21}\frac{{\Gamma }_{22}}{h}+\frac{{\Gamma }_{44}}{h^3}\right){\phi }_{i-2}+</math><math>\left(\frac{1}{21}\frac{{\Gamma }_{22}}{h}+\frac{{\Gamma }_{44}}{h^3}\right){\phi }_{i+2}=</math><math>B_ih-{\frac{3{\Gamma }_{22}}{640}\frac{d^6\phi }{dx^6}}_{x_i}h^5-</math><math>{\frac{{\Gamma }_{44}}{8}\frac{d^6\phi }{dx^6}}_{x_i}h^3+</math><math>\cdots </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (35)
|}

Resolvendo o sistema de equações lineares descritos pelas Eqs. [[#eq0140|(28)]], [[#eq0145|(29)]], [[#eq0155|(31)]], [[#eq0160|(32)]] e [[#eq0175|(35)]], obtém‐se a solução para o modelo completo, ou utilizando‐se as Eqs. [[#eq0060|(12)]], [[#eq0140|(28)]], [[#eq0145|(29)]], [[#eq0155|(31)]] e [[#eq0160|(32)]] obtém‐se a solução para o modelo simplificado.

O sistema de equações lineares foi resolvido utilizando a biblioteca PETSC [[#bib0125|[10]]] ;  [[#bib0130|[11]]]. Foram realizados testes para determinar qual algoritmo para resolver o sistema de equações lineares e qual pré‐condicionador adotar entre as possibilidades existentes, sem que se chegasse a uma única opção. Utilizou‐se CGS [[#bib0135|[12]]] para resolver o sistema de equações lineares e fatoração incompleta (ILU) [[#bib0140|[13]]] como pré‐condicionador.

==4. Solução de referência==

Os resultados apresentados a seguir foram obtidos considerando‐se <math display="inline">B\left(x\right)</math> na Eq. [[#eq0010|(2)]], tal que:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>B\left(x\right)=\left(1+c_0\right)sen\left(n\pi \frac{x}{L}\right){\left(\frac{n\pi }{L^2}\right)}^2\left({\Gamma }_{22}L^2+\right. </math><math>\left. {\left(n\pi \right)}^2{\Gamma }_{44}\right)</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (36)
|}

Para este valor de <math display="inline">B\left(x\right)</math>, a solução analítica da Eq. [[#eq0010|(2)]] é:

<span id='eq0185'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\phi \left(x\right)=\left(1+c_0\right)sen\left(n\pi \frac{x}{L}\right)</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (37)
|}

Na Eq. [[#eq0170|(34)]], ''L'' é o comprimento do domínio, e ''c''<sub>0</sub> e ''n'' são 2 constantes concebidas para se analisar o comportamento da solução numérica, para amplitudes e frequências variadas.

A Eq. [[#eq0185|(37)]] foi utilizada como referência para medir a qualidade da solução numérica através da seguinte equação:

<span id='eq0190'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\epsilon =\frac{\sum_0^{N-1}{\left({\phi }_i-\phi \left(x_i\right)\right)}^2}{\sum_0^{N-1}{\phi }_i^2}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (38)
|}

Na Eq. [[#eq0190|(38)]], ''x''<sub>''i''</sub> e ''ϕ''<sub>''i''</sub> são, respectivamente: a coordenada do centro do volume e o valor da solução numérica para o i‐ésimo volume. Os valores calculados nos volumes fictícios não fazem parte da composição do erro apresentado na Eq. [[#eq0190|(38)]].

Para utilizar a Eq. [[#eq0170|(34)]] como solução de referência, devem‐se ajustar as condições de contorno para que esta equação também satisfaça à modelagem do sistema de equações discretizadas de balanço espacial. Desta forma, para valores quaisquer de <math display="inline">w_{i0}</math> e de <math display="inline">w_{i1}</math>, definidos na Eq. [[#eq0065|(13)]], tem‐se:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>w_{i2}=w_{i1}\left(1+c_0\right)\frac{n\pi }{L};\mbox{para};i=</math><math>0,1</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (39)
|}

e para valores quaisquer de ''e''<sub>''i''0</sub> e de ''e''<sub>''i''1</sub>, definidos na Eq. [[#eq0150|(30)]],

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>e_{i2}=\left(1+c_0\right)\left[e_{i0}sen\left(n\pi \right)+\right. </math><math>\left. e_{i1}\frac{n\pi }{L}cos\left(n\pi \right)\right];\mbox{para};i=</math><math>0,1</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (40)
|}

==5. Resultados==

Três situações diferentes foram elaboradas para testar a metodologia objeto deste trabalho. Cada uma das situações possui as seguintes características:
* Caso 1, ''c''<sub>0</sub> = 1 e ''n'' = 1;
* Caso 2, ''c''<sub>0</sub> = 1 e ''n'' = 10; e
* Caso 3, ''c''<sub>0</sub> = 10 e ''n'' = 1.

Como primeiro teste, investigou‐se a influência da relação <math display="inline">{\Gamma }_{22}/{\Gamma }_{44}</math> para 2 malhas distintas, uma com 40 volumes e outra com 400 volumes. Neste caso, procurou‐se verificar se o denominado modelo completo apresentaria resultados significativamente diferentes do modelo simplificado.

A [[#fig0025|figura 5]] mostra o resultado obtido para a simulação do caso 1. Para a situação em que <math display="inline">{\Gamma }_{22}/{\Gamma }_{44}>1</math>, observa‐se alguma diferença entre a solução obtida com modelo completo e a solução obtida com o modelo simplificado. Esta diferença aumenta com o aumento da relação <math display="inline">{\Gamma }_{22}/{\Gamma }_{44}</math>''.'' Como esperado, a solução com 400 volumes apresentou resultados mais precisos, mais próximos dos esperados quando comparados aos gerados com 40 volumes.

<span id='fig0025'></span>

{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; max-width: 100%;" 
|-
|


[[Image:draft_Content_428379833-1-s2.0-S0213131516300232-gr5.gif|center|282px|Influência de Γ22/Γ44 para o caso 1.]]


|-
| <span style="text-align: center; font-size: 75%;">

Figura 5.

Influência de <math display="inline">{\Gamma }_{22}/{\Gamma }_{44}</math> para o caso 1.

</span>
|}

Os resultados obtidos para o caso 2 são bastante semelhantes aos apresentados na [[#fig0025|figura 5]] e por esta razão não serão apresentados aqui. Devido à semelhança dos resultados obtidos com o caso 1 e com o caso 2, pode‐se concluir que a frequência não desempenha um papel relevante na formação do erro.

Os resultados obtidos para o caso 3, apresentados na [[#fig0030|figura 6]], exibem um comportamento distinto daqueles da [[#fig0025|figura 5]]. Pode‐se observar que o nível do erro para o caso 3 é superior àqueles obtidos nos casos 1 e 2, entretanto, na [[#fig0025|figura 5]], os maiores erros se localizam na região com os maiores valores de ''Γ''<sub>22</sub> > ''Γ''<sub>44</sub>. Na simulação do caso 3 ocorre exatamente o contrário. Por fim, observa‐se que quanto maior a amplitude menor é a diferença entre os erros do modelo completo e o modelo simplificado.

<span id='fig0030'></span>

{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; max-width: 100%;" 
|-
|


[[Image:draft_Content_428379833-1-s2.0-S0213131516300232-gr6.gif|center|281px|Influência de Γ22/Γ44 para o caso 3.]]


|-
| <span style="text-align: center; font-size: 75%;">

Figura 6.

Influência de <math display="inline">{\Gamma }_{22}/{\Gamma }_{44}</math> para o caso 3.

</span>
|}

O próximo teste foi elaborado para observar o comportamento da solução numérica em malhas variáveis para cada um dos casos analisados neste trabalho.

O primeiro resultado é apresentado na [[#fig0035|figura 7]] em que, para o caso 1, procurou‐se verificar se havia diferenças significativas entre as soluções obtidas utilizando‐se o modelo completo e o Simplificado. A conclusão deste teste é similar ao obtido no teste anterior, ou seja, diferenças na solução numérica obtida utilizando‐se os modelos Completo e Simplificado são observadas somente para as situações em que o fluxo primário predomina na solução.

<span id='fig0035'></span>

{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; max-width: 100%;" 
|-
|


[[Image:draft_Content_428379833-1-s2.0-S0213131516300232-gr7.gif|center|282px|Caso 1 – influência da variação da malha.]]


|-
| <span style="text-align: center; font-size: 75%;">

Figura 7.

Caso 1 – influência da variação da malha.

</span>
|}

Contudo, não passa despercebido que para este problema os erros de truncamento passam a ser relevantes em malhas relativamente grossas. No computador em que foi simulado este programa, somente para valores de ''h'' < 10<sup>−6</sup> é que os erros de truncamento passam a ser observados na análise do problema de difusão envolvendo apenas o fluxo primário, i. e., para aquele em que ''Γ''<sub>44</sub> = 0. Ao se considerar ''Γ''<sub>44</sub> ≠ 0, mesmo nos casos em que <math display="inline">{\Gamma }_{22}/{\Gamma }_{44}\gg 1</math> os efeitos de truncamento da solução numérica já são relevantes para ''h'' < 10<sup>−3</sup>.

A análise da variação da malha no comportamento do erro da variável foi realizada também para o caso 2. Novamente, os resultados foram bastante semelhantes aos apresentados na [[#fig0035|figura 7]] e não serão aqui apresentados.

Como já havíamos constatado na primeira sequência de testes, os resultados para o caso 3 são diferentes dos observados para os casos 1 e 2, ou seja, todos os valores de ''ɛ'' apresentados na [[#fig0040|figura 8]] são significativamente maiores dos que os apresentados na [[#fig0035|figura 7]] para os mesmos valores de ''h''. Desta forma, confirmou‐se que o valor da amplitude da solução influencia de forma significativa na qualidade da solução obtida.

<span id='fig0040'></span>

{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; max-width: 100%;" 
|-
|


[[Image:draft_Content_428379833-1-s2.0-S0213131516300232-gr8.gif|center|284px|Caso 3 – influência da variação da malha.]]


|-
| <span style="text-align: center; font-size: 75%;">

Figura 8.

Caso 3 – influência da variação da malha.

</span>
|}

Com os resultados previamente apresentados, podemos concluir que a modelagem proposta foi corretamente implementada no código computacional construído para a simulação deste problema e que estes resultados são bastante razoáveis, independente da relação <math display="inline">{\Gamma }_{22}/{\Gamma }_{44}</math> escolhida.

Utilizando‐se a metodologia detalhada em Roy [[#bib0145|[14]]], obteve‐se a ordem das aproximações numéricas para todos os resultados apresentados nas [[#fig0035|Figura 7]] ;  [[#fig0040|Figura 8]]. Os resultados foram, para todas as curvas, valores em torno de 2,0. Estes valores são esperados para as curvas obtidas com o modelo simplificado, mas esperava‐se uma ordem em torno de 3,0 para os resultados obtidos com o modelo completo. Esta discrepância já foi observada por Jones e Menzies [[#bib0075|[15]]] e serão apresentadas respostas a esta diferença num trabalho que se encontra em fase de elaboração, baseado na metodologia descrita por Jones e Menzies [[#bib0075|[15]]].

Para confirmar se o modelo funcionava adequadamente para condições adversas, foi realizado um teste onde os valores de −10 ≤ ''c''<sub>0</sub> ≤ 10 e de 10<sup>−3</sup> ≤ ''n'' ≤ 10 foram definidos aleatoriamente. Também de forma aleatória, as condições de contorno tanto do lado esquerdo quanto do lado direito do domínio foram definidas, dentro das seguintes possibilidades:
* Dirichlet/Dirichlet
* Dirichlet/Neumann
* Dirichlet/Robin
* Neumann/Robin

Os valores dos coeficientes das condições de contorno, Eqs. [[#eq0065|(13)]] e (30), foram obtidos aleatoriamente dentro do intervalo <math display="inline">\left[-10,10\right]</math>.

Nestas condições, executou‐se o programa elaborado para modelar a Eq. [[#eq0010|(2)]] por 50.000 vezes com os valores citados gerados aleatoriamente dentro do intervalo previamente definido. Os resultados obtidos podem ser vistos na [[#fig0045|figura 9]], onde cada ponto representa o erro da variável. O gráfico à esquerda mostra os resultados obtidos com o modelo simplificado e o da direita os obtidos com o modelo completo. Observa‐se nos 2 gráficos uma tendência da redução do erro para quando a difusão primária governa o problema, o que também já foi observado em resultados anteriores. Mas o que importa, ao se observar a [[#fig0045|figura 9]], é que para diferentes valores de condições de contorno e de ''c''<sub>0</sub> e de ''n  '', para diferentes razões de <math display="inline">{\Gamma }_{22}/{\Gamma }_{44}</math>, os erros obtidos estiveram sempre em níveis adequados.

<span id='fig0045'></span>

{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; max-width: 100%;" 
|-
|


[[Image:draft_Content_428379833-1-s2.0-S0213131516300232-gr9.gif|center|489px|Resultados obtidos para condições aleatórias – malha de 400 volumes.]]


|-
| <span style="text-align: center; font-size: 75%;">

Figura 9.

Resultados obtidos para condições aleatórias – malha de 400 volumes.

</span>
|}

O próximo resultado, apresentado na [[#fig0050|figura 10]], foi obtido de forma semelhante ao obtido no resultado anterior, i.e., os parâmetros que governam o problema foram definidos de forma aleatória. A única diferença da simulação que será apresentada a seguir e a apresentada anteriormente é que nesta simulação o valor de ''c''<sub>0</sub> = 10. Neste caso, fica patente que os resultados obtidos pelo modelo simplificado apresentam erros maiores que os obtidos pelo modelo completo nas situações em que <math display="inline">{\Gamma }_{22}/{\Gamma }_{44}>1</math>. Mas o que mais se depreende observando os resultados da [[#fig0050|figura 10]] é que os erros são bastante sensíveis ao valor de ''c''<sub>0</sub>, mais sensíveis a esta variável do que a qualquer outra que tenhamos observado nos testes elaborados com esse programa. Este comportamento já havia sido observado inclusive no resultado anterior, em que se pode ver claramente linhas com um mesmo erro. Estes resultados foram obtidos com valores de ''c''<sub>0</sub> semelhantes. Resultados semelhantes, e que não serão apresentados aqui, foram obtidos para outros valores de ''c''<sub>0</sub> constantes.

<span id='fig0050'></span>

{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; max-width: 100%;" 
|-
|


[[Image:draft_Content_428379833-1-s2.0-S0213131516300232-gr10.gif|center|288px|Resultados obtidos para condições aleatórias – malha com 400 volumes e c0=10.]]


|-
| <span style="text-align: center; font-size: 75%;">

Figura 10.

Resultados obtidos para condições aleatórias – malha com 400 volumes e ''c''<sub>0</sub> = 10.

</span>
|}

Os resultados apresentados na [[#fig0045|figura 9]] e na [[#fig0050|figura 10]] nos permitem afirmar que a equação objeto deste estudo foi adequadamente modelada pelo MVF para condições de contorno variáveis e com valores bastante heterogêneos.

==6. Conclusões==

Neste trabalho, apresentou‐se uma metodologia baseada no MVF para estudar a equação da difusão com fluxo bimodal em regime permanente, equação esta proposta por Bevilacqua et al. [[#bib0105|[6]]]; [[#bib0110|[7]]] ;  [[#bib0115|[8]]]. Determinaram‐se também alguns caminhos que foram utilizados para se resolver este problema, caminhos estes que nem sempre conduziram a soluções satisfatórias. Cita‐se o exemplo das aproximações que levaram a Eq. [[#eq0090|(18)]] e que, como demonstrado, é uma solução aproximada inconsistente. Dependendo como se faz as aproximações das equações discretizadas de volumes finitos, estas equações podem vir a ser inconsistentes

Os testes elaborados e os resultados apresentados nos permitem dizer que o método descrito neste trabalho funciona de forma adequada para a solução de diferentes casos do problema em análise, que vão desde problemas dominados pelo fluxo primário aos problemas em que o fluxo secundário é bastante superior ao primário.

Destacamos que a amplitude da solução é o parâmetro mais relevante para definir a ordem do erro da solução numérica.

Podem haver diferenças importantes entre a solução obtida utilizando o modelo completo e o modelo simplificado. Por outro lado, não encontramos nada que justificasse a aplicação do modelo simplificado, visto que os pontos extras utilizados na modelagem do modelo completo são os mesmos necessários na discretização do fluxo secundário. Sendo assim, em princípio, o esforço computacional para resolver a equação do modelo completo é da mesma ordem do modelo simplificado.

==Referências==

<ol style='list-style-type: none;margin-left: 0px;'><li><span id='bib0080'></span>
[[#bib0080|[1]]] J. Wu, K.M. Berland; Propagators and Time‐Dependent Diffusion Coefficients for Anomalous Diffusion; Biophys J., 95 (4) (2008), pp. 2049–2052</li>
<li><span id='bib0085'></span>
[[#bib0085|[2]]] M.J. Saxton; Anomalous diffusion due to obstacles: A Monte Carlo study; Biophys J., 66 (2) (1994), pp. 394–401 Part 1</li>
<li><span id='bib0090'></span>
[[#bib0090|[3]]] J.I. Ramos; On the numerical treatment of an ordinary differential equation arising in one‐dimensional non‐Fickian diffusion problems; Comput Phys Commun., 170 (3) (2005), pp. 231–238</li>
<li><span id='bib0095'></span>
[[#bib0095|[4]]] P.F. Nealey, R.E. Cohen, A.S. Argon; Limited‐supply non‐Fickian diffusion in glassy polymers; Polymer (Guildf)., 36 (19) (1995), pp. 3687–3695</li>
<li><span id='bib0100'></span>
[[#bib0100|[5]]] P.M. Smith e, M.M. Fisher; Non‐Fickian diffusion of water in melamine‐formaldehyde resins; Polymer (Guildf)., 25 (1) (1984), pp. 84–90</li>
<li><span id='bib0105'></span>
[[#bib0105|[6]]] L. Bevilacqua, A.C.N.R. Galeão, F.P. Costa; A new analytical formulation of retention effects on particle diffusion processes; An Acad Bras Cienc., 83 (2011), pp. 1443–1464</li>
<li><span id='bib0110'></span>
[[#bib0110|[7]]] L. Bevilacqua, A.C.N.R. Galeão, F.P. Costa; On the significance of higher order differential terms in diffusion processes; J Brazilian Soc Mech Sci Eng., 33 (2011), pp. 166–175</li>
<li><span id='bib0115'></span>
[[#bib0115|[8]]] L. Bevilacqua, A.C.N.R. Galeão, J.G. Simas, A.P.R. Doce; A new theory for anomalous diffusion with a bimodal flux distribution; J Brazilian Soc Mech Sci Eng. (2013), pp. 1–10</li>
<li><span id='bib0120'></span>
[[#bib0120|[9]]] C.R. Maliska; Transferência de Calor e Mecânica dos Fluidos Computacional; LTC Editora (1994)</li>
<li><span id='bib0125'></span>
[[#bib0125|[10]]] S. Balay, J. Brown, K. Buschelman, W.D. Gropp, D. Kaushik, M.G. Knepley, L.C. McInnes, B.F. Smith, e H. Zhang, «PETSc Web page», 2011.</li>
<li><span id='bib0130'></span>
[[#bib0130|[11]]] S. Balay, W.D. Gropp, L.C. McInnes, e B. F. Smith, «Efficient Management of Parallelism in Object Oriented Numerical Software Libraries» in ''Modern Software Tools in Scientific Computing'', 1997, pp. 163‐202.</li>
<li><span id='bib0135'></span>
[[#bib0135|[12]]] P. Sonneveld; CGS, A Fast Lanczos‐Type Solver for Nonsymmetric Linear systems; SIAM J Sci Stat Comput, 10 (1) (1989), pp. 36–52</li>
<li><span id='bib0140'></span>
[[#bib0140|[13]]] T. Dupont, R.P. Kendall, H.H. Rachford Jr.; An Approximate Factorization Procedure for Solving Self‐Adjoint Elliptic Difference Equations; SIAM J Numer Anal., 5 (3) (1968), pp. 559–573</li>
<li><span id='bib0145'></span>
[[#bib0145|[14]]] C.J. Roy; Review of code and solution verification procedures for computational simulation; J Comput Phys., 205 (1) (2005), pp. 131–156</li>
<li><span id='bib0075'></span>
[[#bib0075|[15]]] W.P. Jones e, K.R. Menzies; Analysis of the Cell‐Centred Finite Volume Method for the Diffusion Equation; J Comput Phys., 165 (1) (2000), pp. 45–68</li>
</ol>