Análise numérica da equação da difusão anômala com fluxo bimodal

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 Fick's 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.

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 [1–5]. Bevilacqua et al. [6–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:

O termo de quarta ordem da Eq. (1) foi denominado originalmente de termo de difusão de retenção [6,7]. Em trabalhos posteriores, houve uma mudança na nomenclatura para termo de difusão secundária [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. (1), 1−βk22 representa o coeficiente de difusão de Fick, e que denominaremos de difusão primária, e β1−βk44 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₂₂ e k₄₄ são constantes. Nas situações em que β = 0 ou k₄₄ = 0, a Eq. (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. [6–8] abordaram quase que exclusivamente a apresentação da teoria que levou ao desenvolvimento da Eq. (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. (1) e suas respectivas condições de contorno. Neste ponto, para melhor compreensão das hipóteses que levaram à elaboração da Eq. (1), recomenda‐se a leitura dos artigos supracitados.

O grande desafio na modelagem da Eq. (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. (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₂₂ ou k₄₄. Portanto, desconsiderou‐se o termo transiente da Eq. (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. (2):

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

Na Eq. (3), x_i representa a coordenada do centro do i‐ésimo volume de tamanho h, detalhado na figura 1.

Figura 1.

Esquema dos volumes internos.

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:

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

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

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. (6) e (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. (4) e (5), utilizou‐se as seguintes expansões na definição do fluxo secundário na face esquerda do volume:

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

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

A Eq. (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. (11) não pode ser empregada.

2.1. Equação dos volumes internos

As Eqs. (6), (7), (10) e (11) definem os fluxos em cada uma das faces dos volumes. Substituindo‐as na Eq. (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. (9):

e Bi=Bxi. 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. (10) ou (11). Observe‐se que nas Eqs. (6), (7), (10) e (11) as derivadas foram avaliadas nas interfaces e na Eq. (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. (6), (7), (10) e (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. (2) é necessário utilizar as 4 condições de contorno e desenvolver equações diferentes da Eq. (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:

Na Eq. (13), w10, w11, w12, w20, w21 e w22 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: w10=0 e w11≠0.

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. (2). A primeira forma idealizada para implementar as condições de contorno baseou‐se na utilização da Eq. (13), em conjunto com expansões em série de Taylor utilizando‐se os 2 primeiros volumes, Eqs. (5) e (9). Nessa forma, os fluxos primário e secundário da face esquerda do primeiro volume, conforme mostra a figura 2, assumem a forma:

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 figura 2. Desta forma, para esta face, os fluxos foram assim determinados:

Figura 2.

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

Substituindo as Eqs. (14) a (17) na Eq. (3), chega‐se a equação discretizada para o primeiro volume:

Este modelo não chegou a ser utilizado, pois uma análise prévia demonstrou que ao se escrever os valores de ϕ₀, ϕ₁, ϕ₂ e ϕ₃ em termos de série de Taylor e substituindo‐se estes termos na Eq. (18) para o limite em que h → 0 chega‐se a equação:

Conclui‐se desta forma que a Eq. (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. (16) e (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:

e

Os fluxos primários para este modelo ficam inalterados, Eqs. (14) e (16).

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

A análise que determinou a inconsistência da Eq. (18) foi realizada também na Eq. (22) e demonstrou que a Eq. (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:

e

Substituindo estas equações na Eq. (3), chega‐se a:

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. (22) e (27) serem consistentes, os resultados numéricos obtidos utilizando estas equações foram invariavelmente insatisfatórios para as situações em que Γ₂₂/Γ₄₄ < 1. Os erros, medidos pela Eq. (38), foram da ordem ou superior a 10⁻¹ 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 figura 3 e na figura 4, respectivamente. Para cada um destes volumes foram elaboradas equações especiais que satisfazem apenas as condições de contorno.

Figura 3.

Esquema dos volumes fictícios em x = 0.

Figura 4.

Esquema dos volumes fictícios em x = L.

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

A partir da Eq. (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:

e para o segundo volume teremos:

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. (30):

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

e para o penúltimo volume teremos:

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 ϕ_i, a Eq. (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. (6) e (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 Γ₂₂/Γ₄₄ 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:

e

Ao se descrever o fluxo primário utilizando‐se as Eqs. (33) e (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:

Resolvendo o sistema de equações lineares descritos pelas Eqs. (28), (29), (31), (32) e (35), obtém‐se a solução para o modelo completo, ou utilizando‐se as Eqs. (12), (28), (29), (31) e (32) obtém‐se a solução para o modelo simplificado.

O sistema de equações lineares foi resolvido utilizando a biblioteca PETSC [10,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 [12] para resolver o sistema de equações lineares e fatoração incompleta (ILU) [13] como pré‐condicionador.

4. Solução de referência

Os resultados apresentados a seguir foram obtidos considerando‐se Bx na Eq. (2), tal que:

Para este valor de Bx, a solução analítica da Eq. (2) é:

Na Eq. (34), L é o comprimento do domínio, e c₀ 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. (37) foi utilizada como referência para medir a qualidade da solução numérica através da seguinte equação:

Na Eq. (38), x_i e ϕ_i 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. (38).

Para utilizar a Eq. (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 wi0 e de wi1, definidos na Eq. (13), tem‐se:

e para valores quaisquer de e_i0 e de e_i1, definidos na Eq. (30),

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:

• 1.

  Caso 1, c₀ = 1 e n = 1;

• 2.

  Caso 2, c₀ = 1 e n = 10; e

• 3.

  Caso 3, c₀ = 10 e n = 1.

Como primeiro teste, investigou‐se a influência da relação Γ22/Γ44 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 figura 5 mostra o resultado obtido para a simulação do caso 1. Para a situação em que Γ22/Γ44>1, 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 Γ22/Γ44. Como esperado, a solução com 400 volumes apresentou resultados mais precisos, mais próximos dos esperados quando comparados aos gerados com 40 volumes.

Figura 5.

Influência de Γ22/Γ44 para o caso 1.

Os resultados obtidos para o caso 2 são bastante semelhantes aos apresentados na 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 figura 6, exibem um comportamento distinto daqueles da 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 figura 5, os maiores erros se localizam na região com os maiores valores de Γ₂₂ > Γ₄₄. 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.

Figura 6.

Influência de Γ22/Γ44 para o caso 3.

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

Figura 7.

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

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⁻⁶ é 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 Γ₄₄ = 0. Ao se considerar Γ₄₄ ≠ 0, mesmo nos casos em que Γ22/Γ44≫1 os efeitos de truncamento da solução numérica já são relevantes para h < 10⁻³.

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 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 figura 8 são significativamente maiores dos que os apresentados na 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.

Figura 8.

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

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 Γ22/Γ44 escolhida.

Utilizando‐se a metodologia detalhada em Roy [14], obteve‐se a ordem das aproximações numéricas para todos os resultados apresentados nas figuras 7 e 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 [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 [15].

Para confirmar se o modelo funcionava adequadamente para condições adversas, foi realizado um teste onde os valores de −10 ≤ c₀ ≤ 10 e de 10⁻³ ≤ 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:

• 1.

  Dirichlet/Dirichlet

• 2.

  Dirichlet/Neumann

• 3.

  Dirichlet/Robin

• 4.

  Neumann/Robin

Os valores dos coeficientes das condições de contorno, Eqs. (13) e (30), foram obtidos aleatoriamente dentro do intervalo −10,10.

Nestas condições, executou‐se o programa elaborado para modelar a Eq. (2) por 50.000 vezes com os valores citados gerados aleatoriamente dentro do intervalo previamente definido. Os resultados obtidos podem ser vistos na 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 figura 9, é que para diferentes valores de condições de contorno e de c₀ e de n, para diferentes razões de Γ22/Γ44, os erros obtidos estiveram sempre em níveis adequados.

Figura 9.

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

O próximo resultado, apresentado na 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₀ = 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 Γ22/Γ44>1. Mas o que mais se depreende observando os resultados da figura 10 é que os erros são bastante sensíveis ao valor de c₀, 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₀ semelhantes. Resultados semelhantes, e que não serão apresentados aqui, foram obtidos para outros valores de c₀ constantes.

Figura 10.

Resultados obtidos para condições aleatórias – malha com 400 volumes e c₀ = 10.

Os resultados apresentados na figura 9 e na 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. [6–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. (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.