Uma boa diferença finita para a equação de continuidade

22

O que seria uma boa discretização de diferenças finitas para a seguinte equação:

ρt+(ρu)=0 ?

Podemos pegar o caso 1D:

ρt+ddx(ρu)=0

Por alguma razão, todos os esquemas que encontro são para a formulação em coordenadas lagrangianas. Eu vim com esse esquema por enquanto (desconsidere o índice j ):

ρi,jn+1ρi,jnτ+1hx(ρi+1,jn+1+ρi,jn+12uxi+1/2,jnρi,jn+1+ρi1,jn+12uxi1/2n)=0

Mas parece ser realmente instável ou tem alguma condição de estabilidade horrível. É assim mesmo?

A velocidade é realmente calculada através da lei de darcy . Além disso, temos a equação de estado. O sistema completo também consiste em uma equação de energia e a equação de estado para o gás ideal. As velocidades podem ficar negativas .u=kμp

tiam
fonte
No caso 1D, o problema é essencialmente um pde hiperbólico de 1ª ordem. Você já tentou usar um esquema de diferenças finitas contra o vento de primeira ordem?
Paul
Até agora, estou executando o que escrevi na pergunta. Meu caso é realmente 2D, no entanto. Mas como essa é uma equação clássica, pensei que haveria alguma discretização clássica disponível.
tiam
Você poderia mostrar como seria um esquema a favor do vento? Eu estou familiarizado com o conceito do método do volume finito quando você o usa no termo convectivo, mas não há mais um produto derivado espacial.
tiam
O campo de velocidade é dado ou também satisfaz uma equação de evolução?
David Ketcheson
A velocidade é realmente calculada através da lei de darcy . O sistema completo também consiste em uma equação de energia e a equação de estado para o gás ideal. As velocidades podem ficar negativas. u=kμp
tiam

Respostas:

21

Você está olhando para a equação de conservação de massa:

dmdt=0

Ao considerar a evolução da massa por unidade de volume, isso se resume à equação de advecção da densidade na forma de fluxo:

ρt=(ρu)

O bom disso é que é apenas a equação de advecção de um campo escalar arbitrário (no nosso caso, isso é densidade ) e é (relativamente) fácil de resolver, desde que haja esquemas adequados de diferenciação de tempo e espaço, e inicial e condições de contorno.ρ

Ao projetar um esquema de diferenciação finita, nos preocupamos com convergência, estabilidade e precisão. Um esquema está convergindo se quando . A estabilidade dos esquemas garante que a quantidade permaneça finita quando . A precisão formal do esquema indica onde está o erro de truncamento na série de expansão de Taylor da derivada parcial. Examine um livro CFD para obter mais detalhes sobre essas propriedades fundamentais de um esquema diferencial. Δt0AtΔAΔtAtΔt0At

Agora, a abordagem mais simples é ir diretamente para a diferenciação a montante de 1ª ordem. Esse esquema é positivo-definido, conservador e computacionalmente eficiente. As duas primeiras propriedades são especialmente importantes quando modelamos a evolução de uma quantidade sempre positiva (ou seja, massa ou densidade).

Para simplificar, vejamos o caso 1-D:

ρt=(ρu)x

Agora é conveniente definir o fluxo , para que:Φ=ρu

(ρu)x=ΦxΔΦΔxΦi+1/2Φi1/2Δx

Aqui está um esquema do que estamos simulando:

            u           u
|          -->         -->          |
|    rho    |    rho    |    rho    |
x-----o-----x-----o-----x-----o-----x
     i-1  i-1/2   i   i+1/2  i+1

ρiΦi1/2Φi+1/2. É aqui que começamos a divergir da resposta de Paulo. Na verdadeira diferenciação conservadora a montante, a quantidade no centro da célula está sendo transportada pela velocidade na borda da célula, na direção de seu movimento. Em outras palavras, se você imagina que é a quantidade prevista e está sentado no centro da célula, está sendo carregado na célula à sua frente pela velocidade na borda da célula. A avaliação do fluxo na borda da célula como um produto de densidade e velocidade, tanto na borda da célula, não está correta e não conserva a quantidade prevista.

Os fluxos de entrada e saída são avaliados como:

Φi+1/2=ui+1/2+|ui+1/2|2ρi+ui+1/2|ui+1/2|2ρi+1

Φi1/2=ui1/2+|ui1/2|2ρi1+ui1/2|ui1/2|2ρi

O tratamento acima da diferenciação de fluxo assegura a definição a montante. Em outras palavras, ajusta a direção de diferenciação de acordo com o sinal de velocidade.

O critério de estabilidade de Courant-Friedrichs-Lewy (CFL), ao fazer a diferença de tempo com uma primeira ordem simples, a diferenciação de Euler para a frente é dada como:

μ=uΔtΔx1

Observe que em 2 dimensões, o critério de estabilidade da CFL é mais rigoroso:

μ=cΔtΔx12

cu2+v2

Algumas coisas a considerar. Esse esquema pode ou não ser apropriado para sua aplicação, dependendo do tipo de processo que você está simulando. Esse esquema é altamente difusivo e é apropriado para fluxos muito suaves, sem gradientes acentuados. Também é mais difusivo para etapas mais curtas. No caso 1-D, você obterá uma solução quase exata se os gradientes forem muito pequenos e se . No caso 2-D, isso não é possível e a difusão é anisotrópica.μ=1

Se o seu sistema físico considerar ondas de choque ou altos gradientes de outro tipo, você deve considerar a diferenciação upstream de ordem superior (por exemplo, terceira ou quinta ordem). Além disso, pode valer a pena examinar a família de esquemas de transporte corrigido por fluxo (Zalesak, 1979, JCP); correção anti-difusão para o esquema acima por Smolarkiewicz (1984, JCP); Família de esquemas MPDATA de Smolarkiewicz (1998, JCP).

Para diferenciar o tempo, a diferenciação de Euler para a frente de 1ª ordem pode ser satisfatória para suas necessidades. Caso contrário, procure métodos de ordem superior, como Runge-Kutta (iterativo) ou Adams-Bashforth e Adams-Moulton (multinível).

Vale a pena olhar para alguns livros didáticos de pós-graduação em CFD para um resumo dos esquemas mencionados acima e muito mais.

milancurcic
fonte
u
1
Δt=Δxmax(u)Δt
u=Cρ
uρ
uρu=Cρρt=C[(ρ)2+ρ2ρ]
13

(ddx)

ρik+1ρikΔt+ρikUikρi1kUi1kΔx=0

ΔxΔt

Paulo
fonte
ρk+1Δt
Não tenho certeza ... acho que você precisaria verificar o erro de truncamento para garantir que ele se aproxima do PDE corretamente. Você pode considerar os outros esquemas implícitos neste site: web.mit.edu/dongs/www/publications/projects/…
Paul