假粘性和色散

电热输运

假粘性和色散

考虑线性对流方程 \[ \frac{\partial u}{\partial t} + c \frac{\partial u}{\partial x} = 0, c>0 \] 数值求解该方程时，会出现两种典型的由数值格式的误差引起的非物理现象，假粘性和色散.

假粘性

考虑迎风格式 \(O(\Delta t+\Delta x)\) \[ \frac{u^{n+1}_j - u^n_j}{\Delta t} + c\frac{u^n_j - u^n_{j-1}}{\Delta x} = 0 \] 数值格式本身只是偏微分方程的低阶近似，我们可以通过泰勒展开，写出数值格式的高阶误差项，来分析数值格式本身的特点. 我们把新得到的方程称作数值格式的修正方程.

定义\(\lambda = \Delta t/\Delta x\)，可以写出迎风格式的修正方程为 \[ \frac{\partial u}{\partial t} + c\frac{\partial u}{\partial x} = \frac{c\Delta t}{2\lambda}\left(1-\lambda c\right)\frac{\partial^2 u}{\partial x^2} + O(\Delta t+\Delta x^2) \] 如果迎风格式求解该修正方程，可以得到空间的二阶精度，也就是说看起来迎风格式在近似求解线性对流方程，实际上更近似于在求解这个对流扩散方程，我们将 \[ \mu = \frac{c\Delta t}{2\lambda}\left(1-\lambda c\right) \] 定义为假粘性. 因为空间二阶导数对应着梯度项的散度项，而梯度和流量成正比的过程，比如物质扩散过程、热传导过程，都是倾向于把物理场的差异拉平. 所以采用迎风格式求解对流方程时，会出现不真实的耗散现象. 从\(u\)的表达式可以看出，随着空间网格尺寸的减小，假粘性会减弱，不真实的耗散也就随之降低.

考察下面这个简单的问题，波速\(c=1\)，初始时刻在\(x\in(0.5,1)\)的区域内\(u=2\)，其余位置\(u=1\)，观察一下采用迎风格式数值求解的结果. 和抛物方程不同，由于双曲方程的波速有限，扰动不会传播到全场，因此对于这个简单的问题我们不需要给出边界条件.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
import numpy as np
import matplotlib.pyplot as plt

%matplotlib inline
%config InlineBackend.figure_format = 'svg'
plt.style.use('science')

nx = 41
dx = 3 / (nx - 1)
nt = 60
dt = 0.025
c = 1

u = np.ones(nx)
u[int(0.5 / dx) : int(1 / dx + 1)] = 2  # setting u = 2 between 0.5 and 1.
un = np.ones(nx)

for n in range(nt):
    if n % 20 == 0:
        plt.plot(np.linspace(0, 3, nx), u, label="t={}".format(dt * n))
    un = u.copy()
    for i in range(1, nx):
        u[i] = un[i] - c * dt / dx * (un[i] - un[i - 1])
        
plt.legend()
plt.xlabel(r"$x$")
plt.ylabel(r"$u$")
plt.savefig('wangwang.svg')

可以看到随着求解过程，这个分布在假粘性作用下，越来越向周围扩散了.

减小网格尺寸，将nx设置为121，此时就不会发生耗散现象了.

色散

考虑Lax-Wendroff格式 \(O(\Delta t^2 + \Delta x^2)\) \[ u^{n+1}_j = u_j^n - \frac{c\Delta t}{2\Delta x}\left(u^n_{j+1} - u^n_{j-1}\right) + \frac{c^2\Delta t^2}{2\Delta x^2}\left(u^n_{j+1} - 2u_j^n + 2^n_{j-1}\right) \] 和迎风格式一样，我们也可以得到该格式的修正方程 \[ \frac{\partial u}{\partial t} + c\frac{\partial u}{\partial x} = \mu\frac{\partial^3 u}{\partial x^3} + O(\Delta t^3+\Delta x^3) \] 其中 \[ \mu = \frac{c\Delta x^2}{6}\left(c^2\lambda^2 - 1\right) \] 该方程称为色散方程. 加入我们将解 \[ u=\mathrm{e}^{i\left[kx-a(k)t\right]} \] 代入该修正方程，可以得到 \[ a(k) = ck + \mu k^3 \] 即不同波数\(k\)的波速是不同的，相速度（不同波数的波的传播速度）\(\neq\) 群速度（波的平均传播速度）. 此时随着时间的演化，不同波数\(k\)的运动会分开，这个波会散掉，就像白光通过棱镜一样. 同时，我们可以得到群速度 \[ c_g = \partial a/\partial k = c + 3k^2\mu \] 由于收敛的必要条件（CFL条件）是差分格式的依赖区域必须包含PDE的依赖区域，即一个时间步内不可能流过多个空间网格，即 \[ c\lambda < 1 \] 于是\(\mu < 0\)，此时\(c_g < 0\). 也就是说除了色散，该格式还会导致平均数值波速小于真实波速.

P.S.

如果将色散方程中的恒定流速\(c\)改为非线性对流，就得到了有名的KdV方程 \[ \partial_t u - 6u \partial_x u + \partial_x^3 u = 0 \]

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
nx = 41
dx = 3 / (nx - 1)
nt = 60
dt = 0.025
c = 1

u = np.ones(nx)
u[int(0.5 / dx) : int(1 / dx + 1)] = 2  # setting u = 2 between 0.5 and 1.
un = np.ones(nx)

for n in range(nt):
    if n % 20 == 0:
        plt.plot(np.linspace(0, 3, nx), u, label="t={}".format(dt * n))
    un = u.copy()
    for i in range(nx-1): ## ... uncommenting this line and see what happens!
        u[i] = (
            un[i]
            - c * dt / (2 * dx) * (un[i + 1] - un[i - 1])
            + c ** 2 * dt ** 2 / (2 * dx ** 2) * (un[i + 1] - 2 * un[i] + un[i - 1])
        )
        
plt.legend()
plt.xlabel(r"$x$")
plt.ylabel(r"$u$")

plt.savefig('miaomiao.svg')

可以看到随着求解过程，这个分布发生了色散，不同波数\(k\)的解分开了，体现为局部的振荡.

减小网格尺寸，将nx设置为121，此时就不会发生色散现象了.