适用于非静力大气模式的完全平衡多矩约束有限体积方法.pdf

适用于非静力大气模式的完全平衡多矩约束有限体积方法*张寅钲1,2陈春刚3沈学顺2,4肖锋5李兴良2,4ZHANGYinzheng1,2CHENChungang3SHENXueshun2,4XIAOFeng5LIXingliang2,41.中国气象科学研究院，北京，1000812.中国气象局地球系统数值预报中心，北京，1000813.西安交通大学，西安，7100494.国家气象中心，北京，1000815.东京工业大学机械工程系，东京，152-88501.Chinese Academy of Meteorological Sciences，Beijing 100081，China2.CMA Earth System Modeling and Prediction Centre，Beijing 100081，China3.Xian Jiaotong University，Xian 710049，China4.National Meteorological Centre，Beijing 100081，China5.Department of Mechanical Engineering，Tokyo Institute of Technology，Tokyo 152-8850，Japan2022-10-27 收稿，2023-02-28 改回.张寅钲，陈春刚，沈学顺，肖锋，李兴良.2023.适用于非静力大气模式的完全平衡多矩约束有限体积方法.气象学报，81（4）：619-629Zhang Yinzheng，Chen Chungang，Shen Xueshun，Xiao Feng，Li Xingliang.2023.A well-balanced multi-moment constrained finitevolume method for nonhydrostatic atmospheric models.Acta Meteorologica Sinica，81（4）:619-629AbstractTheexactdiscretehydrostaticbalancecannotbekeptingeneralwhentheverticalmomentumequationisdiscretizedinatmosphericnumericalmodels.Inordertoexactlybalancethediscreteverticalpressuregradientforceandgravityforce,awell-balancedmulti-momentconstrainedfinitevolumemethodfornonhydrostaticatmosphericflowisdevelopedbyintroducingawell-balanced numerical formulation,in which a numerical reconstruction of the gravity source term is conducted in terms of athermodynamicreferencestatethatsatisfiesthehydrostaticbalance.One-dimensionalbenchmarksshowthatthewell-balancedmulti-momentconstrainedfinitevolumemethodcanmaintainthenumericalerrorofthehydrostaticreferencestateuptomachineroundoffwithacoarsegridspacingandcanwellsimulatethepropagationofsmallperturbationsevenincaseofinitialsmallperturbation.Two-dimensionalnonhydrostaticthermalbubbletestfurtherconfirmsitsabilityofsimulatingnonhydrostaticatmosphericmotionaccurately.Theabovenumericalexperimentshaveverifiedthewell-balancedpropertyandapplicabilityofthewell-balancedmulti-momentconstrainedfinitevolumemethod,whichprovidesagoodreferenceforthedevelopmentofnonhydrostaticatmosphericmodels.Key wordsWell-balanced，Multi-moment constrained finite volume method，Nonhydrostatic atmosphere，Pressure gradientforce，Gravitationalsourceterm摘要大气数值模式离散垂直动量方程时，一般而言气压梯度力和重力不易维持严格的静力平衡关系。为精确平衡数值离散的垂直气压梯度力和重力，基于高精度多矩约束有限体积方法引入完全平衡数值公式，即以满足静力平衡关系的热力学参考态*资助课题：国家自然科学基金项目（42275168）和国家重点研发计划项目（2017YFC1501901）。作者简介：张寅钲，主要从事数值天气预报研究。E-mail：通信作者：李兴良，从事数值算法和数值模式研究。E-mail：doi:10.11676/qxxb2023.20220182气象学报1061014对重力源项进行数值离散构造，发展了适用于非静力大气的完全平衡多矩约束有限体积方法。一维标准数值试验表明，完全平衡多矩约束有限体积方法能在较粗糙的计算网格点上保持静力平衡参考态的数值计算误差在计算机的单精度（）和双精度（）水平，在具有小量级扰动的初始条件下，完全平衡多矩约束有限体积方法能较好地模拟扰动的传播，二维非静力热泡试验进一步验证了完全平衡多矩约束有限体积方法对非静力大气运动的模拟能力。数值试验结果验证了所发展方法的完全平衡属性和适用性，这为非静力大气模式发展提供了良好参考价值。关键词 完全平衡，多矩约束有限体积方法，非静力大气，气压梯度力，重力源项中图法分类号P4351引言实际大气中的热力学状态量存在一种满足静力平衡关系的参考态，与大气非静力效应相关的天气系统可以用叠加在这种静力平衡参考态之上的扰动来描述。随着大气模式向中小尺度、对流尺度和微尺度发展，非静力效应变得十分明显，模式方程中垂直气压梯度力的数值离散与重力源项的不平衡可能会带来不期望的垂直加速度，因此如何使用正确的数值方法精确平衡垂直气压梯度力和重力以及解析非静力扰动的变化传播显得十分重要。含源项的双曲守恒律允许存在稳定态，即不为0 的通量导数与源项精确平衡，物理场量不随时间变化。对大气控制方程而言，热力学量存在的静力平衡参考态便可视为一种稳定态形式，数值方法须在离散计算网格下保持这种稳定态至机器精度。LeVeque（1998）指出维持这种微妙平衡的稳定态与数值方法密切相关，且分步法（FactionalStepMethod）对于维持这种静力平衡稳定态容易导致失败，因为通量导数会导致守恒量的较大变化，再由重力源项去平衡这种变化从而导致较大的误差。高分辨率双曲守恒律有限体积方法（Bale，etal，2003）也被应用到解决维持稳定态问题上。对中尺度非静力大气流动而言，可能会带来物理场上的较大梯度，因而传统大气模式中被用于离散大气控制方程的差分格式（如中央差格式、蛙跳格式）都会带来稳定性问题，Ahmad 等（2007）使用高分辨率的 Godunov型方法求解二维非静力欧拉方程，该方案能够在不引入离散化误差的情况下合并重力源项，以确保能更好地维持静力平衡稳定态。Xing 等（2005，2006，2013）基于含重力源项的浅水方程、弹性波方程、气体动力学方程等双曲模型，基于五阶有限差分WENO（WeightedEssentiallyNon-Oscillatory）方法（Jiang，etal，1996）发展了高精度完全平衡 WENO方法。Ghosh 等（2016）将完全平衡方法引入大气流动中，给出广义化的重力加速度重构表达式，来达到不为 0 的通量导数与重力源项精确平衡的目的。多 矩 有 限 体 积 方 法（Multi-moment finiteVolumemethod）是近年来提出的一种局地高阶重构有限体积方法，该方法是基于多矩概念设计出的高阶多矩有限体积公式而来，是一个通用计算流体力学数值方案框架。在多矩有限体积方法中，它定义了不同的矩，如点值（PointValue，PV）、积分平均值（VolumeIntegratedAverage，VIA）和导数值（DerivativeValue，DV），它们可作为直接预报变量（Xiao，2004；Xiao，etal，2006）或者作为约束条件（Ii，etal，2009；Xiao，etal，2013）来推导演变方程，目前在大气数值模式中获得了广泛应用（Li，etal，2008，2013；Chen，etal，2008；Qin，etal，2019；舒谦等，2020）。文中采用的是多矩约束有限体积方法，即矩作为约束条件来导出求解点的演化方程，具体来说，在构造各个物理量以及其导数的空间数值近似时，多矩约束有限体积方法基于单个计算网格单元的矩给出约束条件，然后根据数值计算精度需要局地高阶地构造出高精度多矩约束拉格朗日（MC-Lagrange）插值重构函数。相比于在传统有限体积方法里为了获得高精度的数值近似需要跨越多个网格单元，新的多矩约束有限体积方法计算模板紧致、更加简便和高效。采用完全平衡方法与多矩约束有限体积方法相结合的研究思路，试图通过多矩算法离散垂向气压梯度力和重力源项，以避免不期望的垂直加速度，从而在满足静力参考态初始条件下，使得计算网格点上的计算误差保持在机器精度，以严格保持静力平衡关系。2完全平衡多矩约束有限体积方法基于一维非静力大气控制方程，采用 Ii 等（2009）提出的三阶多矩约束有限体积方法（3-orderMulti-620Acta Meteorologica Sinica气象学报2023，81（4）momentConstrainedfinite-Volumescheme，MCV3）和完全平衡（well-balanced）方法的构建思路，发展完全平衡多矩约束有限体积方法。考虑以重力为源项的一维欧拉方程qt+f(q)z=s(q,z)（1）q=(,w,e)Tf(q)=w,w2+p,w(e+p)Ts(q,z)=(0,g,wg)Tge式中，预报变量为包含大气流体质量、动量和能量的守恒量，通量，重力源项，上标 T表示转置，为大气密度（kg/m3），w 为垂直速度（m/s），p 为大气压（Pa），为重力加速度（m/s2），z 为高度（m），为包含大气流体内能与动能的能量（J），表达为e=p1+12w2（2）=cp/cv式中，其中 cp为大气定压比热（J/(kgK)），cv为大气定容比热（J/(kgK)）。2.1 静力平衡参考态g对大气而言，重力是一种静态的力场，即重力加速度 为常数。热力学参考态满足的静力平衡关系描述为垂直气压梯度力与重力的平衡，表达为dpdz=g（3）考虑满足上述静力平衡关系式（3）的一般参考态表达形式w=0=0(z)p=p0(z)（4）、p0=0RT0式中，为高度 z 的标量函数，下标 0 表示位于参考高度处的值，满足状态方程，T 为温度，R 为气体常数。文中大气流动标准数值试验包含如下两种静力平衡参考态类型。T=T0常数温度静力平衡参考态。假设温度，结合式（3）和（4）可得到常数温度静力平衡参考态类型的标量函数形式为(z)=exp(gzRT)(z)=exp(gzRT)（5）=T0=(pp0)1常数位温静力平衡参考态。假设，引入无量纲气压，由此静力平衡关系可表示为R1ddz=g（6）对式（6）积分后为=1(1)gzR（7）式（7）回代入无量纲气压定义式中，结合式（4）得到常数位温静力平衡参考态类型的标量函数形式为(z)=1(1)gzR11(z)=1(1)gzR1（8）2.2 完全平衡公式式（1）中，动量方程的通量导数与重力源项的平衡关系允许存在如式（4）中描述的一般参考态，常规数值方法的离散误差不易维持静力参考态或淹没叠加在静力参考态上扰动的模拟传播，因而选择合适的数值方法保持静力参考态至机器精度是至关重要的。w=0、在满足静力参考态式（4）下，由于，式（1）中质量方程和能量方程满足精确平衡关系，动量方程转换为垂直气压梯度与重力源项的平衡关系。静力平衡关系式（3）结合一般静力参考态式（4），得到标量函数满足p0(z)=0(z)g（9）p0=0RT0式中，（z）为（z）对 z 的导数，结合参考高度处的状态方程，得到重力加速度表达式g=RT0(z)1(z)（10）将重新构造的重力加速度代入式（1）qt+f(q)z=s(q,z)（11）q=(,w,e)Tf(q)=w,w2+p,w(e+p)Ts(q,z)=0,RT0(z)1(z)z,wRT0(z)1(z)zT式 中，守 恒 量，通 量，重 新 构 造 后 的 重 力 源 项。D1D2完全平衡方法将不为 0 的通量导数与重力源项在数值计算网格上离散后仍能与网格分辨率z 趋于 0 时达到同样的精确平衡效果。如式（1），数值离散是对控制方程中的通量导数项进行离散的数值重构逼近。一般而言，微分方程中空间导数在数值离散后存在截断误差。在满足静力平衡参考态下，对动量方程考虑通量的线性离散算子和重力源项的线性离散算子，数值离散的截断误差张寅钲等：适用于非静力大气模式的完全平衡多矩约束有限体积方法621ET（）可表示为ET=D1(w2+p)RT0(z)1D2(z)（12）在热力学量满足静力平衡参考态式（4）下ET=D1p0(z)0(z)RT0(z)1D2(z)（13）D(ch)=cD(h)ch一般而言，线性离散算子 D 满足，其中 为常数，为任意格点函数。从而截断误差可表示为ET=p0D1(z)0RT0D2(z)（14）p0=0RT0D1=D2=D0结 合 参 考 高 度 状 态 方 程，若 满 足，则动量方程的截断误差为ET=D1(w2+p)0RT0(z)1D2(z)=p0D(z)p0D(z)=0（15）(z)zf(q)z由此，若重新构造后的重力加速度中导数项采用与通量导数相同的数值离散方法，二者的截断误差算子 D 也应保持一致，从而式（11）在满足静力平衡参考态式（4）的初始条件下截断误差为 0，即可保证数值离散后垂直气压梯度力与重力源项的完全平衡。2.3 完全平衡多矩约束有限体积离散通常情况下，多矩约束有限体积方法以单元内的点值作为预报变量、单元积分平均作为约束条件以严格保证数值守恒，本节采用三阶多矩约束有限体积方法（Ii，etal，2009），发展完全平衡多矩约束有限体积方法。Ciill=1,2,3i将计算域划分为非重叠、边界点共用且连续的单元，对任意场量，单个单元内定义点值矩，如图 1 所示，其中，以及单元积分平均值矩i(t)=1ziwzii(z,t)dz（16）kzil(t)=k(zil,t)zkk=0,1（17）zi=zi+1/2zi1/2zilk=00zil=ilk=11zil=zil式中，为单元的大小，为单元两端或者内部的约束点。对式（17），时，表示定义在单元内的点值；时，表示一阶导数。单元内定义的点值采用等距分布，约束点位于单元中心和两端。单元内构造多矩约束拉格朗日（MC-Lagrange）插值函数i(z)=3l=1clil（18）=q=Q=f=F=Mcl当时，；当时，；当时，。为拉格朗日插值基函数cl=3p=1,p,lzzpzlzp（19）对式（11），将预报变量的约束条件施加在单元积分平均以及单元两端界面处的点值qi=1ziwziQi(z)dz=qi1+4qi2+qi36（20）qi1/2=Qi(zi1/2)=qi1（21）qi+1/2=Qi(zi+1/2)=qi3（22）Ci单元内各约束矩的演化方程为ddtqi=1zi(fi+1/2fi1/2)+si（23）ddtqi1/2=fzi1/2+si1/2（24）ddtqi+1/2=fzi+1/2+si+1/2（25）从而导出各点矩或预报量的演化方程dqi1dt=fzi1/2+si1/2（26）dqi2dt=32si1zi(fi+1/2fi1/2)+14(fzi1/2+fzi+1/2 si1/2 si+1/2)（27）dqi3dt=fzi+1/2+si+1/2（28）si单元积分平均中导数项通过与通量相同方式表达si=0，iRT01i(i3i1)zi，wiRT01i(i3i1)ziT（29）Zi 1/2Zi+1/2i1i2Cii3Ciil图1单元中点值（求解点）分布ilCiFig.1Thelocationsof（unknowns）inone-dimensionalcell622Acta Meteorologica Sinica气象学报2023，81（4）iizi1/2式中，与通过单元积分平均表达式给出。位于单元界面处的数值通量导数以及重力源项的数值导数部分可通过导数黎曼问题（DerivativeRiemannproblem）得到。以单元的界面为例fzi1/2=DRiemann(fzi1/2，f+zi1/2)（30）si1/2=DRiemann(si1/2，s+i1/2)（31）单 元 界 面 处 导 数 黎 曼 问 题 采 用 Local Lax-Friedrich（LLF）近似黎曼求解器fzi1/2=12f+zi1/2+fzi1/2?i1/2?max(q+zi1/2qzi1/2)（32）i1/2A=f/q式中，为单元界面处通量雅可比矩阵的特征值向量，其中下标“max”表示特征向量中的最大特征值。值得注意的是，重力是一种体力，作为源项它不应该用迎风格式来计算，故为了避免引入非物理的数值粘性，单元界面处重力数值计算采取算术平均给定 si1/2=12(s+i1/2+si1/2)（33）iqfCizi1/2单元界面处守恒量、通量以及重力源项的左值和右值则通过对应单元的多矩约束拉格朗日插值多项式（z）给定。以任意场量 为例（可以是守恒量、通量 或者重力源项包含的（z），考虑单元的界面，其界面左右点值以及导数值表达为i1/2=(i1)3（34）zi1/2=(i1)14(i1)2+3(i1)3zi1（35）+i1/2=i1（36）+zi1/2=3i1+4i2i3zi（37）qq、qq=(z)1，w(z)1，eTe值得注意的是，采用守恒量 计算单元界面处导数黎曼问题时会破坏静力参考态，因此需将守恒量通过标量函数修正为，其中，修正后的能量 表达为e=p(z)11+12(z)1w2（38）通过修正之后，式（32）变换为下式计算fzi1/2=12f+zi1/2+fzi1/2?ip?max(q+zi1/2qzi1/2)（39）当大气流体满足静力平衡参考态时，数值离散q时的不会破坏完全平衡属性，当大气流体中存在小扰动的模拟传播或有较强梯度的系统时，也能够保证足够的粘性。2.4 二维扩展本节将完全平衡多矩约束有限体积方法扩展至以重力为源项的二维笛卡尔坐标系下守恒型欧拉方程qt+F(q)x+G(q)z=s(q,z)（40）q=(,u,w,e)TF(q)=u,u2+p,uw,u(e+p)TG(q)=w,wu,w2+p,w(e+p)Ts(q,z)=(0,0,g,wg)Te=p1+12(u2+w2)式中，预报变量为包含大气流体质量、动量和能量 的 守 恒 量，水 平 方 向 通 量 项，垂直方向通量项，源 项，其中 u、w 分别为 x、z 方向的速度，能量包含流体内能与二维动能，重力方向为 z 轴负方向。g通过重新构造的重力加速度式（10）替换重力加速度 后的二维重力源项为s(q,z)=0,0,RT0(z)1(z)z，wRT0(z)1(z)zT（41）代入式（40）得到如下形式qt+F(q)x+G(q)z=s(q,z)（42）在笛卡尔均匀网格坐标系下，可直接将一维离散方法通过对应坐标方向的空间导数近似推广至二维（Li，etal，2013）qijlmt=q(x)ijlmt+q(z)ijlmt（43）CijCijq(x)ijlmq(z)ijlm式中，i、j 分别为 x、z 方向上的单元网格索引，l、m 分别为单元中的 x、z 方向的点值索引，l、m 取值范围为 1，2，3，单元内的点值矩分布如图 2 所示。式（43）中，重力源项作用于垂直方向 z，从而和的时间演变方程可表示为q(x)ijlmt=Fijlmx=3m=13l=1Fxijlm（44）q(z)ijlmt=Gijlmz+sijlm=3l=13m=1(Gzijlm+sijlm)（45）张寅钲等：适用于非静力大气模式的完全平衡多矩约束有限体积方法623x方向（z方向）水平方向（）式（44）中通量导数采用三阶多矩约束有限体积方法数值离散，垂直方向式（45）中通量导数与重构重力源项采用完全平衡三阶多矩约束有限体积方法数值离散。具体而言，水平方向的离散格式为3m=1Fxij1m=3m=1Fxij1m（46）3m=1Fxij2m=3m=114(Fxij1m+Fxij3m)3(Fij3mFij1m)2x（47）3m=1Fxij3m=3m=1Fxij3m（48）垂直方向的离散格式为3l=1Gzijl1+sijl1=3l=1Gzijl1+sijl1（49）3l=1Gzijl2+sijl2=3l=114(Gzijl1+Gzijl3 sijl1+sijl3)+32(Gijl3Gijl1z+sijlm)（50）3l=1Gzijl3+sijl3=3l=1Gzijl3+sijl3（51）q同样地，用于导数黎曼问题的静力平衡参考态修正的二维表达为q=(z)1,u(z)1,w(z)1，eT（52）e式中，修正能量 为包含修正的内能和二维速度矢量（u，w）的动能形式e=p(z)11+12(z)1(u2+w2)（53）一般来说，设计非静力大气数值模式时需要兼顾算法的计算效率与计算精度。与传统有限体积方法相比，多矩有限体积方法在进行数值重构时不仅可以采用传统有限体积的积分平均值，而且还可以利用点值及其导数值，属于局地高阶重构一类算法，具有类似谱方法的高精度性质（Ii，etal，2009）。由于其局地重构计算模板紧致，计算强度大，它非常适用于当前及未来众核高性能计算，具有良好的并行可扩展性（Zhang，etal，2017）。文中发展的完全平衡多矩约束有限体积方法仅需在初始时刻计算一次重构重力源项的导数，它对算法计算效率的影响可忽略不计。2.5 时间积分qt=R(q)R(q)通过上述构建的完全平衡多矩方法的数值离散后，方程组写为，其中为空间离散算子。采用三阶 TVD（TotalVariationDiminishing）龙格库塔（Runge-Kutta）时间积分方案（Shu，1988）求解q(1)=q(n)+tRq(n)（54）q(2)=34q(n)+14q(1)+tRq(1)（55）q(n+1)=13q(n)+23q(2)+tRq(2)（56）q(n)nq(1)q(2)q（n+1）n+1式中，为 时刻的守恒量，、为龙格库塔每积分子步的守恒量，为时刻的守恒量。3数值试验与结果 3.1 一维等温平衡解=1.4 g=1T=T0=1Rp0=1 Pa,0=1 kg/m3一维等温平衡解试验用于验证完全平衡多矩约束有限体积方法对保持静力参考态至机器精度以及模拟扰动传播的能力，首次由 LeVeque（1998）提出。计算域设置为0，1m，考虑常数温度静力平衡参考态形式式（5），取，m/s2，参考高度处。完全平衡多矩约束有限体积方法将保持如下形式的常数温度静力平衡参考态0(z)=exp(z)p0(z)=exp(z)w0(z)=0（57）qij13qij12qij11qij21qij31mqij32qij22Cijlqij23qij33zj+1/2zj 1/2xi+1/2xi 1/2xzCijqijlm图2二维单元中点值（求解点）分布qijlmCijFig.2Thelocationsof（unknowns）intwo-dimensionalcell624Acta Meteorologica Sinica气象学报2023，81（4）3.1.1完全平衡测试t=2 sl1=imaxi=1|q0iqni|/imax首先是发展的完全平衡多矩约束有限体积方法的完全平衡属性测试。初始态为如式（57）的等温静力平衡参考态。分别采用 40 和 200 个单元进行计算，积分总时间。为验证静力平衡参考态的计算误差是否能保持到机器精度，分别使用了4 字节单精度浮点数和 8 字节双精度浮点数进行计算，首先采用非平衡的多矩约束有限体积离散格式测试垂直压力梯度和重力的平衡性，此时重力加速度保持常值 g=1，守恒量、w、e 的 l1误差（）计算结果如表 1 所示，保持静力平衡参考态数值试验的 l1误差表明，不论 40 网格数或是 200 网格数，非平衡的多矩约束有限体积离散格式对保持静力参考态的离散数值误差均未达到机器单、双精度的舍入误差，因此采用完全平衡的数值离散显得十分必要。采用完全平衡的多矩约束数值离散格式计算结果的 l1误差如表 2 所示，从误差结果可以看出，在不同的计算机器精度下以及不同计算网格分辨率下，完全平衡多矩约束有限体积方法都保持到对应的机器精度舍入误差，这验证了所发展方法具有的完全平衡属性。3.1.2等温平衡解扰动测试p0(z)本试验用于测试完全平衡多矩约束有限体积方法对于非静力扰动的精确模拟能力。对静力平衡参考态式（57），在初始气压上施加一个扰动p(z,t=0)=p0(z)+exp100(z0.5)2（58）=104=106式中，是一个非 0 常数，用于控制施加在初始气压上的扰动量级大小，试验中分别取和两种情况。=104=106初始扰动将分裂为两道扰动波，分别向z轴正方向和负方向传播。试验使用 100 个单元网格数和无通量边界条件（No-fluxboundaryconditions）计算，积分时间 t=0.25s，数值计算结果如图 3 所示，气压的初始扰动用灰色虚线表示，红色圆圈为完全平衡三阶多矩约束有限体积方法的计算结果，蓝色空心三角形为未结合完全平衡的三阶多矩约束有限体积方法的计算结果。为了进行比较，同时绘制了使用 2000 个单元的完全平衡多矩约束有限体积方法计算结果，以黑色实线表示。可以看出，在较大的 误差量级（）下，完全平衡多矩约束有限体积方法与未采用完全平衡处理的多矩约束有限体积方法计算结果基本趋于一致，当 量级较小（）时，非平衡的多矩约束有限体积方法扰动模拟结果出现较为明显的偏差，然而完全平衡多矩约束有限体积方法的模拟结果合理，这表明完全平衡多矩约束有限体积方法既使在较小量级的以及相对粗糙的网格分辨率z 下也能准确模拟扰动的传播。同时进一步的试验表明，在满足一定计算机精度范围内的扰动量级下，所发展的完全平衡多矩约束有限体积方法同样具有扰动传播的准确捕获和模拟能力。3.2 重力场下 Sod 激波管g=1=1.4标准激波管试验（Sod，1978）用于测试构建的数值算法对强梯度物理场的模拟能力，同时结合重力的影响。计算域范围设置为 0，1m，重力加速度取m/s2，初始化条件如下=1，w=0，p=1z0.5=0.125，w=0，p=0.1z0.5（59）使用完全平衡多矩约束有限体积方法，导数黎曼问题计算选择守恒量 q，采用反射边界条件以及200 个单元进行计算，为抑制高阶数值方法的非物理振荡，应用了 TVB 限制器（Shu，1987；Ii，etal，2009）。积分时间 t=0.2s，密度（）、速度（w）、气压表1非平衡的多矩约束有限体积方法，相对初始态的l1误差Table1Withoutthewellbalancedmulti-momentconstrainedfinite-volumemethod，thel1errorwithrespecttoinitialsolution网格数浮点数精度（kg/m3）w（kg/(m2s)）e（J）40单精度1.071043.831053.30104双精度2.041065.191067.15106200单精度3.901041.481059.26104双精度7.511081.771072.63107表2完全平衡的多矩约束有限体积方法，相对初始态的l1误差Table2Withthewell-balancedmulti-momentconstrainedfinite-volumemethod，thel1errorwithrespecttoinitialsolution网格数浮点数精度（kg/m3）w（kg/(m2s)）e（J）40单精度4.521063.921061.43106双精度6.7610155.2810151.361016200单精度3.781053.651051.40105双精度1.7010142.2810146.141015张寅钲等：适用于非静力大气模式的完全平衡多矩约束有限体积方法625（p）以及能量（e）的单元积分平均值分布结果（红色空心方块）如图 4 所示。由于存在重力场的影响，模拟的流体密度、气压、能量都存在向 z 轴负方向堆积的趋势，并伴随00.500.25Pressure perturbation(Pa)00.250.500.751.001.251.50104Initial stateCell=2000Cell=1000.20.4z0.60.81.000.500.25Pressure perturbation(Pa)00.250.500.751.001.251.50106Initial stateCell=2000Cell=1000.20.4z0.60.81.0MCV3MCV3图3等温平衡解扰动测试Fig.3Perturbationtestoftheequilibriumsolution000.20.4Density(kg/m3)0.60.81.01.20.20.40.60.81.0Cell=200Cell=200000.40.20Velocity(m/s)0.20.60.40.81.00.20.40.60.81.0Cell=200Cell=20000(c)(a)(d)(b)00.20.4Pressure(Pa)0.60.81.01.21.40.20.4z0.60.81.0Cell=200Cell=2000000.51.0Energy(J)1.52.03.02.53.50.20.4z0.60.81.0Cell=200Cell=2000图4重力场下 Sod 激波管（t=0.2s；a.密度，b.速度，c.气压，d.能量）Fig.4Shocktubetestundergravitationalfield（t=0.2s；a.density，b.velocity，c.pressure，d.energy）626Acta Meteorologica Sinica气象学报2023，81（4）有负速度的出现。同时绘制了使用 2000 个单元计算的参考解结果，以黑色实线表示，可以看出在相对粗糙网格下的结果（红色空心方块），完全平衡多矩约束有限体积方法的计算结果与参考解一致。3.3 二维上升热泡=T0=300 Kp0=1105Pag=9.81 m/s2二维上升热泡是中性大气中暖气泡模拟的一个基准算例，被广泛用于验证数值算法对热力学效应引起大气运动的模拟能力（Li，etal，2013；Ghosh，etal，2016）。将二维计算域设置为0，1000m0，1000m，静力平衡参考态为恒定位势温度类型，取（位于参考高度 z=0 处的温度T0），参考高度处的气压，重力加速度，初始场水平以及垂直速度均为 0，物理场处于静止状态。二维热泡通过施加在位势温度上的扰动（）描述(x,z,t=0)=c21+cos(crrc)rrc0rrc（60）r=(xxc)2+(zzc)2ccrcx=z=2.5 m式中，为扰动量级，为圆周率常数，为热泡的半径大小，取 rc=250m，（xc，zc）=（500m，350m）是热泡的初始中心位置，四周边界均采用无通量边界条件，水平和垂直方向分别取 400 个单元，相当于单元大小，时间步长 dt=0.00125s，积分到最终时间 t=700s。图 5 为位温扰动于 0（初始时刻）、150、300、450、600s 以及最终积分时刻 700s 的结果。二维02004006008001000(a)z(m)00.10.20.30.40.502004006008001000(b)z(m)00.10.20.30.40.5002004006008001000(d)z(m)200400600 x(m)800100000.10.20.30.40.5002004006008001000(c)z(m)200400600 x(m)800100000.10.20.30.40.5KKKK图5二维上升热泡不同时刻位温扰动分布（a.0s，b.150s，c.300s，d.450s，e.600s，f.700s）Fig.5Risingthermalbubble:potentialtemperatureperturbationcontoursatdifferenttimesteps（a.0s，b.150s，c.300s，d.450s，e.600s，f.700s）0s，f.700s）张寅钲等：适用于非静力大气模式的完全平衡多矩约束有限体积方法627热泡受到浮力的作用上升，上升过程中由于位温扰动梯度的存在，热泡中心更热从而上升速度更快，导致热泡上层出现较强的位温扰动梯度。积分最终时刻的速度场分布如图 6 所示，由于水平和垂直速度的产生，热泡被剪切从而形成蘑菇云形状。由于计算域上边界存在一定的阻隔作用，二维上升的热泡与其相互作用形成一层隔层，从而导致热泡进一步卷曲以及向两侧扩展。图 5af 各时刻的位温扰动计算结果分布与 Ghosh 等（2016）研究结果保持一致。4结论发展了适用于非静力大气运动的完全平衡多矩约束有限体积方法。具体来说，首先引入静力平衡参考态，重构了控制方程中重力加速度；进而基于高精度多矩约束有限体积方法，对通量导数和重构的重力源项导数项采用相同的空间数值离散，实现垂直气压梯度力与重力精确平衡，静力平衡参考002004006008001000(e)z(m)200400600 x(m)800100000.10.20.30.40.5002004006008001000(f)z(m)200400600 x(m)800100000.10.20.30.40.5KK续图5Fig.5Continued002004006008001000(a)(b)z(m)200400 x(m)60080010001.81.41.00.60.20.20.61.01.41.8m/s002004006008001000z(m)200400 x(m)60080010001.81.41.00.60.20.20.61.01.41.8m/s图6二维上升热泡最终时刻 t=700s 速度场分布（a.水平速度（u），b.垂直速度（w）Fig.6Risingthermalbubble:windspeedcontoursatthefinaltimestept=700s（a.horizontalwind（u），b.verticalwind（w）628Acta Meteorologica Sinica气象学报2023，81（4）态的数值误差达到机器精度。一维数值试验结果表明，即使在较粗的计算网格上，完全平衡多矩约束有限体积方法模拟的静力平衡参考态数值计算误差能保持至机器精度舍入误差（单精度舍入误差为 106，双精度舍入误差为 1014），对非静力扰动具有良好的模拟能力；激波管试验表明完全平衡多矩约束有限体