Problema de Advecção-Difusão

Prévia do material em texto

<p>Problema de Advecção-Difusão Estacionário</p><p>- Versão Preliminar</p><p>Elson Magalhães Toledo</p><p>Luis Paulo S. Barra</p><p>Outubro 2020</p><p>1 Introdução</p><p>Examinamos aqui o problema da solução numérica do problema (unidimensio-</p><p>nal) que inclui apenas os termos da derivada segunda e da derivada primeira.</p><p>Por sua origem este problema é denominado de Advecção-Difusão. Ao termo da</p><p>derivada segunda usualmente denominamos de termo difusivo por sua conexão</p><p>com a modelagem do fenômeno da difusão ou dispersão ou "espalhamento"de</p><p>uma substância em outra(uma sólida e uma líquida ou duas líquidas). O termo</p><p>da derivada primeira está associado a modelagem do fenômeno de transporte,</p><p>de uma substância ou até mesmo do transporte de uma concentração de uma</p><p>substância diluida num fluido pelo movimento(escoamento) deste fluido.</p><p>Como veremos a solução numérica deste problema por diferenças finitas,</p><p>quando se empregam diferenças de segunda ordem podem conduzir a resul-</p><p>tados cuja natureza diferem substancialmente da solução exata do problema</p><p>quando a relação que envolve os coeficientes dos termos da equação e o valor</p><p>de h de uma discretização ultrapassa certo limite. Discutiremos meios de</p><p>evitar esta impropriedade, justificando a aproximação denominada de Upwind.</p><p>Dependendo do método de discretização empregado podemos ou não ter pro-</p><p>blemas da mesma natureza quando tratamos do caso que inclui um termo</p><p>com a função e não com sua derivada. A este problema assim caracterizado</p><p>denominamos de Difusão-Reação.</p><p>1</p><p>2</p><p>2 Advecção-Difusão - Solução Analítica</p><p>Considere o PVC:</p><p>−εu′′(x) + βu′(x) = 0 0 ≤ x ≤ 1 (1)</p><p>u(0) = 0 u(1) = 1 ε, β > 0</p><p>cuja solução analítica pode ser fácilmente determinada considerando-se a</p><p>função u(x) = eλx e suas duas primeiras derivadas:</p><p>u′(x) = λeλx</p><p>u′′(x) = λ2eλx</p><p>A substituição na equação diferencial (1) fornece:</p><p>−ελ2eλx + βλeλx = 0</p><p>Tendo em vista que eλx 6= 0, obtemos:</p><p>− ελ2 + βλ = 0→</p><p>{</p><p>λ1 = 0</p><p>λ2 =</p><p>β</p><p>ε</p><p>(2)</p><p>A partir destes valores, pode se obter a solução u(x) = C1e</p><p>λ1x + C2e</p><p>λ2x fica:</p><p>u(x) = C1 + C2e</p><p>β</p><p>ε</p><p>x</p><p>Para satisfação das condições de contorno temos que a primeira desta condição</p><p>de contorno fornece:</p><p>u(0) = 0→ C1 + C2e</p><p>0 = 0</p><p>Logo:</p><p>C1 = −C2 (3)</p><p>A segunda condição se contorno nos fornece:</p><p>u(1) = 1→ C1 + C2e</p><p>β</p><p>ε = 1</p><p>3</p><p>e usando (3), chega-se a :</p><p>C2</p><p>(</p><p>e</p><p>β</p><p>ε − 1</p><p>)</p><p>= 1</p><p>Desta forma:</p><p>C2 =</p><p>(</p><p>1− e</p><p>β</p><p>ε</p><p>)−1</p><p>E a solução pode ser escrita como:</p><p>u(x) = C2</p><p>(</p><p>e</p><p>β</p><p>ε</p><p>x − 1</p><p>)</p><p>(4)</p><p>Logo:</p><p>u(x) =</p><p>e</p><p>β</p><p>ε</p><p>x − 1</p><p>e</p><p>β</p><p>ε − 1</p><p>(5)</p><p>Se definimos o número de Peclet Global do problema como:</p><p>Pe =</p><p>βL</p><p>2ε</p><p>com L = 1→ Pe =</p><p>β</p><p>2ε</p><p>Logo:</p><p>u(x) =</p><p>ePex − 1</p><p>ePe − 1</p><p>(6)</p><p>que pode ser escrita como:</p><p>u(x) = C1 + C2e</p><p>2Pex (7)</p><p>3 Valores Extremos de Pe</p><p>Examinamos em seguida a influência dos valores do número de Peclet global</p><p>na solução analitica do problema.</p><p>Pe = 0</p><p>Neste caso a equação diferencial se reduz a:</p><p>−εu′′(x) = 0 com u(0) = 0 e u(1) = 1</p><p>Logo:</p><p>u(x) = x</p><p>4</p><p>Pe� 1</p><p>Neste caso tem-se:</p><p>u(x) ≈ e2Pex</p><p>e2Pe</p><p>= e−2Pe(1−x)</p><p>Solução analítica desenvolvendo uma camada limite para valores altos número</p><p>de Peclet.</p><p>5</p><p>Figura 1: Soluções Analíticas</p><p>6</p><p>3.1 Solução Aproximada com Diferenças Finitas Centrais</p><p>Utilizando aproximações de diferenças finitas centradas de segunda ordem</p><p>para a primeira e segunda derivadas da equação (1) obtem-se para a aplicação</p><p>ao ponto i :</p><p>−ε</p><p>(</p><p>ui+1 − 2ui + ui−1</p><p>h2</p><p>)</p><p>+ β</p><p>(</p><p>ui+1 − ui−1</p><p>2h</p><p>)</p><p>= 0 i = 0, 1, 2, · · ·M</p><p>Definindo-se o número de Peclet da discretização:</p><p>Pe =</p><p>βh</p><p>2ε</p><p>onde h = L/M , onde L é o comprimento do dominio, que neste caso, adotamos</p><p>L = 1.</p><p>A equação (8) pode ser escrita como:</p><p>−</p><p>(</p><p>ui+1 − 2ui + ui−1</p><p>h2</p><p>)</p><p>+ Pe</p><p>(</p><p>ui+1 − ui−1</p><p>h2</p><p>)</p><p>= 0</p><p>ou ainda, de forma mais compacta como:</p><p>(Pe− 1)ui+1 + 2ui − (Pe+ 1)ui−1 = 0 (8)</p><p>3.2 Oscilações da Solução Numérica</p><p>Nas figuras a seguir apresentamos as soluções numéricas do sistema de equa-</p><p>ções formado pela aplicação da equação (8) variando-se o número de pontos da</p><p>malha, para valores de Pe = 5 e Pe = 50. Observe que no primeiro caso todos os</p><p>casos analisados possuem Pe 1 → 1 + Pe</p><p>1− Pe</p><p>ε =</p><p>{</p><p>1 + Pe∗</p><p>1− Pe∗</p><p>}</p><p>o que nos leva a:</p><p>Pe</p><p>∗ = tanhPe</p><p>onde Pe é o Numero de Peclet da malha para solução do problema original. A</p><p>conclusão é que um problema com um dado Peclet com uma certa discretização</p><p>h pode ser substituido por outro com um Peclet dado pela expressão acima e a</p><p>solução numérica do problema original ( com Pe > 1 ) será nodalmente exata</p><p>quando alteramos ( ou incluimos uma difusão adicional) de forma que este</p><p>novo problema numérico possua Pe∗ dado por esta última expressão.</p><p>5 Estratégia implementação Difusividade Ótima</p><p>Podemos reescrever a equação modificada(com uma difusividade adicional)</p><p>como:</p><p>−(ε+ ςεh)u</p><p>′′ + βu′ = 0</p><p>12</p><p>Operacionalmente podemos calcular um fator ς que multiplique o termo com ε</p><p>em nossa equação original de forma que tenhamos:</p><p>εh =</p><p>βh</p><p>2</p><p>O valor de Peclet para este problema deve valer:</p><p>Pe∗ = tanhPe</p><p>onde:</p><p>ε∗ = ε+ ςεh</p><p>Logo devemos ter:</p><p>Pe∗ =</p><p>βh</p><p>2(ε+ ςεh)</p><p>=</p><p>sinhPe</p><p>coshPe</p><p>De onde podemos obter ς dado por:</p><p>ς = cothPe− 1</p><p>Pe</p><p>Com esta variável podemos numa implementação computacional considerar</p><p>tres possibilidades a saber:</p><p>Se</p><p>ς = 0</p><p>ς = 1</p><p>ς = cothPe− 1</p><p>Pe</p><p>Teriamos no primeiro caso a solução original- oscilante se Pe > 1, para o</p><p>segundo caso a solução com aproximação upwind, que é sobredifusiva quando</p><p>Pe cresce muito, e no terceiro caso teriamos uma solução nodalmente exata.</p><p>6 Observações Finais</p><p>Para o caso em que o termo forçante - termo do lado direito da equação diferen-</p><p>cial - não for nulo perde-se esta característica de solução nodalmente exata</p><p>quando utilizamos a difusividade ótima como apresentada. Esta característica</p><p>13</p><p>fica mantida apenas para a parte homogênea da solução da equação diferencial</p><p>e não para a solução total da equação. Outro aspecto importante diz respeito</p><p>ao cálculo de ς que, por depender de funções hiperbólicas, requer cuidados</p><p>especiais na sua avaliação.</p><p>Introdução</p><p>Advecção-Difusão - Solução Analítica</p><p>Valores Extremos de Pe</p><p>Solução Aproximada com Diferenças Finitas Centrais</p><p>Oscilações da Solução Numérica</p><p>Solução Analítica da Aproximação por Diferenças Finitas</p><p>Diferenças Finitas - Upwind</p><p>Soluções Nodalmente Exatas: Difusividade Ótima</p><p>Estratégia implementação Difusividade Ótima</p><p>Observações Finais</p>