大气化学模式的数值计算方法

环境科学研究第13卷第1期esearch of environmental sciencesVol.13No.12000大气化学模式的数值计算方法刘峻峰,李金龙,白郁华北京大学技术物理系北京100871)摘要∶大气化学模式计算的一个重要研究热点是如何快速、精确地对相互关联的非线性常微分方程组(OEs进行数值积分求解。由于大气中督物种寿命长短相差很大,导致ωDEs具有了很强的刚性因而在数值求算过程中存在着计算精度和计算效率的平衡问题。笔者介绍了当前几种常见的大气化学模式的计算方法包括 GEAR QSSA( MQSSA)Y&BMY&B),RADM, LSODE SMVGEAR小参数法 Gong and Cho及 TWOSTEP等并结合当前本领域内的部分研究成果对一些重要的方法进行了比较和分析。关键词:大气化学模式;刚性常微分方程组;数值积分中图分类号:X11文献标识码:A文章编号:1001-692×200001-0044-06The Numerical Solvers for Atmospheric Chemistry ModelsLIU Jun-feng, LI Jin-long, BAI Yu-huaDepartment of Technical Physics Peking University Beijing 100871)Abstract: The numerical integration of nonlinear systems of ordinary differential equations( ODEs )is a major problem to atmosphechemistry models. Because of the great disparities in the rates of the atmospheric chemistry reactions it is importantd a properway to accurately and efficiently solve the stiff systems of ODEs. In this paper several solvers including GEAR QSSA( MQSSA)Y&E MY&B) RADM SODE SMVGEAR SPM ong and Cho and TWOSTEP aare presented and compared the based on current research achievements in this fieldKey words: atmospheric chemistry model stiff ODE numerical integration自70年代以来,模式研究已成为大气环境化学领域中或(1)代表了一组数量较大、非线性且互相关联的常微不可缺少的组成部分。一—个性能优异的大气化学模式应具分方程细ODEs)故无法求出其解析解,只能在给定每种化备2个必要条件⑩合理的模式机理。此机理能较准确地反合物初值后用各物种浓度对时间数值积分的方法,求岀各映实际大气中物种数量及彼此间的化学反应同时给出较准物种在特定时刻的浓度。确的动力学数据。②高效精确的计算方法。该方法在满足最早的方法是对具有初值的式(1)在固定时间步长的条定误差范围的前提下应尽可能快速地计算出最终结果件下进行 Euler积分。但由于各反应物种在大气中的寿命长根据研究的对象及目的不同,人们设计出了不同的机短不一,物种间的反应速率相差高达几个数量级,造成了理。常见的几种大气化学模式机理及计算方法见表1。ODEs的刚性。为了保证求解方程组的准确性和稳定性由表1可以看出,般机理包含的化合物种类由几十种须将时间步长定得很小造成计算效率极低。因此在数值积到上百种反应方程式由几百个到上千个。由这些化合物及分过程中存在着积分时间步长同各物种的寿命是否相匹配反应方程式根据质量作用定律可以写出每种化合物的反的问题。若积分时间步长取的较大,对于寿命短的化合物会应速率方程产生较大误差时间步长取得过小,又会使计算过慢造成dh=F(C)=P(Ct)-L(CC(1)对计算资源的浪费中国煤化工中新数值求解OEs的方法式中:P(Cd)第i种化合物的产生速率HCNMHG应用,但是周看人们研究物L(C,)第i种化合物的去除速率。质量模式中所取的网络数成干上万,一般又要模拟较长收稿日期:1999-09-14段时间的反应因而对计算速度的要求也在提高,即需要作者简介刘峻峰(1974-)男北京市人硕士研究生种快速且较为精确的计算方法。GEAR法逐渐成为人们对新研究的快速方法精确度的检验手段。1977年Y&B法161第1期刘峻峰等大气化学模式的数值计算方法45和1978年QSsA"拟稳态假设)法的提出标志着快速计算较18-201,并提出了一批新的或是改进的方法,如法的诞生。此后直到现在又有许多科学家对大气化学模式 LSODEL211,RADM21, Gong andTWOSTEP 241中刚性ODEs数值解法的效率和精度进行了研究和比 SMVGEAR25] MQSSA和MY&B201表1几种常见大气化学模式及其计算方法Table 1 Atmospheric chemistry models and their numerical solvers机理缩写名参考物种数反应方程数计算法RADM2-IFUStockwell et al.[(1990)RADM2-FZKStockwell et al. [( 1990)RADMRADM2-KFAStockwell et al. [I(1990)Implicit EulaEuro-RADMStockwell and Kley?( 1994)CBM-ⅣLOTOGery et al. [3( 1989 )Builties+( 1992)63CBM-ⅣCBM-N-TNOGery, et al. 3(1989 Sterr'( 1994)GearC4.1Gery et al. [3(1989EMEPEMEPSimpson et al. [7( 1993), Simpson&( 1995140MCMPilling et al. [(1998)>2500>7000ADOM-ⅡlADOMLurmann et al. 10( 1986)47114Strand and HoyI( 1994)Andersson-Skold 12(1995RuhnkeRuhnke 13( 1995)145342资料主要来源 Kuhn et al.41998笔者在前面简要介绍了GEAR, LSODE, SMVGEARhF(Cm)Cm-1为前一步的C值,CQSSA, MQSSA,Y&B,MY&B,小参数法, RADM Gong and为本步的C值Cho及 TWOSTEP等计算方法。在下面结合各方面的研究对这些方法进行了分类并对重要方法的精确性和计算效率n——计算步数做了比较和分析每步校正次数1方法介绍积分开始时,Z0=(CC′o),C0为初值,C'o=dC/d1.1GEAR法和改进GEAR法(C)第一步步长取值较小P取1以后根据计算的误GEAR法及与GEAR法类似的 LSODE, LSODES,差大小可按一定的方法自动调整h和P。SMVGEAR等代表了一类较精确的数值计算方1.1.2 LSODE与 LSODES1.1.1GEAR法LSODE( Livermore solver for ordinary differential equaGEAR数值计算法151具有适用范围广、可以自动调节ons521是建立在GEAR法基础上,被广为应用计算ODEs时间步长与误差,以及计算结果准确度高等特点,是一种的较精确方法。此方法允许使用者根据ODEs的刚性程度为理想的计算方法。但其最大的缺点是计算步骤繁琐要频选择Adam多步计算公式(对于非刚性体系)或BDF公式繁计算 Jacobian行列式,当计算的物种比较多时会耗费大 backward difference multistep formulas对于刚性体系)大量的机时和内存。GEAR法的数值计算公式如下气化学模式的ODEs体系刚性较大故多用BDF法。BD)F的a.预测公式计算公式为(2)C)+△tbf(Cb.校正公式式中:C"+1—当前t"1的浓度Zm+1=Zm+1-[(l1-lnh)(xZm+1)](3)中国煤化工C.最终公式(4)CNMHG邇以刊酊炬阵结构稀少的不同来源(是式中zn—(k)×(P+1矩阵来自于用户直接输入,还是来自于调用 Jacobian子程序或是A-(P+1)×(P+1)矩阵l≤P≤6多次调用方程评估子程序)分别采用不同方法简化计算步Ⅰ——单位矩阵骤提高计算效率26。(xm})k×1矩阵,G(Zm+1)=AF(Cm1.1. 3 SMVGEAR46环境科学研究第13卷SMVGEAR( Sparse-Matrix Vectorized Gear Code25是另(△tm≤Mt≤△tmx)(9)种GEAR式的解ODEs的方法。它与 LSODE有一定的相式中:最大允许△t的改变量取值为1%~100%。似性但也有其特点1.3Y&B和MY&B其最大的特点是利用了矢量式计算机在计算内循环时1.3.1Y&B可以显著提高浮点运算效率这一现象设计计算代码的。矢Y&B数值计算法又称 Hybrid法l6,由 Young和Bori量化是指计算机在计算内部DO…LOOP循环时,可以提速于1977年提出其目的是建立一种先预测再迭代校正的积到其最大速度的90%以上。当然若循环中包含一些特殊分方法。此法后被 McRae odman等采用并改进。 Yamani作如递归、矩阵运算时)矢量化也会减慢因而此法也减1o等将其应用于 CALGRID大尺度模式计算中271少了不必要的矩阵操作以提高速度Y&B法开始时按反应物的特征时间与积分时间步长的对于化学模式计算中的ODEs可以认为是二维的即关系将反应物分成2类:①刚性物质(τ<△t),其在大气中维是化学物种数,一维是网格数。 Jacobson等研究认为油于的寿命较短非刚性物质(r;>△t)在大气中的寿命较长般情况下网格数会远远大于物种数故以网格数矢量化的公式如下速度比以化学物种数矢量化的速度快。计算每步的阶和时a.刚性(x;≥△tn)间步长的方法同于 LSODE。减少矩阵运算的措施是引入预测值JSPARSE法。此法先对物种进行排序,将出现频率较少、和其他物种关系简单的物种排在矩阵的顶部而出现频率C:= Ci+(P(10)多、关系复杂的物种排在矩阵的底部这样处理可以最大限修正值度地提高矩阵和零相乘的几率洏而后 JSPARSE将所有零和矩阵相乘的项删除。这是其相对于 LSODE法的主要差别。(11)当然有一点SMVGEAR运行在非矢量化的计算b.非刚性(r<△tn)机上时和 LSODE法在性能上是相差不大的预测值1.2QSsA和 MQSSAC(2r"-4tn)+2P?r"△1.2.1 QSSAQSSA数值计算法17又称拟稳态近似法是 Hesstwedu修正值等于1978年提出的。其主要思想是按照积分时间步长△t△t(P"+P;r"+r;y2与寿命τ的比将计算公式分成3种r+r;-△tna.Δt/τ>10,此时化合物的寿命很短按稳态假设处(13)理即将P;和L;看成是常数1.3.2MY&BPMathur20等人指出倘若对n的物种应用稳态假(6)设则会大大提高计算效率。他们提出了将反应物分成3类b.0.01<△t/r<10,此时化合物的寿命与所取步长的方法即快速物质τ;<0.2△n中速物质0.2△tn≤r;接近5△tn和慢速物质r2>5△tn。对于中速和慢速物质应用公式(7)(10~13)快速物质则应用公式(6) Mathur等同样将自动c,△t/x<0.01此时化合物的寿命较长时间步长控制机制引入到了Y&B法见式(9+1=CLn+△(P-LC)(8)14小参数法(SPM)公式(7是建立在P和L;是常数的基础上的公式6厢(8)小参数法由我国韩天敏提出281适用于中等刚性的是为了增加计算速度。当△t/r>10时用公式(6算所引ODEs。计算公式为起的误差小于公式(7)所引起误差的4.5×10-3%而公式M= p[ef(tn+I Ch+1)+Ch+I-C] (14)(8X简单 Euler积分公式)虽引入了误差项(△t/r;)/2但引+1=Cn+(1-a)xn+ax+1(15)起的误差小于公式(7)所引起误差的5.0×103%。式中1.2. 2 MQSSAp=hh+∈)Mh2等引入矫正因子P和L:P=-2,中国煤化工L;=L;。每次计算出的结果用上式校正后将P和L;代入1.5CNMHG公式(6~8)中继续迭代RADM法22是针对RADM( Regional acid depositon mod此外 Mathur等还引入了时间步长控制机制:el)模式提出。其对18种自由基作了稳态假设,与QSSA法类似对时间步长也规定了上下界。计算公式如下a.OH自由基浓度由解下式方程得到第1期刘峻峰等大气化学模式的数值计算方法47al oh]+ b[ OH]+c=16)各不相同但最终目的都是为了提高计算速度和精确程度对于其他自由基应用稳态假设见公式13然而就象测不准原理中描述的时间和空间不可能同时测准C,H,O,与HO计算公式一样刚性ODEs的数值解法也存在着计算速度和计算精度(L:)(17)间难以同时满足的矛盾。计算速度的提高是建立在牺牲计算精度的基础上的而高精度的计算方法必然伴随着占有大d.其他物种量的计算机资源。在模式的实际应用过程中有时对计算速C(1-Ltn2)+P△t度要求很高而计算精度只要满足一定的误差范围即可故L"△tn/2(18)采用一定的方法简化计算过程可以大大提高计算速度1.6 Gong and Cho笔者介绍的方法基本上可以分成3类:第1类是GEARng and Cho法23将时间步长定为0.5h根据物质的类包括 GEAR LSODE LSODES) SMVGEAR等方法其基寿命将其分为2类:①r≥104s称为慢物质熜τ<104s称本思路来源于GEAR法;第2类是稳态假设类,此类方法的为快物质。计算公式为基本思路是将物种按其寿命大间步长之比分类,对不C=(C1C2)同类的物种分别采用不同的计算公式,包括QSSA(MQSSA),Y&BMY&B) RADM Gong and Cho等方法第3类是(19)d迭代法包括小参数法和 TWOSTEP法慢物质在策(n+1)s浓度用式(20)算在各种数值计算方法中GEAR法是一种较为精确的计Cr= Ci +f(C )At(20)算方法也是研究其他快速方法的基础。在研究模式机理阶对于快物质则用牛顿迭代公式计算出n-(n+1)s的浓度段对计箅结果的精确程度要求较高,般使用该方法这点增量δC可从表1看出。但是该法的计算速度较慢不适于大尺度的C2+-1+C2+f(Cm1+1,C2+-1)气质模式计算。 LSODE和 SMVGEAR两者为了提高计算速(21)度都减少了计算 Jacobian矩阵的量,即只有在收敛失败的情况下才计算。 SMVGEAR则又利用了矢量计算机的特点来其中C2=C21+8C2而送代计算出的C2又可以利用加快计算速度。总体而言第1类的3种方法基本上可以满式(22增加C1的精度足较高的计算精度,其时间步长可根据计算过程动态改变。C+=C+(C1)(22)稳态假设法代表了一类高速计算方法。由于大气体系中物1.7 TWOSTEP种的寿命相差很大是造成ODEs刚性的来源故其核心內容WOSTEP法2428是专为解大气化学中的ODEs设计是对寿命很小和很大的物种采用较为简单的公式而对寿命的。其用 Gauss-Seidel迭代代替了牛顿迭代公式如下适中的物种进行复杂计算这样可以提高计算速度。但是此C+l= Y+ yf( tn+I Cn+I)(23)类方法的精度依赖于时间步长取值大小,当Δt的取值较大式中+5=时计算精度较差而Δt取值较小虽然提高了精度但计算c+2速度较低。QSSA法和Y&B法从形式上可以看到是比较相似的,尤其经过改进后的 MQSSA和MY&B法可近似认为是一种方法。( ong and Chou法介于 LSODE法和Y&B法之间其对寿命较长的物种采用一般的积分方法,而对寿命较利用 Gauss-Seidel迭代式(23)可写成短物质的计算方法类似于 LSODE法。RADM法建立在C=C)Y+mxtt11、(24)RADM模式的基础上,它虽然针对特定模式提出但也有其普遍的适用意义特别是它提出了18种自由基的计算顺序TWOSTEP的迭代计算公式为众所周知在数值计算过程中不同的计算顺序会产生不同的结果,该算法给出的合理计算顺序对后人的研究会带来一定的帮助。第3类方法主要是用迭代法来求解,其中(k=12(25) TWOSTEP法假定了本步所用浓度是通过前2次的浓度计式中:m—物种数。算得汁汗田工内竿刚性的体系当然式(25)中物种的计算顺序对Fk的计算结果有影中国煤化工因而有必要对其进行误响但试验证明采用式25厢(26)对结果的影响是不大的。差和CNMHG交方法是计算各物种的均C=F(C1C2…Ck=12,Ck…,,Cm)(26)方根(RMS)对于第k种物质其RMS定义为2各算法间的比较RISk(27)以上所列方法虽然在基本原理、公式结构及条件参数上式中C",C"—分别代表精确值和计算值环境科学研究第13卷总步数,m=12,,M[5] Stern R M. Entwichlung und Anwendung eines ells mit ver在满足了一定的误差范围前提下看其对相同计算机所schiedenen chemischen Mechanismen[ A ] In: Meteorologische用时间来比较速度的高低。Say19等对 LSODE Gong&Abhandlungen Vol. 8. Institut Fur Meteorologie der Freien UrCho,Y& B QSSA4种方法作了比较。其模式选用了CBMersitat Berlin[ C ] Verlag von Dietrich Reimer, Berlin, ISSNⅣ机理(包含33种化合物和81个反应),对10几种物质采0342-4324,1994用不同方法时的RMS与CPU耗时情况进行了比较。可归[6] Milford J B Gao D Russell A G et al. Use of sensitivity analysis纳如下:①随着计算误差的降低,各种方法的CPU耗时也在to compare chemical mechanisms for air-quality modeling J ]. Environ Sci Technol 1992 26: 1179-1 189成倍增加。② LSODE法的RMS随着CPU耗时的增加降低[ 7] Simpson D Andersson-Skold Y Jenkin M E. Updating the chemi-最快而QSSA法降低最慢。③当对CPU耗时要求较少时4cal scheme for the EMEP MSC-w oxidant model: current status种算法的RMS近似QSA法的性能较高。 ODMAN等研[R]. EMEP MSC-W Note 2/93, Norwegian Meteorological Insti究认为使用Y&B和改进后的QSA法可以使计算精度得到tute Oslo ,Norway 1993保证的前提下大幅提高计算速度。他们同样采用CBM-Ⅳ[8] Simpson D. Biogenic emissions in Europe2: implication for ozone机理(一些条件不同于上者),以Y&B和改进的QSSA法分control strategie J ] Journal of Geophysical Rescarch 1995 .100别与 LSODE法作了比较结果发现在与 LSODE的误差控22891-22906在1%之内时,Y&B和QSSA法分别是 LSODE法的40与20[9] Pilling M Saunders s,Jenkin M,etal. Tropospheric Chemistry低计算过程中对QSSA法使用了质量守恒的控制技术,此Modelling[Eb/oLI.Downloadedfromhttp:/www.chem技术占时不大此外 MATHUR等与 JACOBSON等分别leeds. ac, uk/Atmospheric/MCM/ main. html. 1998对 MQSSA MY& B RADM及 SMVGEAR LSODE等方法作[10] Lurmann F W, Lloyd A C ,Atkinson R, A chemical mechanism了比较结论较为相似。for use in long-range transport/ acid deposition computer modeling结论[J]. J Geophys Res 1986 91: 10905--1093611] Strand A Hov O. A two-dimensional global study of tropospheric通过对以上各种方法的介绍和讨论可以看出大气化学ozone production[ J]. j Geophys Res 1994 99 22877-2289模式中刚性ODFs的计算问题的解决方法是多种多样的但[12] Andersson-Skold Y, Updationg the chemical scheme for the IVL大多都围绕在计算精度和计算时间的取舍上建立的可分为tochemical trajectory model[ R]. IVL B-1151, IVL,Gote-2类,一类在可以接受的误差范围内尽可能地提高计算速borg sweden 1995另一类则追求在可以忍受的计算速度下提高计算精度。显[13] Ruhnke r. Ein Verfabren zur Analyse von chemischen Reak然随着三维模式和大尺度模式的广泛应用,研究前者的人占tionssystemen und zur Erstellung von Chemienodulen fur atmo-多数spharische Modelld Ph-D]. Thesis University Essen germany在各种方法中 GEAR LSODE QSSA和Y&B法一直是人们研究的热点。其中GEAR和 LSODE等方法代表了一类[14] Kuhn M Builtjes P J H Poppe D et al. Intercomparison for the精度较高的算法,是人们研究其他方法时精确度的参照法gas-phase chemistry in several chemistry and transport models而 LSODE因其在计算速度上的优势已被广泛应用。QSSA[J]. Atmos Environ 1998 32 593--709和Y&B法作为快速计算方法其计箅精度与系统稳定性能15 Gear C W. Numerical Initial Value Problems in Ordinary Differ-也在不断提高,些大型模式计算就采用了该方法。其他方ential Equation R ] Prentiedhill Englewood Cliffs , J, 197116] Young T R Boris J P. A numerical technique for solving stiff or法的建立原则也很明显,一种是寻找新的原理,另一种是对dinary differential equations associated with the chemical kineticEAR与QSSA等已有方法进行综合改进其中有许多地方blems[J]. J Phys Chem, 1977, 81: 2424值得后人借鉴。2427参考文献[17] Hesstvedt E Hov O,Isaksen I S A. Quasi-steady-state approxi[ 1] Stockwell WR, Middleton P Chang J S et al. The second generamations in air pollution modeling: Comparison of two numericaltion regional acid deposition model chemical mechanism for regionalschemes for oxidant predictio[ J ] Int J Chem Kinet, 1978, 10air quality modeling J ]. J Geophys Res 1990 95: 16343-16367[2 Stockwell W R, Kley D. The Euro-RADM Mechanism: A Gas. 18 J Odman M T, Kumar N Russell A G. A comparison of fast chemPhase Chemical Mechanism for European Air Quality Studie[ C]eling[j]. Atmos Environ 1992Forschungszentrum Julich gmbH KFA) Julich Germany 1994中国煤化工[ 3] Gery M w Whitten G Z Killus J P et al. A photochemical kinet- [19]CN MH Gmparison of numerical methodsics mechanism for urban and regional scale computer modelinfor the integration of kinetic equations in atmospheric chemistry[J]. I Geophys Res19899412925-12956and transport models[ J ]. Atmos Environ 1995,29: 2585[ 4] Builtjes P J. The LOTOS- Long tem ozone simulation project summary report[ R ] IMW-TNO Report 92/240, TNO, Delft, The 20 Mathur R Young JO Schere K L , et al. A comparison of niNetherlands 1992merical techniques for solution of atmospheric equation[ J ].At第1期刘峻峰等大气化学模式的数值计算方法49mos Environ 1998 32: 1535--1553torized gear code for atmospheric models[ J ]. Atmos Environ[21] Hindmarsh A C. LSODE and LSODI, Two new initial value or199428273-284dinary differential equation solver J]. ACM Newsl 1980, 15: 10 [26] Hindmarsh A C. ODEPACK: a systematized collection of ode22] Chang I S Brosr R A Isaksen I S A, et al. A three-dimensionalNorth-Holland Amsterdam 1983 55-64acid deposition model: physical concepts anion[J]. [27 Yamartino RJ Scire J S Carmichael GR et al. The CALGRIDGeophys Res98792:14681-14700rid model: I. Model formulation[J[ 23] Gong W Cho HR. A numerical scheme for the integration of the中国煤化工1512gas-phase chemical rate equations in three-dimensional atmo- [28 1CNMHG境化学M]北京高等教育spheric model![ J ] Atmos Environ 1993 27A 2147--2160出版社1990.156~158[24] VerwerJ(. Gauss-Seidel iteration for stiff ODEs from chemical[29]关治陆金甫.数值分析基砒[M].北京高等教育出版社kinetic J ]. SIAM J Sci Comput 1994, 15: 1243-12501998.445-45125] Jacobson M Z, Turco R P. SMVGEAR: A sparse-matrix ve