Estudo do Coeficiente de Difusão Secundária em Problema de Difusão com Fluxo Bimodal

RESUMO

Uma formulação recentemente desenvolvida para o problema de difusão anômala, a qual emprega um termo diferencial de quarta ordem, apresentou em determinadas situações observadas na literatura valores fisicamente irreais em sua solução. Neste trabalho é realizado um estudo do efeito coeficiente de difusão secundária visando contribuir para o entendimento do comportamento das soluções nessas situações. É implementada uma função para representar a variação da parcela sujeita a difusão primária e secundária, de acordo com a quantidade da propriedade em difusão. É apresentada a formulação para a solução desse novo modelo pelo Método de Diferenças Finitas. Os resultados obtidos são compatíveis com aqueles apresentados em trabalhos anteriores na literatura.

Palavras-chave:
difusão anômala; Método de Diferenças Finitas; equação diferencial de quarta ordem

ABSTRACT

A recent formulation for anomalous diffusion problems, which envolves a fourth-order differential term, has presented negative values in its solution for certain cases. This paper presents a study of the effect of secondary diffusion coefficient in order to contribute with the understanding of the solutions behavior in these situations. A function was implemented to represent the variation in the parcel subject to primary and secondary diffusion, according to the amount of the property in diffusion. The formulation for the solution of this new model by the Finite Difference Method is presented. The results obtained are compatible with those presented in previous works in the literature.

Keywords:
anomalous difusion; Finite Difference Method; fourth-order differential equation

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 ^12),(15),(25. 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^{\lambda }$, com $\lambda =1$ ²⁵.

Para aquelas situações em que $\lambda \ne 1$ têm-se o que se denomina de difusão anômala e há na literatura diversos modelos, cada qual com características, aplicações e métodos próprios para estudá-los ^{5), (10}.

A difusão com retenção é aquela em que o fluxo é menor do que aquele que seria esperado quando se tem uma escala linear do descolamento com o tempo, isto é, $\lambda =1$. Segundo Simas ²⁰, a difusão com retenção tem sido modelada com três maneiras distintas: incluindo termos fontes e sumidouros na equação da difusão clássica, como o fez Atsunami ¹; introduzindo um mecanismo de retardo na difusão e ajustando os coeficientes deste mecanismo a partir de experimentais do problema, como observado em D’Angelo e colaboradores ⁶; definindo-se um coeficiente de difusão não linear e o ajustando a partir de dados experimentais, como o visto no trabalho de Muhammad ¹⁶.

Na área de Transferência de Calor, uma das formas de se similar a retenção do fluxo em um problema de condução de calor se dá pela inclusão de um termo com derivada segunda no tempo, como na equação a seguir:

$\alpha {\nabla }^{2}T=\frac{\partial T}{\partial t}+\tau \frac{{\partial }^{2}T}{\partial t^2}$ (1.1)

em que α é a difusividade térmica, T a temperatura e τ o tempo de relaxação. Quando $\tau =0$, recupera-se a equação de condução de calor clássica. Há ainda muitas outras formas para se modelar a retenção do fluxo de calor, como observado em Tzou ²¹.

Equações da difusão com derivadas de ordem fracionária no espaço ou no tempo, ou no espaço e no tempo simultaneamente têm sido utilizadas para modelar processos de difusão anômalas ^{2), (25), (26}.

Neste trabalho será investigado o modelo de difusão anômala proposto por Bevilacqua e colaboradores ^{3), (4}, denominado de difusão bimodal ou de difusão bi-fluxo. Diferentemente de muitos outros modelos de difusão anômala, o modelo de difusão bi-fluxo é um modelo fenomenológico, isto é, ele foi obtido a partir de equações básicas de Fenômenos de Transporte.

Bevilacqua e colaboradores ^3),(4 apresentaram a formulação matemática para problemas de dispersão considerando um novo parâmetro β que representa a fração de partículas capaz de se difundir e que é igualmente distribuída às células vizinhas, enquanto a fração $1-\beta$ permanece retida temporariamente em cada passo de tempo ∆t.

Seguindo esta abordagem, e após uma cuidadosa manipulação algébrica, chega-se a seguinte expressão:

$\frac{\partial \varphi }{\partial t}=\beta {\kappa }_2\frac{{\partial }^2\varphi }{\partial x^2}-\beta \left(1-\beta \right){\kappa }_4\frac{{\partial }^4\varphi }{\partial x^4}$ (1.2)

em que $\varphi =\varphi \left(x,t\right)$ é a variável de interesse, x é a coordenada espacial, t é o tempo, com $0<\beta \le 1$. O efeito de retenção, como definido por Bevilacqua et al.³, está representado na Eq. (1.3) pelo termo em que há a derivada de quarta ordem. ${\kappa }_2>0$ é o coeficiente associado ao fluxo primário e ${\kappa }_4\ge 0$ ao fluxo secundário.

Nesta nova abordagem, o processo de difusão está associado a dois fluxos, sejam eles de matéria, momento ou energia, um correspondente à fração β, denominado fluxo primário e caracterizado pela derivada de primeira ordem, e outro à fração $\left(1-\beta \right)$, denominado fluxo secundário e caracterizado pela derivada de terceira ordem. Dessa forma, nesta formulação, a equação da difusão é reescrita, em sua versão unidimensional, como:

$\frac{\partial \varphi }{\partial t}=\frac{\partial }{\partial x}\left({\kappa }_2\beta \frac{\partial \varphi }{\partial x}\right)-\frac{\partial }{\partial x}\left({\kappa }_4\beta \left(1-\beta \right)\frac{{\partial }^3\varphi }{\partial x^3}\right)$ (1.3)

Resultados obtidos da simulação da Eq. (1.3), ou da sua versão em duas dimensões, foram obtidos por diversos autores empregando variadas metodologias para solução numérica de equações diferenciais, tais como o Método de Volumes Finitos ²⁴ para o modelo unidimensional e estacionário, a Técnica da Transformada Integral Generalizada ²³ e o Método de Elementos Finitos ⁹. Contudo, em casos já documentados ^{3), (8), (19}, valores fisicamente irreais foram observados em decorrência de uma aceleração inicial no fluxo associada às condições inicial e de contorno.

No presente trabalho buscou-se contribuir para melhorar o entendimento das situações em que resultados fisicamente inconsistente surgem, reproduzindo alguns casos com tais ocorrências e propondo uma variação nos valores do coeficiente β para que situações fisicamente irreais deixem de ocorrer.

Os detalhes de uma nova formulação para a Eq. (1.3) considerando $\beta =\beta \left(\varphi \right)$ são apresentados na seção 2. Na seção 3 são detalhados os aspectos da discretização dessa nova equação pelo Método de Diferenças Finitas ¹¹. Na seção 4 é apresentada uma análise da metodologia apresentada empregando a técnica de solução manufaturada ^{17), (18} a fim de certificar sua acurácia, bem como os resultados para diferentes conjuntos de condição de contorno e inicial. Por fim na seção 5 são encontradas as conclusões deste trabalho.

2 FLUXO SECUNDÁRIO ADAPTATIVO

Considere o seguinte problema que pode ser modelado pela Eq. (1.3): difusão de um soluto em um meio unidimensional qualquer. Definindo ϕ(x,t) como a concentração deste soluto em um ponto e tempo qualquer, então a simulação deste problema deverá, obrigatoriamente, fornecer valores de $\varphi \left(x,t\right)\ge 0$, pois $\varphi \left(x,t\right)<0$ para este problema não faz sentido físico. Contudo, ao se simular numericamente a Eq. (1.3), dependendo da combinação dos valores de κ ₂, κ ₄, β, das condições iniciais e das condições de contorno pode-se obter valores $\varphi \left(x,t\right)<0$ ³.

Para evitar situações como esta, um ajuste foi elaborado na Eq. (1.3): fazer $\beta \to 1$ sempre que ϕ tender a um valor fora de limites previamente especificados. Para alcançar este objetivo considerou-se que na Eq. (1.3) tanto κ2 como κ ₄ são valores constantes, porém fez-se $\beta =\beta \left(\varphi \right)$ e se reescreveu a Eq. (1.3) considerando a variação de β:

$\frac{\partial \varphi }{\partial t}={\kappa }_2\frac{\partial }{\partial x}\left(\beta \frac{\partial \varphi }{\partial x}\right)-{\kappa }_4\frac{\partial }{\partial x}\left(\beta \left(1-\beta \right)\frac{{\partial }^3\varphi }{\partial x^3}\right)$ (2.1)

a qual, após alguma manipulação algébrica, assume a seguite forma:

$\frac{\partial \varphi }{\partial t}={\kappa }_2\left(\frac{\partial \beta }{\partial x}\frac{\partial \varphi }{\partial x}+\beta \frac{{\partial }^2\phi }{\partial x^2}\right)-{\kappa }_4\left(\left(1-2\beta \right)\frac{\partial \beta }{\partial x}\frac{{\partial }^3\varphi }{\partial x^3}+\beta \left(1-\beta \right)\frac{{\partial }^4\varphi }{\partial x^4}\right)$ (2.2)

Neste trabalho será apresentada a formulação matemática para a solução da Eq. (2.2) empregando o Método de Diferenças Finitas ¹¹. Detalhes sobre a discretização desta equação serão tratados na seção 3. Na próxima seção são discutidas algumas das características da função β(ϕ) empregada nas simulações apresentadas na seção 4.

2.1 Função para o fluxo secundário adaptativo

A função β(ϕ) a ser escolhida deve possuir a seguinte propriedade: ${\beta }_{min}\le \beta \left(\varphi \right)\le {\beta }_{max}$ e $\beta \left(\varphi \right)\to {\beta }_{min}$ quando $\varphi \to {\varphi }_{min}$. Ao fazer $\beta \left(\varphi \right)\to {\beta }_{min}$ se está eliminando o efeito da difusão secundária neste modelo e, desta forma, valores fisicamente irreais de ϕ não deverão surgir na simulação empregando-se a Eq. (2.2).

Após alguns teste escolheu-se a seguinte função:

$\beta \left(\varphi \right)={\beta }_{max}-\frac{{\beta }_{max}-{\beta }_{min}}{1+e^{\left[\gamma \left({\varphi }_{min}-\varphi \right)\right]}}$ (2.3)

com $0<{\beta }_{min}<{\beta }_{max}<1,{\varphi }_{min}$ o limite inferior para ϕ e $\gamma >0$ é o parâmetro que define a sensibilidade de β com a variável ϕ.

Como a Eq. (2.3) fornece $\beta =\beta \left(\varphi \right)$ e a Eq. (2.2) necessita de $\frac{\partial \beta }{\partial x}$ então esta derivada poderá ser calculada a partir da seguinte transformação:

$\frac{\partial \beta }{\partial x}=\frac{\partial \beta }{\partial \varphi }\frac{\partial \varphi }{\partial x}$ (2.4)

com

$\frac{\partial \beta }{\partial \varphi }=\gamma \left({\beta }_{max}-{\beta }_{min}\right)\frac{e^{\gamma \left({\varphi }_{min}-\varphi \right)}}{{\left(1+e^{\gamma \left({\varphi }_{min}-\varphi \right)}\right)}^2}$ (2.5)

A Eq. (2.2) pode ser reescrita então da seguinte maneira:

$\frac{\partial \varphi }{\partial t}=\left({\Lambda }_2\frac{\partial \varphi }{\partial x}+{\Lambda }_{22}\frac{{\partial }^2\varphi }{\partial x^2}\right)-\left({\Lambda }_4\frac{{\partial }^3\varphi }{\partial x^3}+{\Lambda }_{44}\frac{{\partial }^4\varphi }{\partial x^4}\right)$ (2.6)

em que

$$\begin{gathered} {\Lambda }_2={\kappa }_2\frac{\partial \beta }{\partial \varphi }\frac{\partial \varphi }{\partial x},{\Lambda }_{22}={\kappa }_2\beta \\ {\Lambda }_4={\kappa }_4\left(1-2\beta \right)\frac{\partial \beta }{\partial \varphi }\frac{\partial \varphi }{\partial x},{\Lambda }_{44}={\kappa }_4\beta \left(1-\beta \right) \end{gathered}$$ (2.7)

Na Figura 1 está representada a função β(ϕ) considerando ${\beta }_{min}=0.8,{\beta }_{max}=1.0$ e ${\varphi }_{min}=0.001$. Para $\gamma =100$ a região de β variável é bastante extensa, enquanto para valores de $\gamma >800$ esta região passa a ser bastante reduzida. Neste exemplo, sempre que $\beta >0.8$ o peso do fluxo secundário nesta região está sendo reduzido.

Figura 1:
Comportamento de β(ϕ).

3 MODELO NUMÉRICO

Após ter sido definida a equação alvo deste trabalho, Eq. (2.6), a sequência natural é resolvêla. Contudo, se existir uma solução analítica desta equação, considerando β(ϕ) como proposto na Eq. (2.3), muito provavelmente esta solução só será alcançada empregando-se recursos matemáticos bastante sofisticados. Desta forma, para se resolver esta equação deve-se utilizar um método numérico para solução de equações diferencias parciais. O método de Diferenças Finitas ¹¹ foi a escolha dos autores para este propósito.

O domínio de cálculo, $0\le x\le L\in \mathfrak{R}$, foi subdividido em N − 1 intervalos igualmente espaçados de comprimento h e obteve-se um sistema linear de equações com N incógnitas.

Como a Eq. (2.6) demanda quatro condições de contorno, utilizou-se os dois primeiros e os dois últimos nós da malha para montar equações com informações sobre as condições de contorno e com exceção destes quatro nós, todos os demais são considerados nós internos e as equações discretizadas para estes nós serão apresentadas a seguir.

3.1 Equação discretizada para os nós internos

O Método de Diferenças Finitas ¹¹ baseia-se na substituição dos termos da equação diferencial por aproximações obtidas utilizando-se série de Taylor. Sendo assim, os termos $\frac{\partial \varphi }{\partial x}$ da Eq. (2.6) podem ser aproximados para a malha unidimensional igualmente espaçada como:

$\frac{\partial \varphi }{\partial x}{|}_{x_i}=\frac{{\varphi }_{i-2}-8{\varphi }_{i-1}+8{\varphi }_{i+1}-{\varphi }_{i+2}}{12h}+O\left(h^4\right)$ (3.1)

em que $3\le i\le N-2$. Usualmente emprega-se somente três nós na discretização de derivadas de primeira ordem, contudo, como a derivada de quarta ordem exige cinco nós em sua discretização, utilizou-se a mesma quantidade de pontos nas expansões em Séries de Taylor empregadas na discretização das demais derivadas desta equação. A derivada de segunda ordem foi aproximada como:

$\frac{{\partial }^2\varphi }{\partial x^2}{|}_{x_i}=\frac{-{\varphi }_{i-2}+16{\varphi }_{i-1}-30{\varphi }_i+16{\varphi }_{i+1}-{\varphi }_{i+2}}{\overline{)12h^2}}+O\left(h^4\right)$ (3.2)

A derivada de terceira ordem foi aproximada segundo a seguinte expressão:

$\frac{{\partial }^3\varphi }{\partial x^3}{|}_{x_i}=\frac{-{\varphi }_{i-2}+2{\varphi }_{i-1}-2{\varphi }_{i+1}+{\varphi }_{i+2}}{2h^3}+O\left(h^2\right)$ (3.3)

e, por fim, a derivada de quarta ordem fica:

$\frac{{\partial }^4\varphi }{\partial x^4}{|}_{x_i}=\frac{{\varphi }_{i-2}-4{\varphi }_{i-1}+6{\varphi }_i-4{\varphi }_{i+1}+{\varphi }_{i+2}}{h^4}+O\left(h^2\right)$ (3.4)

O termo transiente foi modelado utilizado-se a aproximação tradicionalmente empregada na modelagem deste termo:

$\frac{\partial \varphi }{\partial t}{|}_{x_i,t+\Delta t}=\frac{{\varphi }_i^{t+\Delta t}-{\varphi }_i^t}{\Delta t}+O\left(\Delta t\right)$ (3.5)

A equação na forma discretizada para todos os nós internos é a seguinte:

$a_{i-2}{\varphi }_{i-2}^{t+\Delta t}+a_{i-1}{\varphi }_{i-1}^{t+\Delta t}+a_i{\varphi }_i^{t+\Delta t}{a}_{i+1}{\varphi }_{i+1}^{t+\Delta t}+a_{i+2}{\varphi }_{i+2}^{t+\Delta t}=S_i$ (3.6)

em que os coeficientes desta equação são:

$$\begin{gathered} a_i=\frac{1}{\Delta t}+\frac{5{\lambda }_{22}}{2}+6{\lambda }_{44},S_i=\frac{1}{\Delta t}{\varphi }_i^t \\ {a}_{i-2}=-\frac{{\lambda }_2}{12}+\frac{{\lambda }_{22}}{12}-\frac{{\lambda }_4}{2}+{\lambda }_{44},a_{i-1}=\frac{2{\lambda }_2}{3}-\frac{4{\lambda }_{22}}{3}+{\lambda }_4-4{\lambda }_{44} \\ {a}_{i+2}=\frac{{\lambda }_2}{12}+\frac{{\lambda }_{22}}{12}+\frac{{\lambda }_4}{2}+{\lambda }_{44},a_{i+1}=-\frac{2{\lambda }_2}{3}-\frac{4{\lambda }_{22}}{3}-{\lambda }_4-4{\lambda }_{44} \end{gathered}$$ (3.7)

em que ${\lambda }_2={\Lambda }_2/h$, e ${\lambda }_{22}={\Lambda }_{22}/h^2,{\lambda }_4={\Lambda }_4/h^3$, e ${\lambda }_{44}={\Lambda }_{44}/h^4$.

Optou-se por uma formulação totalmente implícita ¹³ na discretização temporal da Eq. (3.6). De forma a evitar não lineariedades no sistema de equações, a avaliação dos termos termos β(ϕ) e $\frac{\partial \beta }{\partial \varphi }$ presentes nos coeficientes Λ₂ e Λ₄ foi realizada empregando ϕ^t.

3.2 Condições de contorno

No caso do problema em análise, Eq. (2.6), para que este seja bem definido em um domínio finito, quatro condições de contorno são necessárias. Para a primeira condição de contorno em $x=0$ será utilizada a seguinte equação:

$w_{10}\varphi \left(x=0\right)+w_{11}\frac{\partial \varphi }{\partial x}{|}_{x=0}+w_{12}\frac{{\partial }^2\varphi }{\partial x^2}{|}_{x=0}=w_{13}\left(t\right)$ (3.8)

em que w ₁₀, w ₁₁, w ₁₂ são valores conhecidos e constantes. Por exemplo, se a primeira condição de contorno em $x=0\text{for}\varphi \left(x\right)=2$ (Dirichlet) então $w_{10}=1,w_{11}=0,w_{12}=0$ e $w_{13}\left(t\right)=2$.

A segunda condição de contorno em <mml:math><mml:mi>x</mml:mi><mml:mo>=</mml:mo><mml:mn>0</mml:mn></mml:math> é definida de forma similar ao que foi feito para a primeira condição de contorno:

$w_{20}\varphi \left(x=0\right)+w_{21}\frac{\partial \varphi }{\partial x}{|}_{x=0}+w_{22}\frac{{\partial }^2\varphi }{\partial x^2}{|}_{x=0}=w_{23}\left(t\right)$ (3.9)

sendo w ₂₀, w ₂₁, w ₂₂ valores conhecidos e constantes.

As condições de contorno em $x=L$ também são expressas de forma análoga:

$e_{10}\varphi \left(x=L\right)+e_{11}\frac{\partial \varphi }{\partial x}{|}_{x=L}+e_{12}\frac{{\partial }^2\varphi }{\partial x^2}{|}_{x=L}=e_{13}\left(t\right)$ (3.10)

$e_{20}\varphi \left(x=L\right)+e_{21}\frac{\partial \varphi }{\partial x}{|}_{x=L}+e_{22}\frac{{\partial }^2\varphi }{\partial x^2}{|}_{x=L}=e_{23}\left(t\right)$ (3.11)

com e ₁₀, e ₁₁, e ₁₂, e ₂₀, e ₂₁, e ₂₂ sendo valores conhecidos e constantes.

As Eqs. (3.8) e (3.9) serão utilizadas na obtenção das equações discretizadas dos dois primeiros nós e as Eqs. (3.10) e (3.11) na discretização dos dois últimos.

3.2.1 Primeiro nó

A equação para o primeiro nó x ₁ foi obtida utilizando a Eq. (3.8). Foram realizadas novas expansões em Série de Taylor em torno dos nós posteriores ao primeiro nó obtendo-se as seguintes expressões para as derivadas de primeira e de segunda ordem, respectivamente:

$\frac{\partial \varphi }{\partial x}{|}_{x_1}=-\frac{25{\varphi }_1}{12h}+\frac{4{\varphi }_2}{h}-\frac{3{\varphi }_3}{h}+\frac{4{\varphi }_4}{3h}-\frac{{\varphi }_5}{4h}+O\left(h^4\right)$ (3.12)

$\frac{{\partial }^2\varphi }{\partial x^2}{|}_{x_1}=\frac{35{\varphi }_1}{12h^2}-\frac{26{\varphi }_2}{3h^2}+\frac{19{\varphi }_3}{2h^2}-\frac{14{\varphi }_4}{3h^2}+\frac{11{\varphi }_5}{12h^2}+O\left(h^3\right)$ (3.13)

Substituindo tais expressões na Eq. (3.8), a equação discretizada para o primeiro nó é dada por:

$a_1{\varphi }_1^{t+\Delta t}+a_2{\varphi }_2^{t+\Delta t}{a}_3{\varphi }_3^{t+\Delta t}+a_4{\varphi }_4^{t+\Delta t}+a_5{\varphi }_5^{t+\Delta t}=S_1$ (3.14)

na qual:

$$\begin{gathered} a_1=w_{10}-\frac{25w_{11}}{12h}+\frac{35w_{12}}{12h^2},a_2=\frac{4w_{11}}{h}-\frac{26w_{12}}{3h^2} \\ {a}_3=-\frac{3w_{11}}{h}+\frac{19w_{12}}{2h^2},a_4=\frac{4w_{11}}{3h}-\frac{14w_{12}}{3h^2} \\ {a}_5=-\frac{w_{11}}{4h}+\frac{11w_{12}}{12h^2},S_1=w_{13}\left(t\right) \end{gathered}$$ (3.15)

3.2.2 Segundo nó

O segundo nó, x ₂, deve satisfazer a Eq. (2.6) fazendo-se necessário obter novas expressões para as derivadas presentes em tal equação. Para tanto, primeiramente, são realizadas novas expansões em séries de Taylor em torno do ponto $i=2$ considerando três pontos a frente, $i=3,4\text{e}5$. Da mesma forma as das derivadas presentes na Eq. (3.9), que define uma das condições de contorno em $x=0$, são substituídas por expansões em Séries de Taylor. A partir desse conjunto de equações chega-se, então, às expressões para as derivadas de primeira, de segunda, de terceira e de quarta ordens, conforme equações a seguir:

$$\begin{gathered} D_{w2}\frac{\partial \varphi }{\partial x}{|}_{x_2}=-3h{w}_{23}\left(t\right)+\left[-10{\varphi }_2+18{\varphi }_3-6{\varphi }_4+{\varphi }_5\right]h{w}_{20} \\ \left[\frac{197}{6}{\varphi }_2-\frac{93}{2}{\varphi }_3+\frac{33}{2}{\varphi }_4-\frac{17}{6}{\varphi }_5\right]w_{21}+ \\ \left[-\frac{331}{6}{\varphi }_2+81{\varphi }_3-\frac{63}{2}{\varphi }_4+\frac{17}{3}{\varphi }_5\right]\frac{w_{22}}{h}+O\left(h^4\right) \end{gathered}$$ (3.16)

$$\begin{gathered} D_{w2}\frac{{\partial }^2\varphi }{\partial x^2}{|}_{x_2}=11h{w}_{23}\left(t\right)+\left[-20{\varphi }_2+6{\varphi }_3+4\varphi -{\varphi }_5\right]w_{20} \\ \left[-\frac{7}{3}{\varphi }_2+\frac{41}{2}{\varphi }_3-23{\varphi }_4+\frac{29}{6}{\varphi }_5\right]\frac{w_{21}}{h}+ \\ \left[37{\varphi }_2-87{\varphi }_3+63{\varphi }_4-13{\varphi }_5\right]\frac{w_{22}}{h^2}+O\left(h^3\right) \end{gathered}$$ (3.17)

$$\begin{gathered} D_{w2}\frac{{\partial }^3\varphi }{\partial x^3}{|}_{x_2}=\frac{18w_{23}\left(t\right)}{h}+\left[60{\varphi }_2-72{\varphi }_3+36{\varphi }_4-6{\varphi }_5\right]\frac{w_{20}}{h} \\ \left[-53{\varphi }_2+96{\varphi }_3-51{\varphi }_4+8{\varphi }_5\right]\frac{w_{21}}{h^2}+ \\ \left[19{\varphi }_2-39{\varphi }_3+21{\varphi }_4-{\varphi }_5\right]\frac{w_{22}}{h^3}+O\left(h^2\right) \end{gathered}$$ (3.18)

$$\begin{gathered} D_{w2}\frac{{\partial }^4\varphi }{\partial x^4}{|}_{x_2}=\frac{12w_{23}\left(t\right)}{h^2}+\left[-48{\varphi }_2+72{\varphi }_3-48{\varphi }_4+12{\varphi }_5\right]\frac{w_{20}}{h^2} \\ \left[52{\varphi }_2-114{\varphi }_3+84{\varphi }_4-22{\varphi }_5\right]\frac{w_{21}}{h^3}+ \\ \left[-36{\varphi }_2+96{\varphi }_3-84{\varphi }_4+24{\varphi }_5\right]\frac{w_{22}}{h^4}+O\left(h\right) \end{gathered}$$ (3.19)

Nas equações anteriores, considerou-se $D_{w2}=10w_{20}{h}^2-25w_{21}h+35w_{22}\ne 0$. Substituindo-se estas expressões na Eq. (2.6), obtém-se a equação discretizada para o segundo nó:

$a_2{\varphi }_2^{t+\Delta t}+a_3{\varphi }_3^{t+\Delta t}+a_4{\varphi }_4^{t+\Delta t}+a_5{\varphi }_5^{t+\Delta t}=S_2$ (3.20)

na qual os coeficientes a _i são definidos por:

$$\begin{gathered} a_2=\frac{1}{\Delta t}+\left[10{\Lambda }_{2}h+20{\Lambda }_{22}+\frac{60{\Lambda }_4}{h}-\frac{48{\Lambda }_{44}}{h^2}\right]\frac{w_{20}}{D_{w2}}+\left[-\frac{197{\Lambda }_2}{6}+\frac{7{\Lambda }_{22}}{3h}-\frac{53{\Lambda }_4}{h^2}+\frac{52{\Lambda }_{44}}{h^3}\right] \\ \frac{w_{21}}{D_{w2}}+\left[\frac{331{\Lambda }_2}{6h}-\frac{27{\Lambda }_{22}}{h^2}+\frac{19{\Lambda }_4}{h^3}-\frac{36{\Lambda }_{44}}{h^4}\right]\frac{w_{22}}{D_{w2}} \end{gathered}$$ (3.21)

$$\begin{gathered} a_3=\left[-18{\Lambda }_{2}h-6{\Lambda }_{22}-\frac{72{\Lambda }_4}{h}+\frac{72{\Lambda }_{44}}{h^2}\right]\frac{w_{20}}{D_{w2}}+\left[\frac{93{\Lambda }_2}{2}-\frac{41{\Lambda }_{22}}{2h}+\frac{96{\Lambda }_4}{h^2}-\frac{114{\Lambda }_{44}}{h^3}\right] \\ \frac{w_{21}}{D_{w2}}+\left[\frac{-81{\Lambda }_2}{h}+\frac{87{\Lambda }_{22}}{h^2}-\frac{39{\Lambda }_4}{h^3}+\frac{96{\Lambda }_{44}}{h^4}\right]\frac{w_{22}}{D_{w2}} \end{gathered}$$ (3.22)

$$\begin{gathered} a_4=\left[6{\Lambda }_{2}h-4{\Lambda }_{22}+\frac{36{\Lambda }_4}{h}-\frac{48{\Lambda }_{44}}{h^2}\right]\frac{w_{20}}{D_{w2}}+\left[-\frac{33{\Lambda }_2}{2}+\frac{23{\Lambda }_{22}}{h}-\frac{51{\Lambda }_4}{h^2}+\frac{84{\Lambda }_{44}}{h^3}\right] \\ \frac{w_{21}}{D_{w2}}+\left[+\frac{63{\Lambda }_2}{2h}-\frac{63{\Lambda }_{22}}{h^2}+\frac{21{\Lambda }_4}{h^3}-\frac{84{\Lambda }_{44}}{h^4}\right]\frac{w_{22}}{D_{w2}} \end{gathered}$$ (3.23)

$$\begin{gathered} a_5=\left[-{\Lambda }_{2}h+{\Lambda }_{22}-\frac{6{\Lambda }_4}{h}+\frac{12{\Lambda }_{44}}{h^2}\right]\frac{w_{20}}{D_{w2}}+\left[+\frac{17{\Lambda }_2}{6}-\frac{29{\Lambda }_{22}}{6h}+\frac{8{\Lambda }_4}{h^2}-\frac{22{\Lambda }_{44}}{h^3}\right] \\ \frac{w_{21}}{D_{w2}}+\left[-\frac{17}{3}\frac{{\Lambda }_2}{2h}+13\frac{{\Lambda }_{22}}{h^2}-\frac{{\Lambda }_4}{h^3}+\frac{24{\Lambda }_{44}}{h^4}\right]\frac{w_{22}}{D_{w2}} \end{gathered}$$ (3.24)

$S_2=\left[-3h{\Lambda }_2+11{\Lambda }_{22}+18\frac{{\Lambda }_4}{h}-12\frac{{\Lambda }_{44}}{h^2}\right]\frac{w_{23}\left(t\right)}{D_{w2}}$ (3.25)

Neste caso os coeficientes Λ₂ e Λ₄ foram calculados empregando a seguinte discretização para $\frac{\partial \varphi }{\partial x}$

$\frac{\partial \varphi }{\partial x}{|}_{x_i}=\frac{{\varphi }_{i+1}-{\varphi }_{i-1}}{2h}+O\left(h^2\right)$ (3.26)

3.2.3 Penúltimo nó

As equações empregadas na discretização do penúltimo nó, x _N−1 , foram obtidas de forma semelhante as do segundo nó e, por esta razão, são bastante semelhantes. Utilizando a condição de contorno definida pela Eq. (3.11) foram obtidas as seguintes expressões para as derivadas:

$$\begin{gathered} D_{e2}\frac{\partial \varphi }{\partial x}{|}_{x_{N-1}}=3h{e}_{23}\left(t\right)+\left[10{\varphi }_{N-1}-18{\varphi }_{N-2}+6{\varphi }_{N-3}-{\varphi }_{N-4}\right]h{e}_{20} \\ \left[\frac{197}{6}{\varphi }_{N-1}-\frac{93}{2}{\varphi }_{N-2}+\frac{33}{2}{\varphi }_{N-3}-\frac{17}{6}{\varphi }_{N-4}\right]e_{21}+ \\ \left[\frac{331}{6}{\varphi }_{N-1}-81{\varphi }_{N-2}+\frac{63}{2}{\varphi }_{N-3}-\frac{17}{3}{\varphi }_{N-4}\right]\frac{e_{22}}{h}+O\left(h^4\right) \end{gathered}$$ (3.27)

$$\begin{gathered} D_{e2}\frac{{\partial }^2\varphi }{\partial x^2}{|}_{x_{N-1}}=11h{e}_{23}\left(t\right)+\left[-20{\varphi }_{N-1}+6{\varphi }_{N-2}+4{\varphi }_{N-3}-{\varphi }_{N-4}\right]e_{20} \\ \left[\frac{7}{3}{\varphi }_{N-1}-\frac{41}{2}{\varphi }_{N-2}+23{\varphi }_{N-3}-\frac{29}{6}{\varphi }_{N-4}\right]\frac{e_{21}}{h}+ \\ \left[37{\varphi }_{N-1}-87{\varphi }_{N-2}+63{\varphi }_{N-3}-13{\varphi }_{N-4}\right]\frac{e_{22}}{h^2}+O\left(h^3\right) \end{gathered}$$ (3.28)

$$\begin{gathered} D_{e2}\frac{{\partial }^3\varphi }{\partial x^3}{|}_{x_{N-1}}=\frac{18h{e}_{23}\left(t\right)}{h}+\left[-60{\varphi }_{N-1}+72{\varphi }_{N-2}-36{\varphi }_{N-3}+6{\varphi }_{N-4}\right]\frac{e_{20}}{h} \\ \left[-53{\varphi }_{N-1}+96{\varphi }_{N-2}-51{\varphi }_{N-3}+8{\varphi }_{N-4}\right]\frac{e_{21}}{h^2}+ \\ \left[-19{\varphi }_{N-1}+39{\varphi }_{N-2}-21{\varphi }_{N-3}+{\varphi }_{N-4}\right]\frac{e_{22}}{h^3}+O\left(h^2\right) \end{gathered}$$ (3.29)

$$\begin{gathered} D_{e2}\frac{{\partial }^4\varphi }{\partial x^4}{|}_{x_{N-1}}=\frac{12e_{23}\left(t\right)}{h^2}+\left[-48{\varphi }_{N-1}+72{\varphi }_{N-2}-48{\varphi }_{N-3}+12{\varphi }_{N-4}\right]\frac{e_{20}}{h^2} \\ \left[-52{\varphi }_{N-1}+114{\varphi }_{N-2}-84{\varphi }_{N-3}+22{\varphi }_{N-4}\right]\frac{e_{21}}{h^3}+ \\ \left[-36{\varphi }_{N-1}+96{\varphi }_{N-2}-84{\varphi }_{N-3}+24{\varphi }_{N-4}\right]\frac{e_{22}}{h^4}+O\left(h\right) \end{gathered}$$ (3.30)

Nas equações acima $D_{e2}=10e_{20}{h}^2+25e_{21}h+35e_{22}\ne 0$. Substituindo-se estas expressões na Eq. (2.6) e rearranjando os termos segundo, obtém-se a equação discretizada para o segundo nó:

$a_{N-4}{\varphi }_{N-4}^{t+\Delta t}+a_{N-3}{\varphi }_{N-3}^{t+\Delta t}+a_{N-2}{\varphi }_{N-2}^{t+\Delta t}+a_{N-1}{\varphi }_{N-1}^{t+\Delta t}=S_{N-1}$ (3.31)

em que os coeficientes a _i são definidos por:

$$\begin{gathered} a_{N-1}=\frac{1}{\Delta t}+\left[20{\Lambda }_{22}-10{\Lambda }_{2}h-\frac{60{\Lambda }_4}{h}-\frac{48{\Lambda }_{44}}{h^2}\right]\frac{e_{20}}{D_{e2}}+\left[-\frac{197{\Lambda }_2}{6}-\frac{7{\Lambda }_{22}}{3h}-\frac{53{\Lambda }_4}{h^2}-\frac{52{\Lambda }_{44}}{h^3}\right]\frac{e_{21}}{D_{e2}}+ \\ \left[-\frac{331{\Lambda }_2}{6h}-\frac{37{\Lambda }_{22}}{h^2}-\frac{19{\Lambda }_4}{h^3}-\frac{36{\Lambda }_{44}}{h^4}\right]\frac{e_{22}}{D_{e2}} \end{gathered}$$ (3.32)

$$\begin{gathered} a_{N-2}=\left[18{\Lambda }_{2}h-6{\Lambda }_{22}+\frac{72{\Lambda }_4}{h}+\frac{72{\Lambda }_{44}}{h^2}\right]\frac{e_{20}}{D_{e2}}+\left[\frac{93{\Lambda }_2}{2}+\frac{41{\Lambda }_{22}}{2h}+\frac{96{\Lambda }_4}{h^2}+\frac{114{\Lambda }_{44}}{h^3}\right]\frac{e_{21}}{D_{e2}}+ \\ \left[\frac{81{\Lambda }_2}{h}+\frac{87{\Lambda }_{22}}{h^2}+\frac{39{\Lambda }_4}{h^3}+\frac{96{\Lambda }_{44}}{h^4}\right]\frac{e_{22}}{D_{e2}} \end{gathered}$$ (3.33)

$$\begin{gathered} a_{N-3}=\left[-6{\Lambda }_{2}h-4{\Lambda }_{22}-\frac{36{\Lambda }_4}{h}-\frac{48{\Lambda }_{44}}{h^2}\right]\frac{e_{20}}{D_{e2}}+\left[-\frac{33{\Lambda }_2}{2}-\frac{23{\Lambda }_{22}}{h}-\frac{51{\Lambda }_4}{h^2}-\frac{84{\Lambda }_{44}}{h^3}\right]\frac{e_{21}}{D_{e2}}+ \\ \left[-\frac{63{\Lambda }_2}{2h}-\frac{63{\Lambda }_{22}}{h^2}-\frac{21{\Lambda }_4}{h^3}-\frac{84{\Lambda }_{44}}{h^4}\right]\frac{e_{22}}{D_{e2}} \end{gathered}$$ (3.34)

$$\begin{gathered} a_{N-4}=\left[{\Lambda }_{2}h+{\Lambda }_{22}-\frac{6{\Lambda }_4}{h}+\frac{12{\Lambda }_{44}}{h^2}\right]\frac{e_{20}}{D_{e2}}+\left[+\frac{17{\Lambda }_2}{6}+\frac{29{\Lambda }_{22}}{6h}+\frac{8{\Lambda }_4}{h^2}+\frac{22{\Lambda }_{44}}{h^3}\right]\frac{e_{21}}{D_{e2}}+ \\ \left[\frac{17}{3}\frac{{\Lambda }_2}{h}+13\frac{{\Lambda }_{22}}{h^2}+\frac{{\Lambda }_4}{h^3}+\frac{24{\Lambda }_{44}}{h^4}\right]\frac{e_{22}}{D_{e2}} \end{gathered}$$ (3.35)

$S_{N-1}=\left[3h{\Lambda }_2+11{\Lambda }_{22}-18\frac{{\Lambda }_4}{h}-12\frac{{\Lambda }_{44}}{h^2}\right]\frac{e_{23}\left(t\right)}{D_{e2}}$ (3.36)

Novamente, os coeficientes Λ₂ e Λ₄ foram calculados empregando a Eq. (3.26).

3.2.4 Último nó

A equação para o último nó é obtida utilizando procedimento análogo ao empregado para obter a equação do primeiro nó. Neste caso, são feitas expansões em Série de Taylor em torno dos quatro pontos anteriores ao último nó, obtendo-se as seguintes expressões para as derivadas de primeira e de segunda ordem:

$\frac{\partial \varphi }{\partial x}{|}_{x_N}=\frac{25{\varphi }_N}{12h}-\frac{4{\varphi }_{N-1}}{h}+\frac{3{\varphi }_{N-2}}{h}-\frac{4{\varphi }_{N-3}}{3h}+\frac{{\varphi }_{N-4}}{4h}+O\left(h^4\right)$ (3.37)

$\frac{{\partial }^2\varphi }{\partial x^2}{|}_{x_N}=\frac{35{\varphi }_N}{12h^2}-\frac{26{\varphi }_{N-1}}{3h^2}+\frac{19{\varphi }_{N-2}}{2h^2}-\frac{14{\varphi }_{N-3}}{3h^2}+\frac{11{\varphi }_{N-4}}{12h^2}+O\left(h^3\right)$ (3.38)

Substituindo estas duas derivadas na Eq. (3.8) obtém-se a equação discretizada para o primeiro:

$a_{N-4}{\varphi }_{N-4}^{t+\Delta t}+a_{N-3}{\varphi }_{N-3}^{t+\Delta t}+a_{N-2}{\varphi }_{N-2}^{t+\Delta t}+a_{N-1}{\varphi }_{N-1}^{t+\Delta t}{a}_N{\varphi }_N^{t+\Delta t}=S_N$ (3.39)

na qual os coeficientes a _i são definidos por:

$$\begin{gathered} a_N=e_{10}+\frac{25e_{11}}{12h}+\frac{35e_{12}}{12h^2},a_{N-1}=-\frac{4e_{11}}{h}-\frac{26e_{12}}{3h^2} \\ {a}_{N-2}=\frac{3e_{11}}{h}+\frac{19e_{12}}{2h^2},a_{N-3}=-\frac{4e_{11}}{3h}-\frac{14e_{12}}{3h^2} \\ a-4=\frac{e_{11}}{4h}+\frac{11e_{12}}{12h^2},S_N=e_{13}\left(t\right) \end{gathered}$$ (3.40)

Ao incorporar as condições de contorno o sistema de equações discretizadas pode ser reescrito na forma matricial:

$\mathit{A}{\varphi }^{t+\Delta t}=\mathit{b}$ (3.41)

na qual b contém as informações das condições de contorno e do tempo anterior, coeficientes S _i . A é uma matriz esparsa de dimensão N × N e, por tanto, técnicas de armazenamento para tais matrizes foram utilizadas quando da confecção do código elaborado para resolver a equação objeto deste estudo.

Em cada etapa de discretização foi realizada a análise da consistência ⁷ das equações obtidas, sendo verificado que tais equações são consistentes com a Eq. (2.6) e de ordem O(h ² , ∆t).

4 RESULTADOS

Nesta seção é apresentada a verificação numérica da rotina computacional desenvolvida em linguagem C++ para a solução do modelo bi-fluxo pelo Método de Diferenças Finitas ¹¹ conforme discretização apresentada na seção anterior.

São tratados dois casos considerando condições inicias e de contorno distintas descritas conforme a seguir.

• Problema A: No primeiro problema foi considerada como condição inicial uma distribuição normal expressa por:

$\varphi \left(x\right)=\frac{1}{\sigma \sqrt{2\pi }}{e}^{\frac{-{\left(x-\mu \right)}^2}{2{\sigma }^2}}$ (4.1)

• em que µ é a média e σ o desvio padrão. Neste caso as condições de contorno são definidas pelo valor da concentração e da derivada de segunda ordem em ambos os contornos.

• Problema B: Para o segundo problema optou-se pela utilização de uma condição de contorno já tratada por diferentes métodos em situações com β constante ^{14), (19}, tendo sido verificados valores negativos da concentração.

$\varphi \left(x\right)={\sin }^{100}\left(\frac{x}{L}\pi \right)$ (4.2)

Para este caso, considera-se que os valores da concentração e do fluxo primário nos dois contornos sejam nulos.

O método empregado para a solução do sistema de equações lineares Eq. (3.41) foi o método do Gradiente bi-conjugado estabilizado ²². Todas as simulações foram realizadas em um computador com 4GB de memória RAM e um processador Intel ^®Core i5-5200U, 2.20GHz CPU.

Os resultados apresentados nesta seção estão divididos em duas partes. Inicialmente é realizada a análise de convergência do método para alguns casos a fim de certificar a qualidade do esquema numérico adotado e do código computacional desenvolvido. Em um segundo momento são realizados alguns testes com o intuito de investigar a ocorrência de valores negativos na solução da equação de difusão desenvolvida.

4.1 Análise de convergência

Para verificar o desempenho da discretização proposta na solução da Eq. (2.2) foram utilizados problemas testes com soluções analíticas construídas utilizando a técnica de solução manufaturada ^{17), (18}. Com este propósito inclui-se mais um termo na equação original resultando em:

$\frac{\partial \varphi }{\partial t}=\frac{\partial }{\partial x}\left({\kappa }_2\beta \frac{\partial \varphi }{\partial x}\right)-\frac{\partial }{\partial x}\left({\kappa }_4\beta \left(1-\beta \right)\frac{{\partial }^3\varphi }{\partial x^3}\right)+g\left(x,t\right)$ (4.3)

A inclusão de g(x,t) na discretização pode ser feita sem grandes alterações nas equações existentes: ao termo independente de todos os nós, com exceção do primeiro e do último, acrescenta-se o valor $g_i^{t+\Delta t}$.

A ordem de convergência do MDF foi avaliada através da seguinte expressão:

$ordem={\log }_2\left(\frac{{||\varphi -{\varphi }_{\delta }||}_j}{{||\varphi -{\varphi }_{\delta /2}||}_j}\right),j=2,\infty$ (4.4)

em que ϕ_δ é a solução numérica obtida com passo δ, ou seja, δ será igual a h quando feita a análise em relação ao refinamento espacial da malha ou igual a ∆t para a análise em relação ao refinamento do passo no tempo empregado. A ordem de convergência considerando as normas euclidiana, $j=2$, e máxima, $j=\infty$.

Nesta etapa considerou-se um caso com solução analítica definida por:

$\varphi \left(x,t\right)=\frac{1}{\sigma \sqrt{2\pi }}{e}^{\frac{-{\left(x-\mu \right)}^2}{2{\sigma }^2}}{e}^{\left(-t\right)}$ (4.5)

cujas condições de contorno e inicial são definidas conforme o problema A. A expressão para g(x,t) é obtida substituindo as Eqs. (4.5) e (2.5) a (2.3) em:

$g\left(x,t\right)=\frac{\partial \varphi }{\partial t}-{\kappa }_2\left(\frac{\partial \beta }{\partial x}\frac{\partial \varphi }{\partial x}+\beta \frac{{\partial }^2\varphi }{\partial x^2}\right)+{\kappa }_4\left(\left(1-2\beta \right)\frac{\partial \beta }{\partial x}\frac{{\partial }^3\varphi }{\partial x^3}+\beta \left(1-\beta \right)\frac{{\partial }^4\varphi }{\partial x^4}\right)$ (4.6)

A soluções foram obtidas no domínio espacial $0\le x\le L$, com <mml:math><mml:mi>L</mml:mi><mml:mo>=</mml:mo><mml:mn>4</mml:mn></mml:math> e no intervalo de tempo $0<t\le t_f$, com $t_f=2$. Os valores dos parâmetros empregados nas simulações foram ${\kappa }_2=1\times {10}^{-2},{\kappa }_4=5\times {10}^{-3},\mu =2,\sigma =0.2$.

Inicialmente apresenta-se os resultados para o caso em que β é constante, $\beta =0.8$. As tabelas a seguir apresentam a norma máxima e a norma euclidiana do erro, ε, entre a solução analítica e a obtida pelo MDF, o custo computacional, bem como a ordem de convergência da solução calculada com ambas as normas.

Verifica-se que as ordens de convergência encontradas estão de acordo com os resultados teóricos: na Tabela 1 toma-se h fixo e observa-se que a ordem de convergência no tempo é aproximadamente igual a 1, enquanto que na Tabela 2 é feita a análise com ∆t fixo e a ordem de convergência obtida é aproximadamente 2.

Tabela 1:
Análise dos erros com β = 0.8 e h = 1 × 10⁻²

∆t	||ε||∞	Ordem	||ε||2	Ordem	CPU (s)
4 × 10−2	2,12 × 10−2	−	4,24 × 10−1	−	9,09
2 × 10−2	1,09 × 10−2	0.9499	2.24 × 10−1	0.9195	14,09
1 × 10−2	5.52 × 10−3	0.9892	1.12 × 10−1	1.0050	22.95
5 × 10−3	2.88 × 10−3	0.9401	5.81 × 10−2	0.9437	26.82

Tabela 2:
Análise dos erros com β = 0.8 e ∆t = 2 × 10⁻³

h	||ε||∞	Ordem	||ε||2	Ordem	CPU (s)
2 × 10−1	2.03 × 10−1	−	1.78522	−	0.079
1 × 10−1	3.83 × 10−2	2.40866	3.43 × 10−1	2.37931	0.172
5 × 10−2	9.45 × 10−3	2.01759	8, 38 × 10−2	2.03382	0.478
2.5 × 10−2	2.35 × 10−3	2.01004	2.24 × 10−2	1.94384	4.365

É realizada agora, uma segunda análise com $\beta =\beta \left(\varphi \right)$ definido pela Eq. (2.3). Para este caso empregou-se ${\beta }_{max}=1.0,{\beta }_{min}=0.2,\gamma =2500\text{e}{\varphi }_{min}=1\times {10}^{-3}$. Novamente os erros e as ordens de convergência encontradas utilizando a Eq. (4.4) estão em acordo com a ordem esperada, como pode ser verificado nas Tabelas 3 e 4.

Tabela 3:
Análise dos erros com β(ϕ) e h = 1 × 10⁻²

∆t	||ε||∞	Ordem	||ε||2	Ordem	CPU (s)
4 × 10−2	1.91 × 10−2	−	−3.90 × 10−1	−	12,05
2 × 10−2	9.71 × 10−3	0.9805	2.03 × 10−1	0.9401	22,43
1 × 10−2	5.00 × 10−3	0.9584	1.13 × 10−1	0.8443	32,3
5 × 10−3	2.57 × 10−3	0.9585	5.81 × 10−2	0.9229	43,24

Tabela 4:
Análise dos erros com β(ϕ) e ∆t = 2 × 10⁻³

h	||ε||∞	Ordem	||ε||2	Ordem	CPU (s)
1 × 10−1	2.17 × 10−1	−	2.1754	−	0.062
1 × 10−1	4.04 × 10−2	2.4204	4.17 × 10−1	2.3811	0.114
5 × 10−2	9.98 × 10−3	2.0188	1.05 × 10−1	1.9875	0.453
2.5 × 10−2	2.50 × 10−3	1.9983	2.50 × 10−2	2.0036	6.641

4.2 Resultados numéricos

Nesta seção são apresentados os resultados numéricos obtidos considerando as condições de contorno e inciais definidas na seção anterior. Os casos analisados estão dispostos na Tabela 5. Em todos os casos o valor de β _max utilizado foi ${\beta }_{max}=1$.

Tabela 5:
Descrição dos casos analisados

Variável	Caso 1	Caso 2	Caso 3	Caso 4	Caso 5	Caso 6
Problema	A	A	A	B	B	B
κ2	10−1	10−1	10−1	10−3	10−3	10−3
κ4	0	5 × 10−2	5 × 10−2	10−5	10−5	10−5
β	0.8	0.8	β(ϕ)	0.2	β(ϕ)	β(ϕ)
βmin	−	−	0.8	−	0.2	0.2
γ	−	−	2500	−	1000	50
ϕmin	−	−	10−3	−	10−2	0.5
L	4.0	4.0	4.0	1.0	1.0	1.0
t f	1.0	1.0	1.0	10.0	10.0	10.0
h	8 × 10−3	8 × 10−3	8 × 10−3	10−2	10−2	10−2
∆t	10−3	10−3	10−3	10−2	10−2	10−2

As Figuras de 2 e 3 ilustram o comportamento da solução numérica para o problema A. Os valores dos parâmetros que definem a condição inicial, dada pela Eq. (4.5), empregados nas simulações numéricas foram $\sigma =0,03$ e $\mu =2$. Na Figura 2 observa-se o comportamento dos casos 1 e 2. No primeiro caso tem-se ${\kappa }_4=0$, ou seja, somente a difusão primária é considerada e a solução se comporta conforme esperado, uma distribuição simétrica ao longo do eixo x e com $\varphi \left(x,t\right)>0$. No segundo caso é o incluído o termo de difusão secundária e ϕ assume valores negativos. Se esta equação representar algum problema de base física então resultados negativos podem ser fisicamente inconsistente.

Figura 2:
Resultados para os casos 1 e 2.

Figura 3:
Resultados para o caso 3 e comparação entre os casos 1 a 3.

A Figura 3 apresenta os resultados para o caso 3, bem como uma comparação entre os três casos num mesmo instante de tempo. Os resultados para o caso 3 foram obtidos definindo β como a função sigmoide, Eq. (2.3). Ao utilizar esta abordagem, não se observam valores negativos para ϕ(x,t). Na comparação entre os casos, esperaria-se que os valores de ϕ(x,t) obtidos no caso 2 fossem mais próximos daqueles obtidos no caso 3.

Pode-se argumentar que talvez as contantes que geraram β(ϕ) no caso 3 não são as mais adequadas para reproduzir o caso 2. Uma possível solução para essa situação consiste em resolver o problema inverso associado a este problema, isto é, determinar γ e ϕ_min que produzam uma solução mais próxima da observada no caso 2.

Nas Figuras 4 e 5 são observados os resultados do problema B, casos 4 a 6. No caso 4, considerou-se β constante e novamente foram observados valores negativos da concentração ϕ(x,t). Já no caso 5, em que β varia de acordo com a função sigmoide, tais valores não foram observados.

Figura 4:
Resultados para os casos 4 e 5.

Figura 5:
Resultados para o caso 6 e comparação entre os casos 4 a 6.

O mesmo se repetiu para o caso 6 apresentado na Figura 5 onde foram empregados outros valores dos coeficientes da função β(ϕ). Ao comparar os resultados de ϕ(x,t) nos casos 4, 5 e 6, observa-se que as curvas dos casos 4 e 5 são bem mais próximas do que aquelas obtidas no caso 6. Isto indica que os coeficientes escolhidos, arbitrariamente, para o caso 5 são mais adequadas para representar o fenômeno definido no caso 4.

5 CONCLUSÕES

Este trabalho realizou uma análise da equação de difusão bi-fluxo proposta por Bevilacqua e coautores ³, a qual utiliza uma equação diferencial com um termo de quarta ordem para modelar problemas de difusão anômala.

Os resultados demonstram que tal formulação, em determinadas circunstâncias, apresenta resultados fisicamente inconsistentes, como valores negativos, incompatíveis com o fenômeno simulado.

Foi proposto uma nova formulação para o modelo bi-fluxo na qual o parâmetro de controle β varia de acordo com a quantidade da propriedade em difusão. A formulação matemática para a solução do problema em análise pelo Método de Diferenças Finitas foi apresentada e os resultados demonstram a convergência das soluções, bem como a possibilidade de eliminar a ocorrência de valores negativos na solução. Os resultados apresentados mostraram que a metodologia proposta funciona como esperado.

Observou-se ainda que, a escolha dos parâmetros estabelecidos para representar a variação de β , leva a diferentes efeitos na difusão anômala, sendo o estudo de tais efeitos e dos valores ótimos para os coeficientes que definem a função β(ϕ) objetos de trabalhos futuros.

AGRADECIMENTOS

Os autores agradecem o 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.

REFERÊNCIAS

• ¹
  H. Atsumi. Hydrogen bulk retention in graphite and kinetics of diffusion. Journal of Nuclear Materials, 307-311 (2002), 1466-1470. doi: 10.1016/S0022-3115(02)01069-3. URL https://linkinghub.elsevier.com/retrieve/pii/S0022311502010693
  » https://doi.org/10.1016/S0022-3115(02)01069-3» https://linkinghub.elsevier.com/retrieve/pii/S0022311502010693
• ²
  D.A. Benson, S.W. Wheatcraft & M.M. Meerschaert. Application of a fractional advection-dispersion equation. Water Resources Research, 36(6) (2000), 1403-1412. doi: 10.1029/2000WR900031. URL http://dx.doi.org/10.1029/2000WR900031
  » https://doi.org/10.1029/2000WR900031» http://dx.doi.org/10.1029/2000WR900031
• ³
  L. Bevilacqua, A.C.N.R. Galeão & 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
• ⁴
  L. Bevilacqua, A.C.N.R. Galeão & 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
  » https://doi.org/10.1590/S1678-58782011000200007» http://www.scielo.br/scielo.php?script=sci_arttext&pid=S1678-58782011000200007& lng=en&nrm=iso&tlng=en
• ⁵
  J. Crank. “The Mathematics of Diffusion”. Oxford Univ. Press, 2 ed. (1976).
• ⁶
  M. D’Angelo, E. Fontana, R. Chertcoff & M. Rosen. Retention phenomena in non-Newtonian fluids flow. Physica A: Statistical Mechanics and its Applications, 327(1-2) (2003), 44-48. doi: 10.1016/S0378-4371(03)00436-9. URL http://linkinghub.elsevier.com/retrieve/ pii/S0378437103004369
  » https://doi.org/10.1016/S0378-4371(03)00436-9» http://linkinghub.elsevier.com/retrieve/ pii/S0378437103004369
• ⁷
  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
  » https://doi.org/10.1016/B978-075066594-0/50049-7» http://linkinghub.elsevier.com/retrieve/pii/B9780750665940500497
• ⁸
  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
  » www.coc.ufrj.br/pt/documents2/doutorado/2017/2922-jiang-m-td-17/file
• ⁹
  M. Jiang, L. Bevilacqua, A.J. Silva Neto, A.C.N.R. Galeão & J. Zhu. Bi-flux theory applied to the dispersion of particles in anisotropic substratum. Applied Mathematical Modelling, 64 (2018), 121-134. doi: 10.1016/j.apm.2018.07.022. URL https://linkinghub.elsevier.com/ retrieve/pii/S0307904X1830338X
  » https://doi.org/10.1016/j.apm.2018.07.022» https://linkinghub.elsevier.com/ retrieve/pii/S0307904X1830338X
• ¹⁰
  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
  » https://doi.org/10.1088/2058-7058/18/8/33» http://stacks.iop.org/2058-7058/18/i=8/a=33?key=crossref.49d9670701f966532fc6c0805ce1fde2
• ¹¹
  R.J. LeVeque. “Finite difference methods for ordinary and partial differential equations: steady-state and time-dependent problems”. SIAM, Philadelphia (2007), 339 pp.
• ¹²
  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
  » https://doi.org/10.1533/9781845692544.195» https://linkinghub.elsevier.com/retrieve/pii/B9781845691752500040
• ¹³
  C.R. Maliska. “Transferência de calor e mecânica dos fluidos computacional”. LTC Editora (1994), 472 pp.
• ¹⁴
  G.M. Marinho, D.C. Knupp, J.F. Vasconcellos & A.J. Silva Neto. Aplicação do Método da Transformada Integral Generalizada para Solução da Equação de Difusão Anômala com Distribuição de Fluxo Bimodal. In “XVIII Encontro Nacional de Modelagem Computacional”. Salvador (2015).
• ¹⁵
  R. Metzler, J.H. Jeon, A.G. Cherstvy & E. Barkai. Anomalous diffusion models and their properties: non-stationarity, non-ergodicity, and ageing at the centenary of single particle tracking. Phys. Chem. Chem. Phys., 16(44) (2014), 24128-24164. doi: 10.1039/C4CP03465A. URL http://xlink.rsc.org/?DOI=C4CP03465A
  » https://doi.org/10.1039/C4CP03465A» http://xlink.rsc.org/?DOI=C4CP03465A
• ¹⁶
  N. Muhammad. “Hydraulic, diffusion, and retention characteristics of inorganic chemicals in bentonite”. Ph.D. thesis, University of South Florida (2004). doi: E14-SFE0000383. URL https://digital.lib.usf.edu/SFS0025076/00001
  » https://doi.org/E14-SFE0000383» https://digital.lib.usf.edu/SFS0025076/00001
• ¹⁷
  W.L. Oberkampf & T.G. Trucano. Verification and validation benchmarks. Nuclear Engineering and Design, 238(3) (2008), 716-743. doi: 10.1016/j.nucengdes.2007.02.032. URL http://dx.doi.org/10.1016/j.nucengdes.2007.02.032http://linkinghub.elsevier.com/retrieve/pii/S0029549307003548
  » https://doi.org/10.1016/j.nucengdes.2007.02.032» http://dx.doi.org/10.1016/j.nucengdes.2007.02.032http://linkinghub.elsevier.com/retrieve/pii/S0029549307003548
• ¹⁸
  C.J. Roy. Review of code and solution verification procedures for computational simulation. Journal of Computational Physics, 205(1) (2005), 131-156. doi: 10.1016/j.jcp.2004.10.036. URL http://www.sciencedirect.com/science/article/pii/S0021999104004619
  » https://doi.org/10.1016/j.jcp.2004.10.036» http://www.sciencedirect.com/science/article/pii/S0021999104004619
• ¹⁹
  L.G. Silva. “Problemas inversos em processos difusivos com retenção”. Ph.D. thesis, Universidade do Estado do Rio de Janeiro (2013).
• ²⁰
  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). URL https://tede.lncc.br/handle/tede/177
  » https://tede.lncc.br/handle/tede/177
• ²¹
  D.Y. Tzou. “Macro to microscale heat transfer”. Wiley (2014). doi: 10.1002/9781118818275. URL https://onlinelibrary.wiley.com/doi/book/10.1002/9781118818275
  » https://doi.org/10.1002/9781118818275» https://onlinelibrary.wiley.com/doi/book/10.1002/9781118818275
• ²²
  H.A. van der Vorst. Bi-CGSTAB: A Fast and smoothly converging variant of Bi-CG for the solution of nonsymmetric linear systems. SIAM Journal on Scientific and Statistical Computing, 13(2) (1992), 631-644. doi: 10.1137/0913035. URL http://epubs.siam.org/doi/10.1137/0913035
  » https://doi.org/10.1137/0913035» http://epubs.siam.org/doi/10.1137/0913035
• ²³
  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ão, 4(3) (2019). doi: 10.21575/ 25254782rmetg2019vol4n3850. URL http://periodicos.ifpr.edu.br/index.php?journal=MundiETG&page=article&op=view&path[]=850
  » https://doi.org/10.21575/ 25254782rmetg2019vol4n3850» http://periodicos.ifpr.edu.br/index.php?journal=MundiETG&page=article&op=view&path[]=850
• ²⁴
  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
  » https://doi.org/10.1016/j.rimni.2016.05.001» https://www.scipedia.com/public/Vasconcellos_et_al_2016a
• ²⁵
  L. Vlahos, H. Isliker, Y. Kominis & K. Hizanidis. Normal and Anomalous Diffusion: A Tutorial (2008).
• ²⁶
  T. Zhang & Q. Guo. The finite difference/finite volume method for solving the fractional diffusion equation. Journal of Computational Physics , 375 (2018), 120-134. doi: 10.1016/j.jcp.2018.08.033. URL https://linkinghub.elsevier.com/retrieve/pii/S0021999118305564
  » https://doi.org/10.1016/j.jcp.2018.08.033» https://linkinghub.elsevier.com/retrieve/pii/S0021999118305564

Datas de Publicação

• Publicação nesta coleção
  03 Ago 2020
• Data do Fascículo
  May-Aug 2020

Histórico

• Recebido
  10 Dez 2018
• Aceito
  19 Jan 2020