1 引 言

在气象数值模式中，平流作为一种最基本的动力过程对大气环境场和中小尺度运动场的演变都起着关键作用，因此，高精度物质平流方案在数值模式中的有效应用对提高数值预报质量至关重要。对于大气示踪物(如水汽)平流，平流计算算法的正定性更为重要，而在气候数值模式中，计算守恒性也是对水汽平流计算的重要要求。半拉格朗日平流方案由于具有良好的计算稳定性和高效的特点，在计算流体力学和大气数值模式中得到广泛应用(沈学顺等，2011)。在三次约束曲线插值法(Constrained Interpolation Profile，简称CIP)原理的基础上，Xiao等(2002)提出了一种基于有理函数插值的CIP-CSLR(守恒型有理函数插值半拉格朗日法，Conservative Semi-Lagrangian with Rational function，简称CSLR)物质平流方案，利用单元格上的积分平均值和界面值，通过半拉格朗日后向追踪方法和分段有理插值函数计算上游点值，进而计算单元格积分值，有理函数的正定性保证了物质平流计算的非负要求，消除了高阶多项式插值在不连续界面处带来的虚假数值振荡。同时采用通量形式平流方程，利用单元网格界面的通量收支平衡计算单元格内积分平均值，保证了平流计算的守恒性。Xiao等(2002)、Peng等(2005)、常燕等(2008)曾对CSLR物质平流方案在一些理想试验和气象数值模式实际应用中进行了检验，并与其他同类算法(如CIP，PPM(Piecewise Parabolic Method))进行了比较，发现CSLR物质平流方案在保证守恒型的同时，获得了正定的结果，取得了较理想的预报和模拟效果，但是，都没有明确指出该方法的计算精度和收敛速度阶数，而收敛速度是平流数值算法应用于高分辨率数值模式时所要考虑的。另外，由于多维有理函数构建时面临困难，CSLR物质平流方案向多维空间扩展时采用简单的分维计算技术，这种分维算法易于实现，容易编程，但是，当流场在x方向或y方向存在辐合、辐散时，就会产生虚假的梯度，那么这种分维算法多大程度上影响多维CSLR物质平流方案的计算精度呢？为了加深对CSLR物质平流方案计算性能的进一步了解，本研究着重分析讨论CSLR物质平流方案在平面直角坐标和球面坐标中的计算性能，并评估分维技术的计算精度和数值收敛问题，为其在复杂数值模式中的实际应用提供理想试验资料和参考。2 CSLR物质平流方案及其误差评估

一维无辐散平流问题的通量型平流方程可写为

式中，t为时间，x为空间坐标，u为特征速度，φ是被平流物理量，此方程是数值模式动力框架的重要组成部分。CSLR物质平流方案就是采用CIP方法的原理在有限体积法概念基础上计算单元网格的通量收支积分来实现平流计算，在单元格界面点上采用半拉格朗日后向追踪方法确定界面处的质量通量，上游出发点采用一维平流速度的三阶泰勒展开计算，单元格内积分平均值则按有限体积单元网格两侧界面的通量收支得到更新。CSLR物质平流方案中上游出发点处物理量的空间重构由一维有理函数(Xiao et al，2002)来确定。系数a_i、b_i、c_i、β_i可由单元格积分平均值和界面值计算，其中，式中，φ是单元格积分平均值，f是单元格界面值，ε是接近于0的正数，以防止分母为0。参数β_i是由网格界面处f的连续性确定的。对于光滑连续函数，β_i趋近于0，则有理函数R(x)表现为二次曲线函数(CSL2)，保证较高的计算精度；而对于不连续分布，单元格界面值f和积分平均值φ相差较大，β_i的绝对值增大，有理函数R(x)趋于线性函数，保证数值解的正定性和单调无振荡。因此，CSLR物质平流方案利用有理函数本身特性就能很好地刻画不连续或强梯度分布的计算，而不需要额外的单调性限制器。为便于比较，本文同时给出CSLR物质平流方案和CSL2方案的计算结果，在平面直角坐标系和球面坐标系中分析CSLR物质平流方案的计算效果和数值收敛性。

关于数值计算收敛速度的评估方法较多(马利斌等，2010; Lauritzen et al，2012)，本研究采用分辨率倍增后误差减小倍数的对数和分辨率提高倍数的对数之比作为收敛阶数，即在柯朗数(Courant-Friedrichs-Lewy，CFL)(n_CFL=|u|Δt/Δx)一定情况下，空间网格距为Δx时的计算误差为E₁，当分辨率提高1倍(格距为Δx/2，时间步长为Δt/2)后，积分计算相同时间后的计算误差为E₂，则定义其收敛阶数

高维计算类似，这种收敛率计算方式同时考虑了空间分辨率和时间积分步长的影响。同时为了更加直观地表示收敛速度，还采用了Lauritzen等(2012)的方法，即利用最小二乘法计算误差收敛阶数

式中，L_i是误差范数L₁、L₂、Linf、ERMS，A_i是某一常数，K_i是误差收敛阶数，Δλ表示某一空间网格分辨率。本文中计算区域的总积分量数值为则计算误差(Williamson et al，1992; Zhang et al，2012)定义为均方根误差为其中，φ_n表示节点的数值解，φ_e则表示节点的精确解，Ω表示计算区域，N表示计算节点总数。3 数值试验

为全面考察CSLR物质平流方案的数值收敛性，评估引入空间分维算法及其在曲线坐标系中的计算效果，本研究按照平面二维和球面坐标的顺序给出数值试验结果。在理想试验的选取和设计上，参考了Williamson等(1992)、 Xiao等(2002)、沈学顺等(2011)、Zhang等(2012)常用的理想试验设定，研究平流函数连续性、尺度大小和水平风场辐散、辐合等对模拟结果的影响，试验设计由简单到复杂，以求全面了解CSLR物质平流方案的数值计算特点。

首先在平面直角坐标系下检验了一维CSLR物质平流方案的计算性能(图略)。对于光滑平缓的物理量分布(如正弦波)，CSLR物质平流方案获得了二阶收敛速度和严格的正定单调数值解，CSL2方案(β=0)则获得三阶收敛精度，但数值解有微小数值振荡；对不连续分布(如方波)的平流计算，CSLR物质平流方案仍然给出了很好的保形和单调无振荡特性数值结果，而CSL2方案在间断处给出明显的振荡数值解，表现为“上冲”和“下冲”。此外，两者均能准确获得数值解相位，体现了半拉格朗日方法具有很强的位相描述能力。

3.1 平面直角坐标系下二维理想试验

在实际应用中，面临二维、三维平流计算，现行方法主要通过分维技术来实现(Xiao et al，2002; Peng et al，2005)，即将多维问题进行空间和时间的分离，分解成几个连续的一维问题来求解。分维技术简单方便，计算量小，所以被广泛应用，但是分维技术将多维平流过程人为划分成多个一维平流过程，造成平流过程在时间和空间上的分离，势必引入一定的计算误差，并且，当流场在x方向或y方向存在辐合、辐散时，也会产生虚假的梯度，针对分维技术的误差订正方法已有研究(Clappier，1998)，本研究主要就二维问题中应用简单分维技术后分析CSLR物质平流方案的收敛速度和计算效果。对于二维通量型平流方程

式中，t为时间，x、y为空间坐标，u、v为特征速度，单元格积分平均值记为φ。将上述二维平流方程分解成x和y方向两个相对独立的一维平流方程分别计算，计算顺序和具体步骤详见Xiao等(2002)。对于二维计算过程，采用先x方向后y方向的固定顺序进行，在每一个方向完成积分平均值的更新之后，再利用单元格积分平均值的增量对另一个方向单元格上的界面值进行一次中央差分修正。

3.1.1 二维正弦波试验

首先针对平缓光滑的二维正弦波平流进行分析，平流物理量场的初始条件为二维正弦波动

计算区域为(x，y)∈[-1，1]×[-1，1]，其解析解可写为采用周期性边界条件且水平均匀速度场为u=v=1.0 m/s，计算总时间为5个周期即10 s，n_CFL=1.04，n_CFL大于1，体现出半拉格朗日方法计算绝对稳定的优点。由于采用了分维技术，柯朗数定义为式中，Δx、Δy分别表示x和y方向上均匀网格格距，Δx=L/N_x，Δy=L/N_y。图 1 给出CSLR物质平流方案和CSL2方案的各种计算误差收敛率，与一维正弦波试验结果不同，在相同计算条件下CSL2方案的计算误差要略大于CSLR物质平流方案，即对于光滑平缓的二维正弦波动平流计算，CSLR物质平流方案的计算精度要优于CSL2方案，说明在二维平流计算中，CSL2方案受分维计算影响较大，CSL2方案的收敛速度也由一维情况下严格取得三阶降为二维计算中的二阶左右。同样，CSLR物质平流方案在二维情况下的收敛速度也有所下降，可见分维技术对数值收敛速度的影响比较明显，同时也反映出分维技术在多维平流计算中仍需进一步改进。

3.1.2 凹槽圆柱平流试验

在二维计算中选用凹槽圆柱平流试验(Zhang et al，2012)来考察CSLR物质平流方案对于不连续分布的模拟能力，同时检验分维技术对有理函数保证计算结果正定性的影响。凹槽圆柱体在角速度为ω的旋转风场中运动，初始场为

式中，圆柱高度φ₀和圆柱半径σ分别是某一常数，W_s和L_s分别是凹槽的宽度和长度。ξ和ζ分别是相对于圆柱中心((x_c-γcos(ωt)，y_c-γsin(ωt)))的相对坐标式中，r是圆柱半径，γ为流场中心(x_c，y_c)到圆柱中心的距离。

选择角速度ω=1.0的逆时针旋转速度场，圆柱高度φ₀=1，圆柱最大半径σ=20，W_s=10，L_s=25，γ=25，风场中心为(50，50)，n_CFL最大值为1。图 2给出旋转风场示意图、初始场以及用CSLR物质平流方案和CSL2方案计算凹槽圆柱旋转5圈(时间步长是0.02 s，积分步数为1571)后回到出发点的模拟结果，CSLR物质平流方案能够保证不连续处的正定分布，不会产生新的极值，而CSL2方案由于采用二次插值函数则在不连续处产生明显的数值振荡，即带来了虚假上冲和负值，产生了较大的频散误差，破坏了标量场固有的空间分布特点和守恒特性。由此，本试验说明CSLR物质平流方案在向二维问题扩展中，针对这种空间变化幅度大、存在大梯度和强间断的标量场分布，通过分维技术，CSLR物质平流方案仍能够保持一维有理函数的基本特性，可以保证计算的正定和守恒性，这一点对流体物质平流输送是非常重要的特性。

3.1.3 余弦钟平流试验

采用孤立的二维平缓余弦钟形物做进一步的平流试验，相对凹槽圆柱的不连续分布，这种连续平缓余弦钟平流试验更适合对数值算法的收敛速度评价。初始场为

式中，余弦钟的半径R=，首先选择均匀定常速度场，这里R₀=0.2，x₀=0，y₀=0，n_CFL=0.5，给定周期性边界条件，计算区域为[-0.2，0.2]×[-0.2，0.2]，积分时间为余弦钟运动一个周期T=0.4 s。图 3给出不同分辨率下CSLR物质平流方案和CSL2方案的计算误差及网格收敛率，可以看出，在这个试验中CSLR物质平流方案的计算误差要比CSL2方案略大，其收敛阶数不到2阶，而CSL2方案的收敛阶数平均能达到2阶；从实际模拟结果上来看(图略)，在同等分辨率下CSL2方案的模拟效果优于CSLR物质平流方案，两种算法数值解的位相和精确解保持一致，说明分维算法保持了半拉格朗日方法精确位相求解的特点。

在相同初始分布条件下，还进行了无辐散旋转风场下余弦钟平流试验，即将二维风场改为，其中，角速度ω=1，n_CFL=0.5。图 4给出余弦钟旋转一个周期T=2π s后的计算误差和对应的网格收敛速度，可以看出，在这种旋转风场下，CSLR物质平流方案和CSL2方案的计算误差基本达到相当水平，而两者的收敛阶数相对于均匀风场下也都有下降，且随分辨率提高收敛阶数都有所降低，这与在分维计算中二维风场存在切变梯度(如u在y方向的梯度)造成上游点后向追踪计算误差增大有关。

3.1.4 平面变形流场试验

上述二维试验所用的二维均匀风场无论一维还是二维方向都是无散度流场，平流过程中没有形变，其中，一维问题的半拉格朗日后向追踪计算都可以精确得到上游点位置，而在旋转凹槽圆柱试验和余弦钟旋转平流试验中，虽然二维散度为0，但风场存在切变梯度，在利用分维技术计算二维平流问题时就会带来一定的计算误差。本试验条件更为苛刻，进一步考察CSLR物质平流方案在单一方向风场散度不为0的变形流场中的模拟情况(Zhang et al，2012)，评估CSLR物质平流方案在这种相当于有源汇项存在的流场中的计算性能。数值模拟计算区域为[1，100]×[1，100]，同样设置周期边界条件，初始场计算格点均匀为1，速度场是一种二维方向无辐散而一维方向存在风速梯度的强变形流场

解析解(真值)在这种无辐散流场下应该是均匀、定常初始值，不随时间和空间变化。图 5给出CSLR物质平流方案和CSL2方案获得平流100步后的模拟结果，可以看出，对于这种强烈的变形流场，两者利用半拉格朗日方法计算上游点均可取得较理想的模拟效果，整体计算误差和单个网格点的绝对误差都比较小，对于这种均匀连续分布的物理场模拟，分维技术和有理(二次)插值函数都没有产生严重的计算误差。进一步对比两个算法的计算结果，可以看出CSL2方案的绝对误差要比CSLR物质平流方案大，CSLR物质平流方案的计算结果中格点最大值为1.10852，最小值为0.89195，而CSL2方案计算结果中的格点最大值为1.11088，最小值为0.75698，说明在本试验中CSL2方案计算能力稍逊于CSLR物质平流方案，在平流计算中产生了较大的极值，而CSLR物质平流方案借助有理插值函数能够更好地处理这种变形流场中的平流模拟。

3.1.5 理想锋生试验

平面锋生试验(Zhang et al，2012)包含一个二维定常切线速度为

的圆形涡旋，r是到涡旋中心(x_c，y_c)的径向距离，V₀是常数，平流标量解析解为式中，ω=V(r)/r是角速度，σ是锋区的特征宽度。锋生试验针对的是实际天气现象中水汽输送以及温度、湿度和风场等不连续分布情况(如梅雨锋)，这种试验可以验证平流方案利用较少信息对细微结构的捕捉能力，平流输送作为数值模式动力框架部分的重要过程之一，较好的平流方案是数值预报和模拟成功的关键因素。取V₀=4，σ=5，计算区域为[1，10]×[1，10]，图 6给出初始物理场分布和流场以及两种平流方案在分辨率为100×100下积分3.0 s(时间步长是0.05 s，积分60步)后的数值解及其对应的解析解，对于这种平流物理量反对称分布的锋生试验来说，两者都能模拟出锋生的演变过程，只是对于涡旋中心锋面强梯度之处，CSLR物质平流方案的模拟效果要比CSL2方案的模拟效果更接近解析解，更能描述平流物理量的空间分布特点，再次说明利用有理插值函数的CSLR物质平流方案对于空间变化幅度大的物理量具有很好的平流模拟能力，这也是CSLR物质平流方案的优势所在。

3.2 阴阳网格上球面二维波动试验

前面在平面直角坐标上进行了一维和二维理想试验，可以看出在不同维数计算和针对不同平流问题时，CSLR物质平流方案和CSL2方案的计算性能及收敛速度均有不同表现。最后利用球面准均匀网格——阴阳网格(Kageyama et al，2004; Peng et al，2006; 李江浩等，2013)分析球面坐标系下CSLR物质平流方案的计算性能和收敛性，采用球面阴阳网格是因为其避免了传统经纬度网格所具有的“极点问题”，选择初始物理场为光滑平缓的高度场(Li et al，2012)

并选择散度为0的速度场其中，λ是经度，φ是纬度，u₀=2πa/(12 d)是最大的速度值，约为40 m/s，a为地球半径，α代表旋转轴与z轴(南北极轴)的夹角。计算时间为12 d，波动能够回到初始位置，n_CFL=0.5。

选择α=0，采用3种网格分辨率，它们相应的时间步长分别是360、180、90 s，平流积分12 d回到初始位置。图 7给出在本试验中CSLR物质平流方案和CSL2方案的计算误差与网格收敛率，可以看出，CSLR物质平流方案的网格收敛率大概在1—2阶，范数误差ERMS相对另外3种误差的收敛率较高，而CSL2方案的网格收敛率相比前者要略小，因此，两者在相对复杂的球面阴阳网格上的平流计算中，网格收敛速度都维持在1—2阶；进一步比较两者的计算误差，在同等空间分辨率和积分时间步长下，经过相同时间积分后CSL2方案的计算误差略小于CSLR物质平流方案的计算误差，从两者的实际模拟结果(图 8)上也可看出，CSL2的模拟效果要优于CSLR物质平流方案，这与一维波动试验得出的结论一致；另外，与直角平面平流试验相比，可以发现球面曲率的存在影响了算法的计算精度和收敛速度，这可能与球面上半拉格朗日出发点的计算有关。

4 结 论

本文讨论和分析了平面直角坐标系和阴阳网格球面坐标系下CSLR物质平流方案的计算性能，与CSL2方案对比给出了不同理想数值试验下的网格收敛速度和计算精度，得到以下结论：

(1)采用一维有理函数插值的CSLR物质平流方案在一维计算中，能够保持正定和单调，并且在采用分维技术的多维计算中仍能够保持一维有理函数的性质，对物质平流计算具有较好的刻画。

(2)CSLR物质平流方案在一维光滑平缓的物理场计算中可以获得约二阶的网格收敛速度，而CSL2方案能够严格取得三阶收敛速度，有理函数的正定单调性一定程度上降低了CSLR物质平流方案本身的计算精度；但在有不连续物理场情况下，CSLR物质平流方案具有正定性，可以更好地计算物质平流，得到更好的计算质量。在多维计算中采用了简单的分维技术，而这种空间/时间平流分离计算降低了平流方案的计算精度和数值收敛速度，但两者仍能取得1—2阶的网格收敛速度。

(3)对于均匀无辐散流场中平缓波动和余弦钟等光滑连续函数平流计算时，CSL2方案的网格收敛速度和模拟效果要优于CSLR物质平流方案，而对于理想锋生试验和有源汇项存在的变形流场等存在不连续分布的平流试验，CSLR物质平流方案表现要优于CSL2方案，因此CSLR物质平流方案在增强不连续分布平流计算能力的同时，对平滑函数的计算略有负面影响。

(4)球面曲率影响半拉格朗日出发点的计算，CSLR物质平流方案和CSL2方案在球坐标中的网格收敛速度相对平面直角坐标下都有所降低。