Conservação de uma quantidade física ao usar condições de contorno de Neumann aplicadas à equação de advecção-difusão

25

Não entendo o comportamento diferente da equação de advecção-difusão quando aplico diferentes condições de contorno. Minha motivação é a simulação de uma quantidade física real (densidade de partículas) sob difusão e advecção. A densidade de partículas deve ser conservada no interior, a menos que saia pelas bordas. Por essa lógica, se eu aplicar as condições de fronteira de Neumann, os fins do sistema, como ϕx=0 0(nos lados esquerdo e direito), o sistema deve estar"fechado",isto é, se ofluxono limite for zero, nenhuma partícula poderá escapar.

Para todas as simulações abaixo, apliquei a discretização de Crank-Nicolson à equação de difusão de advecção e todas as simulações têm condições de contorno. No entanto, para a primeira e a última linha da matriz (as linhas de condição de contorno), permito que seja alterado independentemente do valor interno. Isso permite que os pontos finais sejam totalmente implícitos.ϕx=0 0β

Abaixo, discuto 4 configurações diferentes, apenas uma delas é o que eu esperava. No final, discuto minha implementação.

Limite apenas de difusão

Aqui, os termos de advecção são desativados, definindo a velocidade como zero.

Somente difusão, com β = 0,5 (Crank-Niscolson) em todos os pontos

Somente difusão (limites de Neumann com beta = 0,5)

A quantidade não é conservada, como pode ser visto pela redução da área de pulso.

Somente difusão, com = 0,5 (Crank-Niscolson) nos pontos internos e = 1 (totalmente implícito) nos limitesβββ

Somente difusão (limites de Neumann com beta = 0,5 para interior, beta = 1 totalmente implícito) os limites

Usando uma equação totalmente implícita nos limites, alcanço o que espero: nenhuma partícula escapa . Você pode ver isso pela área sendo conservada à medida que a partícula se difunde. Por que a escolha de nos pontos de fronteira influencia a física da situação? Isso é um bug ou é esperado?β

Difusão e advecção

Quando o termo de advecção é incluído, o valor de nos limites parece não influenciar a solução. No entanto, para todos os casos em que os limites parecem estar "abertos", isto é, partículas podem escapar dos limites. Por que esse é o caso?β

Orientação e difusão com = 0,5 (Crank-Niscolson) em todos os pontosβ

Difusão de Advecção (limites de Neumann com beta = 0,5)

Orientação e difusão com = 0,5 (Crank-Niscolson) nos pontos internos e = 1 (totalmente implícito) nos limitesβββ

Advecção e Difusão (limites de Neumann com beta = 0,5 para interior, beta = 1 totalmente implícito) os limites

Implementação da equação advecção-difusão

Começando com a equação advecção-difusão,

ϕt=D2ϕx2+vϕx

Escrever usando Crank-Nicolson fornece,

ϕjn+1ϕjnΔt=D[1β(Δx)2(ϕj1n2ϕjn+ϕj+1n)+β(Δx)2(ϕj1n+12ϕjn+1+ϕj+1n+1)]+v[1β2Δx(ϕj+1nϕj1n)+β2Δx(ϕj+1n+1ϕj1n+1)]

Observe que = 0,5 para Crank-Nicolson, = 1 para totalmente implícito e = 0 para totalmente explícito.β ββββ

Para simplificar a notação, vamos fazer a substituição,

s=DΔt(Δx)2r=vΔt2Δx

e mova o valor conhecido da derivada de tempo para o lado direito,ϕjn

ϕjn+1=ϕjn+s(1β)(ϕj1n2ϕjn+ϕj+1n)+sβ(ϕj1n+12ϕjn+1+ϕj+1n+1)+r(1β)(ϕj+1nϕj1n)+rβ(ϕj+1n+1ϕj1n+1)

A consideração dos termos fornece,ϕ

β(r-s)ϕj-1n+1+(1+2sβ)ϕjn+1-β(s+r)ϕj+1n+1UMAϕn+1=(1-β)(s-r)ϕj-1n+(1-2s[1-β])ϕjn+(1-β)(s+r)ϕj+1nMϕn

que podemos escrever na forma de matriz como onde,UMAϕn+1=Mϕn

A=(1+2sββ(s+r)0β(rs)1+2sββ(s+r)β(rs)1+2sββ(s+r)0β(rs)1+2sβ)

M=(12s(1β)(1β)(s+r)0(1β)(sr)12s(1β)(1β)(s+r)(1β)(sr)12s(1β)(1β)(s+r)0(1β)(sr)12s(1β))

Aplicando condições de contorno de Neumann

NB está trabalhando com a derivação novamente, acho que percebi o erro. Eu assumi um esquema totalmente implícito ( = 1) ao escrever a diferença finita da condição de contorno. Se você assumir um esquema Crank-Niscolson aqui, a complexidade se tornará muito grande e eu não poderia resolver as equações resultantes para eliminar os nós que estão fora do domínio. No entanto, parece possível, existem duas equações com duas incógnitas, mas não consegui. Isso provavelmente explica a diferença entre o primeiro e o segundo gráfico acima. Penso que podemos concluir que apenas os gráficos com = 0,5 nos pontos de limite são válidos.βββ

Assumindo que o fluxo no lado esquerdo seja conhecido (assumindo uma forma totalmente implícita),

ϕ1n+1x=σL

Escrever isso como uma diferença centralizada dá,

ϕ1n+1xϕ2n+1ϕ0n+12Δx=σL

portanto, ϕ0n+1=ϕ2n+12ΔxσL

Observe que isso introduz um nó que está fora do domínio do problema. Este nó pode ser eliminado usando uma segunda equação. Podemos escrever o nó como,ϕ0n+1j=1

β(rs)ϕ0n+1+(1+2sβ)ϕ1n+1β(s+r)ϕ2n+1=(1β)(sr)ϕj1n+(12s[1β])ϕjn+(1β)(s+r)ϕj+1n

Substituindo o valor de encontrado na condição de limite, obtém o seguinte resultado para a linha = 1,ϕ0n+1j

(1+2sβ)ϕ1n+12sβϕ2n+1=(1β)(sr)ϕj1n+(12s[1β])ϕjn+(1β)(s+r)ϕj+1n+2β(rs)ΔxσL

A execução do mesmo procedimento para a linha final (em = ) produz,jJ

2sβϕJ1n+1+(1+2sβ)ϕJn+1=(1β)(sr)ϕJ1n+(12s(1β))ϕJn+2β(s+r)ΔxσR

Por fim, tornar as linhas de limite implícitas (configuração = 1) fornece,β

(1+2s)ϕ1n+12sϕ2n+1=ϕj1n+1ϕjn+2(rs)ΔxσL

2sϕJ1n+1+(1+2s)ϕJn+1=ϕJn+2(s+r)ΔxσR

Portanto, com condições de contorno de Neumann, podemos escrever a equação da matriz, ,Aϕn+1=Mϕn+bN

Onde,

A=(1+2s2s0β(rs)1+2sββ(s+r)β(rs)1+2sββ(s+r)02s1+2s)

M=(10 00 0(1-β)(s-r)1-2s(1-β)(1-β)(s+r)(1-β)(s-r)1-2s(1-β)(1-β)(s+r)0 00 01)

bN=(2(r-s)Δxσeu0 0...0 02(s+r)ΔxσR)T

Meu entendimento atual

  • Penso que a diferença entre o primeiro e o segundo gráficos é explicada pela observação do erro descrito acima.

  • Em relação à conservação da quantidade física. Acredito que a causa é que, como apontado aqui , a equação de advecção na forma que escrevi não permite propagação na direção reversa, de modo que a onda passa apenas mesmo com condições de contorno de fluxo zero . Minha intuição inicial em relação à conservação só se aplica quando o termo de advecção é zero (esta é a solução no gráfico 2, onde a área é conservada).

  • Mesmo com condições de limite de fluxo zero de Neumann a massa ainda pode sair do sistema, isso ocorre porque as condições de contorno corretas nesse caso são as condições de Robin nas quais o fluxo total é especificado . Além disso, a condição de Neunmann especifica que a massa não pode deixar o domínio por difusão , mas não diz nada sobre advecção. Em essência, o que ouvimos são condições de contorno fechado para difusão e condições de contorno aberto para advecção. Para obter mais informações, consulte a resposta aqui: Implementação de condições de gradiente de limite zero na equação de advecção-difusãoϕx=0 0j=Dϕx+vϕ=0 0.

Você concordaria?

boyfarrell
fonte
Parece que as condições de contorno não foram implementadas corretamente. Você poderia nos mostrar como impôs as condições de contorno?
precisa saber é o seguinte
OK. Atualizei com a implementação e acho que percebi o erro ao aplicar = 0,5 apenas nas linhas de limite. Atualizei meu "entendimento atual" na parte inferior da pergunta. Você tem algum comentário? β
boyfarrell
Então ... como é a discretização dos limites no caso dos limites de Robin? Você mostrou isso para os limites de Neumann, mas não para os limites de Robin.

Respostas:

15

Penso que um dos seus problemas é que (como você observou em seus comentários) as condições de Neumann não são as condições que você procura , no sentido de que elas não implicam a conservação de sua quantidade. Para encontrar a condição correta, reescreva seu PDE como

ϕt=x(Dϕx+vϕ)+S(x,t).

Agora, o termo que aparece entre parênteses, é o fluxo total e essa é a quantidade que você deve colocar em zero nos limites para conservar . (Adicionei por uma questão de generalidade e de seus comentários.) As condições de contorno que você deve impor são então (supondo que seu domínio espacial seja )Dϕx+vϕ=0 0ϕS(x,t)(-10,10)

Dϕx(-10)+vϕ(-10)=0 0

para o lado esquerdo e

Dϕx(10)+vϕ(10)=0 0

para o lado direito. Essas são as chamadas condições de limite de Robin (observe que a Wikipedia diz explicitamente que essas são as condições de isolamento para as equações de difusão de advecção).

Se você configurar essas condições de contorno, obterá as propriedades de conservação que estava procurando. De fato, integrando o domínio espacial, temos

ϕtdx=x(Dϕx+vϕ)dx+S(x,t)dx

Usando a integração de peças no lado direito, temos

ϕtdx=(Dϕx+vϕ)(10)-(Dϕx+vϕ)(-10)+S(x,t)dx

Agora, os dois termos centrais desaparecem graças às condições de contorno. Integrando com o tempo, obtemos

0 0Tϕtdxdt=0 0TS(x,t)dxdt

e se pudermos mudar as duas primeiras integrais ,

ϕ(x,T)dx-ϕ(x,0 0)dx=0 0TS(x,t)dx

Isso mostra que o domínio é isolado do exterior. Em particular, se , obtemos a conservação de .S=0 0ϕ

Dr_Sam
fonte
Percebo agora porque o que só funcionou quando = 0; porque isso implicaria conservação seguindo sua abordagem acima. Qual seria a conseqüência do uso dessa condição de contorno no exposto acima, a onda refletiria? Eu pensei que isso não seria possível porque não há nada na equação que me daria velocidade negativa? v
precisa saber é o seguinte
A melhor maneira de saber é provavelmente tentar! Mas se isso se comportar corretamente (e na IMO), você verá uma certa quantidade de que começa a acumular-se no lado esquerdo do domínio: o aviso empurra nessa direção, mas o limite está fechado. A acumulação para quando a difusão é suficientemente grande para equilibrá-la. Portanto, não, não deve haver onda refletida. ϕϕ
precisa saber é
@DrSam Apenas uma pergunta sobre a implementação. Entendo seu ponto de vista de como fazer a quantidade zero no lado esquerdo. Mas quando você diz "à direita, apenas um termo de fronteira", o que isso significa? Eu pensei que as condições de contorno deveriam ser Neumann ou Dirichlet (ou uma mistura de ambos)?
boyfarrell
@boyfarrell A esquerda / direita na resposta referia-se a uma derivação das condições de contorno corretas, não da maneira como é implementada (editada para maior clareza). As condições de Robin são condições clássicas, embora sejam menos conhecidas que Dirichlet e Neumann.
Dr_Sam
Então, com relação à implementação, você acha que eu deveria derivar as condições de limite de Robin para os dois limites? Além disso, se a equação tiver um termo de reação (por exemplo, faz a necessidade condição de contorno para incluir também este termo?
ϕt=x(Dϕx+vϕ)+S(x,t)
boyfarrell