在 Advection-Diffusion 方程中实现零通量条件的代码

Question

我想知道如何为 avdection-diffusion 定义的方程实现零通量条件：

$Difussion$

分析以上可知满足零磁通条件时：

$BC$ .

所以，我用有限差分方案写了一个代码：

import numpy as np
import matplotlib.pyplot as plt

nx = 101
dx = 0.01
nt = 200
c = -1.
D = .1
dt = 0.0001

u = np.zeros(nx)
u[0:50] = 2
un = u.copy()

for n in range(nt):
    for i in range(1, nx-1):
    un = u.copy()
    u[i] = un[i] + D * dt /dx**2* (un[i+1]-2*un[i]+un[i-1])-\
           c*dt/dx*(un[i] - un[i-1])

    u[0] =  u[1]*(c*dx+D)/D #BC zero flux at left side

plt.plot(np.linspace(0, 1, 101), u)
plt.ylim(0, 4)
plt.xlim(0, 1.)
plt.show()

哪里

u[1]*(c*dx+D)/D

表示零通量条件，其结果来自：

$D\frac{u^{t}_{nx+1} - u^{t}_{nx} }{\Delta t} + c u^{t}_{nx+1} =0 \ u^{t}_{nx} =u^t_{nx+1} \frac{c*\Delta x + D}{D}$

但是结果不满足零通量条件，所以质量随时间不守恒

任何人都可以帮助我发现我的错误吗？

提前致谢，

Answer 1

解决此问题的最简单方法是 explicit upwind finite volume 方案。设 u[i] 表示 u 在单元格 i 上的平均值，q[i] 表示单元格 i 外的通量，则两个方程式为

保护：du[i]/dt + (q[i]-q[i-1])/dx = 0
流量：q[i] = c*u[i] - D*(u[i+1]-u[i])/dx

在一个时间步中，您使用 u[0:nx] 的已知值来计算 interior 通量 q[0:nx]。然后使用这些计算变化 du[0:nx]，以及 flux 边界条件 q[-1]=q[nx]=0（即不需要幽灵节点）。

（请注意，对于此类问题，您可能应该在 Computational Science StackExchange 处提问。）

Answer 2

我已经找到问题的答案了，我会在这里分享，以防有人觉得有用。

我想为平流扩散方程建立zero-flux保守边界条件，这代表罗宾边界条件。鉴于在 SO 中几乎不可能编写等式，您可以在此处查看一个很好的解释：https://scicomp.stackexchange.com/questions/5434/conservation-of-a-physical-quantity-when-using-neumann-boundary-conditions-appli

下面的代码将RBC加入到平流扩散方程中，解决了我的问题。

# 1. Import libraries
import numpy as np
import matplotlib.pyplot as plt

# 2. Set up parameters
nx = 101
dx = 0.01
nt = 7000
c = .5
D = .1
dt = 0.0001

# 3. Initial conditions
u = np.zeros(nx)
u[10:35] = 4
un = u.copy()

# 4. Solving the PDE
for n in range(nt):
    un = u.copy()
    for i in range(1, nx-1):
        u[i] = un[i] + D * dt / dx ** 2 * (un[i+1]-2*un[i]+un[i-1]) -\
           c*dt/dx*(un[i] - un[i-1])
    u[0] = u[1] - dx * c / D * u[1]    # Robin Boundary Condition LEFT
    u[-1] = u[-2] + dx * c / D * u[-1] # Robin Boundary Condition RIGHT

# 5. Plot results
plt.plot(np.linspace(0, 1, 101), u)
plt.ylim(0, 4)
plt.xlim(0, 1.)
plt.show()

在 Advection-Diffusion 方程中实现零通量条件的代码

Implementing code for zero flux condition in Advection-Diffusion equation

python

differential-equations