我正在尝试使用迎风方案对一维平流对流进行数值建模。我使用以下等式来计算内部单元格的值：

$$C_x^{t+1} = C_x^{t} + D\frac{\Delta t}{\Delta x^2} ( C_{x-1}^t-2 C_{x}^t+C_{x+1}^t) + u\frac{\Delta t}{\Delta x}(C_{x+1}^t-C_{x}^t)$$

我正在使用 Robin 边界条件（又名绝缘边界）。换句话说，我试图保证没有质量会离开模型环境。通过将最左边和最右边的单元格设置为：

$C_1^{t+1} = D \frac{C_2^{t}}{u+D}$和$C_n^{t+1} = C_{n-1}^t\frac{(u\Delta x+D)}{D}$

在哪里$C$是浓度，$D$扩散系数，$u$是速度和$x$是当前单元格的 x 坐标，从$1$到$n$.

我理解当满足以下两个条件时，该方案在数值上应该是稳定和准确的：

$u\frac{\Delta t}{\Delta x} + 2D\frac{\Delta t}{\Delta x^2} \leq1$和$\frac{u\Delta x}{D} < \frac{2}{1- u\frac{\Delta t}{\Delta x}}$

但是，如果我使用参数运行此模型$D=0.05$,$u=0.1$,$\Delta t=0.001$,$\Delta x=1$, 10,000 个时间步后有一些质量损失，即使满足上述条件。可能看起来很多，但由于我的时间步长非常小，这意味着更少$t<10$当这种情况开始发生时。

任何想法，什么会导致这种质量损失？

下面是在 R 中进行此模拟的代码。

transfun = function(time,state,parms) {
    D = parms$D
u = parms$u
    dt = parms$dt
dx = parms$dx
    n  = length(state)


    new = numeric(n)

    #Boundary cells
    state[1]= state[2]* D/(u+D)
    state[n]= state[n-1]* (u*dx+D)/D


    # Internal cells
    for(a  in 2:(n-1)) {
        state[a] = state[a] +D*dt/dx^2 * ( state[a-1] + -2*state[a] + state[a+1]) + 
                 u * dt/dx * (-state[a] + state[a+1])
    }

    #Boundary cells
    state[1]= state[2]* D/(u+D)
    state[n]= state[n-1]* (u*dx+D)/D


    return(list(state))
}

## Hand made loop 
# Ensure parameters follow the stability rule: 0<= C^2 <= 2s <= 1 
# where C = u*dt/dx and s = D*dt/dx^2
D=0.05
u=0.1
dt=0.001
dx=1
C = u*dt/dx
s = D*dt/(dx^2)
n = 15
T = 10000
print(C+2*s < 1)
print(u*dx/D<2/(1-C))



resul = matrix(NA,T,n)
resul[1,] = rep(0,n)
resul[1,ceiling(15/2)] = 1
for(a in seq(2,T,1)) {
    resul[a,] = transfun(1,state=resul[a-1,],parms=list(D=D,u=u,dt=dt,dx=dx))[[1]]
}

,