高阶正定守恒的中心约束多矩有限体积平流模式

奥鹏网院作业

高阶正定守恒的中心约束多矩有限体积平流模式*
舒 谦1 唐 杰2,3 陈春刚4 肖 锋5 沈学顺2,3周立隆2,3 李泽椿3 朱克云1 李兴良2,3

1. 成都信息工程大学，成都，610225

2. 中国气象局数值预报中心，北京，100081

3. 国家气象中心，北京，100081

4. 西安交通大学，西安，710049

5. 日本东京工业大学机械工程系，东京，152-8850

摘 要 采用新的均匀三点中心约束多矩有限体积方法（3-point Multi-moment Constrained finite-Volume scheme for Uniform Points with Center Constraints， MCV3_UPCC），发展了一个三阶正定守恒的平流模式。三点多矩有限体积方法在单网格内定义等距的3 个自由度，采用多矩约束条件并通过控制方程获得时间演变方程。新的三点中心约束多矩方法能在单网格内采用等距的3 个点值及中心一阶、二阶导数作为约束条件进行空间4 次多项式数值重构，获得3 个自由度的时间演变方程；所构建的新数值方案具有三阶精度，边界通量连续性保证了其数值严格守恒。为了抑制该方法的非物理数值振荡，引入了边界保型限制器技术，它能够把数值解控制在既定物理场最小值（最小值为0 时则保持数值正定）与最大值之间。数值试验表明新发展的三阶平流模式具有良好的计算精度，能够严格保持数值解的正定性和守恒性，同其他高精度平流模式相当，在实际大气模式水汽等平流输送应用中具备良好的发展潜力。

关键词 平流方案， 多矩有限体积法， 正定守恒， 高精度方案

1 引 言
大气中水物质（如水汽、云水等）分布复杂，具有梯度大、不连续、多相变等特点，使得数值模拟水汽等输送过程容易产生非物理的负值及数值振荡。因而，采用良好的数值方法精确地模拟水物质的分布与输送，对数值天气预报输送模式的改进，特别是准确模拟降水过程，具有重要意义。

低精度的一阶迎风格式不产生数值振荡，但是数值耗散大，满足不了实际高精度模式的需求。高阶（大于一阶）平流算法精度较高，但是会产生非物理数值振荡，导致水物质输送产生较大误差（沈学顺等，2011），例如基于欧拉方法的Lax-Friedrich方法、Lax-Wendroff 方法、MacCormach 方法等。为了抑制数值振荡，中外学者提出了一系列的方法，例如通量修正（Book，et al，1975）方案，基于Godunov 方 法 的 MUSCL （ Monotone Upstreamcentered Schemes for Conservation Laws;van Leer，1977）方案，以及基于分段抛物线函数的PPM（Piecewise Parabolic Method ；Colella，et al，1984）方案，两步保型平流方案（TSPAS）（Yu，1994）等。

与欧拉方法相比，半拉格朗日方案（Robert，1981）可取较大的时间步长，效率较高，数值耗散低，被广泛应用于大气数值平流模式。Bermejo 等（1992）设计出一个混合算法使传统的半拉格朗日方法转化为准单调的半拉格朗日（QMSL）方法，该方法在强梯度、不连续区域低阶插值不能保证守恒。Zerroukat 等（2002）提出的守恒高效型半拉格朗日方案（ SLICE）和 Lauritzen 等（2009）提出的守恒型半拉格朗日多示踪物方案（CSLAM），Xiao等（2002）提出的分段有理数插值（PRM）方案和CSLR （Conservative Semi-Lagrangian with Rational function）方案，都能保证平流方案的守恒性。为了获得良好的计算精度，Guo 等（2014）把半拉格朗日方法和不连续伽辽金算法结合，设计出了SLDG（Semi-Lagrangian Discontinuous Galerkin）方案，该方案能够保证守恒，精度大于二阶，并且能通过限制器去除数值振荡。

多矩有限体积方法（ Multi-moment finite volume scheme）是近年提出来的一种新的有限体积方法（Yabe，et al，1991，2001；Xiao，et al，2001，2006，2013；Xiao，2004；Ii，et al，2007，2009；Chen，et al，2008），该方法引入多矩概念，如体积积分平均值和点值、导数值进行局地高阶重构，通过控制方程更新预报变量时间向前积分（Li，et al，2015；Chen，et al，2014）；其积分平均值约束条件能严格保证数值守恒，是传统有限体积方法的拓展。文中采用的新三点中心约束多矩有限体积方法（3-point Multimoment Constrained finite-Volume scheme for Uniform Points with Center Constraints，MCV3_UPCC；Xiao，2012）能够基于单个网格构造高阶插值多项式，获得高阶计算精度，并且该格式计算简洁、计算模板紧致。该方案与原三阶多矩有限体积方案（MCV3）（Ii，et al，2009）相比，具有网格边界不连续，计算精度高，不用求解边界导数黎曼问题，计算效率高等特点。和其他高阶方法一样，多矩约束有限体积方法在非光滑解附近仍会出现数值振荡，并可能导致负值。文中采用边界保型（BP）限制器对其进行修正，保持该方案的正定性，以期设计一种可以用于大气模式的高精度平流输送模式。

2 均匀三点中心约束多矩有限体积方法
文中采用三点均匀中心约束多矩有限体积法，该方案在单个网格内进行局地高阶通量重构，约束条件施加在单元网格中心点值及其一阶和二阶导数上，保证单元边界处通量函数的连续性，通过控制方程得到点值矩的时间演变方程。该方法计算精度高，严格保证数值守恒。MCV3 方案局地高阶通量重构约束限制发生在单元网格界面（Ii，et al，2009），与之相比，MCV3_UPCC 方案约束主要体现在单元网格中心上，具有较高的精度和允许更大的CFL（Courant-Friedrich-Lewy）条件，并且计算模板紧致。针对格式产生的数值振荡，引入边界保型限制器，该限制器优势在于不影响试验误差和精度的同时，抑制间断处的数值振荡并保持正定。

2.1 均匀三点中心约束多矩有限体积方法数值离散
考虑一维守恒平流方程

pagenumber_ebook=120,pagenumber_book=667
pagenumber_ebook=120,pagenumber_book=667
图 1 一维单元上的矩和节点值
Fig. 1 The locations of unknowns in one dimension

pagenumber_ebook=120,pagenumber_book=667
首先，通过单元内的3 个点值，利用拉格朗日多项式构造出在单个网格内不连续的主重构函数fi(ξ)， 表达 式 为 f i(ξ)=uipi(ξ)， ui为 速 度， p i(ξ)为 关于索引单元i内点值的拉格朗日插值多项式，主重构函数 fi (ξ)表达式如下

pagenumber_ebook=120,pagenumber_book=667
pagenumber_ebook=120,pagenumber_book=667
pagenumber_ebook=121,pagenumber_book=668
为了验证该方案的守恒性，积分守恒方程的平流项得到

pagenumber_ebook=121,pagenumber_book=668
由此可见，式（1）单元内守恒量的变化等于单元边界通量的净变化，亦即单元网格界面通量是连续的，因而均匀三点中心约束多矩有限体积方法能够保证数值严格守恒。

2.2 边界保型限制器
边界保型限制器（Zhang，et al，2010）是通过重新构造局部插值函数来起到限定的效果，该限制器计算精简，应用方便，能够把数值解控制在初始物理场的极大、极小值范围内，并且不损失守恒性（Katta，et al，2015；Guo，et al，2016）。边界保型限制器是基于Liu 等（1998）所设计的限制器，其构造原理如下

pagenumber_ebook=121,pagenumber_book=668
pagenumber_ebook=121,pagenumber_book=668
可以看到，新重构的插值函数pagenumber_ebook=121,pagenumber_book=668能够保持原有积分平均值不变，因而不损失数值守恒性。

2.3 二维数值离散推广
二维笛卡尔坐标下的通量形式平流方程为

pagenumber_ebook=121,pagenumber_book=668
式 中 ， e =uq， f =vq ， q 为 守 恒 变 量 ， u 和 v 分 别 为 x和y方向的速度。

在笛卡尔坐标结构网格下（图2），可以直接将一维数值离散推广到二维

pagenumber_ebook=121,pagenumber_book=668
pagenumber_ebook=122,pagenumber_book=669
图 2 二维单元内点值 （求解点） 分布
Fig. 2 The locations of solution points （unknowns） in two dimensions

2.4 时间积分
pagenumber_ebook=122,pagenumber_book=669
3 数值试验
pagenumber_ebook=122,pagenumber_book=669
式中，pagenumber_ebook=122,pagenumber_book=669和pagenumber_ebook=122,pagenumber_book=669分别为索引单元i 内的精确解的积分平均值和单元内变量的积分平均值。

3.1 一维收敛性试验
一维收敛性试验初始场为

pagenumber_ebook=122,pagenumber_book=669
式（1）中的平流速度 u =1，边界条件为周期性边界条件，试验设置的库朗数为0.1，周期时间t=2。

从表1 标准误差中可以看出，MCV3_UPCC 方案具有三阶收敛。考虑到MCV3_UPCC 方案单网格内自由度和高次多项式构建的优势，跟Ii 等（2009）MCV3 方案数值误差相比（文献中的 L1、L∞误差和文中误差算法不同，根据文献中误差算法），MCV3_UPCC 方案的标准误差要比MCV3 方案小，例如在40 个网格数分辨率下MCV3_UPCC方案的积分平均值的标准误差比MCV3 方案减小45%。由表2 可知，采用边界保型限制器后的MCV3_UPCC 方案 L1 、 L2误差和阶数几乎没有变化，仅 L∞误差有些增大与阶数略降。

3.2 正弦叠加波试验
为了严格测试数值格式模拟强梯度和保持极值的能力，采用Blossey 等（2008）的试验设置，初始场由两个等振幅正弦函数叠加

pagenumber_ebook=122,pagenumber_book=669
试验区域 x ∈[0,1]，网格数为30，平流速度u=1， 边界条件为周期边界。在 t =1时刻的数值模拟的误差及其极大、极小值如表3 和图3 所示。

从表3 中可以看出，MCV3_UPCC 方法的 E2误差比MCV3 方法减少了45%，并且数值模拟的极值较MCV3 更接近解析解。对比Blossey 等（2008）中的通量修正方案 E2误差0.05，MCV3_UPCC 方案的误差为0.036，误差较小。在MCV3_UPCC 方案加上边界保型限制器后，MCV3_UPCC_BP 方案的误差几乎没有改变，且最大值、最小值和MCV3_UPCC方案差别不大。对比Blossey 等（2008）的通量修正方案，由于其使用了单调限制器，数值解在 x =0.5附近耗散较强，而图3c 中MCV3_UPCC_BP 限制器能够很好地保持局地峰值，效果更佳。同时，相比于MCV3 方法，采用MCV3_UPCC 方法对于保持极值效果更好，一个周期后的数值模拟积分平均值（VIA）更接近初值，即使引入了边界保型限制器（图3c），几乎对波峰处的极值也没有影响。

表 1 MCV3_UPCC 标准误差和收敛阶数
Table 1 Numerical errors and convergence rates of the 1D linear scalar equation of MCV3_UPCC scheme

pagenumber_ebook=123,pagenumber_book=670
表 2 MCV3_UPCC_BP 标准误差和收敛阶数
Table 2 Numerical errors and convergence rates of the 1D linear scalar equation of MCV3_UPCC scheme with BP filter

pagenumber_ebook=123,pagenumber_book=670
表 3 正弦叠加波标准误差和极值 （按照Blossey 等（2008）试验中的误差算法）
Table 3 Errors and maximum/minimum values of advection of sum of two sine waves（According to the method of Blossey，et al，（2008））

pagenumber_ebook=123,pagenumber_book=670
pagenumber_ebook=123,pagenumber_book=670
图 3 正弦叠加波试验 （a. MCV3，b. MCV3_UPCC，c. MCV3_UPCC_BP）
Fig. 3 Sum of 7.5 and 10 Δ x sine waves experiments（a. MCV3，b. MCV3_UPCC，c. MCV3_UPCC_BP ）

为了测试平流方案的保正定性，取式（29）中函数的正值部分，区间和网格同正弦叠加波试验，即

pagenumber_ebook=123,pagenumber_book=670
图 4 为MCV3、MCV3_UPCC、MCV3_UPCC_BP 3 种方案模拟正弦叠加波数值结果。可以看到，MCV3_UPCC 方案对波峰的模拟极值高于MCV3，但它们在间断附近有明显的数值下冲，即存在非物理数值振荡，表4 中两方案的极小值（q m in）小于0。然而，引入BP 限制器后，MCV3_UPCC_BP 方案能够很好地抑制下冲，同时不影响波峰处的极值。从表4 可以看到，MCV3_UPCC_BP 方案的极小值（-3.4694E-018）为机器误差，即保持数值解的正定性。

pagenumber_ebook=124,pagenumber_book=671
图 4 正值部分的正弦叠加波试验 （a. MCV3，b. MCV3_UPCC，c. MCV3_UPCC_BP）
Fig. 4 The positive parts of sum of 7.5 and 10 Δ x sine waves experiments（a. MCV3，b. MCV3_UPCC，c.MCV3_UPCC_BP）

表 4 正值部分正弦叠加波标准误差和极值 （按照Blossey 等（2008）试验中的误差算法）
Table 4 Errors and maximum/minimum values of advection of sum of 7.5 and 10 Δ x positive sine waves（According to the method of Blossey，et al，（2008））

pagenumber_ebook=124,pagenumber_book=671
同样，对比Blossey 等（2008）的通量修正方案在间断处的模拟结果，文中MCV3、MCV3_UPCC和MCV3_UPCC_BP3 种方案计算精度高，它们的数值模拟误差 E2 和 E∞均低于文献结果。

3.3 一维方波平流
为了验证数值格式在强间断的附近是否能够有效地抑制非物理数值振荡，设计的方波试验初值分布如下

pagenumber_ebook=124,pagenumber_book=671
一维的计算区域为[-1，1]，其内网格数为200，网格距0.01，水平速度场 u =1。 在 t =2时刻MCV3、MCV3_UPCC 和 MCV3_UPCC_BP 3 种方案数值模拟结果如图5 所示。从图中可以看到，在x=±0.4处 MCV3 和 MCV3_UPCC 的点值（PV）和积分平均值（VIA）存在明显的数值上冲和下冲；这一点也可以从表5 中的极大值（q m ax ） 大于1 和极小值（q m in）小于0 看出。此外，从表5 还可以看到，MCV3_UPCC和MCV3_UPCC_BP 方案的标准误差要比MCV3方案小；同时，MCV3_UPCC 算法引入边界保型限制器能够有效抑制振荡，并能在平流输送过程中把点值严格控制在初始场的最大值和最小值之间。

pagenumber_ebook=125,pagenumber_book=672
图 5 方波试验结果 （a. MCV3，b. MCV3_UPCC，c. MCV3_UPCC_BP）
Fig. 5 The square wave experiments （a. MCV3，b. MCV3_UPCC，c. MCV3_UPCC_BP）

表 5 一维方波标准误差和极值
Table 5 Errors and maximum/minimum values of advection of a square wave in 1D

pagenumber_ebook=125,pagenumber_book=672
3.4 二维收敛性试验
二维收敛性试验初始场分布如下

pagenumber_ebook=125,pagenumber_book=672
试验速度 u =1, v=1， 周期 t =2，边界条件为周期性边界条件，库朗数为0.1。

数值试验结果如表6 所示，从表中可知，MCV3_UPCC 方案获得了三阶收敛阶数，并且比MCV3 方案误差减小了一半。MCV3_UPCC 方法引入边界保型限制器后，对试验误差几乎没有影响， L1阶数没有改变， L2 和 L∞有所减小，与一维试验结果完全一致。

3.5 二维复杂波试验
为了验证数值格式同时捕获光滑和间断处的解的能力，一维复杂波试验最早由Jiang 等（1996）提出，Ii 等（2009）将它推广至二维，该试验由高斯波、方波、尖三角波和半椭圆波组合而成，速度场分布为pagenumber_ebook=125,pagenumber_book=672计算区域为 [- 1,1]×[-1,1]，试验时间为一个时间周期，网格数为100×100，时间步长为1500 步，初始场分布为

pagenumber_ebook=125,pagenumber_book=672
pagenumber_ebook=125,pagenumber_book=672
pagenumber_ebook=125,pagenumber_book=672
从图6 中可以看出，MCV3_UPCC 方法和MCV3方法一样，都在强间断处出现了数值振荡，并且两种方法对于波峰处的极值保持都不错，MCV3_UPCC方法引入边界保型限制器后没有影响到波峰处的极值。由表7 可知，加上边界保型限制器后，整个波的点值在精确解的最大值和最小值之间，除去了数值振荡，而且得到的 L1误差较MCV3_UPCC 方法的误差更小。综上所述，边界保型限制器在保持守恒和正定的同时，能够很好地模拟复杂波的平流过程。

表 6 二维收敛性试验标准误差和收敛阶数
Table 6 Numerical errors and convergence rates of the 2D linear scalar equation

pagenumber_ebook=125,pagenumber_book=672
pagenumber_ebook=126,pagenumber_book=673
图 6 二维复杂波试验结果 （a. 初始场，b、c、d 分别为MCV3、MCV3_UPCC 和MCV3_UPCC_BP 旋转一周后的结果）
Fig. 6 The complex waves experiment （a. the initial field of complex waves，b. MCV3，c. MCV3_UPCC，d. MCV3_UPCC_BP）

3.6 凹槽圆柱旋转试验
该试验主要是验证平流的保型能力，在二维网格数为100×100，时间积分总步数为1500 步，时间周期t=1 的圆柱凹槽试验中，初始场半径为0.5，中心位于原点，高为1 的凹槽圆柱，速度场为pagenumber_ebook=126,pagenumber_book=673试验初始场分布为

pagenumber_ebook=126,pagenumber_book=673
式中， σ为圆柱半径， Ws 和 Ls分别为凹槽的宽和长，ξ 和 ζ分别为相对于圆柱中心的坐标，长和宽分别为0.4 和 0.74。

从图 7b、c 可知，MCV3_UPCC 方案具有良好的保型能力，在旋转了10 周后凹槽圆柱的形状保持较好，只是在强间断处会产生小的数值振荡。从表8 可知，应用边界保型限制器后，在对凹槽圆柱有很好保型效果的同时，能有效去除数值振荡，并把最大、最小值严格控制在初始场的最大、最小值，即0 到1 之间。对比苏勇等（2013）的凹槽旋转试验，采用准单调半拉格朗日方案在旋转10 周后极值被削弱到了0.953，而MCV3_UPCC 方法引入边界保型限制器旋转10 周后极值还能维持初始场的最大值1，说明MCV3_UPCC 方案的计算精度高，保型能力强。

3.7 二维变形场试验
数值试验采取Durran（1999）的试验，该试验是具有挑战性的复杂变形试验，其速度场随时间变化，该试验初始场为

表 7 复杂波试验误差和最大值、最小值
Table 7 Errors and maximum/minimum values of advection of complex waves

pagenumber_ebook=127,pagenumber_book=674
pagenumber_ebook=127,pagenumber_book=674
图 7 凹槽圆柱模拟试验结果（a. 初始时刻，b 和 c 分别为 MCV3_UPCC 旋转 1 周和 10 周后结果，d 和 e 分别为 MCV3_UPCC_BP 旋转 1 周和 10 周后结果）
Fig. 7 The cut-cylinder rotation experiment：（a） the initial distribution of cut-cylinder rotation；（b），（c） result of MCV3_UPCC scheme after 1 and 10 period，respectively；（d），（e） results of MCV3_UPCC scheme with BP filter after 1 and 10 period，respectively

pagenumber_ebook=128,pagenumber_book=675
续图 7
Fig. 7 Continued

表 8 凹槽圆柱模拟试验误差和最大值、最小值
Table 8 Errors and maximum/minimum values of advection of cut-cylinder rotation experiment

pagenumber_ebook=128,pagenumber_book=675
pagenumber_ebook=128,pagenumber_book=675
速度场为

pagenumber_ebook=128,pagenumber_book=675
速度场的定义使得初始场从平滑的状态在t=2.5的时刻演变成高度变形状态，然后速度反向运动，在 t =5时候又回到初始状态，对于一个好的正定保型平流方案，在这种复杂的变形流动过程中，对于保持每一步的正定是具有挑战的。试验区域在[0,1]×[0,1]范围内，结果如图8、9 所示。

pagenumber_ebook=128,pagenumber_book=675
图 8 变形场试验结果 （a. 初始场，b 和c 分别为50×50 和100×100 网格的MCV3_UPCC 方案运动1 个周期的结果）
Fig. 8 （a） The initial field，（b） the result of MCV3_UPCC scheme after 1 period in 50×50 mesh，（c） the result of MCV3_UPCC scheme after 1 period in 100×100 mesh

pagenumber_ebook=129,pagenumber_book=676
图 9 二维变形流模拟试验结果 （a. MCV3_UPCC 方案 t/2 时刻，b. MCV3_UPCC 方案 t 时刻，c. MCV3_UPCC_BP 方案 t/2 时刻，d. MCV3_UPCC_BP 方案 t 时刻）
Fig. 9 The deformational flow experiment：（a），（b） the result of MCV3_UPCC scheme after t/2 and t ，respectively；（c），（d）the result of MCV3_UPCC scheme with BP filter after t/2 and t， respectively

由图8 可知，在网格数为50 和100 的数值试验中 ， MCV3_UPCC 方 案 E2 误 差 分 别 为 0 .0406和0.0102，对比 Skamarock（2006）中 PPM 方法所得误差 0 .0492 和 0 .0140，MCV3_UPCC 方案误差更小，并且在极小值处的维持上也比Skamarock（2006）中-0.073 和-0.039 更好。在保持正定试验中（图9），采用50×50 的试验网格，MCV3_UPCC 方案引用了边界保型限制器后，在准确地模拟强二维变形流试验的同时又能有效消除数值振荡。从表9 可知，时间积分每一步的最大、最小值都限制在初始场的最大、最小值之间，并且与Skamarock（2006）相同试验设置的限制器结果相比，MCV3_UPCC 方法引入边界保型限制器在一个周期后保持的极大值是0.8150，较前者的0.776 极值保持得更好，所以边界保型限制器即使在变形流场中也能严格保持正定和高精度。

表 9 变形流试验误差和最大值、最小值 （按照Skamarock（2006）试验中误差算法）
Table 9 Errors and maximum/minimum values of advection of deformational flow experiment（According to the method of Skamarock（2006））

pagenumber_ebook=130,pagenumber_book=677
4 结 论
基于新的三点均匀中心约束多矩有限体积方法，发展了一个三阶正定守恒的平流模式。与传统有限体积方法相比，新发展的三点均匀中心约束多矩有限体积平流格式（MCV3_UPCC）能够在单个网格内进行局地高阶重构，计算模板紧致。由于采用积分平均值矩积分通量形式控制方程，并保持边界通量连续，因此能够保证数值解严格守恒。数值收敛试验表明，新发展的平流模式具有三阶计算精度。通过分析标准误差，发现MCV3_UPCC 方案相比于MCV3 方案，计算误差减小了约一半，模拟的不连续函数峰值也更接近数值真解，计算精度高。为了有效控制高阶格式产生的数值噪音（振荡），文中引入边界保型限制器，它能够有效抑制非物理数值振荡，具有良好的保型能力。数值试验表明采用边界保型限制器并不降低方案的计算精度，并且能够严格控制每一步的数值解在初始场的最大、最小值之间。综上所述，文中发展的二维平流模式具有良好的发展潜力，下一步计划将新发展的三点均匀中心约束多矩有限体积方案推广到二维球面。

参考文献

沈学顺，王明欢，肖锋等. 2011. GRAPES 模式中高精度正定保形物质平流方案的研究Ⅰ：理论方案设计与理想试验. 气象学报，39（1）：1-15.Shen X S，Wang M H，Xiao F. 2011. A study of the high-order accuracy and positive-definite conformal advection scheme in GRAPES model Ⅰ：Scientific design and idealized tests. Acta Meteor Sinica，39（1）：1-15 （in Chinese）

苏勇，沈学顺，彭新东等. 2013. PRM 标量平流方案在GRAPES 全球预报系统中的应用. 大气科学，37（6）：1309-1325. Su Y，Shen X S，Peng X D，et al. 2013. Application of PRM scalar advection scheme in GRAPES global forecast system. Chinese J Atmos Sci， 37（6）： 1309-1325 （in Chinese）

Bermejo R，Staniforth A. 1992. The conversion of semi-Lagrangian advection schemes to quasi-monotone schemes. Mon Wea Rev， 120（11）： 2622-2632

Blossey P N，Durran D R. 2008. Selective monotonicity preservation in scalar advection. J Comput Phys，227（10）：5160-5183

Book D L，Boris J P，Hain K. 1975. Flux-corrected transport Ⅱ：Generalizations of the method. J Comput Phys，18（3）：248-283

Chen C G，Xiao F. 2008. Shallow water model on cubed-sphere by multimoment finite volume method. J Comput Phys，227（10）：5019-5044

Chen C G，Li X L，Shen X S，et al. 2014. Global shallow water models based on multi-moment constrained finite volume method and three quasiuniform spherical grids. J Comput Phys，271：191-223

Colella P，Woodward P R. 1984. The piecewise parabolic method （PPM）for gas-dynamical simulations. J Comput Phys，54（1）：174-201

Durran D R. 1999. Numerical Methods for Wave Equations in Geophysical Fluid Dynamics. New York：Springer-Verlag，465

Guo W， Nair R D， Qiu J M. 2014. A conservative semi-Lagrangian discontinuous Galerkin scheme on the cubed sphere. Mon Wea Rev，142（1）：457-475

Guo W， Nair R D， Zhong X H. 2016. An efficient WENO limiter for discontinuous Galerkin transport scheme on the cubed sphere. Int J Numer Methods Fluids，81（1）：3-21

Ii S， Xiao F. 2007. CIP/multi-moment finite volume method for Euler equations：A semi-Lagrangian characteristic formulation. J Comput Phys，222（2）：849-871

Ii S， Xiao F. 2009. High order multi-moment constrained finite volume method. Part Ⅰ：Basic formulation. J Comput Phys，228（10）：3669-3707

Jiang G S， Shu C W. 1996. Efficient implementation of weighted ENO schemes. J Comput Phys，126（1）：202-228

Katta K K，Nair R D，Kumar V. 2015. High-order finite-volume transport on the cubed sphere： Comparison between 1D and 2D reconstruction schemes. Mon Wea Rev，143（7）：2937-2954

Lauritzen P H，Nair R D，Ullrich P A. 2009. A conservative semi-Lagrangian multi-tracer transport scheme （CSLAM）on the cubed-sphere grid. J Comput Phys，229（5）：1401-1424

Li X L，Chen C G，Xiao F，et al. 2015. A high-order multi-moment constrained finite-volume global shallow-water model on the Yin-Yang grid. Quart J Roy Meteor Soc，141（691）：2090-2102

Liu X D，Tadmor E. 1998. Third order nonoscillatory central scheme for hyperbolic conservation laws. J Numer Math，79（3）：397-425

Robert A. 1981. A stable numerical integration scheme for the primitive meteorological equations. Atmos Ocean，19（1）：35-46

Shu C W. 1988. Total-variation-diminishing time discretizations. SIAM J Sci Stat Comput，9（6）：1073-1084

Skamarock W C. 2006. Positive-definite and monotonic limiters for unrestricted-time-step transport schemes. J Mon Wea Rev，134（8）：2241-2250

Van Leer B. 1977. Towards the ultimate conservative difference scheme. J Comput Phys，135（2）：229-248

Xiao F，Yabe T. 2001. Completely conservative and oscillationless semilagrangian schemes for advection transportation. J Comput Phys，170（2）：498-522

Xiao F， Yabe T， Peng X， et al. 2002. Conservative and oscillation-less atmospheric transport schemes based on rational functions. J Geophys Res Atmos，107（D22）：ACL 2-1-ACL 2-11

Xiao F. 2004. Unified formulation for compressible and incompressible flows by using multi-integrated moments Ⅰ：One-dimensional inviscid compressible flow. J Comput Phys，195（2）：629-654

Xiao F， Akoh R， Ii S. 2006. Unified formulation for compressible and incompressible flows by using multi-integrated moments Ⅱ：Multidimensional version for compressible and incompressible flows. J Comput Phys，213（1）：31-56

Xiao F. 2012. Two variants of the MCV3 scheme. ArXiv：1207.6844

Xiao F，Ii S，Chen C G，et al. 2013. A note on the general multi-moment constrained flux reconstruction formulation for high order schemes. Appl Math Model，37（7）：5092-5108

Yabe T，T Aoki. 1991. A universal solver for hyperbolic equations by cubicpolynomial interpolation Ⅰ ： One-dimensional solver. Comput Phys Commun，66（2-3）：219-232

Yabe T，Tanaka R，Nakamura T，et al. 2001. An exactly conservative semi-Lagrangian scheme（CIP-CSL）in one dimension. Mon Wea Rev，129（2）：332-344

Yu R C. 1994. A two-step shape-preserving advection scheme. Adv Atmos Sci，11（4）：479-490

Zerroukat M， Wood N， Staniforth A. 2002. SLICE： A semi-Lagrangian inherently conserving and efficient scheme for transport problems. Quart J Roy Meteor Soc，128（586）：2801-2820

Zhang X X，Shu C W. 2010. On maximum-principle-satisfying high order schemes for scalar conservation laws. J Comput Phys，229（9）：3091-3120

SHU Qian1 TANG Jie2,3 CHEN Chungang4 XIAO Feng5 SHEN Xueshun2,3 ZHOU Lilong2,3 LI Zechun3 ZHU Keyun1 LI Xingliang2,3

1. Chengdu University of Information Technology，Chengdu 610225，China
2. Numerical Weather Prediction Center of CMA，Beijing 100081，China
3. National Meteorological Center，Beijing 100081，China
4. Xi'an Jiaotong University，Xi'an 710049，China
5. Department of Mechanical engineering，Tokyo Institute of Technology，Tokyo 152-8850，Japan

Abstract A two dimensional 3rd order positive-definite conservative advection model is developed in this study based on the novel 3-point multi-moment constrained finite-volume scheme for uniform points with center constraints (MCV3_UPCC). In the context of 3-point multi-moment finite volume method, three equidistant solution points are defined within a single cell and the time evolution equations can be obtained by flux form formulation. The multi-moment constraints in this novel scheme are imposed at the cell center on the point value, the first and second order derivatives, and a polynomial of 4th degree can then be reconstructed in a single cell.The resulting MCV3_UPCC scheme has 3rd order accuracy and ensures the exact numerical conservation due to the continuity of the flux function at cell interfaces. To suppress the numerical oscillations and the positivity of certain physical quantities, the boundpreserving limiting projection is introduced into the new MCV3_UPCC scheme, which satisfies the minimum and maximum principle. The new positive-definite conservative advection model is validated by widely used benchmarks. The presented transport model has a good accuracy in comparison with other existing high-order models, and it has the potential to be applied in real moist transport model.

Key words Advection scheme，Multi-moment finite volume method，Positive-definite and conservation，High order scheme

2019-11-22 收稿，2020-03-12 改回.

舒谦，唐杰，陈春刚，肖锋，沈学顺，周立隆，李泽椿，朱克云，李兴良. 2020. 高阶正定守恒的中心约束多矩有限体积平流模式. 气象学报，78（4）：665-678

Shu Qian， Tang Jie， Chen Chungang， Xiao Feng， Shen Xueshun， Zhou Lilong， Li Zechun， Zhu Keyun， Li Xingliang. 2020. A high order positive-definite conservative multi-moment center constrained finite volume transport model. Acta Meteorologica Sinica，78（4）:665-678

* 资助课题：国家重点研发计划项目（2017YFA0603901、2017YFC1501901）、北京市杰出青年科学家基金项目（JQ18001）。

作者简介：舒谦，主要从事数值天气预报研究。E-mail：516189471@qq.com

通信作者：李兴良，从事数值预报模式研究和数值算法应用研发。E-mail：lixliang@cma.gov.cn

中图法分类号 P435