Conservación de una cantidad física cuando se utilizan condiciones de contorno de Neumann aplicadas a la ecuación de advección-difusión

25

No entiendo el comportamiento diferente de la ecuación de advección-difusión cuando aplico diferentes condiciones de contorno. Mi motivación es la simulación de una cantidad física real (densidad de partículas) bajo difusión y advección. La densidad de partículas debe conservarse en el interior a menos que fluya desde los bordes. Según esta lógica, si impongo condiciones de contorno de Neumann, los extremos del sistema, como ϕx=0(en los lados izquierdo y derecho), entonces el sistema debe estar"cerrado",es decir, si elflujoen el límite es cero, entonces no pueden escapar partículas.

Para todas las simulaciones a continuación, he aplicado la discretización de Crank-Nicolson a la ecuación de advección-difusión y toda la simulación tiene ϕx=0condiciones de contorno. Sin embargo, para la primera y la última fila de la matriz (las filas de condiciones de contorno), permito queβse cambie independientemente del valor interior. Esto permite que los puntos finales sean totalmente implícitos.

A continuación, analizo 4 configuraciones diferentes, solo una de ellas es lo que esperaba. Al final hablo de mi implementación.

Límite de difusión solamente

Aquí los términos de advección se desactivan al establecer la velocidad en cero.

Difusión solamente, con = 0.5 (Crank-Niscolson) en todos los puntosβ

Difusión solamente (límites de Neumann con beta = 0.5)

La cantidad no se conserva como puede verse por la reducción del área del pulso.

Difusión solamente, con = 0.5 (Crank-Niscolson) en los puntos interiores, y = 1 (totalmente implícito) en los límitesβββ

Difusión solamente (límites de Neumann con beta = 0.5 para interior, beta = 1 totalmente implícito) los límites

Al usar una ecuación totalmente implícita en los límites, logro lo que espero: no se escapan partículas . Puede ver esto por el área que se conserva como la partícula difusa. ¿Por qué la elección de en los puntos límite influye en la física de la situación? ¿Es esto un error o esperado?β

Difusión y advección

Cuando se incluye el término de advección, el valor de en los límites no parece influir en la solución. Sin embargo, para todos los casos en que los límites parecen estar "abiertos", es decir, las partículas pueden escapar de los límites. ¿Por qué es este el caso?β

Advección y difusión con = 0.5 (Crank-Niscolson) en todos los puntosβ

Difusión-advección (límites de Neumann con beta = 0.5)

Advección y difusión con = 0.5 (Crank-Niscolson) en los puntos interiores, y = 1 (totalmente implícito) en los límitesβββ

Advección y difusión (límites de Neumann con beta = 0.5 para interior, beta = 1 totalmente implícito) los límites

Implementación de la ecuación de advección-difusión

Comenzando con la ecuación de advección-difusión,

ϕt=D2ϕx2+vϕx

Escribir usando Crank-Nicolson da,

ϕ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)]

Tenga en cuenta que = 0.5 para Crank-Nicolson, = 1 para totalmente implícito y = 0 para totalmente explícito.β ββββ

Para simplificar la notación hagamos la sustitución,

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

y mueva el valor conocido de la derivada del tiempo al lado derecho,ϕ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)

Factorizar los términos da,ϕ

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

que podemos escribir en forma de matriz como donde,Aϕ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β))

Aplicación de condiciones de contorno de Neumann

NB está trabajando a través de la derivación nuevamente, creo que he descubierto el error. Asumí un esquema totalmente implícito ( = 1) al escribir la diferencia finita de la condición de contorno. Si asume un esquema de Crank-Niscolson aquí, la complejidad se vuelve demasiado grande y no podría resolver las ecuaciones resultantes para eliminar los nodos que están fuera del dominio. Sin embargo, parece posible, hay dos ecuaciones con dos incógnitas, pero no pude manejarlo. Esto probablemente explica la diferencia entre la primera y la segunda gráfica de arriba. Creo que podemos concluir que solo las gráficas con = 0.5 en los puntos límite son válidas.βββ

Suponiendo que se conoce el flujo en el lado izquierdo (asumiendo una forma totalmente implícita),

ϕ1n+1x=σL

Escribir esto como una diferencia centrada da,

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

por lo tanto, ϕ0n+1=ϕ2n+12ΔxσL

Tenga en cuenta que esto introduce un nodo que está fuera del dominio del problema. Este nodo puede eliminarse utilizando una segunda ecuación. Podemos escribir el nodo 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

Sustituyendo en el valor de encuentra en la condición límite, se obtiene el siguiente resultado para la fila = 1,ϕ0n+1j

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

Realizando el mismo procedimiento para la fila final (en = ) produce,jJ

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

Finalmente, hacer implícitas las filas límite (configuración = 1) da,β

(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

Por lo tanto, con las condiciones de contorno de Neumann podemos escribir la ecuación matricial, ,Aϕn+1=Mϕn+bN

dónde,

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

M=(100(1β)(sr)12s(1β)(1β)(s+r)(1β)(sr)12s(1β)(1β)(s+r)001)

bN=(2(rs)ΔxσL002(s+r)ΔxσR)T

Mi comprensión actual

  • Creo que la diferencia entre la primera y la segunda gráfica se explica al observar el error descrito anteriormente.

  • En cuanto a la conservación de la cantidad física. Creo que la causa es que, como se señaló aquí , la ecuación de advección en la forma que he escrito no permite la propagación en la dirección inversa, por lo que la onda simplemente pasa incluso con condiciones de límite de flujo cero . Mi intuición inicial con respecto a la conservación solo se aplica cuando el término de advección es cero (esta es la solución en la parcela número 2 donde se conserva el área).

  • Incluso con las condiciones de límite de flujo cero de Neumann la masa aún puede abandonar el sistema, esto se debe a que las condiciones de límite correctas en este caso son condiciones de límite de Robin en las que el flujo total se especifica . Además, la condición de Neunmann especifica que la masa no puede abandonar el dominio por difusión , no dice nada acerca de la advección. En esencia, lo que hemos escuchado son condiciones de límite cerrado a la difusión y condiciones de límite abierto a la advección. Para obtener más información, consulte la respuesta aquí, Implementación de la condición de límite de gradiente cero en la ecuación de advección-difusiónϕx=0j=Dϕx+vϕ=0.

¿Estarías de acuerdo?

boyfarrell
fuente
Parece que las condiciones de contorno no se implementan correctamente. ¿Podría mostrarnos cómo ha impuesto las condiciones de contorno?
David Ketcheson
OK, actualicé con la implementación y creo que detecté el error relacionado con la aplicación de = 0.5 solo en las filas de límite. He actualizado mi "comprensión actual" al final de la pregunta. ¿Tienes algún comentario? β
boyfarrell
Entonces ... ¿cómo se ve la discretización en los límites en el caso de los límites de Robin? Lo has mostrado para los límites de Neumann, pero no para los límites de Robin.

Respuestas:

15

Creo que uno de sus problemas es que (como observó en sus comentarios) las condiciones de Neumann no son las condiciones que está buscando , en el sentido de que no implican la conservación de su cantidad. Para encontrar la condición correcta, reescriba su PDE como

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

Ahora, el término que aparece entre paréntesis, es el flujo total y esta es la cantidad que debe poner a cero en los límites para conservar . (He agregado en aras de la generalidad y por sus comentarios.) Las condiciones de contorno que debe imponer son (suponiendo que su dominio espacial sea )Dϕx+vϕ=0ϕS(x,t)(10,10)

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

para el lado izquierdo y

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

por el lado derecho Estas son las llamadas condiciones límite de Robin (tenga en cuenta que Wikipedia dice explícitamente que estas son las condiciones aislantes para las ecuaciones de advección-difusión).

Si configura estas condiciones de contorno, obtendrá las propiedades de conservación que estaba buscando. De hecho, al integrarnos en el dominio espacial, tenemos

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

Usando la integración por partes en el lado derecho, tenemos

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

Ahora, los dos términos centrales desaparecen gracias a las condiciones de contorno. Integrando en el tiempo, obtenemos

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

y si se nos permite cambiar las dos primeras integrales ,

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

Esto muestra que el dominio está aislado del exterior. En particular, si , obtenemos la conservación de .ϕS=0ϕ

Dr_Sam
fuente
Ahora me doy cuenta de por qué solo funcionaba cuando = 0; porque esto implicaría conservación siguiendo su enfoque anterior. ¿Cuál sería la consecuencia de usar esta condición límite en lo anterior, la onda se reflejaría? Pensé que esto no sería posible porque no hay nada en la ecuación que me diera velocidad negativa. v
boyfarrell
¡La mejor manera de saber es probablemente intentarlo! Pero si esto se comporta correctamente (e IMO lo hace), debería ver una cierta cantidad de que comienza a acumularse en el lado izquierdo del dominio: la advección empuja en esa dirección pero el límite está cerrado. La acumulación se detiene cuando la difusión es lo suficientemente grande como para equilibrarla. Entonces no, no debería haber onda reflejada. ϕϕϕ
Dr_Sam
@DrSam Solo una pregunta sobre la implementación. Entiendo su punto sobre cómo hacer que la cantidad sea cero en el lado izquierdo. Pero cuando dices "a la derecha solo un término límite", ¿qué significa eso? ¿Pensé que las condiciones de contorno deberían ser Neumann o Dirichlet (o una combinación de ambos)?
boyfarrell
@boyfarrell La izquierda / derecha en la respuesta se refería a una derivación de las condiciones de contorno correctas, no a la forma en que se implementa (editado para mayor claridad). Las condiciones de Robin son condiciones clásicas, aunque hay menos conocidas que Dirichlet y Neumann.
Dr_Sam
Entonces, con respecto a la implementación, ¿cree que debería derivar las condiciones de límite de Robin para ambos límites? Además, si la ecuación tiene un término de reacción (por ejemplo, ¿la condición de frontera debe incluir también este término?
ϕt=x(Dϕx+vϕ)+S(x,t)
boyfarrell