数值方法

前言：想要写出一篇令人眼前一亮的文章吗？我们特意为您整理了5篇数值方法范文，相信会为您的写作带来帮助，发现更多的写作思路和灵感。

数值方法范文第1篇

关键字：内域；数值计算；有限元

中图分类号：TU315.3

1.引言

虽然目前技术和计算设备发展十分迅速，计算能力不断提高，一些大型通用有限元软件已经具备十分强大的分析功能，它解决了地震反应的许多实际工程结构分析的问题。但对于一些大型复杂体系而言，空间离散的自由度数目非常庞大，数值稳定性的限制要求时间离散的步距也不能过大。这样，开展结构地震反应分析时所需要完成的时空四维数值计算的工作量将变的很大。在工程设计中，需要分析各种工况下结构的地震反应行为，对比不同的设计方案并做出优化决策，从而要求在设计期内多次完成结构的地震反应计算。在数值精度的基础上，保证系统的稳定和提高结构地震反应分析方法的核心计算效率。研究高效率的大型复杂体系地震反应数值分析方法，减少计算费用仍然是非常现实的考虑，它具有重要理论意义和实用价值。

地震反应分析的数学物理模型就是波动方程。波动数值模拟包含两个部分，一是对人工边界的数值模拟，二是对内域的数值模拟。这里仅讨论对内域的运动节点的数值模拟问题。

现有的内域波动的数值模拟方法，按能量等效式，可划分为三类。一类是空间有限元方法。所谓空间有限元方法是指的空间用有限元法进行离散，而使用一个直接的方式离散时间。第二类是时空有限元方法。这里是指对时间使用和对空间一样的有限元法方式离散。时空有限元的时空域是空间有限元的空间域在时间域上增加了一个时间维度，区别只是处理时间域的先后上，空间有限元是先处理空间问题，然后处理时间问题，而时空有限元是同时处理时间和空间的问题。第三类是微分求积方法，是直接的方式对空间和时间的离散。下面分别阐述这三类方法的应用发展过程――

2.空间有限元方法

在1956年，有限元的概念首次被Turner等人提出，最早应用于弹性力学平面应力问题上。1963年，Besseling、Melosh和Jones等人发现有限元法和基于变分原理的里兹法是等效的。有限元法在处理连续介质问题上比普通里兹法更有优势。随后几十年，在解决复杂工程问题上，有限元法得到广泛的应用。

波动方程是时空耦合的，基于广义Hamilton原理的波动有限元方法通常也是时空解耦的数值过程。传统有限元方法的离散过程通常包含两步，先进行空间离散，将微分方程转化为常微分方程，然后对时间进行离散，即在时域对常微分方程进行数值积分。由于时空耦合的数值过程包含过多的自由度，求解这类方程在实际工程中很难实现，建立时空解耦的波动数值分析方法是这方面重要的工作。最直接的做法是实现空间及时间域的解耦，通常是只建立空间的有限元离散方法，而时间采用直接的假设，最常用的是采用逐步积分的方式进行离散。

逐步积分法简单来讲就是把最终速度和位移由它们的初始值和一个积分表达式来表示。加速度历程的积分决定速度的变化，速度的积分决定位移的变化。换句话说，加速度控制了速度的变化，因而可以由这一步向前获得下一时间步。解答这类问题，第一步先考虑时间步内的加速度问题，假设加速度是如何变化的，依据加速度和位移的关系，得到关于时间步的递推公式。所谓的Euler-Gauss法就是假设在时间间隔内加速度为常数。而Newmark[1]法是加入系数从而可以改变初始和最终加速度的权重从而得到加速度的一种方法。Wilson- [2]方法是假定在时间步距内加速度为线性加速度的一种数值方法，用内插公式得到体系在下一刻的运动。α方法[3]是在Newmark方法的基础上，通过修改结构动力方程的时间离散形式得到的。Chung 和 Hulbert 发展了一种无条件稳定的隐式广义α法，它由三个参数控制数值损耗。Runge-Kutta[4， 5]方法是在一个时间步距中内插若干计算点，利用这些计算点上函数值线性组合来代替函数的泰勒展开中的高阶导数，从而提高精度阶。不同于两步信息预测，线性多步法[6]发展了多步信息来预测下一步，从而获得了更高的精度。

逐步积分是最主要的时域积分方法，而它最常规的做法是差分法。时域有限差分法（Finite Difference Method， 简称FDM），是地震波传播模拟最广泛被使用的一类方法[kelly―Marfurt，1990][7-10]。有限元差分是将微分方程中的微分项用相应的差商代替，从而将微分方程转化为代数形式的差分方程。

由微商和差商的定义可以知道，微分的有限形式是差分，而导数的有限形式是差商。而微分和导数是以极限形式表示。数值计算方法导数可以用差商的自变量趋近于零来代替，换句话说，位移对时间的求导可以用有限差分的方式得到，位移的一阶导数是速度，二阶导数加速度。当世界步距为等步长，得到中心差分。差商代替微分方程中的导数，就可以得到微分方程的有限差分形式。。较之传统有限元法，虽然在定义几何结构上不够灵活，但时域有限差分法具有显式计算的优势，计算效率高，计算精度高于显示有限元法。这些方法的缺点是精度不高，只有一阶或者二阶精度，难以模拟高频问题，这类无法避免算法阻尼，从而形成较大的误差。

为了克服低精度的问题，很多高精度的数值积分方法相继提出。不仅仅四阶，五阶精度，甚至更高精度的数值积分方法处在发展之中。Golley[11]为了得到三阶精度的格式，对时间域采用高斯点作为配置点。在哈密顿变分原理的基础上，Riff和Barch[12] 采用3次多项式构造插值函数，得到了有条件稳定四阶精度数值方法。Argyris[13]等在前人的基础上，采用Hermitian插值，得到无条件稳定四阶精度数值方法。Kujawski和Gallagher[14]从另外一个角度，利用广义最小二乘法，在无阻尼结构中，构造了一种四阶精度的无条件稳定积分格式。Tarnow和Simo[15]利用二阶精度算法的结果，在此基础上缩短时间间隔，构造了一种四阶精度算法。钟万勰[16-19]在1994年提出了精细时程积分法。在保守系统下，积分结果保持系统守恒量不变。1995年在之前工作之上，钟万勰提出了子域精细时程积分法，提高了计算效率，使工作量大大减小，存储量大幅减小，为精细法应用提供了基础。2000年，顾元宪[20]提出了增维精细积分法，改进了矩阵的运算，提高了精度。但局限条件较多。2004年，汪梦甫[21]利用精细积分方法基本原理，采用高斯积分方法，建立了更新精细积分方法，这种精细方法适应度高，为得到了广泛的应用提供了条件，并且提高了精度。理论上汪梦甫分析了精细方法得到任意精度的可能。2009年，富明慧[22-25]在汪梦甫研究的基础上提出了高效高精度广义精细积分法。

这类方法的困难在于不容易构造精度较高的时间离散模式，并且空间有限元在时域上每个时间步逐步推进，因而会产生误差累积[26-32]。

3.时空有限元方法

时空有限元最早由J.TOden，I.Fried和J.H.Argyris等人提出。利用哈密顿原理建立关于时间边界的变分原理。几十年来，在各个领域得到全面的发展。在波动问题，动力问题，流体问题等方向得到广泛的应用。

传统的数值方法假定时间和空间是相互独立的，这样的假定广泛应用在实践当中并且在数学上很好理解。但是，使用上述技术同时也导致了数学上的困难。因为有用的信息可能通过速度在结构传播，而没有分离的时空格式能规避这种类型的数值困难。这就要求对时间的离散和对空间离散一样使用有限元。例如，当结构是三维时，这种格式需要四维的维度来表示。从而需要对时间和对空间使用相同的离散方式。空间有限元对时间和空间分别离散，空间节点形成的网格只能处在同一时间下，形成的是规则网格。时空有限元同时对时间和空间进行离散，理论上可以把网格划成任意形状，不必考虑相同的时间值，可以灵活的划分。有限元方法区别于其它方法在于它利用能量等效原理将偏微分方程进行积分。得到方程的弱形式。恰当的变分形式是有限元是否能够利用能量等效的关键。而时空有限元能否成功取决于能否找到对应的变分原理。

R.Riff和M.Baruch等人建立了一种能同时求解所有变量的时间有限元，这种有限元格式借鉴了空间有限元，把整个时间域看作是空间的延续，从而能够求解不同时刻的变量。冯康是提出 罗恩在冯康的基础上，完善了哈密顿型变分原理，发展了非传统哈密顿变分原理。罗恩对比等价的正则方程和相空间变分原理，认为即使是等价的，但由于形式的不同，产生的算法并不一定会有相同的效果。其结果就是相空间变分原理计算效率更高，也更接近物理问题的本质特征。沿着这个思路，罗恩建立了非传统相空间哈密顿变分原理。刘世奎对哈密顿原理进行了推广，构建时间边界条件，建立了有两个参数的广义哈密顿变分原理，完成了从单一变量向多种变量的转变。卓家寿等人在哈密顿体系下，推导了几种变分原理的等价形式。罗恩在此基础上推导了类变量广义哈密顿变分原理，它包含了所有的条件。基于这种变分原理提出在时间子域上进行五次样条插值的时间子域法。 2007年钟万勰[38]发现时空有限元可以更高效的解决边界问题和多尺度刚度问题。2013年朱宝[39]在钟万勰的基础上进一步研究多尺度和边界问题，对其稳定性进行了研究。

获得高阶精度，只需要提高多项式基函数的次数，从理论上来说，时空有限元可以获得任意精度。空间域增加时间域之后，单元的几何性质简单，避免了空间有限元的复杂边界。让传统的有限元得到更广的应用。是一种有很大发展空间的数值方法。

4.微分求积方法

Bellman和Casti[40]在1971年提出微分求积法的基本原理。此后，微分求积法因为原理简单，广泛应用在工程问题中，微分求积法得到快速发展。

微分求积法即DQM方法，本质上来说，函数和它的导数在给点节点的值用全部节点的函数值乘以系数并求和来代替。从而让微分方程可以转化成关于节点函数值的一组代数方程组。由此可知，DQM是一种数值技术，它通常被用来解决初值和边界问题。从本质上来说，DQM是特殊的一种加权残值法，而且是高阶的有限差分法。DQM方法相对有限元方法而言，并不需要变分原理就可以求解微分方法。

从微分积分法的原理出发，可以发现影响数值精度主要由两个因素构成。一方面是权系数的值，另一个方面是选取合适的网格离散点。其中网格离散点的选取和假设试函数的模式可以确定权系数，从另一个角度来说，试函数的假设和网格点的选取是决定精度的关键因素。从而研究人员也沿着这个思路对微分求积方法进行了探索。利用多项式是有限元试函数选取的基本思路。Bert和Wang[41]为试函数来求权系数，此时构成线性方程组的系数矩阵是勒得蒙矩阵。但由于当离散点数目增多以后，勒得蒙矩阵会出现病态。所以出现很大的误差。后来。Quan[42， 43]用采用了Lagrange插值，得到了微分积分法一阶和二阶精度的公式。Bert和Striz[41]建立了HDQ方法，采用用不同于多项式的谐函数作为试函数，开阔了研究思路。由上可知，试函数的选取并不是单一的，可以从多个角度来选取。不仅是谐函数或者多项式，甚至是指数函数，对数函数等初始函数都可以作为试函数进行研究。根据需要选择混合初始函数来得到试函数是值得探索的方向。网格点的选择方面，研究发现一些问题对节点选取是很敏感的，等距网点因为使用方便而被先采用，但是结果发现得到的解不够理想。其实真实的状况让均匀网格模拟显得不够合理，发展非均匀网格更可能得到高精度方法。 Bellman[40]就用勒让得多项式零点进行了研究，发现用其作为网格点提高了精度。在这启发之下，Quan等研究了切比雪夫多项式零点，并且用之与其它正交多项式作了对比，发现切比雪夫多项式零点更有优势。此外，微分求积的研究的方向是更加具体的研究，Bellman[44]用微分求积法应用到初值非线性偏微分方法得到高效精确的解法。 在多维问题上，微分积分法也得到了应用，Civan[45]得到了多维积分微分方程。Bert[46]首先将这种方法用到结构力学问题的求解。2001年在DQM法则的基础上，Fung[47， 48]建立了一种不同于边值问题的动力微分方程解法，解决了初值问题的动力微分方程。并研究时间网点选择方式对数值精度和稳定性的影响，提出了一种高精度的时间网格离散方法。

微分求积法虽然发展的历史比较短，但是由于这种方法原理简单，精度高，计算效率高，处理数据方便等优点。在工程领域有广泛的应用。

5.结语

有限元方法一直在数值模拟中很占有重要地位，这种思想在理论研究和实际应用中发挥着很重要的作用。以有限元原理为基础，发展的新方法让数值计算展现出新的活力。但是数值模拟是一门深奥的学问，在理论上和实际应用中还有很多不完善的地方，需要克服的缺点还有很多，本文作者仅仅就自己所涉及的研究领域做了一些简单的论述。

参 考 文 献

[1]. Newmark， N.M. A method of computation for structural dynamics[C]. in Proc. ASCE. 1959.

[2]. Wilson E L. A computer program for the dynamic stress analysis of underground structures[R]. DTIC Document， 1968.

[3]. Trefethen L N. Finite difference and spectral methods for ordinary and partial differential equations[M]. Cornell University [Department of Computer Science and Center for Applied Mathematics]， 1996.

[4]. 易大义，陈道琦. 数值分析引论[M].杭州：浙江大学出版社，1998.

[5]. 胡健伟，汤怀民. 微分方程数值方法[M].北京：科学出版社，1999.

[6]. Chew W C. Waves and fields in inhomogeneous media[M]. IEEE press New York， 1995.

[7].Wilson E L， Bathe K. Numerical methods in finite element analysis[M]. Prentice-Hall， 1976.

[8]. Clough R W， Penzien J. Dynamics of structures[R]. 1975.

[9]. Wilson E L， Farhoomand I， Bathe K J. Nonlinear dynamic analysis of complex structures[J]. Earthquake Engineering & Structural Dynamics. 1972， 1（3）： 241-252.

[10]. 刘晶波，杜修力.结构动力学[M]. 北京：机械工业出版社，2005.

[11]. Golley B W. A TIMESTEPPING PROCEDURE FOR STRUCTURAL DYNAMICS USING GAUSS POINT COLLOCATION[J]. International Journal for Numerical Methods in Engineering. 1996， 39（23）： 3985-3998.

[12]. Riff R， Baruch M. Time finite element discretization of Hamilton's law of varying action[J]. AIAA journal. 1984， 22（9）： 1310-1318.

[13]. Argyris J H， Vaz L E， Willam K J. Higher order methods for transient diffusion analysis[J]. Computer Methods in Applied Mechanics and Engineering. 1977， 12（2）： 243-278.

[14]. Kujawski J， Gallagher R H. A generalized leastsquares family of algorithms for transient dynamic analysis[J]. Earthquake engineering & structural dynamics. 1989， 18（4）： 539-550.

[15]. Tarnow N， Simo J C. How to render second order accurate time-stepping algorithms fourth order accurate while retaining the stability and conservation properties[J]. Computer Methods in Applied Mechanics and Engineering. 1994， 115（3）： 233-252.

[16]. 钟万勰. 结构动力方程的精细时程积分法[J].大连理工大学学报.1994，34（2）：131-136.

[17]. 钟万勰. 计算结构力学与最优控制[M].大连：大连理工大学出版社，1993.

[18].钟万勰. 应用力学的辛数学方法[M].北京：高等教育出版社，2006.

[19]. Zhong W X， Williams F W. A precise time step integration method[J]. Proceedings of the Institution of Mechanical Engineers， Part C： Journal of Mechanical Engineering Science. 1994， 208（6）： 427-430.

[20]. 顾元宪，陈飚松，张洪武.结构动力方程的增维精细积分法[J].力学学报.2000，32（4）.

[21].汪梦甫，周锡元. 结构动力方程的更新精细积分方法[J].力学学报.2004，36（2）：191-195.

[22]. 富明慧，刘祚秋. 固体力学中的变分差分方法[J].中山大学学报：自然科学版.2001，40（2）：13-15.

[23]. 富明慧，林敬华，刘祚秋.结构动力分析的广义精细积分法[Z].第九届全国振动理论及应用学术会议论文摘要集.2007.

[24]. 富明慧，林敬华. 精细积分法在非线性动力学问题中的应用[J].中山大学学报：自然科学版.2008，47（3）： 1-5.

[25]. 富明慧，廖子菊，刘祚秋. 结构动力方程的样条精细积分法[J].计算力学学报.2009， 26（3）： 379-384.

[26]. 裘春航，吕和祥，钟万勰. 求解非线性动力学方程的分段直接积分法[J].力学学报.2002， 34（3）： 369-378.

[27]. 郑兆昌，沈松，苏志霄.非线性动力学常微分方程组高精度数值积分方法[J].力学学报.2003，35（3）： 284-295.

[28]. 胡海昌.弹性力学的变分原理及其应用[M].北京：科学出版社，1981.

[29]. 赵秋玲.非线性动力学方程的精细积分算法[J].力学与实践.1998，20（6）：24-26.

[30]. 王金东，高鹏. 波动方程的精细逐步积分法[J].力学季刊.2000，21（3）：316-321.

[31]. 张洪武.关于动力分析精细积分算法精度的讨 ）[J]. 力学学报. 2001，33.

[32]. 梁立孚胡海昌.一般力学中三类变量的广义变分原理[J].中国科学： A 辑. 2000，30（12）：1130-1135.

[33]. 罗恩.几何非线性弹性动力学中广义 Hamilton 型拟变分原理[J].中山大学学报（自然科学版）.1990， 29（2）：15-19.

[34]. 罗恩，黄伟江.相空间非传统 Hamilton 型变分原理与辛算法[J].中国科学：A辑.2002，32（12）： 1119-1126.

[35].罗恩，梁立孚，李纬华.分析力学的非传统 Hamilton 型变分原理[J].中国科学：G 辑.2007，36（6）： 633-643.

[36]. 梁立孚，罗恩，冯晓九.分析力学初值问题的一种变分原理形式[J].力学学报. 2007， 23（1）： 106-111.

[37].罗恩，朱慧坚，黄伟江.Hamilton 弹性动力学及其辛算法[J].中山大学学报（自然科学版）. 2003，42（5）： 131-132.

[38]. 钟万勰，高强. 时间-空间混和有限元[J].动力学与控制学报，2007，5（1）：1-7.

[39]. 朱宝，高强，钟万勰. 三维时间-空间混和有限元[J].计算力学学报，2013，30（3）：331-335.

[40]. Bellman R， Casti J. Differential quadrature and long-term integration[J]. Journal of Mathematical Analysis and Applications. 1971， 34（2）： 235-238.

[41]. Bert C W， Xinwei W， Striz A G. Differential quadrature for static and free vibration analyses of anisotropic plates[J]. International Journal of Solids and Structures. 1993， 30（13）： 1737-1744.

[42]. Quan J R， Chang C T. New insights in solving distributed system equations by the quadrature method―I. Analysis[J]. Computers & chemical engineering. 1989， 13（7）： 779-788.

[43]. Quan J R， Chang C. New insights in solving distributed system equations by the quadrature method―II. Numerical experiments[J]. Computers & chemical engineering. 1989， 13（9）： 1017-1024.

[44]. Bellman R， Kashef B G， Casti J. Differential quadrature： a technique for the rapid solution of nonlinear partial differential equations[J]. Journal of Computational Physics. 1972， 10（1）： 40-52.

[45]. Civan F， Sliepcevich C M. Differential quadrature for multi-dimensional problems[J]. Journal of Mathematical Analysis and Applications. 1984， 101（2）： 423-443.

[46]. Bert C， Jang S， Striz A. New methods for analyzing vibration of structural components[Z]. Structures， Structural Dynamics and Materials Conference， 28 th， Monterey， CA， Apr. 6-8， 1987 and AIAA Dynamics Specialists Conference， Monterey， CA. 1987，936-943.

[47]. Fung T C. Solving initial value problems by differential quadrature method―part 1： firstorder equations[J]. International Journal for Numerical Methods in Engineering. 2001， 50（6）： 1411-1427.

[48]. Fung T C. Solving initial value problems by differential quadrature method―part 2： secondand higherorder equations[J]. International Journal for Numerical Methods in Engineering. 2001， 50（6）： 1429-1454.

数值方法范文第2篇

Lighthill利用保角变换的方法首先提出了二维翼型的反设计方法，Hicks，Murman和Henne等人将此方法发展为可应用机设计的工程设计方法。后Campbell等提出过一种带约束的直接迭代的表面曲率（CDISC）方法，Yu将其与N-S解算器耦合形成了一种翼型和机翼的设计方法。波音公司则将此方法发展成工程应用的设计方法，并广泛地应用于波音的B757，B777和B737NG等型号的设计过程，取得了很好的效果。例如在B777研制中由于使用了反设计方法，仅经过三轮机翼的设计便取得了满意的结果，使风洞实验的机翼模型大大少于过去B757和B767设计时的数目，充分表明了该设计工具的作用。可以说，反设计方法曾对民机设计起过革新性的推动作用；但反设计方法也有其固有的弱点（参见文献[13]的附录D）：首先，对于高度三维的流动要找到“好”的压强分布很困难；其次，不能保证所得结果为最优，即既具有高速巡航低阻的特性又在非设计条件下具有可接受的性能；最后，其他学科的约束会导致反复迭代。

低可信度CFD模型的数值优化方法

随着计算能力和数值优化方法的快速发展，应用基于CFD的数值优化方法于民机设计得到了很大的发展。这一方法的应用也从低可信度CFD模型开始，逐渐发展到采用先进的N-S方程解算器。波音公司发展了一种耦合TRANAIR[16]（一种全速势方程的有限元方法，可参见文献[13]附录B）和梯度优化方法的数值优化气动力设计方法，并在1992年形成了TRANAIR优化器的雏形[17]。经过近十年的改进，得到了一个适用于位势流/边界层耦合飞行条件的气动力优化设计工具[18-20]，具有多点优化设计能力，可处理高达600个几何自由度和45000个非线性不等式的约束条件（图1表示了TRANAIR优化过程示意图）。作为一个例子，图2给出了采用该软件对机翼/发动机短舱设计计算前后压强分布的对比，图a和图b分别表示了设计前后等马赫数线的分布。可以看出图a中挂架处出现激波；图b中短舱附近的机翼表面上消除了由于短舱干扰形成的激波。算例结果表明该设计软件可以处理很复杂的飞机/发动机综合设计问题。

高可信度CFD模型的数值优化方法现代优化算法可以分为依赖和不依赖梯度的方法两大类。

1.依赖梯度的优化算法

目前可用的大多数依赖梯度的数值优化方法都是从控制理论出发的，Jameson是此类方法的先驱者之一。尽管最初是由Pironneau提出利用控制理论进行椭圆方程系主控的外形优化的[21-22]，但Jameson首先提出了通过控制理论自动进行外形优化的伴随方程方法[23]并应用于跨声速流动。后来，Jameson和他的合作者，还有其他研究者，大力发展此方法，从全位势方程到Euler/N-S方程，从无粘设计到有粘设计，甚至从气动设计到气动/结构的耦合设计，形成了大量文献[24-36]。此方法不同于一般梯度优化方法之处在于它将外形作为一个自由表面，促使流动解和最终优化的外形同时趋于收敛，因而使优化方法具有很高的效率（其基本思想可参见文献[13]附录D）。

2.不依赖梯度的优化算法

最早无需梯度的优化算法有Powell（共轭方向法）[37]和Nolder-Mead的单纯形法[38]。最近Sturdza还应用后者于空气动力的设计[39]。近二十多年来人们更多地使用诸如模拟退火法[40]和遗传算法(GeneticAlgorithm-GA)等的搜索方法，特别后者更为人们所关注。Holland利用进化理论创造了遗传算法[41]（可参阅文献[13]附录D），即模仿生物的自然选择进行搜索以寻求最优解。与传统的搜索和优化方法相比，遗传算法具有下述4个特点[42-45]：1）不是直接作用于参变量集本身，而是对参变量集的某种编码运算。2）不是对单个点而是对多个点构成的群体进行搜索。3）直接计算适应值（函数），无需导数和其他辅助信息。4）利用概率转移原则，而非传统优化方法中的确定性原则。已有愈来愈多的研究和民机研制机构表现出了对这种随机寻优方法的浓厚兴趣，也已出现了不少利用遗传算法进行翼型或机翼优化计算的文献[46-56]。

3.对高可信度CFD模型数值优化方法的要求

分析最近十余年中出现的大量基于Euler/N-S方程的数值优化方法和文献，可以看出多数仍表现为学院式的探讨，提供可直接用于工程设计的方法和工具显得尚很有限，尽管已开始向这方面努力。这可能是因为：1）只是近几年来随DPW研讨会等的进行，数值模拟才可以比过去更正确地估算阻力值。2）工程界的空气动力外形优化需要在高维搜索空间中进行并存在大量的非线性约束，使优化问题十分复杂且计算开销巨大；3）巨大的计算量要求很丰富的计算资源和很长的计算时间，这与工程问题要求的迅速反馈相悖。

因此要使基于CFD的空气动力优化方法和软件成为日常的工程设计手段和工具需解决如下技术关键：1）具有建立准确计算诸如升力、阻力、力矩等敏感气动特性的正确流动模型的能力。比较现有的气动力优化方法可知，大多数方法还在使用不完善的流动模型，如基于Euler方程，甚至全位势方程等。虽然它们在一定条件下，如巡航小迎角飞行状态，可以提供合理的结果，但工程应用常要求准确地估算出阻力、俯仰力矩等敏感的气动特性，要求可计算整个飞行包线的飞行状态以及不同的复杂的几何外形等，这只能通过求解N-S方程来实现。顺便指出，有些文献（如文献[28]）虽以N-S方程为主控方程，但优化时的伴随运算子却是在没有考虑粘性流动的假设下得出的（参见文献[28]第6节）。为了提高计算准确度,最好在离散N-S方程时使用高阶的差分算子[53-54]。2）具有寻求全局最优的能力。通常基于梯度的算法容易陷入局部最优，而遗传算法等随机搜索的方法则具有取得总体最优的优点。3）能有效地处理大量几何和气动力的非线性约束。优化问题的最优解常常是位于不同维超曲面（hyper-surface）的交汇处，遗传算法不同于基于梯度的方法，不限于目标函数的光滑扩展，可应用于多重约束的情况[53-54]。4）可应用于不同的几何外形和设计条件。5）扫描高维搜索空间的计算有效性高，以满足设计周期和研制成本的要求。遗憾的是这正是遗传算法的主要缺点，即估算适应函数的高代价。可以采用多处理器上的有效并行计算来大大减少计算时间[57]，或在估算适应函数值时采用近似模型，如降阶模型[54,58]或响应面模型[50]等。

数值优化方法的发展现状和验证研究#p#分页标题#e#

1.空气动力优化设计计算的系列研讨会

近年来CFD学术界和航空业界都十分关注计算阻力的精度问题，这也是CFD应用于工程设计时所面临的第一个具有挑战性的计算。AIAA的应用空气动力学专业委员会在各方支持下，自2001年开始举行了DPW（DragPredictionWorkshop）系列会议[59]，参与者都用N-S方程求解相同的几何外形（翼/身组合体，翼/身/短舱/挂架的复杂组合体等），得到了一个巨大的计算结果数据集，可与已有的已经过修正的风洞试验值比较。由于参与的计算者所采用的数值方法、湍流模型、计算网格形式及数目等各不相同，此数据集可用作分析和讨论各种因素对CFD计算结果的影响。该系列会议至今已举行了5届，对推动和提高CFD计算阻力的精度很有意义。文献[13]的附录C中给出了前3届结果的分析和讨论。鉴于DPW系列会议的成功，AIAA应用空气动力学专业委员会针对CFD面临的第二个挑战---计算三维高升力外形的最大升力CLmax，于2009年发起并组织了类似的高升力计算研讨会，其第一次会议(HiLiftPW-I)已于2010年6月在美国举行,文献[60]是该次会议的总结。在上述工作的基础上，2013年1月AIAA又进一步在其ASM会议过程中形成了以加拿大McHill大学Nadarajah教授为首的空气动力优化设计讨论组，作为空气动力优化设计计算系列研讨会实际的组委会。讨论组讨论了：1）建立可供在一个有约束的设计空间中测试气动优化方法的一组标准算例。2）举行研讨会的时间。与会者一致认为，由于工业界对基于CFD的气动外形数值优化方法有强烈的需求，优化方法和工具的研制也已有了相当的发展，可以以类似于DPW的研讨会形式，通过对一系列复杂气动外形的优化，来评估现有的寻求最小阻力外形的各种优化方法的能力，并将结果向工业界/研究机构公布。与会者还认为第一次会议从二维和三维机翼外形开始是合适的，并请加拿大的与会者准备标准算例。第一次会议拟于2013或2014年的AIAA应用空气动力会议期间举行。

2.先导性的研究

事实上为准备此研讨会，波音的Vassberg，斯坦福的Jameson，以色列的Epstein及Peigin等三方从2007年起即开始了先导性的研究（pilotproject）,以积累经验和发现问题。三方用各自己开发的优化软件（MDOPT，SYN107，OPTIMAS）对第三届DPW会议的测试机翼DPW-W1独立地作优化计算[61,62]。波音研制的MDOPT[63]（也可参见文献[13]的1.7.3节）可使用响应面模型（InterpolatedRe-sponseSurface—IRS）的数值优化格式[64]，也可直接从流场解计算设计变量的灵敏度代替IRS模型完成优化。其流场解软件为TLNS3D[65]，计算网格点为3582225。Jameson开发的SYN107采用基于梯度的“连续”伴随方程方法[23,31]，其流场解软件为FLO107，计算网格点为818,547。

以色列航空公司开发的OPTIMAS采用降阶模型的GA算法，流场解软件为NES[66-68]，计算网格点为250,000。对三方独立优化后所得的外形再用不参与优化的流场解软件OVERFLOW[69]作评估计算，计算网格点数为4,000,000，以便能准确地计算阻力。结果表明，4个分析软件计算得到的阻力增量值的分散度在Ma=0.76时为5counts(1count=0.0001)，Ma=0.78时为10counts，因此很难确定哪个优化后外形最优。但从Ma=0.76，C=0.5单设计点的阻力改进结果（表1）[61]看，OPTIMAS优化后的04外形明显优于MDOPT优化后的M5和SYN107优化后的S4。文献[61]还讨论了从比较中可吸取的经验和教训。

一种基于高可信度CFD模型的数值

优化方法的构造本节将以OPTIMAS为例对如何满足可应用于工程实践的高可信度CFD模型数值优化方法的要求做一说明。

1.优化方法的构造及其特点

OPTIMAS是将遗传优化算法和求解全N-S方程的分析算法相结合的一种有效并鲁棒的三维机翼优化方法。1)其全N-S方程的流场并行解算器NES[66-67]基于高阶低耗散的ENO概念（适用于在多区点对点对接网格中的多重网格计算）[66,71]和通量插值技术相结合的数值格式，采用SA湍流模型，可快速准确地完成气动力计算，因此具有计算大量不同流动和几何条件的鲁棒性。作为例子图3和4给出了ARA翼身组合体Ma=0.80，Re=13.110时的升阻极线和CL=0.40时的阻力发散曲线[68]，使用的网格点数分别为，细网格（3lev）：900,000，中等网格（2lev）：115,000。可见升阻极线直到大升力状态的计算与实验都很一致。对比图中还给出的TLNS3D在细网格（2,000,000）中的计算值可见，无论升阻极线或阻力发散曲线NES的都更优。作为数值优化软件的特点之一是其在流场解算器中首次使用了高精度格式。2)优化计算的遗传算法中采用了十进制编码、联赛选择算子[42]、算术交叉算子、非均匀实数编码变异算子[72]和最佳保留机制。为解决搜索时总体寻优耗时大和求解N-S方程估算适应函数代价高的问题，在寻优过程中估算适应函数时采用当地数据库中的降阶模型[54,58]获取流场解（当地数据库是在搜索空间中离散的基本点处求解全N-S方程建立的），并以多区预测-修正方法来弥补这种近似带来的误差。多区预测-修正方法即在搜索空间的多个区域并行搜索得到各区的优化点，再通过求解全N-S方程的验证取得最优点。为保证优化的收敛，寻优过程采用了迭代方法。3)在整个空间构筑寻优路径（图5），扩大了搜索空间和估算适应函数的区间[54]。4）为提高计算效率，OPTIMAS包含了五重并行计算：Level1并行地求解N-S方程Level2并行地扫描多个几何区域，提供多个外形的适应函数的计算（level1隐于level2中）。Level3并行的GA优化过程（level3隐于level4中）。Level4并行地GA搜索多个空间。Level5并行地生成网格。5）采用单参数或双参数的BezierSpline函数对搜索空间参数化;并基于优化外形与原始外形的拓扑相似自动地实现空间网格的快速变换。

2.优化设计的典型结果

文献[53]~文献[58]给出的大量算例充分表明了OPTIMAS优化软件的优异性能。本文5.2中给出了其优化三维机翼的性能，这里再补充两例。1）翼身组合体整流（fairing）外形的优化文献[73]讨论了某公务机翼身组合体机翼外形优化的单点和多点设计两者性能的比较。结果表明，多点优化设计能同时保证设计的巡航状态时，和高Ma数飞行，起飞等非设计状态时的良好性能。文献[74]进一步讨论了翼身组合体整流外形的优化设计。流动的复杂性（三维粘流/无粘流强相互作用的流动区域）和几何的复杂性（三维非线性表面）使整流外形的设计经历了传统的试凑法，基于Euler解的试凑法等，最后才发展为现代完全N-S解的数值优化方法。文献[74]采用了这种方法，先作机翼外形优化，再作整流外形优化，然后再作机翼优化，整流外形优化，……依次迭代，直至收敛。优化中用双参数的BezierSpline函数将整流外形参数化，所得搜索空间的维数ND=（2N-2）*（M-1）决定的参数化整流外形与实际外形的差别在M=10，N=4，ND=54时可准确到0.3mm（满足工程需求）。计算网格数为90万。表3给出了设计条件和约束，表4给出了设计点的阻力值比较。由表4可知，GBJ2的减阻为16.7，50%DC，GBJFR1的减阻为10.7，32.1%DC，GBJFR2的减阻为5.9，两次优化机翼的减阻总计为22.6，67.9%DC，优化机翼和优化整流外形减阻作用分别约占2/3和1/3，可见整流外形的优化也是十分重要的。约束则使减阻损失4.6（如GBJFR3-GBJFR1）。图6至图9分别为原始外形，GBJ2，GBJFR2和GBJFR4的整流处等压线分布，可见整流外形的优化消除了原始外形和GBJ2中存在的激波。图10和图11分别给出了Ma=0.8时升阻极线和CL=0.4时阻力发散曲线的比较，可见优化设计不仅对设计点，对非设计状态也都有好处。2）翼身融合体飞机气动外形的优化[75]优化设计以英国克朗菲尔德大学设计的BWB外形[76]为出发外形，该外形的主要设计点为，。在数值优化计算中还考虑了，的第二个设计点和，（起飞状态）的第三个设计点。几何约束有剖面相对厚度，前缘半径，后缘角，每个剖面的樑处还附加两个厚度约束，其中上标b表示出发外形，*表示优化外形，下标i表示第i个剖面。附加的空气动力约束为对俯仰力矩的规定。采用Bezier样条描述几何外形，总设计变量为93个。表5给出了设计计算各状态的条件和约束，其中是权系数。表6给出了优化计算结果。#p#分页标题#e#

单点优化的BWB-1结果与文献[77]的结果相比较可见，文献[77]采用Euler方程的无黏优化使阻力降低了26counts；而这里的BWB-1全N-S方程优化使阻力降低了52counts，显示了此黏性优化方法的优点。比较有、无俯仰力矩约束时优化得到的BWB-2和BWB-1表明，尽管BWB-1阻力降低的效果突出，但其值过大，出于稳定性考虑而不能接受；BWB-2的阻力虽比BWB-1大了1.9counts，却满足了力矩的要求。表6中的双点优化设计（BWB-4），使第三设计点（低速状态）的达到1.671（消除了BWB-2达不到设计要求1.63的缺点），且基本保持了主设计点的阻力收益，为196.6。然而BWB-4在时的阻力达216.6，高于BWB-2的213.4，表明需要三个设计点的优化设计（BWB-3）。BWB-3在时，为202.5（比两点设计值减小了14.1），同时满足了其它两个设计点的性能要求。图12至图15给出了所有设计状态和时的极曲线，时的阻力发散曲线和时的随迎角α变化的曲线。由图可见，时所有优化设计的极曲线都非常接近，相比于原始外形的极曲线，性能有了很大改进；时也一样，特别是三点优化设计的BWB-3，优点更明显。阻力发散曲线也都有了很大改进，在前所有的总阻力基本保持常值，单点与两点优化的阻力发散点接近，而三点优化的可达附近。由图15可知，没有考虑低速目标的BWB-1和BWB-2具有较低的，将低速目标计入设计状态的BWB-3和BWB-4所得的皆优于原始外形的。上述结果表明三点优化设计具有最佳的优化效果和总体最好的气动性能。Fig.15LiftcoefficientCLvsangleofattackatMa=0.2最后，上述各优化结果在（主设计点）时的阻力值基本相同，但几何外形却差别不小，由此可见，外形阻力优化问题没有唯一解[75]。上述计算是在具有456GBRAM，114MB二级高速缓存的机群环境下通过“过夜”方式完成单点优化设计，在1.5-2天的计算时间内完成三点优化设计的，计算时间可满足应用于工程设计的要求[75]。

结束语

数值方法范文第3篇

1概述

随着科学技术水平的不断提高和工程建设规模的不断扩大，在土木建筑、水利工程和路桥工程中桩基承载力、沉降量大小和堤坝稳定性等力学问题变得十分复杂。而事实上，软土路基呈现出空间非线性沉降变形规律，土体的变形协调条件和应力平衡条件也十分复杂。这些问题已经很少能用数学方法求得精确解或通过模拟试验得到定量解，大多数课题必须借助于计算机和计算数学用数值分析的方法求出近似解。目前用于地基沉降量分析的数值方法主要是差分法、有限元法和边界元法，其发展趋势是有限元法与差分法或与边界元法相结合解决课题，以发挥各种方法的优越性。

2分析方法解析

2.1差分法

差分法的基本精神就是将研究区域用差分网格离散，对每一个节点通过差商代替导数把固结微分方程化成差分方程，然后结合初始条件和边界条件，求解线性方程组得到数值解。

以平面问题为例，差分法可得到所研究平面内在各个时间的孔隙压力的分布，因此可以导出初始沉降Si与任何时间t的总沉降S或固结沉降量Sc。由于土的竖向应变为

（2.1）

地基中沿着某一铅垂线的沉降为：

（2.2）

H为有效压缩层厚度。由于上式中采用不排水指标，，，所以孔隙压力，应用总应力来计算。总沉降为：

（2.3）

式中，E和随有效应力而变化，因此上式可算得任意时间的总沉降。任何时间的固结沉降为。当计算最终总沉降时，式中的孔隙压力。

2.1有限元法

有限元是地基和结构作为一个整体来分析，将其划分网格，形成离散体结构，形成有限数目的区域单元，这些单元体只在结点处有力的联系。材料的应力－应变关系可表示为

（2.4）

由虚位移原理可建立单元体的结点力与结点位移之间的关系，进而写出总体平衡方程：

（2.5）

式中，分别为劲度矩阵，结点位移列阵和结点荷载列阵。然后结合初始和边界条件求解线性方程组，在荷载作用下算得任一时刻地基和结构各点的位移和应力，得到问题的数值解。有限元法可以将地基作为二维甚至三维问题来考虑，反映了侧向变形的影响[1]。

（1）弹性有限元法

土的弹性应力——应变数学模型包括线性和非线性弹性模型两种。用线性弹性模型计算地基的位移和沉降，只适用于不排水加荷情况，并且对破坏要有较大的安全系数，一般不发生屈服的情况。实际上土体中的应力状态都可能发生屈服，其应力——应变的关系是非线性的。此外，除了理论建模还可通过试验拟合的途径，即根据土体的应力应变试验曲线，用弹性系数的连续变化来逼近真实的试验曲线建立起各种不同形式的非线性弹性模型[49]，第一种是以E（弹性模量）和μ（泊松比）两个弹性常数表达的，称E-μ弹性模型，第二种是以K（弹性体积模量）和G（剪切模量）两个弹性常数表达的，称K－G模型。另外一种是南京水利科学研究所采用的非线性的变弹性体模型，它的特点是不用常规的弹性常数，改用两个非线性函数，来表达应力与应变之间的关系。

典型的E-μ弹性模型是Duncan-Chang的双曲线模型。Duncan-Chang的双曲线模型可以考虑应力历史对变形的影响，若应力低于或高于前期固结应力，则采用不同的弹性模量计算公式。它还可考虑土与结构共同工作，考虑复杂的边界条件，考虑施工逐级加荷，考虑土层的各向异性等。

土体在路堤荷载作用下的变形过程，伴随着主应力大小的不断变化及主应力方向的不断偏转，土的弹性模量及泊松比也随之改变。路基真实应力场为初始应力场叠加上由路堤填筑荷载引起的附加应力场。假定路堤填筑荷载引起的附加应力场可近似采用同样荷载作用下在线弹性半空间无限体所产生的附加应力场。采用邓肯－张非线性模型描述土体本构关系，由于路基纵向应变为零，则其增量形式的本构关系可表示为[2]：

（2.6）

式中：，分别为切线弹性模量和切线泊松比，可表示为：

邓肯－张模型含有8个参数须三轴排水试验确定。

K-G弹性模型用三轴等向压缩试验测p和体积应变，由此建立K的公式。Naylor取切线体积模量为p的线性函数，即：

（3.34）

由三轴剪切试验可建立切线模量的公式：

（3.35）

参数均由三轴试验确定。

（2）弹塑性有限元法

土的弹塑性模型将土的应变分为可以恢复的弹性应变和不可恢复的塑性应变两部分。弹性应变增量可以用弹性理论计算，塑性应变增量可以用增量塑性理论求解。土的弹塑性计算模型一般分为理想塑性和硬（软）化塑性模型2种。

弹性非线性模型是假定全部变形都是弹性的，用改变弹性常数的方法来反映非线性；弹塑性模型则把总的变形分成弹性变形和塑性变形两部分，用虎克定律计算弹性变形部分，用塑性理论来求解塑性变形部分。典型的弹塑性应力－应变模型有：剑桥模型、修正剑桥模型、拉德模型、椭圆－抛物双屈服面模型、沈珠江双屈服面弹塑性模型、“空间准滑面”模型（SMP模型）等[3]。

3结语

土的应力－应变关系非常复杂，任何模型都有它的局限性。过去利用E-μ弹性模型和K-G弹性模型，虽然对坝和地基的应力应变分析都曾作出有意义的贡献，但它们都忽视了土的剪胀性和应力路径的影响。而剑桥模型只适用于只有剪缩而没有剪胀的正常固结粘土与松砂。计算时要找出一个数学模型来全面正确地表达土的这种特性是难以想象的。

参 考 文 献

1 钱家欢，殷宗泽. 土工原理与计算（第二版）. 北京：中国水利水电出版社，1996

2 费正华，邓水明. 应用邓肯－张非线性模型近似计算路基沉降. 中南公路工程，2001

3 河海大学，江苏宁沪高速公路股份有限公司. 交通土建软土地基工程手册. 北京：人民交通出版社，2001

4 Casagrande,A. Classification and Identification of soils,Trans.ASCE,113,1948:901-991

数值方法范文第4篇

关键词： 特殊角三角函数值 数形结合 函数图像 函数单调性

高一下学期一开始，教学内容就进入了三角函数。这一节公式很多，需要记忆的东西很多，但是只要学生能够每天定时定量地练习题目，公式自然能够熟练应用，而且烂熟于心。而且学生本身对公式也比较重视，因为公式的各种灵活运用，能够激发学生的兴趣。他们做完一道题目之后，会互相讨论，看还有没有其他方法。这源于笔者平时在教学过程中不断地鼓励学生去思考、去总结，不但要学会，而且要会学；把新课标强调的“提高学生自主学习能力和探究学习能力”这一思想。尽管公式学生已经很熟悉了，但是仍有学生会在三角函数的题目上卡住。为什么呢？因为这一节还出现了大量的特殊角，如30°，45°，120°，甚至还有75°。学生觉得特殊角不如公式灵活，只能去死记硬背。因为对特殊角不熟悉，导致他们看到，却不知道这就是tan60°；看到cos120°，还要苦想该用哪个诱导公式来诱导。虽然他们不止一次地体会到特殊角的重要性，但是他们仍不能接受硬背这样传统的学习方法。随着高一课程的结束，高二的解析几何、立体几何中仍旧会出现这些特殊角。现在学生若是没有记住，到了高二的时候怎么办？

针对这个问题，笔者查阅了很多资料，大概是这个问题基本都是靠硬背来解决，因此所能找到的资料甚少。一个偶然的机会，笔者看到学生在算sin30°的时候，画了一个30°的直角三角形，很显然这个方法不能解决sin210°，但是笔者还是表扬了这个学生，因为他在想办法解决问题。这个发现使笔者体会到，通过高一上函数部分的强化，学生现在已经有了画图解决问题的思想，能不能用数形结合的办法来解决这个一直让学生比较头痛的问题呢？其实学生在特殊角这部分暴露的问题很明显，对［0，90°］范围内的角度接触时间较久，比较熟悉，只是对高中阶段才推广的“大”角比较陌生。通过跟学生共同探讨，笔者发现以下几个方法比较适用。

一、利用三角函数图像

y=sinx， y=cosx， y=tanx的图像，在教材里面有三节内容，对它们的图像和性质研究是三角函数部分的重点内容。因此，若学生产能够画出它们的图像，不要说cos150°，哪怕是sin225°，或者是更大的角，也能够一眼看出。但这种方法的前提条件是，学生必须得记住这三个三角函数图像。

二、 利用直角坐标系

以sin225°为例，在平面直角坐标系中，画出225°所在的终边，再做出它的延长线，这样在第一象限内就出现了一个以它的延长线为终边的角，而此时学生就可选取非常熟悉的45°为此延长线的代表角。接下来做的事情，只需在延长线上取点P（x，y）。由图像可知，两线关于原点对称，故此，两线上的点的纵横坐标均互为相反数，则在原线上可以作出P点关于原点的对称点P′（x，y），由任意角三角函数的定义即可得出sin225°==- =- sin45°=- 。通过刚才的推导过程，也可以得出这样的结论：若两角的终边关于原点对称，则它们的正弦值（或余弦值）互为相反数。这个结论的得出，让学生知道了第三象限的角与他们熟悉的［0，90°］的角的关系，自然他们想到了第二象限和第四象限。

通过刚才的推导可知：若两角终边关于y轴对称，则它们的余弦值互为相反数，正弦值相等。从另一个方面来看，这两个角为互补的关系，所以刚才得出的结论也可叙述为：互补的两角正弦值相等，余弦值相反。第四象限的角推导过程与上述过程类似。通过作出其关于x轴的对称轴可知：若两角终边关于x轴对称，则它们的余弦值相等，正弦值互为相反数。

综上可知，若要解决特殊角的三角函数值，只需要在坐标系中，画出它的终边所在的位置，通过做它关于原点（或x轴、或y轴）对称线，找出第一象限我们非常熟悉的角，判断出符号即可。

三、利用函数单调性

对于某些连［0，90°］都记不住的学生，除了用本文一开始提出的画特殊三角形以外，还可以利用三角函数本身的单调性。由于特殊角的三角函数值只有几个数值：0，，，，1，联系y=sinx在［0，90°］内单调递增，故对号入座，sin0°=0，sin30°=，sin45°=，sin60°=，sin90°=1，相应余弦值则可通过直角三角形里得出的结论，互余的两角正余弦值互换得到。对于数值比较麻烦的15°和75°，我们可以通过构造成两角和或两角差的方法，快速算出它们对应的三角函数值。

笔者提出了这几个方法后，很多学生都不再觉得特殊角三角函数值很难背了。究其原因，是在笔者提出的方法的基础上，他们非常熟练地运用函数图像、函数单调性等函数知识。这些方法中所涉及的数形结合思想，锻炼了他们的数学思维能力，记忆的过程也成为了他们思考问题的过程。学生觉得学有所得，学有所用，这些特殊角三角函数值的记忆过程不再是枯燥无趣的几个数字，而是生动形象的图像、函数知识。而且，在这一过程中所涉及的数形结合方法，其本身就是高中数学阶段重要方法之一。

参考文献：

数值方法范文第5篇

【关键词】地铁；火灾；数值模拟

随着城市的发展，地铁已经成为城市交通的命脉。人们大多喜欢乘坐地铁出行，高峰时期地铁中人员密度非常大。由于地铁火灾的特征不同于地面建筑，当地铁发生火灾时及易造成严重的财产损失和人员伤亡。地铁火灾实验很难进行，为了便于问题的探讨，常选择数值模拟的方法来对其进行研究。[1]

1.地铁车站模型选择

为了便于问题的研究，地铁站火灾模拟中选取沈阳地区典型浅埋岛式地铁车站为研究对象。车站公共区域通风空调系统按站厅、站台均匀送-回（排）风设计（闭式系统）。站台层设置轨底排风风道，排风口与刹车电阻箱对齐；轨顶排风风道，排风口与列车空调冷凝器对齐；站台排风风道。车站空调回（排）风机兼作车站公共区消防排烟风机，回（排）风风道兼作排烟风道。站厅层送风风量为66000 m3/h，站台层排烟量为132000m3/h。地铁车站公共区域剖面如图1所示，车站主体两层结构均在地面以下。

2.地铁车站计算区域

合理简化计算区域，在不影响整体研究效果的情况下，适当的降低计算成本。站厅与站台层之间的结构不作为计算区域（即忽略两层的导热作用）。站台层长120m，候车区计算层高3.5m，宽11.8m；隧道出入口4个，尺寸为3.75m×5.26m；轨顶排烟口24个，尺寸为1m×0.6m；轨底排烟口48个，尺寸为0.8m×0.3m；送回风口72个，尺寸为0.5m×0.5m。站厅层长88m，宽19.3m，计算层高3.5m；送风口30个，尺寸为0.5m×0.5m；车站出入口4个，尺寸为4.5m×2.5m。

3.火源边界条件设置

世界各国对于各种可燃物的释热速率并没有明确的界定，根据香港地铁工程技术人员的保守估计，地铁车站公共区火灾释热速率为2MW[2]；而英国Building Research Establish出版的报告中显示，在人员聚集的公共场所可能的火灾释热速率为2～2.5MW。根据文献[3]，本次研究地铁车站公共区域火灾释热速率5分钟后稳定在5MW的火灾情景。

4.数值实验方法验证

数值模拟方法是否正确常常需要通过相应的实验来进行验证，但是地铁站火灾实验很难进行。下面验证的过程是用地铁站火灾数值模拟实验的方法，针对已有隧道火灾实验数据来进行的对比，以表明火灾数值模拟实验方法的正确。

4.1模型建立

Kumar曾对隧道火灾进行了实验测试，隧道尺寸长390m、宽5m、高4m[4，5]。火源释热速率在10s内，从0线性增长到7.225MW。根据上述条件建立相应物理模型，在此基础上采用与地铁站火灾数值模拟实验相同的方法来进行对比。总计算时间为120s，截取部分计算区域，对数值实验结果进行分析。

4.2数值实验结果

自然通风情况下，120s时火源中心X=2.5m截面处，可以看到强烈对流产生的温度分布和速度矢量分布，基本以火源为中心成对称分布；热烟气向上运动，并在隧道顶部形成热烟气层向隧道两端扩散，最大扩散速度为7m/s；由于烟羽流的作用，烟气层下方的空气向着火源方向运动，速度较慢小于2m/s，与文献[5]的计算值相似。

4.3数值模拟实验与实验测试对比

图2为距火源10m处，隧道中心线上的温度分布和速度分布与文献[5]的计算值和Kumar实测值的对比。速度分布相近，趋势相同；温度分布稍有偏差，主要是没有考虑辐射传热而引起的温度稍高，但总体趋势相同。综上所述，数值模拟方法可以较好的描述火灾的发展过程。

5.总结与展望

在地铁站火灾数值模拟研究中，地铁站物理模型的选择，模拟时边界条件的设置，数值模拟方法是否正确，都对实验结果造成一定的影响。物理模型的选择可以根据已有或在建的项目来构建，边界条件的设置可以根据相关研究来确定，数值模拟方法是否正确就要和实际燃烧实验来进行对比。

研究地铁站发生火灾时，不同防排烟系统模式下，温度场、速度场、以及火灾烟气蔓延和排烟效果，以及列车活塞风对地铁站火灾温度场、速度场以及烟气蔓延的影响，对火灾中人员安全疏散起到了指导作用。

参考文献：

[1]徐硕.地铁火灾烟气蔓延研究初探.中国科技纵横，2011年7月（下）总第122期.

[2]李启荣，黎少其.地铁火灾系统保障研究.香港：城市轨道交通研究，2001.

[3]杨昀，曹丽英.地铁火灾场景设计探讨.自然灾害学报，2006，Vol.15（4）：122～125.