数值方法范文

数值方法范文第1篇

关键词：数值方法，岩石力学

中图分类号：O3文献标识码： A 文章编号：

1 引言

当前岩石力学数值计算方法得到迅猛发展，出现了有限差分、有限元、边界元、离散元、块体元、无限元、流形元及其混合应用等各种数值模拟技术，使复杂岩石力学工程问题的设计发生了根本性的变化[1]。不同数值计算方法的结合,更能发挥各种数值方法优势互补的作用，如有限元—边界元的混合、有限元—离散元的混合、有限元—无限元和有限元—块体元的混合采用等。然而，由于岩体具有非连续、非均质、各向异性、天然初始地应力影响、地下水影响及复杂边界条件处理等诸多复杂性使得当前岩石力学数值计算仍然是一个值得探讨的问题。

2 常用岩石力学数值计算方法应用分析

2. 1 有限元法

1966 年,布理克[2] (W. Blake)最先应用有限元法解决地下工程岩石力学问题。目前，在岩石力学数值计算方面，有限元法主要用来求解线弹性、弹塑性、粘弹塑性、粘塑性等问题，是地下工程岩体应力-应变分析最常用的方法。其优点是可以部分地考虑地下结构岩体的非均质和不连续性，可以给出岩体的应力、变形大小和分布，并可近似地依据应力、应变规律去分析地下结构的变形破坏机制。

2. 2 边界元法

边界元法在20 世纪70年代中期得到迅速发展，在处理半无限域或无限域问题方面非常方便，适用于解决岩石力学问题尤其是岩石力学中地下开采的有关问题。该法只在求解区域的边界上进行离散(剖分单元)，这样就把考虑问题的维数降低了一维，这也是边界元法的优点。例如,在线弹性区域或无限域、半无限域采用边界元法，在非线性的区域采用有限元法，充分发挥各自的优势，使计算效率、计算精度得到提高和改进，这对工程实际应用是很有意义的[2 ,4]。王泳嘉[4]等人讨论了边界元的应力不连续法和直接边界积分法，并用应力不连续法给出了位移不连续时的解。

2. 3 离散单元法

离散单元法(Distinct Element Method)是20世纪70年代后发展起来的一种用于节理岩体应力分析的数值方法，特别适用于节理岩体及其与锚杆(索)的应力分析。该方法以结构面切割而成的离散体为基本单元，其几何形状取决于岩体结构中不连续面的空间位置及其产状，应用牛顿运动定律描述各块体的运动过程，块体可以发生有限移动与转动，体现了变形和应力的不连续性。沈宝堂等学者通过两个模型实验的结果，并与离散元法数值模拟的结果相比较，验证了离散元法用于边坡稳定性分析的可行性[5]。郭爱民[6]等人利用离散单元法研究矿山边坡的破坏机理，笪盍等人[7]利用离散单元法对盘石镍矿边坡进行稳定性分析和计算，计算结果与现场岩移进行比较表明两种结果吻合较好。

2. 4 关键块理论

关键块理论KB T ( Key Block Theory)是在1985年首先由Goodman 教授和石根华博士提出并用于工程稳定性分析的另一种数值计算方法。块体理论认为，在开挖面上所揭露的块体，可以分为可能产生向开挖面运动的块体和不可能向开挖面运动的块体。不可能向开挖面运动的块体即为稳定块体。石根华[8]等人提出搜寻关键块体一般方法，并介绍了传统关键块体理论的近期发展、应用和局限。

2. 5DDA法

不连续变形分析方法DDA (Discontinuous Deformation Analysis)是由石根华博士首创，基于岩体介质非连续性发展起来的，以模拟复杂加载条件下离散块体系统不连续大变形的力学行为为目的的平行于有限元法的一种数值方法，与有限元法不同之处是可以计算不连续面的位错、滑移、开裂和旋转等大位移的静力和动力问题。孙亚东[9]等人利用DDA法分析岩质边坡的倾倒破坏，并与Goodman 和Bray提出的基于极限平衡原理的分析方法进行比较。邬爱清、丁秀丽等学者应用DDA法研究了复杂地质条件下地下厂房围岩的变形与破坏特征[10]。

2. 6 FLAC方法

Cundall根据有限差分法原理,提出了FLAC( FastLagrangian Analysis of Continuum)分析方法。该方法采用了混合离散方法、动态松弛方法和显示差分方法，不形成刚度矩阵。它的求解方法虽同离散元法的显式按时步迭代求解，但是结点的位移连续，本质上仍属于求解连续介质范畴的方法。

2. 7 数值流形法

数值流形法是利用现代数学—“流形”的有限覆盖技术建立起来的一种新的数值方法，将有限元、不连续变形分析(DDA)和解析方法统一到一种计算方法中，它吸收了有限元、DDA和解析法等的优点，通过采用分片光滑的覆盖函数，对连续和非连续问题建立了统一的计算格式，是一种十分适合于岩土工程分析的数值方法[12] 。它最早由石根华博士在1991年提出并率先应用在块体与节理岩体的变形模拟中，在国内外学者的共同努力下，二维数值流形方法已经拥有一套比较完善的理论，而且应用方面的探讨也日渐增多[13]。郑榕明、张勇慧等也曾经对DDA法作了大量的研究工作[14 ]，在此基础上发展了二维流行元程序，并应用在地下洞室开挖的模拟中[15]，刘红岩等人[16]利用数值流形方法对一层状岩石边坡的倾倒破坏过程进行了模拟。

3岩石力学数值方法应用中的几个瓶颈与展望

3. 1 计算参数的取值问题

由于岩体性态与环境的复杂性,准确确定这些参数并非易事,因而数值分析手段至今仍不能最终为工程设计和工程决策提供可靠依据,幸运的是至今仍火热的反演分析方法有望在这方面为原始参数的取值提供一种新的途径[17]。中国近年来在反分析方面进行了大量研究工作，已由简单的线弹性反演问题发展到非线性、粘弹塑性反分析，从单一的毛洞围岩到考虑支护结构体系的反演,从有限元位移反演到边界元位移反演，从确定性反演到非确定性随机反演等等。而且反分析的目的已不仅仅是得到模型参数，更重要的是应用这些参数进行相应的时间序列值分析以及从参数估计发展到模型识别进而建立新的模型，以便对工程效果做出更合理的评价和有依据的预测。联邦德国对计算参数的取值问题也十分重视，他们认为计算输入的参数必须源于客观实际,地质勘探、岩石力学和数值计算必须紧密地结合起来，甚至从事计算的人需要自己动手到现场去取得第一手资料，而不只是单纯依靠委托单位提供的参数，这一点同主观臆断假定参数或依靠委托单位提供参数而不深入实际研究参数可靠性等脱离实际的作法形成鲜明的对比[18]，是我们应该值得学习和重视的。

3. 2 本构关系的选取问题

从事数值模拟计算的人都知道岩体性质比迄今为止人们所熟知的任何工程材料都要复杂得多，它既非连续介质,又非松散介质，既非弹性体，又非塑性、粘性体，从而导致计算中采用的本构关系很难准确把握。尽管用数值模拟对岩体结构进行力学分析的方法得到了广泛的应用，并且取得了许多成果，但是不敢断言这种方法在将来是否会对这样一类问题的研究有新的突破，至少目前还不可能将这一类问题的研究提高到一个全新的高度[19]。基于这样的原因，人们也在力求寻找其它的弥补途径，有学者改变单一的确定性正向分析方法，适应岩体的非确定性特征，将数值模拟方法与反分析方法、随机方法、模糊方法、人工智能等结合起来，这或许是数值模拟方法的一个正确发展方向[19]。

3. 3 单一数值方法的局限性问题

为了吸取各种计算方法的长处以弥补其不足,近年来涌现出大量的各种数值方法的耦合计算,这种思路进一步反映了岩体工程的计算特点，清华大学的周维恒先生在1993年就断言这种思路对岩石力学数值计算的发展是十分有益的[20] 。离散元、块体元和有限元、边界元、无穷元之间的结合又可提出若于种耦合方法，以发挥离散元和块体元在模拟不连续岩体方面的长处，并利用有限元、边界元、无穷元在模拟连续介质方面成熟的理论和计算技术，使应力分析、破坏、垮落及运动整个过程的数值模拟有可能实现。数值方法同解析方法或半解析方法的结合则是另一条可行的途径，这种结合的特点是用相应的数学推导得到更精确的(也更复杂了)单元模式，再通过离散化求解，解题效率及精度提高，不足之处是应用的局限性也随之增加。何翔[19]等人在Biot固结方程的基础上,引入介质渗透系数张量随应力张量的变化函数,建立能反映介质中应力场与渗流场非线性耦合作用的微分方程组,并在此基础上进行渗流场与应力场耦合问题的有限元求解,采用了精细积分方法进行时间离散。

参考文献:

[1] 谢和平,刘夕才,王金安,等. 关于21 世纪岩石力学发展战略的思考[J].岩土工程学报,1996 , 18 (4) : 98 - 102.

[2] 王泳嘉. 边界元法及其在岩石力学中的应用[J].东北工学院学报, 1984 (1) : 1 - 12.

[3] 刘红岩,秦四清. 层状岩石边坡倾倒破坏过程的数值流形方法模拟[J].水文地质工程地质, 2006 (5) :22 - 24.

[4] 王泳嘉, 冯夏庭. 关于计算岩石力学发展的几点思考[J].岩土工程学报. 1996 ,18 (4) :103 - 104.

[5] 刁心宏, 冯夏庭, 杨成祥,等. 岩石工程中数值模拟的关键问题及其发展[J].金属矿山, 1999 (6) :5 - 7.

[6] 周维坦. 岩体力学数值计算方法的现状与展望[J] . 岩石力学与工程学报, 1993 , 12 (1) :84 - 88.

数值方法范文第2篇

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

中图分类号：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.

[49]. 李鸿晶，王通，廖旭. 关于 Newmark-β 法机理的一种解释[J]. 地震工程与工程振动. 2011， 31（2）： 55-62.

数值方法范文第3篇

关键词:函数;值域;方法

中图分类号:G623.5文献标识码:A

函数是中学数学重要的基本概念之一,它与代数式、方程、不等式、三角函数等内容密切联系,应用十分广泛。函数的值域是函数的一个重要组成部分,值域是由定义域和对应法则所确定的。在研究函数值域时,不但要重视对应法则作用,而且要特别注意定义域对值域的制约作用,在初等数学的范围内,求函数的值域是没有通用方法的,它不象定义域有一定可依据的法则和程序,因此要根据问题的不同特点,综合而灵活地运用各种方法求之。下面列举几种函数值域求解的常用方法。

一、观察法

对于一些简单函数,可通过对于函数定义域及对应法则的观察分析求值域。

例1:y=1-|x|

解:|x|≥0 y≤1 故所求的值域为(-∞,1〕

例2:求函数 的值域

解:

故值域为(-∞,1)∪(1,+∞)

注:若观察出函数是单调函数,则可利用函数单调性求值域。

二、判别式法

用判别式求(a1,b1,c1,a2,b2,c2皆为常数,a1、a2不同时为零)型函数值域时分两种具体情况:

1.当分子与分母有公因式时,这时,可先约去公因式后再求值域,但须除去使公因式为零的x所对应的y值。

例3:求函数的值域

解:由

可知y≠1/2,同时因x≠2,故y≠3/2

所求值域为{y|y∈R且y≠1/2,y≠3/2｝

2.当分子与分母无公因式时,使用判别式求值域时,先转化为含参量y的一元二次方程,但要从判别式求出的结果中除去关于x的一元二次方程的二次项系数为零且又使方程无实数解的y值。同时要注意弄准函数定义域,并要检验边界点能否达到,否则可能得到错解。

例4:求函数的值域

解:经检验分子、分母无公因式,把原式变形:

(y-2)x2+(y-2)x+y-3=0①

显然y≠2(若y=2,则方程①为-1=0不成立,y≠2)

y∈R方程①有实根的充要条件

y≠2

{=(y-2)2-4(y-2)(y-3)≥0,得2y≤10/3

故所求值域为(2,10/3〕

三、配方法

适用于求二次函数和与二次函数有关的函数值域

例5:求函数y=x2+4x+3(X∈〔-1,0〕的值域

解:y=(x+2)2-1-1≤x≤0

当x=-1时,ymin=0,当x=0时,ymax=3

故所求值域为〔0,3〕

注:用配方法可求二次函数在指定区间上值域时,切勿直接套二次函数的最值公式,因为这时最值未必在顶点处取得。

例6:求函数y=sin2x+cosx+1=1-cos2x+cosx+1

=-cos2x+cosx+2=-(cosx-1/2)2+9/4

-1≤cosx≤1当cosx=1/2时,ymax=9/4

当cosx=-1时,ymin=0

y∈〔0,9/4〕

四、反函数法

利用反函数定义域求原函数值域

例7:求函数 值域

解:函数的反函数为y=log2x/(1-x)且定义域为

(0,1)故函数值域为y∈(0,1)。

五、换元法

1.形如y=x+b± ,(、b、c、d皆为常数,且、c不为零)。

例8:求函数y=2x-3+值域

解:由4x-13≥0得x∈〔13/4,+∞〕

令t= (t≥0)则x=

于是y=2 -3+t=2/1(t+1)2+3≥7/2

即所求函数值域为〔7/2,+∞)

2.三角代换

3.如(,b,c,d均为常数,且c≠0)

①>0,c

③>0,c>0,与④

例9:求函数 值域

解:设

则:2s2-t2=11①

s-t=y②

函数值域转化为方程①与②在直角坐标平面内有公共点时,y取值范围。

方程①可化为( t≥0)

此方程图象是双曲线,在第一象限部分,方程②的图象是斜率为1,在t轴上的截距为-y的直线,由图可知:

故所求函数值域为(-∞,〕

注:本题例12方法也可以解决函数

f(x)=的值域问题,只是转化为椭圆与直线有公共点时,直线在t轴上的截距问题。

六、不等式法

利用某些重要不等式,结合等号成立条件,得出函数值域。

例10:求函数值域

解:此函数定义域为(1,+∞)设

μ=

=(x>1,x-1>0)

上式当且仅当 即x=3取等号,即y=log2μ≥2

故所求函数值域为〔2,+∞)

七、最值法

对闭区间〔a,b〕上的连续函数,可求出y=f(x)在区间〔a,b〕内极值,并与边界值f(a)、f(b)作比较,求出函数最值,可得到函数值域

八、单调法

确定函数在定义域(或某个定义域的子集)上的单调性,求出函数值域。

九、数形结合法

画出函数图象,利用图象直观得出函数值域。

例11:求函数y=|x-3|-|x+1|值域

解:设y=|x-3|-|x+1|

-4(x≥3)

=2-2x(-1≤x≤3)

4(x≤-1)

的图象如右图

故值域为〔-4,4〕

总之,在具体求某个函数值域时,首先要仔细认真观察其题型特征,然后再选择适当方法,做题就会得心应手,取得事半功倍的效果。

参考文献:

[1]厉以宁,秦宛顺.现代西方经济学概论[M].北京:北京大学出版社,1983

[2]黎诣远.西方经济学(上册)[M].北京:清华大学出版社

[3]李种,现代西方微观经济学概论[J].北京:高等教育出版社,1995

[4][美]H•范里安.[M].上海:上海三联书店,上海人民出版社,1994

[5]朱建中,高汝熹.数理经济学[Ml武汉:武汉大学出版社,1992

作者简介:

数值方法范文第4篇

关键词： 函数 函数值域 方法

1.观察法

对于一些简单的函数，可在定义域及函数对应关系基础上确定函数的值域，这叫观察法。

由于函数值域是对应于函数定义域的函数值集合，因此首先要考察函数结构。在此基础上，从定义域出发，逐步推断出函数的值域。

例1：求函数y=（x-3）的值域。

解：函数定义域为-1≤x＜1，又≥0，x-3＜0，y≤0，即函数值域y∈（-∞，0］。

2.反函数法

如果函数在定义域内存在反函数，而求函数值域又不易求解时，可在通过求反函数的定义域的过程中而使问题获解，叫反函数求函数值域的方法。

即由y=f（x），反解出求函数x=f（x），原函数值域包含在f（y）的定义域中。然后分析二者的关系以确定函数值域。此法的成功取决于反解成立，分析正确，并注意在反解过程中保持同解性。

例2：求函数y=+，x∈（0，1］的值域。

错解一：y=+≥2，函数值域y∈［2，+∞）。

剖析：当x=（0，+∞］时，结论x=［2，+∞）才是正确的。但当x∈（0，1），这个结论就不可靠了。

错解二：y=+?圳x-2yx+4=0，

x∈R，4y-16≥0，解得y≤-2或y≥2。

函数值域为（-∞，-2］∪［2，+∞）。

剖析：以上求出的结果，只能是x∈（-∞，0）∪（0，+∞）时函数的值域，解法二同样忽略了0≤x≤1了这一限制条件，而x∈（0，1］的值域用“判别式法”是无法解决的。

正解：（反函数法）y=+?圳x-2yx+4=0，

x∈（0，1］，y≥2，y+≥2（1），方程（1）的根只能是x=y-，由0＜y-≤1，解得y≥，函数值域为［，+∞）。

3.转化法

利用已知值域的函数或所给函数的定义域，作为“媒介”，将待求值域的函数式变形。通过适当的运算，求得所给函数的值域。将所求函数值域问题转化为熟知的基本初等函数的值域问题，常能化难为易。

例3：求函数y=的值域。

解：由函数表达式得：2sinx+ycosx=3-y?圳sin（x+θ）=3-y，其中θ由sinθ和cosθ=确定。

|sin（x+θ）|≤1，（）≥（3-y）?圳y≥，即原函数值域y∈［，+∞）。

4.不等式法

运用不等式的性质，特别是含等量的不等式，分析等号成立的条件，以确定函数值域，叫不等式求函数值域的方法。

例4：已知α∈（0，π），求函数y=sinα+的值域。

错解：α∈（0，π），sinα＞0，＞0，sinα+≥2=2，函数值域为［2，+∞）。

剖析：由于忽略了“当且仅当sinα+时上式才能取等号”，但因|sinα|≤1故sinα≠，因此上式不能取等号，至少应有y≠2。

正解：α∈（0，π），sinα＞0，＞0，sinα+=sinα++≥3≥3。

当且仅当sinα=，即sinα=1时，上式能全取等号。

小结：用“不等式法”求函数值域，主要是利用“几个正数的算术平均值不小于其几何平均值”，但须注意取等号时条件是否能得到满足。

5.最值法

由于初等函数在其定义域内是连续的，所以我们可以通过求函数在定义区间内的最大值，最小值的办法，并求函数的值域。

例5：求函数y=的值域。

解：由函数定义域知，cosx∈［-1，-）∪（-，1］。

（1）当cosx∈［-1，-）时，y=x+=1-（-1），（）=-1，注意到cosx?邛（-），y?邛-∞-∞＜y≤-1。

（2）当cosx∈（-，1］时，（1+2cosx））=-1，（）=，注意到cosx?邛（-），y?邛+∞，≤y＜+∞。

故函数值域为（-∞，-1］∪［，+∞）.

一般二次函数的值域常用此法求解。有些高次整函数也可用此法。

6.判别式

根据一元二次方程ax+by+c=0有实根时，=b-4ac≥0。的性质，求函数值域的方法叫做判别式法。

例6：求函数y=2x-7x+3的值域。

解：2x-7x+3-y=0，且x∈R，=b-4ac=49-8（3-y）≥0，y≥，该函数值域为［，+∞）.

此法可用于行如：y=（A，P不同时为零，分子分母无公因式）的函数的值域。但必须强调：（1）是既约公式；（2）验证端点值是否能取到；（3）整理成行如一元二次方程的形式后，若平方项系数含字母要讨论；（4）若定义域人为受限，则判别式法失效。

7.换元法

通过代数换元法或者三角函数换元法，把无理函数、指数函数、对数函数等超越函数转化为代数函数求函数值域的方法叫换元法。

例7：已知函数f（x）的值域是［，］，求y=f（x）+的值域。

解：f（x）∈［，］，≤f（x）≤，故≤≤。令t=，则t∈［，］。有f（x）=（1-t），y=g（t）=（1-t）+t=-（t-1）+1，由于g（t）在t∈［，］时单调递增

当t=，y=，当t=，y=，

y=f（x）+的值域是［，］.

8.图像法（数行结合法）

通过分析函数式的结构、定义域、单调性、奇偶性、极值等。确定若干有代表性的点，勾画出函数的大致图形，从而确定函数的值域。

例8：求函数y=|x-1|+x的值域。

解：原函数可以表达成：当x≤-1或x≥1，y=|x-1|+x=（x+2）-；当-1≤x≤1，y=|x-1|+x=-（x+）+。

作出函数图像（见图1）

由图像知函数值域为［-1，+∞）。

9.单调性法

利用函数单调性，先求出函数的单调区间，再求每个区间上函数的值域，最后取其并集即得函数值域。

例9：求y=x-的值域。

解：y=x和y=-均为单调增函数，

y=y+y=x-为增函数，由定义域x≤知y=，故y≤.

10.配方法

如果给定一个复合函数，y=f［g（x）］，若g（x）或f（x）可以视为一元二次多项式，则要用配方法求其函数值域。

例10：求y=x+的值域。

解：y=x+=1-（-1），在定义域x≤内，显然有（-1）≥0，y≤1，函数值域为（-∞，1］。

本文仅从求函数值域的十种常用方法谈起，在不同的文献中可能会有与本文有出入的其它不同的方法，但解法大致相同，如构造法、极限法、解析法、复数换元法、三角代换法、恒等变换法、有理化法等。当然，本论文求函数值域的方法不是一成不变的，应在多次解题过程中综合并灵活应用这几种方法。

参考文献：

［1］董艳梅，吴武琴.求函数值域的常用方法［J］.昆明冶金高等专科学校学报，1999，15，（2）：19-23.

［2］王英.求函数值域的技巧方法探讨［J］.南都学坛（自然科学报），2001，21，（3）：115-117.

［3］侯剑方.求函数值域的几种方法［J］.中学数学，2002，（3）：28-30.

［4］谭廷经.求函数值域的几种初等方法与常见错误剖析［J］.中学数学教学，1995，（3）：28-30.

［5］张秦.求函数值域的方法与技巧［J］.榆林高等专科学校学报，1997，7，（4）：46-49.

［6］林如恺，江杰.求函数值域的几种方法［J］.乐山师范高等师范专科学校学报，1999，（3）：100-103.

［7］王慧贤，张莉.求函数值域的几种方法［J］.白城师范高等师范专科学校，2001，15，（4）：40-42.

［8］纯刚.求函数值域的方法与技巧［J］.安顺师专学报（自然科学报），1996，（4）：53-60.

［9］赵振威.中学数学方法指导［M］.北京：科学出版社，1999：71-75.

［10］谭光宙.中学数学解题方法［M］.北京：北京师范大学出版社，2001：149-151.

数值方法范文第5篇

关键词：燃油箱；Fluent；VOF方法；数值仿真

中图分类号：U463 文献标识码：A 文章编号：1005-2550（2016）06-0093-05

Abstract： The car fuel sloshing in the process of a continuous shaking impact on the fuel tank， the fuel tank under the condition of insufficient oil， liquid sloshing of freedom will affect the stability of the oil supply， especially racing for high speed and acceleration of the pursuit， if the oil supply is not stable， will affect the car driving， will ultimately affect the result of the match.Based on the VOF method， using the Fluent software to a small fuel tank car， using two different kinds of analysis model， namely the standard turbulence model simulation， Laminar layers model simulation. To lower oil amount limit to numerical simulation of stability evaluation of oil and fuel tank design， also provide guidance for future fuel tank design improvement.

Key Words： fuel tank； Fluent； VOF method； numerical simulation

引 言

燃油箱是汽车的油液储备装置，是汽车的动力源泉，为汽车的行驶持续供油，保证汽车正常行驶。在赛车比赛中燃油箱的作用十分重要，燃油箱的稳定供油为赛车提供源源不断动力，保证赛车的高速度和正常行驶。在高速行驶中，燃油箱中的油液并不能一直维持充满的状态，即油液会出现自由液面。由于赛车不能始终保持匀速直线行驶，所以当赛车出现加速度时，其内部的油液会随着赛车速度的变化而晃荡。由于赛车速度变化的不稳定性，油液的运动也会变得十分复杂，而且还会与容器相互冲撞产生作用力[1]。液体处于晃动状态下时，进行的是一种非线性的运动模式，到目前为止还没有理论对其做出具体的解析，所以对其研究的方法主要有实验研究和数值仿真。

本文应用Fluent流体分析软件，对某款小型赛车两种结构不同的燃油箱，进行极限转弯工况下的油液晃动情况的数值仿真。在做数值仿真分析中，采用VOF（Volume of Fluid）[2] 方法处理自由液面，无论是油液小幅度波动产生比较光滑的液面，还是在赛车突然加速后油液与油箱壁碰撞后产生破碎形成的液滴，VOF 方法都能很好地将这些自由液面的动态变化情况展示出来。通过仿真后处理，最终选出结构更合理的燃油箱。

1 VOF理论

早在1974年Debar[3]就采用VOF 方法来解决自由液面问题。VOF模型是一个流体函数，这个函数是目标流体的体积和网格的比值。这种方法占内存小，是一种简单而有效的方法[8]。多种不能混合的流体可以通过VOF模型对目标流体的动量方程进行求解和计算出目标流体通过某一区域的体积分数进行模拟。具体的应用包括预测、射流破碎和气液界面的稳态和瞬态处理[4]。

当VOF模型中包含有两种或两种以上的流体互相没有穿插，那么模型中每增加一个相就会新增一个变量，这个变量用来计算单元里的相的体积分数。当然在每个模型中所有相的体积分数和必定为1。各相共同享有所有变量和所有变量的属性区域，这些变量和它的属性区域代表体积平均值。如果某一相流体的体积分数值是已知的，那么，这些体积分数值就直接决定了，给定单元内的变量和它的属性是其中一相的还是多相混合之后的，这取决于体积分数值。换种说法就是：在给定的单元中，如果第q相流体的体积分数为αq，那么就可能出现三个情况：

① αq=0：单元中没有第q相流体；

② 0

③ αq=1：单元被第q相流体充满。

1.1 体积分数方程

在VOF模型中，通过求解一相或多相的体积分数的连续方程来确定跟踪相与相之间的界面。对于第q相，这个方程如式（1）：

1.2 动量方程

将整个区域内的单一的动量方程所求出的速度场作为各相共享数据。属性 和 的所有相的体积分数决定了这个动量方程表达式。

1.3 能量方程

式中Eq是基于每一相各自的共享温度和比热。属性 和 （有效热传导）被各相共享，源项 包含辐射的影响，同时也存在其他物体的热源。它和速度场一样，当相间存在较的大温度差时，接近界面的温度的精确度也会受到一定的限制。

2 仿真模型理论

1972年，Spalding和Launder提出标准湍流模型[9]，模型为半经验公式，其主要是求湍流耗散率 输运方程和解湍流动能 方程， 方程是由经验公式导出的方程，方程是精确方程，建立其它们与湍流涡粘系数 的关系。 模型要求流场是完全发展的湍流， 忽略流体分子之间的粘性，因而标准 模型只对完全湍流的流场有效。

标准 k-ε 模型同时考虑湍动速度比尺和湍动长度比尺的输运，因此通过标准 k-ε 模型能够确定各种复杂水流的长度比尺分布，比零方程模型和一般方程模型有了很大改进。同时，标准 k-ε 模型基本形式比较简单，能成功地预测许多剪切层型水流和回流[10] 。

3 仿真分析与结果

3.1 结构模型

该小型赛车两款不同结构A和结构B燃油箱透视图如图（1）所示，对该款小型赛车燃油箱进行结构设计，主要包括燃油箱容积选定、燃油箱外形设计、燃油箱出油口位置选取，防波板的设计等。由于油箱布置空间的限定，设计油箱容积约为6L，两款油箱外形相同；为了保证在出油口有较多的油量，同时在结构A和结构B燃油箱出油口处设计了台阶；在防波板的设计中，结构B从轻量化的角度考虑，在降低防波板高度的同时在防波板上开了较多的孔，且提高了出油隔间出油孔的高度。

3.2 数值仿真前处理

对该两款燃油箱进行仿真分析时，考虑到计算中求解的收敛性和稳定性，在保证体积不变，不影响分析结果的情况下，对实际模型进行了相

应的简化，去掉了注油管、小圆角等结构。将燃油箱模型导入ANSYS DM模块进行流体的抽取，在ANSYS Meshing进行网格划分，网格划分利用高级网格划分Use Advanced Size Function 为Proximity and Curvature；Relevance Center设为Fine，Smoothing 设为High，Min Size 设为0.1mm，此外在有些地方进行局部加密；边界层算法采用Program Controlled，边界层形式为Smooth Transition，Transition Ratio为0.272，一共3层，增长率为1.2；Meshing生成的网格数大致在70万左右，网格质量良好，网格划分模型如图（2）所示：

应用流体分析软件Fluent采用VOF方法，其中主相设为汽油相，副相为空气相，界面采用VOF方法捕捉，分析模型分别选择标准 湍流模型，外壁采用壁面边界条件，其速度-压力耦合方式选择Simple，离散方法选择Standard，松弛因子默认，设置剩余燃油量为1L，添加侧向加速度，大小为1g。设置收敛误差为 ，设置时间步为0.001s，时间步长为300s，开始迭代求解计算。

3.3 数值仿真后处理

求解计算完成后，进行数值仿真结果后处理，得到从0s到0.6s时刻的油液晃动云图，结构A燃油箱仿真后处理油液晃动云图如图（3）所示，结构B燃油箱仿真后处理油液晃动云图如图（4），从图上可以得到在低油量1L情况下，极限转弯工况下燃油箱油液运动状态。

对比图（3）和图（4）可以发现，结构B燃油箱出油口处保持了较充足的油量，保证了供油稳定性；同时降低防波板高度，在防波板上多开口，保证了小型赛车行驶过程中，燃油动量的快速转移，能够使得晃动的油箱迅速稳定，且结构B燃油箱质量轻于结构A燃油箱。

4 结论

本文应用计算流体力学软件Fluent，采用VOF方法，分别对小型赛车两款结构不同的燃油箱模型，在低油量极限转弯工况下燃油箱油液晃动过程进行了数值模拟，对比两种分析模型仿真结果，选出结构更合理的燃油箱。从分析后处理可以得到，结构B燃油箱结构在出油口处的燃油较多，优于结构A，且质量较轻于结构A。因此选择结构B燃油箱为该小型赛车燃油箱形式。

参考文献：

[1]龚国毅.基于 VOF 和浸入边界法的液舱晃荡的数值模拟[D].广东：华南理工大学，2013.

[2]Youngs D. Time-dependent multi-material flow with large fluid distortion [J]. Numerical methods for fluid dynamics. 1982， 24： 273-285.

[3]Debar R.Fundamentals of the Kraekn code.Tech.Rep.UCIR-760.Lawrence Liver more Nat.Lab.1974.

[4]江帆.Fluent高级应用与实例分析[M].北京：清华大学出版社，2008.

[5]Zhu X，Sui P C，Djilali N，Three-dimensional numerical simulation of water droplet dynamics in a PEMFC gas channel[J].Journal of power dynamics，2008，181（1）：1-15.

[6]ABE K，KORO K.A topology optimization approach using VOF method[J].Structural and multidisciplinary optimization，2006，31（6）：470-479.

[7]童亮，余罡，彭政，等.基于VOF模型与动网格技术的两相流耦合模拟[J].武汉理工大学：信息与管理工程版，2008，30（4）：525-528.

[8]张国军，闫云聚.基于VOF模型的导弹低速入水数值模拟方法[J].空军工程大学学报：自然科学版，2013，14（6）：24-26.

[9]S. Fu， B. E. Launder， M. A. Leschziner. Modeling strongly swirling recirculating jet flow with Reynolds-stress transport closures. In Sixth Symposium on Turbulent Shear Flows， Toulouse， France， 1987.

数值方法范文第6篇

关键词：数值计算方法；教学；实验；多媒体

数值计算方法是高等学校信息与计算科学专业的一门主要的专业课程，主要研究如何运用计算机近似求解数学问题的方法，逐渐成为数学与计算机科学的交叉性学科，既有纯数学的系统性和严密性特点，又与纯数学的研究重点不同。通过学习本课程，使学生理解并掌握科学计算中常用数值计算的基本原理和方法，并具备建立数学模型的基本技巧；训练学生熟练运用计算机编程语言实现各种数值计算方法；培养学生自行处理常规计算问题的能力和综合运用知识分析、解决问题的能力，达到理论与实践相统一。与数学专业不同，信息与计算科学专业应以实用性为出发点。结合近几年该门课程的教学情况，觉得有必要对这门课的教学内容进行更新、对教学方法进行改革和教学手段进行研制开发。

一、教学内容的更新

数值计算方法中数据的复杂性、算法的抽象性使初学学生感到无所适从，畏难情绪从第一堂课就开始了。如何充分利用有限的学时，在系统讲授数值计算方法的同时，让学生学会应用所学的算法去解决实际遇到的问题，从而使理论成果在实践中得以应用，并在实践中丰富理论算法，是这门课教学的基本方向。每一章基本算法的讲授都应先给出实例，从实例中引出问题，引导学生思考如何运用数学理论去解决问题，最后再给出算法，这样能够激发学生的学习兴趣。同时，数值计算方法该门课程的结构表面上感觉比较松散，实际上各个章节之间都有着密切的联系，所以要重视课程体系结构的讲授。如果没有一条主线贯穿始终，学生无法深层次地理解知识结构之间的联系。因此，在教学中要使各章之间保持一种紧密的联系，这样，学生的思路就会比较清晰，对知识的掌握更加扎实，实现由学习-应用-创新的进阶，并最终掌握科学计算的精髓。

二、教学方法的改革

在课堂讲授中应该遵循教学的主体是学生、主导是教师的原则，采取课堂讲解和提问题相结合，引导学生从实际问题出发，建立数学模型，在对模型的分析中结合以前章节中学习的数学思想，自己思考并动手推出相应的计算公式，而不是机械地记忆所学的数学公式，学生就不会觉得数值计算方法很枯燥。另外，针对数值计算方法课程内容过于抽象、难以理解的特点，采用直观式教学方法，将课堂板书和多媒体相结合。为了加深学生对基本概念和理论的理解，这部分内容以传统板书为主。而对于实例以及复杂公式的计算，应采用作图对比、幻灯片和动画进行演示和练习，来突出所要学的知识和已学知识的联系，以及所要学知识的几何直观性，从而节省时间，有利于培养学生的数学直觉，提高学习的积极性和主动性，提高学习效率。其次，针对数值计算课程抽象理论证明多的特点，尽可能多地从相关资料上收集最新的学科信息，寻找一些本学科在其他学科中应用的实例，引导学生思考问题，活跃课堂气氛，通过分组讨论的形式自由解决问题。同时，在讨论过程中，让学生深刻体会到数值计算方法的结果“没有最好，只有更好”，任何一个问题都没有现成的答案和方法，只有通过独立思考，反复实验比较，才能得到更好的计算结果。而且，不同讨论组所得到的结果会相差甚远，这样可以激发学生互相交流，比较方法的优劣，从而改进问题的求解。这种互动式的教学方法，注重课堂气氛的培养，既能激发学生学习的兴趣，又能使其对课堂内容实践化。

三、教学手段的研制开发

由于数值计算方法属于基础理论课程，在黑板上进行数学推理的过程同时也是学生消化理解知识的一个过程，因此内容的讲授还是主要以黑板为媒介。但是随着现代科学技术的发展，网络和多媒体技术在日常教学中的作用日益显现。在教学中，充分利用计算机和网络资源，通过计算机演示各种数值计算方法的运行结果，并对各种结果进行图形的比较，使得课堂教学环境更加形象生动，不仅大大增加了教学信息量，而且有效地激发学生的形象思维。为适应时代的发展，在教学中应精心制作相应的教学课件，提前准备课程中部分复杂的数学推导过程和计算框图，这样大大节省了在课堂上书写繁琐公式的时间，并且可以将主要精力集中在讲透基本概念、原理、技巧、算法设计与程序实现方面。同时，将一些重要步骤制作成多媒体动画，并配有清楚的文字说明和与图形变化对应的动态数据显示。此外，制作每次课的电子教案，突出教学的重难点，并且在复习时可以将课程内容贯穿在一起，更好地帮助学生理解课程的体系结构。

参考文献：

[1]王能超.数值分析简明教程[M].北京：高等教育出版社，2006.

[2]关治.数值分析学习指导[M].北京：清华大学出版社，2008.

数值方法范文第7篇

关键词：海啸 桥梁 多相流模型

2004年印度洋海啸和2011年东日本大地震引发海啸过后，海啸中桥梁的破坏引起大家的关注。关于2004年印度洋海啸的调查显示，仅在印度尼西亚的Sumatra岛，186座桥梁中的81座就被冲毁或严重破坏。2011年东日本大地震中，日本东北6县就有30座桥梁损坏，其中仅受灾最重的宫城县桥梁损坏就达23座，有些甚至被海啸带到上游的数十米处。海啸波为长周期，常见的海啸波周期为2-40min，波长达几千米至几百千米，当桥梁遭遇洪水或者海啸时，河流的水面高度远高于常规设计水位，此时桥梁极易遭到破坏，尤其对长跨桥梁破坏严重。Wardhana et al. （2003） 分析了自1989年至2000年美国出现的桥梁破坏原因，发现在这个期间大多数的桥梁破坏都是由洪水造成。如1933年，因为密西西比河和密苏里河洪水，爱荷华州有85座桥梁冲毁。

桥梁在遭遇洪水或者海啸时，桥梁的受力及破坏过程，可以由数值模型来研究。本文采用了可以模拟结构物运动的模型，研究桥梁在海啸波作用下的响应，讨论海啸波和桥梁破坏原因之间的关系。

f为追踪界面的体积分数，由VOF方法计算。利用SMAC（Simplified Marker And Cell）法求解N-S方程和连续性方程。

文章中海啸由溃坝生成，海啸冲击桥梁，在水平向和竖向都产生巨大冲击力，造成桥梁移动、破坏。文本利用浸入边界IB（Immersed Boundary）方法分析追踪桥梁面板的运动，并计算作用在桥梁上的水平拖曳力和竖向力。

计算区域概图及尺寸如图1所示。重力加速度g取9.81 m/s2，水的密度取9.97×102 kg/m3，空气的密度取1.18 kg/m3，水的运动粘度取8.93×10-7 m2/s，空气的运动粘度取1.54×10-5 m2/s，表面张力系数取7.20×10-2 N/m。

分析

图2给出海啸冲击桥梁的流速分布各时刻截图。图3给出水平力Fx、竖向力Fz、水平位移Δx、竖向位移Δz、倾角θ的时间历程曲线。为了验证数值模拟的可靠性，本文用经验公式计算的解析解作为对比。

结论

为了研究作用于桥梁上海啸力的特点，建立海啸-桥梁数值模型，模拟了海啸作用下桥梁的受力历时以及运动状态。数值模拟结果显示，长时间作用于桥梁上的水平拖曳力直接导致了桥梁的破坏，这与前人研究结果吻合。2004年印度洋海啸的研究结果显示冲刷桥梁面板的主要力为水平拖曳力，而上部结构与下部墩台结构间带有抗剪键和良好连接的桥梁在灾难中很少遭到损毁 。因此，可加强上部结构与下部墩台结构的构件连接以提高水平抗剪能力，可以有效的防止海啸对桥梁破坏。

参考文献：

[1] Unjoh S.，Endoh K.. Damage Investigation and the Preliminary Analyses of Bridge Damage Caused by the 2004 Indian Ocean Tsunami[C]//Proceedings of the 38th UJNR Jiont Panel Meeting， .

[2] 王培涛，于福江，赵联大，等. 2011年3月11日日本地震海啸越洋传播及对中国影响的数值分析_王培涛[J]. 地球物理学报， 2012， 55（9）： 3088-3096.

[3] 程远兵，李国斌，李靖. 海啸对沿海地区桥梁结构的破坏及抵御对策[J]. 华北水利水电学院学报， 2012， 33（5）： 21-23.

[4] Wardhana K.，Hadipriono F.-C.. Analysis of Recent Bridge Failures in the United States[J]. Journal of Performance of Constructed Facilities， 2003， 17（3）： 144-150.

[5] Cw Hirt，Bd Nichols. Nichols BD.Volume of fluid （VOF） method for the dynamics of free boundaries[J]. Journal of Computational Physics， 1981， 39： 201-225.

[6] Nakamura T，Yim SC. A Nonlinear Three-dimensional Coupled Fluid-sediment. Interaciton Model for Large Seabed Deformation[J]. Journal of Offshore Mechanics and Arctic Engineering， 2011， 133（3）： 14.

[7] Shoji G.，Hirakib Y.，Fujima K.，等. Evaluation of Tsunami Fluid Force Acting on a Bridge Deck Subjected to Breaker Bores[J]. Procedia Engineering， 2011， 14： 1079-1088.

数值方法范文第8篇

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]。

结束语

本文概要地叙述了数值优化方法在民机发展中的应用历史和现状；介绍了即将举行的空气动力优化设计计算系列研讨会；重点讨论了可应用于民机气动设计的基于高可信度CFD模型的数值优化方法，对其要求及其构造方法。以算例表明了这种方法不仅可用于传统圆筒机身+机翼的民机外形，也可用于非常规布局的民机外形（BWB飞机），且从所需计算机资源和计算时间看，可用于民机日常的工程设计。我们如能尽快掌握并应用这种方法无疑将大大缩短我国民机的研制周期和节约研制成本，有利于我们迅速赶超世界先进水平。相信世界范围的空气动力优化设计计算系列研讨会必将进一步推动数值优化方法的发展和应用。（本文图、表略）

数值方法范文第9篇

关键词： 函数 定义域 值域 值域的求解方法

设 是非空的数集，如果按照某种确定的对应关系 ，使对于集合 中的任意一个数 ，在集合 中都有唯一确定的数 和它对应，那么就称 为从集合 到集合 的一个函数，记作 ，其中 叫做自变量。 的取值范围 叫做函数的定义域；与 的值相对应的 的值叫做函数值，函数值的集合 叫做函数的值域

由函数的定义可知，一个函数的构成要素为：定义域、对应关系和值域。其中函数的值域是一个较复杂的问题，又是高中数学中的一个难点。总体来讲，求函数的至于要注意以下几点：（1）值域的概念，即与 的值相对应的函数值的集合 ；（2）函数的定义域。当题目中未明确给出函数的定义域时，应先求出函数的定义域，在定义域的范围内求函数的值域；（3）函数的单调性。求函数的值域时，常常借助函数的最值来求解，而求函数的最值时，对函数的单调性的讨论往往是必不可少的；（4）函数的解析式。在求函数的值域时，往往要根据所给函数的解析式的形式，使用相应的方法。具体常用的求函数值域的方法如下：

（1）观察法

对于一些简单的常见的函数，通过观察就可以求出其值域。例如我们熟悉的一次函数的定义域是 ，值域也是 ；反比例函数 的定义域为 ，值域为 。

（2）配方法（或公式法）

（3）换元法

（4）分离常数法

（5）利用函数的单调性求值域

例5. 求函数 的值域

解：由题可知函数的定义域为 ，因为 和 在 上均为增函数，故原函数为 上的增函数.所以 ，故原函数的值域为

（6）利用函数的最值求值域

对于区间上的连续函数，利用求函数最大值和最小值来求函数的值域。

总之，同学们在学习的过程中应多注意积累，善于总结，从而在求解函数值域的问题中，才能迅速找到求解此类问题的比较简单且合适的方法。

作者简介：

数值方法范文第10篇

论文摘要：在水利水电工程中，存在许多有自由面的无压渗流问题，自由面是渗流场特有的一个待定边界，这使得应用有限元法求解渗流场问题时，较之求解温度场和结构应力等问题更为复杂。归纳总结了无压渗流分析的各种数值计算方法，分析比较了其优缺点和适用条件，提出了无压渗流数值分析方法的发展趋势。

1引言

在许多水利工程中(如土石坝渗流、混凝土坝渗流、拱坝绕流、地下结构渗流等等)，都存在着无压渗流问题，这类问题的关键在于求解渗流场的边界，即确定事先不知道其位置的自由面和溢出面，属于非线性边界问题。求解该问题的有限元法以往采用移动网格法。虽然取得了许多成功的经验，但也表现出方法本身的缺陷。为解决上述问题，国内外学者致力于寻找有自由面渗流分析的新方法。其研究核心就是计算中不变网格，自Neumann于1973年提出用不变网格分析有自由面渗流的Galerkin法以来，出现了多种固定网格法，如剩余流量法、单元渗透矩阵调整法、初流量法、虚单元法和虚节点法等。

2无压渗流的数值分析方法

2.1调整网格法

调整网格法先根据经验假定渗流自由面的位置，然后把它作为一个计算边界，按照vn=0的边界条件进行分析，得出各结点水头H值后，再校核H=z是否已满足。如不满足，调整自由面和渗出点的位置，一般可令自由面的新坐标z等于刚才求出的H，然后再求解。

该方法原理简单，渗流自由面可以随着求解渗流场的迭代过程逐步稳定而自行形成，并且迭代是收敛的。但是，当初始自由面与最终自由面相差较大时，容易造成迭代中的网格畸形，甚至交错重叠；当渗流区内介质的渗流系数不均匀时，特别是有水平分层介质时，程序处理困难；对复杂结构问题，由计算机自动识别和执行网格移动几乎是不现实的。

2.2剩余流量[1]

剩余流量法通过不断求解流过自由面的法向流量（称为剩余流量）建立求解水头增量的线性代数方程组，达到修正全场水头和调整新的自由面位置的目的。迭代过程中只需一次形成总体渗透矩阵，但需要判断自由面被单元分割的各种情形，要求算出穿过单元的自由面被单元切割的面积及流过自由面的法向流速，计算工作量很大，难以推广到三维问题中。剩余流量法的全部调整均基于第一次有限元计算的结果，因而计算精度较差。

2.3单元渗透矩阵调整法[2]

单元渗透矩阵调整法利用对渗流场有限元计算的结果，根据单元结点水头与结点位置势的比较，把渗流场进行分区，各区的渗透系数给不同的值，通过不断调整单元渗透矩阵，模拟渗流不饱和区的作用，来确定出真实的渗流饱和区及渗流场。

该算法实际上是把边界不确定的非线性问题转化成了材料非线性问题来考虑。但是，单元渗透矩阵调整法对三维而言其计算效率是很低的，不能真实反映渗透区域的透水特性，计算精度和收敛稳定性都受到影响。

2.4初流量法[3]

初流量法利用高斯点的水头求出结点的初流量作为求解水头增量的右端项，避免了求自由面被切割的面积，同时避免了每次迭代中确定自由面的位置的做法，大大简化了剩余流量法的计算工作量。由于初流量法在计算跨自由面单元的结点初流量时，自由面以下的高斯点未予计算，计算精度受到影响。初流量法其收敛性不尽人意，解的稳定性不好。

2.5虚单元法[4]

虚单元法以上一次有限元计算的结点水头值为基础，求出自由面与单元边线的交点，移动跨自由面单元的某些结点，使之落于交点处，自由面将单元分成渗流实区和虚区。渗流虚区在下一次计算中退出计算区域，随着渗流计算区域向渗流实区逼近，结果也逼近问题的真解。该方法对三维复杂问题不适用，易产生结果收敛不稳定的现象。同时，虚单元法在处理有自由面穿越的单元时，结点移动路径的确定是比较困难的。

2.6虚节点法[5]

虚节点法以上一次有限元分析求得的节点势为基础，求出自由面和单元节线的交点，根据交点确定单元的积分区域，形成下一次分析的渗透矩阵。不同于虚单元法，虚节点法无需移动任何节点，因此不会出现网格畸形；虚节点法对网格不作改动，并能精确地描述跨越自由面单元的渗透矩阵，具有很好的精度和数值稳定性。

此外，无压渗流的数值分析方法还有边界单元法、流形单元法、无单元法等。

3无压渗流数值分析方法的比较

调整网格法计算原理简单，迭代过程稳定而自行形成，迭代过程收敛，但该算法对有复杂夹层和复杂排水系统的水工结构处理起来太困难，几乎不可能实现；另外对初始渗流自由面位置的假定要求也较高，如果初始位置与最终自由面位置相距甚远，则极易造成单元严重畸变，影响计算的精度；剩余流量法计算工作量很大，难以推广到三维问题中。初流量法在剩余流量法的基础上作了重大改进，大大简化了剩余流量法的计算工作量，但是收敛稳定性较差，而且由于两种算法的整个迭代过程依赖于第一次有限元计算的结果，精度受到一定的影响。单元渗透矩阵调整法对跨自由面单元按复合材料单元处理，复合材料单元渗透系数在复合面突变，其单元渗透矩阵不能代表这一特性，且矩阵主系数常不占优，因而计算精度和计算稳定性均受到影响。虚单元法对三维复杂问题不适用，易产生结果收敛不稳定的现象。虚节点法具有很好的精度和数值稳定性。

结论

本文归纳总结了各种无压渗流数值计算方法的原理及其优缺点，得到如下结论：

传统的调整网格法虽仍被使用，但由于自身的缺陷给应用带来诸多不便，因而正在逐渐被固定网格法所取代。具体选择计算方法时，应从问题的复杂度、收敛性及精度要求等方面加以考虑。现有的大型商用软件如ANSYS提供了良好的二次开发环境，用户可以通过二次开发，来实现无压渗流的数值分析。

参考文献

[1]DESAICS.Finiteelementresidualschemesforunconfirmedflow[J].IntNumMethodEng.1976,10(6):1415~1418.

[2]BATHEJN.Transmitmatrixmethodforseepagewithfreesurfaceproblem[J].IntJNumMethEngng,1983,(7):41~53.

[3]张有天,陈平,王镭.有自由面渗流分析的初流量法[J].水利学报,1988,(8):18~26.

[4]吴梦喜,张学勤.有自由面渗流分析的虚单元法[J].水利学报，1994,(8):67~71.