1. 平流方案

1.1 什么是平流输送

空气污染物作为大气中微量物质参与到整个大气的时空演变中。描述大气演变的基本方程组包含动量方程（牛顿第二定律）、连续方程（质量守恒）、热力学方程（能量守恒）、状态方程（理想气体定律）、以及水汽守恒方程。通过对大气基本动力方程组的求解，可以得到一系列描述大气状态的属性值，如风、温、压、湿、空气密度等。

其中连续方程描述了大气流场带动物质运动的过程。与空气密度一样，空气污染物作为大气组成成分，同样受到连续方程的支配，连续性方程的偏微分形式如下，

$$\frac{\partial c}{\partial t} = -\nabla \cdot (\mathbf{v} c) \tag {1-1}$$

其中 $c$ 表示污染物的浓度（质量浓度，如 $\mu g/m^3$），$\mathbf{v}$ 为速度向量场，$\mathbf{v} c$ 为通量。

可以通过数值分析方法求解连续方程的数值解（大气运动方程组没有解析解）。

为了推导数值解的方程形式，先把散度符号展开。

$$\frac{\partial c}{\partial t} = - (\frac{\partial (u c)}{\partial x} + \frac{\partial (v c)}{\partial y} + \frac{\partial (w c)}{\partial z}) \tag {1-2}$$

为了简化推导过程，只考虑一维连续方程。

$$\frac{\partial c}{\partial t} = - \frac{\partial (u c)}{\partial x} \tag {1-3}$$

有两个需要离散化处理项，分别为 ${\partial x}$ 和 ${\partial t}$ 。先进行空间离散化处理，假定空间离散式均匀的（既等网格）。

$$\frac{\partial c_{i}}{\partial t} = \frac{f_{i-1/2}-f_{i+1/2}}{\Delta x} \tag {1-4}$$

其中 $f_{i-1/2}$ 为格子 $i-1$ 和格子 $i$ 的中间 $i-1/2$ 处的通量（单位时间单位面积通过的物质量的：$\mu g/m^2s$）。

$$f_{i-1/2} = u_{x_{i-1/2}} \times c_{x_{i-1/2}} \tag {1-5}$$

将公式 1.5 带入到公式 1.4 中，可以得到

$$\frac{\partial c_{i}}{\partial t} = \frac {u_{x_{i-1/2}} \times c_{x_{i-1/2}} - u_{x_{i+1/2}} \times c_{x_{i+1/2}}} {\Delta x} \tag {1-6}$$

再对**公式 1.6 ** 时间进行离散化，采用前向欧拉方法，即可得到

$$\frac{c_{x_{i}}^{n+1}-c_{x_{i}}^{n}}{\Delta t} = \frac {u_{x_{i-1/2}}^n \times c_{x_{i-1/2}}^n - u_{x_{i+1/2}}^n \times c_{x_{i+1/2}}^n} {\Delta x} \tag {1-7}$$

其中 $n$ 为时间维度的离散索引。整理一下可得（Godunov's scheme）。

$$c_{x_{i}}^{n+1} = c_{x_{i}}^{n} + (\frac {u_{x_{i-1/2}}^n \times c_{x_{i-1/2}}^n - u_{x_{i+1/2}}^n \times c_{x_{i+1/2}}^n} {\Delta x}) \times {\Delta t} \tag {1-8}$$

从公式 1-8 式分析可知，要求解下一个时刻的浓度（等式左边），需要

1. 当前时刻该网格浓度（$c_{x_{i}}^{n}$）
2. 当前时刻该网格边界处的风速（$u_{x_{i-1/2}}^n$）
3. 当前时刻该网格边界处的浓度（$c_{x_{i-1/2}}^n$）

其中第一项是已知的（在初始时刻，给定初始浓度场，不停迭代）。当前时刻该网格边界处的风速可以由气象数值模式给出。对于网格边界处的浓度，不同的数值方案给出的解法不一样（已知网格中心浓度，如何得到网格边界处浓度）。最简单的方案叫做迎风方案，边界处浓度根据风向来确定。

• 当 $u_{x_{i-1/2}}^n > 0$ （流向为正方向）： $c_{x_{i-1/2}}^n = c_{x_{i-1}}^n$

• 当 $u_{x_{i-1/2}}^n < 0$ （流向为负方向）： $c_{x_{i-1/2}}^n = c_{x_{i}}^n$

等等，为什么能够直接用气象数值模式给出的风场？空气污染物作为大气的一部分，其浓度的时空演变是否对大气运动本身产生影响（离线模拟和在线模拟问题）。

空气污染物浓度演变是否会对大气运动本身产生影响，答案是肯定的。

首先，空气污染物浓度的变化会直接影响空气密度（从而影响下个时次的风场、气压场的求解），但空气污染物浓度是微量的，其浓度的变化对空气密度往往可以忽略不计（PPM量级）。

其次，在能量守恒方程中，污染物化学反应所产生的热量可以作为一个能量来源，该贡献在对流层中基本可以忽略不计。

最后，空气污染物与太阳以及地球辐射存在相互作用，在全球尺度的长期气候研究中，这种污染物演变对其影响显然是不可忽略（耦合模式），但在区域尺度的空气质量研究中或是短期空气污染浓度预测中，则是可以近似忽略，一些对辐射影响较大的空气污染物在气象数值模式中通常会作为常量处理（如气溶胶、臭氧等，可能会考虑月变化），而忽略短期的波动。

气象场演变与空气污染物浓度演变可以解耦的，被称为稀释假设（dilution hypothesis）。

在区域空气质量模式中（CMAQ、CAMx 等），通常都是假定稀释假设成立，通过气象数值模式预处理得到风速、温度、密度、相对湿度等变量，随后将这些气象场提供给污染物的扩散方程。这种方法通常被称为离线耦合（off-line coupling）。

区域气象数值模式（如 WRF ），常采用 Arakawa Staggered C-Grid 计算风场， 因此 $u_{x_{i-1/2}}$ 和 $u_{x_{i+1/2}}$ 为已知量。

再等等，区域气象数值模式给出的风场是中尺度的，污染物网格浓度也是代表整个网格的平均浓度值。因此，目前离散化后的连续性方程并不闭合（公式1.8）。

通过雷诺分解 (Reynolds Decomposition)，将变量分解为平均值与扰动量：

$$c = \bar{c} + c'$$

$$u = \bar{u} + u'$$

其中：

• $\bar{c}$ 和 $\bar{u}$ 分别是污染物网格浓度和区域气象数值模式给出的速度，
• $c'$ 和 $u'$ 是对应的扰动量，该项被视为不成规律或具有随机性，但是对于给定的分辨率平均风场，也不完全成立，比如 1 度的风场没办法描述中尺度的波动（这部分应该被认为是有规律的）。

将分解代入连续方程，并对其进行时间平均操作，忽略二阶小量后可得：

$$\frac{\partial \bar{c}}{\partial t} = - \frac{\partial (\bar{c} \bar{u})}{\partial x} - \frac{\partial (\overline{c' u'})}{\partial x} \tag {1-9}$$

其中公式 1-9 右边第一项就是平流输送的数学表示，而公式 1-8 的右边项则是平流输送的离散表示。

那么平流输送是什么？

平流输送是指在大尺度的风场作用下，引起的该尺度下浓度网格平均变化。

对于该尺度无法描述的由于流场引起的浓度再分配则为扰动项，在区域尺度的语境下，通常该项为湍流扩散项

根据公式 1-8 给出的数值求解方案，似乎很容易求解出平流对浓度的影响。如果仔细想想，我们会面对三个额外的问题。

第一个问题是解稳定性问题，观察公式 1-8 等号右边的中间项（把括号展开, 并稍加整理）

$$(\frac{u_{x_{i-1/2}}^n \times {\Delta t}}{\Delta x}) \times c_{x_{i-1/2}}^n \tag {1-10}$$

我们需要意识到，公式 1-8 中存在一个假设，该公式计算的浓度增量是从格子 $i-1$ 输送到格子 $i$ ， 注意这个增量会最终叠加在格子 $i$ 中，公式1-8 中只涉及对 $c_{x_{i}}$ 的更新， 那么有没有一种可能。风速太大了，吹出了格子 $i$ 呢！

当然有可能，所以我们要对数值求解方案施加限制。

$$\frac{u_{x_{i-1/2}}^n \times {\Delta t}}{\Delta x}<1 \tag {1-11}$$

这就是 Courant-Friedrichs-Lewy（CFL） 条件。为了积分的稳定，在设置${\Delta x}$ 和${\Delta t}$ 的时候，应该尽量保证 CFL 条件。

第二个问题是质量守恒问题，平流只是把物质从一个地方输送到另一个地方，不能因为平流过程的计算增加或者减少了质量，

质量不守恒一般是怎么引发的呢？

其中一个非常重要的因数就是网格的大小不一样大（体积、面积、长度）。 在实际的数值模式中，网格不一致的原因有很多种，包括 1）非均一网格；2）垂直层不为高度坐标；3）地图投影的形变。

那么为什么网格不一致，会导致质量不守恒呢，仔细观察公式 1.8 就可以发现， 物质从左边的格子输送到右边，是以浓度增量的方式更新，而浓度是强度量，需要小心处理。

如果左边的格子远小于右边，比如左边为 $1 m^3$ ，右边为 $10 m^3$ ，如果左边网格损失 $1 \mu g/m^3$ ，

物质的量损失了 $1 \mu g$ 。按照公式，右边的网格增加 $1 \mu g/m^3$ ，但是物种的量却增加了 $10 \mu g$ 。

最后要关注的问题是数值扩散，采用迎风方案，很容易使得平流方案计算的扩散速度大于物理上的实际扩散速度。

想像一下，如果一个格子很大，比如 $81km \times 81km$ ，存在三个格子，左边格子浓度很高（比如 $200 \mu g$ ），中间格子比较低（比如 $50 \mu g$ ），右边格子最低（比如 $10 \mu g$ ），风从左边吹向右边，那么会有一个高浓度气团（来至左边格子）进入中间的格子，由于浓度是直接叠加上去的（注意公式 1.8）。这个高浓度值会瞬间在整个格子里面混合，也就是说这个增量一下就会吹到中间格子的东边界（接下来，里面就会对右边格子产生影响），这肯定比实际扩散的要快。

迎风方案会导致这么明显的数值扩散原因是该方案假定了边界处的浓度与格子中心浓度相等（和瞬间充分的条件等价）。

而一个好的平流方案，会精心处理边界上的浓度值，从而大大缓解数值扩散的影响。

1.2 PPM

分段多项式平流方案（Piecewise Parabolic Method，PPM）是一种用于求解流体动力学方程的数值方法，特别是在处理高分辨率、具有激波等间断的流体问题时表现出色。

PPM具有三阶精度，广泛用于主流的空气质量模式中（CMAQ、CAMx 等）。该方案适用于非结构化网格（结构化网格只是非结构化网格的特例），为了思路简单清楚，在本文档的推导过程中，采用等网格间距的设定，且只考虑一维的情况。

浓度的离散代表网格平均浓度。根据定义，可以得到以下公式

$$c_i^n = \frac{1}{\Delta x} \int_{x_{i-1/2}}^{x_{i+1/2}} c(x, t^n) dx \tag {1-12}$$

其中，$c_i^n$ 表示第 $i$ 个格子在 $t^n$ 时刻的浓度。$c(x, t^n)$ 为 $t^n$ 时刻浓度在空间（当前格子）的分布函数，当然 $c(x, t^n)$ 函数是未知的，也是无法精确给出的。

可以针对某个时刻

$$c_i^n = \frac{1}{\Delta x} \int_{x_{i-1/2}}^{x_{i+1/2}} c^n(x) dx \tag {1-13}$$

注意 $c(x, t^n)$ 变为 $c^n(x)$ 暗含了一层意思，不同时间的浓度空间分布可以用不同的函数表示 ($c^n$) 。

如果知道每个时刻的 $c^n$ ，那么就可以得到边界处的浓度，不同方案对 $c^n$ 的函数形式给出了不同的假设。

比如零阶精度假定格子内的浓度均匀分布（分布函数为常数，格子边界处浓度等于平均浓度），一阶精度假定函数形式为线性的（两个格子的线性插值），而 PPM 假定该函数为多项式。

PPM 方案采用分段连续的二次多项式方案来表示每个格子在不同时间的 $c^n$。为了方便，在后续过程中，省略时间的角标，也就是只考虑某个时次的情况，在每个网格中的二次多项式为（假定左右边界的浓度值已知，该值在后续的推导过程中确认）

$$\begin{align*} c(x) &= c_{x_{i-1/2}} + \xi(\Delta c_i + (1-\xi)c_{6,i}) \quad 0 \le \xi \le 1 \\ c(x) &= c_{x_{i-1/2}} + \frac{x-x_{i-1/2}}{\Delta x}(\Delta c_i + \frac{x_{i+1/2}-x}{\Delta x}{c_{6,i}}), \quad x_{i-1/2} \le x \le x_{i+1/2} \tag {1-14} \end{align*}$$

其中

$$\begin{align*} \xi &= \frac{x-x_{i-1/2}}{\Delta x} \\ \Delta c_i &= c_{x_{i+1/2}} - c_{x_{i-1/2}} \\ \end{align*}$$

已知左右边界的浓度值，已知函数形式为二次多项式，如何得到该函数的参数?

从公式 1-14 中很容易看出，构造的多项式在左右边界是连续（带入$x=x_{i-1/2}$，或者带入 $x=x_{i+1/2}$）。

公式 1-14 中只有一个未知参数 $c_{6,i}$。如何计算 $c_{6,i}$，

考虑其他已知条件，我们知道该网格的平均浓度 $c_{i}$。

对公式 1-14 两边进行积分求，并带入公式 1-12 ，可到 $c_{6,i}$ 的表达式

$$c_{6,i} = 6(c_i- \frac{1}{2}(c_{x_{i+1/2}} + c_{x_{i-1/2}}) ) \tag {1-16}$$

可以采用sympy 进行公式推导，并且进行验证。

from sympy import *
init_printing()

ci = symbols('c_i', real=true, constant=true)  # 网格平均浓度
cL = symbols('cL_i', real=true, constant=true) # 左边界浓度
cR = symbols('cR_i', real=true, constant=true) # 右边界浓度

x = symbols('x')
xL = symbols('xL_i', real=true, constant=true) # 左边界位置
xR = symbols('xR_i', real=true, constant=true) # 右边界位置

c6 = symbols('c6_i', real=true, constant=true) # 参数c6

f = cL + (x-xL)/(xR-xL) * (cR-cL + c6*((xR-x)/(xR-xL))) # 公式1-14
f R = integrate(f, (x, xL, xR))/(xR-xL) # 公式2.2的左边, factor(fR)

solveset(Eq(fR, ci), c6) # 公式2.2

在推导公式 1-14 的过程中，假定了网格之间边界处的浓度已知，但目前并不知道该值。那么，该如何获得 $c_{x_{i+1/2}}$ ？

在PPM中，通过四次多项式（quartic polynomial）描述浓度在空间的积分函数，从而来求解 $c_{x_{i+1/2}}$ 的分布。

为什么是四次多项式，而不是三次多项式呢？

因为 $c_x$ 在 $x$ 上的分布是分段抛物线，抛物线的积分是三次，而分段抛物线的积分是四次。

首先定义浓度在空间的积分函数：

$$C(x_{i+1/2}) = \int^{x_{i+1/2}} c(x)dx = C_{i+1/2}= \sum_{k \le i} c_i \Delta x \tag {1-16}$$

而网格边界处的浓度 $c_{x_{i+1/2}}=\frac{dC}{dx}|x=x_{i+1/2}$

从公式 1-16 可知 $C(x)$ 在边界处的值是已知，现在要求边界处的浓度值，我们需要知道 $C(x)$ 的函数形式，在PPM中，用5个点（$C_{i+1/2}, x_{i+1/2}$）来构建一个四次多项式（5个未知数，5个方程）。

from sympy import *
init_printing() # 更好的打印
# 四次多项式
x, y = symbols('x y', real=true)
a1, a2, a3, a4, a5 = symbols('a1:6') # 参数
f = a1*x**4 + a2*x**3 + a3*x**2  + a4*x + a5 - y
# 线性方程: 已知5个点，确定多项式的参数
dx = symbols('dx') #
x3 = symbols('x_3') # 中心边界点的位置
y3 = symbols('y_3') # 中心边界点处C的值
c1, c2, c3, c4 = symbols('c1:5') # 网格平均浓度
# 方程组
f1 = f.subs([(x, x3 - 2*dx), (y, y3 - c1*dx - c2*dx)])
f2 = f.subs([(x, x3 -   dx), (y, y3         - c2*dx)])
f3 = f.subs([(x, x3       ), (y, y3                )])
f4 = f.subs([(x, x3 +   dx), (y, y3         + c3*dx)])
f5 = f.subs([(x, x3 + 2*dx), (y, y3 - c4*dx + c3*dx)])
# 解线性方程组
args = linsolve([f1, f2, f3, f4, f5], [a1, a2, a3, a4, a5])

# 带入多项式
a1, a2, a3, a4, a5 = tuple(*args)
f = a1*x**4 + a2*x**3 + a3*x**2  + a4*x + a5 - y

df = Derivative(f, x).doit()
simplify(df.subs(x, x3))

通过代数求解，可以获得

$$c_{x_{i+1/2}}=\frac{dC}{dx}|x=x_{i+1/2} = \frac{7}{12}(c_i + c_{i+1}) - \frac{1}{12}(c_{i-1} + c_{i+2}) \tag {1-17}$$

通过该公式，基于网格平均浓度即可计算边界处浓度的值。理论上基于以上公式，就可以进行通量计算了。

对于公式 1-17 ，需要稍加整理。

$$\begin{align*} c_{x_{i+1/2}} &= 0.5(c_{i+1}+c_i) + \frac{1}{6}(\delta c_{i} - \delta c_{i+1}) \\ &= c_i + 0.5(c_{i+1}-c_i) + \frac{1}{6}(\delta c_{i} - \delta c_{i+1}) \\ &\simeq c_i + 0.5(c_{i+1}-c_i) + \frac{1}{6}(\delta_m c_{i} - \delta_m c_{i+1}) \tag {1-18} \end{align*}$$

其中 $\delta c_{i} = 0.5(c_{i+1} - c_{i-1})$ 为第 $i$ 个格子的平均斜率（三个网格连续变化的斜率？）。

为了让浓度不连续的地方模拟得更好（steeper representation），同时保证了边界浓度在两个格子平均浓度之间（条件1，不人为创造极值，比如可能出现负数），用 $\delta_m c_{i}$ 来替换 $\delta c_{i}$ 。

$$\delta_m c_{i} = \begin{cases} min(|\delta c_{i}|, 2|c_{i+1} - c_{i}|, 2|c_{i} - c_{i-1}|)*sign(\delta c_{i}) &\quad {\text{if } (c_{i+1} - c_{i})(c_{i} - c_{i-1}) \gt 0 }\\ \text{0,} &\quad \text{else} \\ \end{cases} \tag {1-19}$$

该替换标准是如何获得的呢？

我们把 $c_{x_{i+1/2}}$ 的约束条件（浓度值在左右两个格子的平均值之间）带入公式 1-18，可以得到一个约束 $|\delta c_{i}-\delta c_{i+1}|<6|c_{i+1}-c_{i}|$。

更强的一个约束是 $|\delta c_{i}|+|\delta c_{i+1}|<6|c_{i+1}-c_{i}|$，那么如果$|\delta c_{i}|<|3|c_{i+1}-c_{i}|$，且 $|\delta c_{i+1}|<|3|c_{i}-c_{i-1}|$，约束条件一定会被满足（为什么PPM的条件会更强一些？）。

注意，替换之后，网格按照 $c_{x}$ 积分计算的平均值就与给定的平均值不相等了（破坏了边界处的连续性，不需要一定相等）。

为了避免在间断处出现非物理的振荡（过冲和欠冲，出现这种情况，通量的计算就会变得比较极端，不知道是对是错），在PPM方案中，还有一个约束条件，在一个网格内，浓度的分布是单调的，避免人造极值，破坏物理的合理性，比如负数（网格平均浓度很低，边界浓度差异很大，为了保持抛物线的形态，并且平均值固定，那么就可能出现负值）。

按照之前的公式，得到的 $c_{x}$ 不一定在 $c_{x_{i-1/2}}$ 和 $c_{x_{i+1/2}}$之间。

• 情况 1: 关注格子的浓度 $c_{x_{i}}$ 本身就是局地极值。插值函数设置为常数（就是在格子 $i$ 中浓度均为 $c_{x_{i}}$ ）。

• 情况 2: 格子 $i$ 位于浓度梯度很大的地方，比如左边界比右边界浓度大很多，但是该格子的平均浓度与左边界接近，那么二次多项式的曲线会形成比边界浓度更大的值（为了抵消右边比平均值低很多的情况），破坏了在一个格子里面的单调性。这种现象被称为overshoot（过充，或者超调，是指信号或者函数超过了预期值），发生overshoot时，网格左边界和右边界的值就会被重新设定，其判定条件为 $|\Delta c_i| < |c_{6,i}|$，该条件通过对公式 1-14 求导，判定恒大于0 或者恒小于0 得到。

$$\begin{align*} c_{x_{i-1/2}} &= c_{x_{i+1/2}} = c_{i} \quad \text{if } (c_{x_{i+1/2}}-c_{i}) (c_{i}-c_{x_{i-1/2}}) \le 0 \\ c_{x_{i-1/2}} &= 3c_{i} - 2c_{x_{i+1/2}} \quad \text{if } \Delta c_i c_{6,i} > (\Delta c_i)^2 \\ c_{x_{i+1/2}} &= 3c_{i} - 2c_{x_{i-1/2}} \quad \text{if } - \Delta c_i c_{6,i} > (\Delta c_i)^2 \tag {1-20} \end{align*}$$

其中，$c_{6,i} > 0$ 时，$c(x)$ 为凸函数， $c_{6,i} < 0$ 时，$c(x)$ 为凹函数函数， $c_{6,i} = 0$ 时，$c(x)$ 为线性函数，

$c_{6,i}$ 与 $\Delta c_i$ 同符号时，且判定条件满足，$c_{i}$ 更靠近右边界浓度值。按照 $c_{6,i}=\Delta c_i$，带入定义公式，求得左边界值，这时在 $x_{i+1/2}$ 处(右边界)，$c(x)$ 的导数为 0。

$c_{6,i}$ 与 $\Delta c_i$ 符号相反时，且判定条件满足，$c_{i}$ 更靠近左边界浓度值。按照 $c_{6,i}=-\Delta c_i$，带入定义公式，求得右边界值，这是在$x_{i-1/2}$ 处(左边界)，$c(x)$ 的导数为0

想象一下，如此处理之后，$c(x)$ 在该网格中表现为抛物线的极值在边界处。

注意通过该步骤进行替换之后，会破坏边界处的连续性，也就是左右极值不相等（似乎边界处的连续性不是很重要？）!

得到每个网格边界处的浓度值，就能较为容易地计算网格边界处的通量。在边界处传输的平均浓度为

$$\begin{align*} c_{i+1/2, L}(y) &= \frac{1}{y} \int_{x_{i+1/2}-y}^{x_{i+1/2}}c(x)dx \\ c_{i+1/2, R}(y) &= \frac{1}{y} \int_{x_{i+1/2}}^{x_{i+1/2}+y}c(x)dx \end{align*} \tag {1-21}$$

其中 $y$ 假设为正数，表示 $c(x)$ 在一次积分时间中，在空间上移动的距离（$y=|u\Delta t|$）。

将公式 1-14 带入，可以得到

$$\begin{align*} c_{i+1/2, L}(y) &= c_{x_{i+1/2}} - \frac{y}{2\Delta x} (c_{x_{i+1/2}} - c_{x_{i-1/2}} -(1-\frac{2y}{3\Delta x})c_{6,i}) \\ c_{i+1/2, R}(y) &= c_{x_{i+1-1/2}} - \frac{y}{2\Delta x} (c_{x_{i+1+1/2}} - c_{x_{i+1-1/2}} -(1-\frac{2y}{3\Delta x})c_{6,i+1}) \tag {1-22} \end{align*}$$

对于第 $i$ 个格子，公式中的浓度变化代表了从右边传出（$i$ 损失浓度），或者从右边传入的量（$i$ 损失浓度），二者只有一种可能，这取决于风的方向。

from sympy import *
init_printing()
ci = symbols("c_i")
cLi = symbols("cL_i")
cRi = symbols("cR_i")

x = symbols("x")
xLi = symbols("xL_i")
xRi = symbols("xR_i")

dx = symbols("dx")
c6i = symbols("c6_i")

y = symbols("y")

f = cLi + (x-xLi)/dx*(cRi-cLi+(xRi-x)/dx*c6i ) # 公式1-14
f R = integrate(f, (x, xRi-y, xRi))/y # 在边界处进行积分
# 对积分后和方程进行整理
ff = simplify(factor(simplify(fR)).subs(xRi, xLi+dx))
ff.subs(y, dx*x)

需要注意的是，该通量数值格式需要左右三个格子的浓度来计算边界处的通量，并行交换时的 halo 应该为 3。

在每个网格边界处的单位时间的平均浓度之后，针对某个网格，考虑左右边界的通量，并对时间进行积分，最终就可以获得平流的影响。

$$c_i^{n+1} = c_i^{n} + \frac{\Delta t}{\Delta x}(f_{i-1/2}-f_{i+1/2}) \tag {1-24}$$

其中$f_{i+1/2}$ 为边界处的通量。

$$f_{i+1/2} = \begin{cases} u_{x_{i+1/2}} \times c_{i+1/2, L} &\quad {\text{if } u_{x_{i+1/2}} \ge 0 }\\ u_{x_{i+1/2}} \times c_{i+1/2, R} &\quad {\text{if } u_{x_{i+1/2}} \le 0 } \\ \end{cases} \tag {1-25}$$

总结一下编程实现步骤

Step 1: 计算 2~n-1 个网格的平均坡度 (slopC)

$$\delta c_{i} = 0.5(c_{i+1} - c_{i-1})$$

为了让浓度不连续的地方模拟得更好， 同时，确保网格边界处浓度在左右格子平均浓度之间， 对 $\delta c_{i}$ 进行约束。 $$\delta_m c_{i} = \begin{cases} min(|\delta c_{i}|, 2|c_{i+1} - c_{i}|, 2|c_{i} - c_{i-1}|)*sign(\delta c_{i}) &\quad {\text{if } (c_{i+1} - c_{i})(c_{i} - c_{i-1}) \gt 0 }\\ \text{0,} &\quad \text{else} \\ \end{cases}$$

Step 2: 计算每个网格边界处的浓度 $$c_{x_{i+1/2}} = c_i + 0.5(c_{i+1}-c_i) + \frac{1}{6}(\delta_m c_{i} - \delta_m c_{i+1})$$ Step 3: 计算抛物线参数 $$c_{6,i} = 6(c_i- \frac{1}{2}(c_{x_{i+1/2}} + c_{x_{i-1/2}}) )$$

并对边界处浓度进行约束

$$\begin{align*} c_{x_{i-1/2}} &= c_{x_{i+1/2}} = c_{i} \quad \text{if } (c_{x_{i+1/2}}-c_{i}) (c_{i}-c_{x_{i-1/2}}) \le 0 \\ c_{x_{i-1/2}} &= 3c_{i} - 2c_{x_{i+1/2}} \quad \text{if } \Delta c_i c_{6,i} > (\Delta c_i)^2 \\ c_{x_{i+1/2}} &= 3c_{i} - 2c_{x_{i-1/2}} \quad \text{if } - \Delta c_i c_{6,i} > (\Delta c_i)^2 \end{align*}$$

Step 4: 计算边界处传输的平均浓度，及其通量 $$\begin{align*} c_{i+1/2, L}(y) &= c_{x_{i+1/2}} - \frac{y}{2\Delta x} (c_{x_{i+1/2}} - c_{x_{i-1/2}} -(1-\frac{2y}{3\Delta x})c_{6,i}) \\ c_{i+1/2, R}(y) &= c_{x_{i+1-1/2}} + \frac{y}{2\Delta x} (c_{x_{i+1+1/2}} + c_{x_{i+1-1/2}} -(1-\frac{2y}{3\Delta x})c_{6,i+1}) \end{align*}$$

其中 $y=|u\Delta t|$

$$f_{i+1/2} = \begin{cases} u_{x_{i+1/2}} \times c_{i+1/2, L} &\quad {\text{if } u_{x_{i+1/2}} \ge 0 }\\ u_{x_{i+1/2}} \times c_{i+1/2, R} &\quad {\text{if } u_{x_{i+1/2}} \le 0 } \\ \end{cases}$$ Step 5: 用前向欧拉法进行时间积分 $$c_i^{n+1} = c_i^{n} + \frac{\Delta t}{\Delta x}(f_{i-1/2}-f_{i+1/2}) \tag {2.13} \label {2.13}$$

1.4 边界条件

水平平流采用混合边界条件。 当边界流入区域内时，为定值边界条件( Dirichlet 边界条件)。 $$c(t, b) = c_b$$ $c_b$ 为边界处的固定值，通过并行通信获得或者外部文件读入。

当边界流出区域内时，采用齐次纽曼边界条件。

$$\frac{\partial c}{\partial x} |_{x=0} = 0$$ $$\frac{\partial c}{\partial x} |_{x=b} = 0$$

流体携带的量自然离开计算区域，不加任何人为干预（避免浓度在边界处堆积）。

混合边界条件的处理在库外部进行，只需要判断边界的风速，确定是否重新赋值边界处的浓度值。

1.5 单元测试

公式 1-12 暗含了一种检查方案，我们假定已知 $c^0(x)=p(x)$ (函数形式可以随意给定，比如 $p(x)=sin(x)$ )，且速度确定并已知，那么可以算出每个网格在任意时间的精确值(函数的形态不变，只是沿着x轴随时间平移)

$$c_i^n = \frac{1}{\Delta x} \int_{x_{i-1/2}}^{x_{i+1/2}} p(x-ut^n) dx$$

假定 $p(x)=sin(x) + 1$。

网格的平均浓度值为：

$$c_i^n = \frac{1}{\Delta x} \int_{x_{i-1/2}}^{x_{i+1/2}} (sin(x-ut^n) + 1) dx$$

积分后可得：

$$c_i^n = 1 + \frac{1}{\Delta x} \left[ \cos(x_{i-1/2} - ut^n) - \cos(x_{i+1/2} - ut^n) \right]$$

1.6 参考文献

1. Sportisse B. Fundamentals in air pollution: from processes to modelling[M]. Springer Science & Business Media, 2009.

2. Colella P, Woodward P R. The piecewise parabolic method (PPM) for gas-dynamical simulations[J]. Journal of computational physics, 1984, 54(1): 174-201.