数值计算方法（精选12篇）
数值计算方法 篇1
准确预估空气动力升力和阻力是飞行器设计中一个至关重要的要求。在用N-S方程求解多段翼型流场后计算阻力时,即使是小迎角下的附体流动,且数值模拟的物面压力和摩擦力与实验数据吻合很好,但用物面压力和摩擦力分布积分求得的阻力却偏离实验值超过100%[1]。文献[1]对N-S方程的流场求解结果,分别用物面压力、摩擦力积分法、远场边界面积分法和尾迹动量损失积分法计算了单个翼型和多段翼型的阻力,发现尾迹积分法计算的单个翼型的阻力和实验吻合最好,计算的多段翼型阻力要比另两种方法准确得多,与实验数据的误差在5%以内,而另两种方法把阻力高估了50%多。
一般在二元风洞中测量翼型阻力都采用动量法,即在翼型下游某横截面安装尾迹排架(又名尾迹排管,简称尾排、尾耙),通过测量该截面的总压、静压获得该截面的动量损失,再对动量损失进行积分获得翼型阻力。这样,尾排安装位置不同可能会对阻力测量精度有一定的影响。如果尾排安装位置过于靠前,总压分布(或动量损失分布)曲线形状就会过于“尖瘦”,对其进行积分就会带来较大误差;反之,如果尾排安装位置过于靠后,虽然动量损失分布曲线足够“肥胖”光滑,但因损失范围扩大,尾迹区动量损失曲线的一部分可能会跑出尾排测点区域以外(尤其在较大迎角时更是如此)而不能测到,也对积分带来较大误差,影响阻力测量精度。所以尾迹排管的安装位置直接影响阻力测量的精度。而要通过在风洞中反复更换位置找到一个最佳安装位置,会对风洞洞壁的光滑性和结构完整性造成一定的损害,且反复试验也很费时费力,同时因为风洞中还安装有其它机构和测试设备,所以想通过这种方式确定一个最佳安装位置的思路是很不现实的。另一方面,计算流体力学和CFD软件的迅猛发展使得通过数值模拟获得整个风洞试验段流场所有物理量的全部信息成为可能,从而可以在任意不同的尾迹横截面进行积分求得阻力,然后进行分析对比,找到合适的安装位置。另外,对在大气中飞行的飞行器绕流场进行数值模拟时,也可以用尾迹横截面动量损失积分法计算阻力,并可与用物面载荷积分所得的阻力进行比较,为数值模拟的阻力计算方法提供更多的选择。文献[2]就低速风洞中多段翼型尾迹阻力测量方法的改进进行了数值计算研究。本文用数值模拟手段研究跨音速风洞中和数值解结果中翼型阻力的尾迹积分法的积分截面选取问题。
1 网格生成与数值方法
本文用Hilgenstock方法[3,4,5]和无限插值法生成计算用网格,用CFX软件进行流场计算。图1是NACA0012翼型在大气中飞行时的C网格。图2是翼型安装在试验段时风洞收缩段、试验段、驻室和扩散段的整体网格和局部网格。流场计算采用雷诺平均N-S方程及Menter的SST k-ω二方程湍流模型。
2 阻力计算方法
2.1 用翼型表面载荷积分法计算阻力
就一般三维情形来说,对物体表面力积分可获得流体对物体的作用力合力,合力在风轴系x轴上投影就可得到飞行器所受到的阻力[6],即
式(1)中,为物面单位外法向矢量,为应力张量,p为静压,为黏性应力张量。对翼型来讲,S=c×1=c,阻力系数为
式(2)中q∞为来流动压,c为翼型弦长。
2.2 用翼型尾迹动量损失积分法计算阻力
根据动量法,对尾迹横截面动量损失进行积分就得到翼型的阻力[7],即
式(3)中h为尾迹排架的高度,下标w表示尾迹横截面(wake)。阻力系数则为
式(4)中,cd wt为空风洞时风洞本身流动损失引起的阻力系数,Cf(y)为尾迹横截面无量纲单位长度上的无量纲动量损失,即[8]
其中,P0,w为尾迹横截面测量的总压值,Pw为尾迹横截面上的静压管测出的静压平均值,P∞为来流静压,P0为来流总压。γ为气体比热比。
3 数值模拟结果与分析
3.1 翼型在大气中飞行时的阻力系数计算
选NACA0012翼型(弦长c=1 m)在大气中飞行的情形,用图1所示的C网格,取静温T=288 K,Re=1.1×107,Ma=0.75,0.8,迎角α=0°,2°,4°,6°,8°进行计算。图3—图5给出了Re=1.1×107,Ma=0.75,α=0°,4°,8°三个典型迎角下计算得到的不同尾迹横截面上的无量纲动量损失分布曲线。从图可以看出,同一迎角下,随着截面往下游推进,动量损失的峰值逐渐减小,但损失波及的范围变宽;随迎角增大,尾迹范围增大。
图6～图10给出了Ma=0.75,Re=1.1×107,α=0°,2°,4°,6°,8°各迎角下不同横截面尾迹积分法得到的阻力系数。从图可以看出,大部分迎角下,随着积分截面向下游推进,阻力系数先减小,大约在后缘后5～6倍弦长后逐渐趋于一个稳定值(图中x为从翼型后缘点算起的平行于来流方向的坐标,c为翼型弦长),α=4°,6°时在x=7c以后才可能达到一个稳定值,这与文献[1]的结果一致,说明在外流场计算时要在后缘后7倍弦长以后的尾迹横截面进行积分才可能得到比较可靠的阻力值。
从图11尾迹积分得到的翼型阻力系数与表面载荷积分得到的阻力系数的比较。由图可知,尾迹积分阻力系数均小于表面载荷积分阻力系数,这与文献[1]的结果一致。从图还可看出,后缘后4倍弦长后尾迹积分阻力系数变化较小。
Ma=0.8,Re=1.1×107各迎角下不同尾迹横截面上的无量纲动量损失分布曲线与Ma=0.75的结果类似,α=0°,6°的典型结果如图12、图13所示。
图14～图18给出了Ma=0.8,Re=1.1×107,α=0°,2°,4°,6°,8°各迎角下不同横截面尾迹积分法得到的阻力系数。阻力系数随积分截面向下游推进变化的趋势与Ma=0.75类似,只是达到稳定值的位置略有差别。大部分迎角下,x=4c～5c后阻力系数达到一个稳定值,α=6°时x=6c以后才可能达到一个稳定值。另外,和Ma=0.75类似,尾迹积分的阻力系数均小于表面载荷积分阻力系数(见图19)。
3.2 翼型在风洞中的阻力系数计算
选NACA0012翼型(弦长c=0.2 m)安装在跨音速风洞试验段情形的流场,用图2所示的多块网格进行计算。计算状态为Re=1.23×107,Ma=0.75,0.8,迎角α=0°,2°,4°,6°,8°。计算结果说明各状态下尾迹横截面动量损失分布曲线形状及其向下游推进的变化趋势,横截面尾迹积分阻力系数往下游推进的变化趋势与翼型在自由大气中飞行的结果很类似。
图20—图22给出了Re=1.23×107,Ma=0.75,α=0°,4°,8°三个典型迎角下计算的不同尾迹横截面上的无量纲动量损失分布曲线。类似地,同一迎角下,越往下游,动量损失的峰值越小,损失波及的范围越宽;随迎角增大,尾迹范围增大。
图23～图27给出了Ma=0.75,Re=1.23×107,α=0°,2°,4°,6°,8°各迎角下不同横截面尾迹积分法得到的阻力系数。图中显示,在各迎角下,随着积分截面向下游推进,阻力系数逐渐减小,在x=2c后达到一个稳定值,说明尾迹排架放在这个截面后测量阻力是比较可靠的。
从图28为尾迹积分阻力系数与表面载荷积分阻力系数的比较。与大气中的情形类似,尾迹积分阻力系数均小于表面载荷积分阻力系数。
Ma=0.8各迎角下的横截面无量纲动量损失分布曲线与Ma=0.75类似,图29,30为Ma=0.8,Re=1.23×107,α=0°,8°两典型迎角下横截面无量纲动量损失分布曲线。
图31～图35为Ma=0.8,Re=1.23×107,α=0°,2°,4°,6°,8°各迎角下不同横截面尾迹积分得到的阻力系数。与Ma=0.75类似,随着积分截面向下游推进,阻力系数逐渐减小,在x=2c后达到一个稳定值。图36为尾迹积分翼型阻力系数与表面载荷积分阻力系数的比较。与Ma=0.75类似,尾迹积分阻力系数均小于表面载荷积分阻力系数。
4 结论
通过对翼型在大气中飞行绕流和在跨音速风洞中绕流的数值模拟结果的分析可以得出初步的结论:(1)尾迹积分法得到的阻力系数比表面载荷积分法得到的阻力系数小;(2)如果在数值模拟翼型在大气中自由飞行情形的绕流时拟用尾迹横截面动量损失积分计算阻力值,则建议选择在翼型后缘下游7c以后某截面进行积分;(3)在跨音速风洞中用尾迹排架测量翼型阻力时,建议尾迹排架安装在翼型后缘2c以后的截面。
参考文献
[1] Vinh H,van Dam C P,Yen D,et al.Drag prediction algorithms fornavier-stokes solutions about airfoils.AIAA paper 1995-1788-CP
[2]支真莉,焦予秦.翼型试验阻力测量方法的数据计算研究.科学技术与工程,2010;10(14):3384—3388
[3] Hilgenstock A.A fast method for the elliptic generation of three-di-mensional grids with full boundary control.Numerical Grid Genera-tion in Computational Fluid Mechanics'88,Pineridge Press Limited,1988:137—146
[4]张正科,罗时钧,李凤蔚.一种生成二维贴体与边界正交网格的方法.第七届全国计算流体力学会议论文集,1994:92—94
[5]张正科,庄逢甘,朱自强,等.两种椭圆型方程求源项方法在喷管内流场网格生成中的应用.推进技术,1997;18(2):95—97
[6] Batchelor G K.An introduction to fluid dynamics.Cambridge:Cam-bridge University Press,England,2000.沈清,贾复,译.北京:科学出版社,1997
[7] Anderson Jr D.Fundamentals of aerodynamics,4th ed.,New York:McGraw-Hill,Inc.,2007
[8]范洁川.风洞试验手册.北京:航空工业出版社,2002
数值计算方法 篇2
一、引言
数学是科学之母。一门学科之是否成为科学，决定于该学科的问题描述是否能化归为数学。工程技术属于应用科学范畴，工程技术问题通过建立数学模型与数学产生直接联系。数学问题的分析解通常是极难得到的，因此必须归结为数值计算问题。例如：人造飞船的轨道研究、汽车耐撞性问题研究、大型桥梁设计、天气预报等都必须数值求解。数值计算方法作为研究数学问题的近似求解方法的课程，既有一般类数学课程理论上的抽象性和严谨性，又有工科类课程的实用性和实验性特征，是一门理论性和实践性都很强的学科。该课程理论涉及面广、方法应用性强、内容丰富，再加上随着计算机技术的飞速发展，优秀数学软件层出不穷，数值计算方法更能与计算机相结合，适应科学发展的需要，现已成为各高校大部分理工科专业的必修课程。在数值计算方法的教学过程中，笔者发现了很多问题。本文对其中的部分问题进行了分析，并提出了几点教学改革建议。
二、教学过程中存在的问题
以笔者所在的机械工程专业为例，起初该课程为学科选修课，选课学生少，且其中大部分是为了凑学分而来的，学习兴趣不高在所难免。后来学科培养计划改变，该课程归入专业必修课，选课学生数量增加了，但是学习热情还是不高。究其原因，主要有以下几点：
1.课程对数学基础要求较高。本课程主要解决以下几大类问题：非线性方程求根、线性代数方程组求解、矩阵特征值与特征向量的数值解法、插值与拟合、函数最佳逼近、数值微分与积分、常微分方程初值问题的求解等。需要先修的数学课程包括高等数学、线性代数等。学生只有掌握这些课程中的基本内容，才能学好数值计算方法课程。而这几门课程均是难度较大的数学课程，学生的掌握程度本来就不好，甚至学过后已经忘记。由于同时要学习其他机械专业课程，学生不愿再花大量的时间和精力去学习或复习相关的数学知识，特别是本来就对数学不感兴趣的学生。所以在课程学习中，学生就会陷入“听不懂，听不懂就没兴趣，没兴趣就不想听课，不听课就不懂”这样一个死循环。
2.教学课时的限制。该课程的内容覆盖很广，如“插值方法”这一章，就算法而言就有Lagrange插值、Aitken插值、Nevile插值、差分与差商形式Newton插值、Hermite插值、分段低次插值、三次样条插值、B-样条插值等。然而，总学时设置仅为32学时。即使不面面俱到，挑选几种典型的插值方法讲解，也需要花费不少时间。因此，教师的课堂时间主要用来讲解各问题相关算法的理论推导和算法设计，几乎没有帮学生回顾相关数学知识的时间，且在课堂内也没有时间及时将理论运用于具体问题。学生觉得这是一门纯数学课，枯燥无味又难懂，没有学习兴趣。
3.没有与计算机很好结合。数值计算方法的一大特点是面向计算机。一个好的数值算法要通过程序设计在计算机上实现，要求用最简练的语言、最快的速度、最少的存储空间来实现某种计算结果。要达到上述要求，要求教师和学生既要掌握数值计算算法，又要能熟悉并熟练使用计算机语言。而现在的课堂教学重点大都侧重在理论讲解上，没有很好地结合计算机编程，学生把这门课当成了数学课来上;且学生在课外也没有将课堂上学到的算法付诸于计算机上实现。导致该门课程理论与实践严重脱节，达不到预期的教学效果和教学目的。
三、如何提高课程的趣味性
综合上述教学中出现的问题，要想教好这门课、使学生学到知识，最重要的是要提高课程本身的趣味性。“兴趣是最好的老师”，学生有了兴趣，才会有学习的热情，才会把精力付诸于课程学习上。那么关键问题是如何提高该课程的趣味性，主要从下面几点出发。
1.结合专业特点，从实际出发，合理安排课时和教学内容。由于课时有限，而且授课对象主要是机械这样的工科类专业的本科生，课程的主要目的是培养学生具有机械工程工作所需的数学能力，培养学生用数学的思想方法分析问题、解决问题的意识和能力，并为后续的工作和学习打下良好基础。那么教师在安排课时要懂得取舍，选择与机械专业紧密相关的内容讲解，课程主要浓缩为如下几个主要内容：非线性方程的求解、线性方程组的求解、插值与拟合、数值微分与积分、常微分方程数值求解。而每个内容应该针对其中的经典算法进行讲解。如非线性方程的求解，只需就二分法、简单迭代法、Newton迭代法进行详细讲解，其他算法如弦割法、简单Newton法等只需简单提及即可;常微分方程的数值解法，只需对Euler法和Runge-Kutta方法详细讲解，其他内容略讲即可。例如非线性方程求解中，判断迭代法收敛的充分条件，复杂的证明过程可以略去不讲。这样一来，教学课时和内容安排合理，整堂课就不会全在枯燥无味的数学定理推导中度过，即使数学基础不是很好的学生也能掌握并运用。而且学生能运用定理判断一种迭代法是否收敛，本身也会获得一定程度的满足感和自信心，而学习兴趣也可以在这之上建立起来。
2.对学生的计算机编程能力要求。该课程研究的是几大类数学问题的数值算法，懂得算法之后，一定要结合计算机进行编程实现。但本门课程又不是专门的计算机编程课程，且由于学时限制，课堂上不可能有多余的时间教授学生编程知识，因此该课程的先修课程还需要掌握一门计算机编程语言。现今的主流商用数学软件主要有如下几种：Matlab、Mathematica、Maple、MathCAD、C/C++、Fortran等，应选用一种熟悉或较易掌握的软件将各种算法进行计算机实现。另外，也可选用如Mathematica这类商用软件进行编程，该类软件界面简洁，语言简单，且功能也比较强大，自学便能很容易上手。
3.将数学理论与计算机相结合。在课堂上利用数学软件，绘制出直观的可视化图片，这样可以把课程中涉及的抽象原理、方法以及复杂的计算过程直观地呈现出来，使学生对相关算法有更直观和清楚的认识，更容易理解，同时也可加强学生运用数学软件进行编程计算的能力。如对非线性方程求根之前，先要找出有根区间，这时可以运用数学软件先画出函数曲线图，找出有根区间的大概位置，然后选择某一算法编程，观察根在迭代过程中的收敛性特征;又例如讲解最小二乘法拟合曲线时，可以运用数学软件将拟合出来的函数图与原函数表对比，可更加直观地理解插值和拟合函数中存在的误差。
4.课程中穿插实践环节，让学生参与到课堂中来。某一算法或某类问题讲解完后，应举出一些算例，让学生在课堂上分组讨论解决的办法，选择怎样的算法合适，怎么编程实现等。对于一些相对较简单的问题，可以请学生直接在课堂上编程求解并运行结果，然后一起讨论该结果的可靠性，或者对编程和运算过程中出现的问题怎么改正等。让所有的学生都参与到课堂中来，以提高学生学习的兴趣，而且同时能提高学生当场解决问题的能力、语言表达能力、计算机编程能力等。
5.课堂教学生动多样化。教学时应充分利用多媒体提高教学效果。如在PPT中增加声音、图像、动画等多种形式，在教学过程中形成多种感观刺激，使原来学生误解的枯燥、抽象的数学课直观化、形象化、生动化，充分激发学生的学习热情，从而大大地提高学生汲取知识的效率。另外，还可以将教学方式多样化，避免教师“满堂灌”、“唱独角戏”的尴尬局面出现。除教师讲解外，还可让学生一起参与到课堂中来，如分成小组讨论某个算法的优缺点，某个具体问题的解法，或采用小组竞赛模式，针对某一问题看谁的算法简洁、效率高、结果可靠等。
6.选择学生感兴趣的.算例。算例的选择应有特点，或与学生专业相关，或与学生感兴趣的事物相关，而不应该是单纯的数学习题，应联系相关的背景或出处。如对于车辆专业的学生，讲述曲线拟合这部分内容时，可以计算汽车车身外形曲线轮廓线为例讲述曲线拟合的过程，那么可先给出一些典型车型的外形轮廓图，然后针对某款车型，给出其轮廓线上某些型值点的数据表。学生在看到丰富多彩的汽车图时，首先会感到眼前一亮，兴趣马上会提高，课堂气氛也会得到活跃，而曲线拟合的知识也能很容易领会。
四、总结
数值计算方法 篇3
关键词：MATLAB软件 数值计算方法 辅助教学
中图分类号：G642 文献标识码：A 文章编号：1674-098X（2014）11（a）-0131-01
随着科技的飞速发展，各工程领域与数学的关系愈加密切，数学应用的广度和深度在现代科技发展中体现的愈加明显。数值计算方法作为利用计算机求解数学问题的学科，是实现实际工程问题的一种重要基础手段。因此，在大学教育阶段开设数值计算方法课程是非常必要的，而这不仅要求学生理解相关的数值计算的理论知识，还要会利用这些理论知识解决实际问题。基于长期的教学实践体会，在数值计算方法课程中做好理论传授和实践能力培养这两个环节变得异常重要。同时，随着科技的不断进步，与数值计算方法相关的软件层出不穷，如何合理的加以利用，是该课程教学过程中必须探讨的课题。该文以具体教学过程为例，介绍了数学软件MATLAB在提高课堂教学质量中的具体操作。
1 MATLAB介绍
MATLAB是由MathWorks公司1976年出品的软件系统，包含科学计算、可视化以及交互式程序设计等计算环境。它将数值分析、矩阵计算、科学数据可视化等诸多强大功能集成在一个易于使用的视窗环境中，为科学研究、工程设计等领域提供全套解决方案，代表了当今国际科学计算软件的先进水平。MATLAB的语法简单，编程易于实现，其强大的数值计算功能，基本涵盖了高等数学中的所有运算。经过多年发展，MATLAB已成为最优化理论，神经网络，计算机模拟仿真等现代学科的基本教学软件，是众多科研工作者的必备工具。
2 数值计算方法课程教学特点与难点分析
2.1 涉及范围广
数值计算方法是面向理工科各专业的基础课程，包括误差分析，插值法，数值微积分，矩阵计算，数值代数，微分方程数值解法等领域，涵盖大学数学的各分支，内容广泛。该课程具有知识结构分散、知识面跨度大、知识要点繁多等特点。因此，本门课程的讲授面临诸多困难，要想对每一种数值解法都做深入研究是不现实的，只能介绍部分经典方法的相关理论。如何在讲授完主要理论后将其应用于实践，是个大难题。
2.2 公式推导多
任意一本数值计算方法教材上的理论都过于复杂，给人的感觉就是这门课一直讲算法，传统的课堂上也以理论推导为主，如此很难有效的调动学生主动学习的积极性。加上课时有限，教师如果对课程不能宏观掌控，常常会在教学内容、方法、节奏等方面出现问题，在强调理论证明的同时，忽略学生对问题实际背景的理解以及数学思想的把握，造成教师对知识讲解的不透彻，学生消化不良。
2.3 计算量大
在解决实际问题时，个别简单问题可以进行少量手工计算。但是，为了很好的说明解决实际问题的效果，本课程一般都需要进行大量的重复计算，而在课堂上进行这种工作会严重影响课堂教学中的互动性。进而造成学生的抵触情绪，教学效果及学习效果差强人意。
3 基于MATLAB软件的数值计算方法课程教学
针对上述数值计算方法课程教学的特点和难点，我们考虑结合MATLAB软件的特点来改进现有的教学方法，将MATLAB软件应用于数值计算方法的教与学，必将会有良好的教学效果。主要做法如下。
3.1 基于MATLAB软件，分析与计算并重
整个教学内容既注重算法的理論分析，也注重算法的实现。对基础概念、基本理论、基本方法注重阐述来源和应用，删减不必要的、繁琐冗长的推导论证和复杂的运算技巧，确保课程内容通俗易懂，算法实用，够用。以具体案例和工程应用实例驱动学生运用数学方法解决实际问题，在此过程中确保理解数值计算方法的相关概念和方法、理论等。
3.2 基于MATLAB软件，经典与现代交融
教学内容在保持经典知识的基础上，加强内容的现代性。用现代数学的观点阐述一些数学概念，延伸数学结论。将现代信息技术和数值实验融入教学，并贯彻于教学全过程。例如，传统的微分方程数值解基本上都是采用差分法来完成，这种方法原理简单，学生容易接受，但数值解的精度较低或者需要较多的迭代次数。MATLAB软件中提供了全新的微分方程工具箱，对于常见的经典偏微分方程如热传导方程、扩散方程等都能给出精度足够的数值解，这对学生理解微分方程数值求解部分的理论是有很好助益的。
3.2 基于MATLAB软件，理论与实践结合
理论联系实际，课内课外相结合，利用习题课，给学生足够的可供选择的实用性较强的习题和数学建模问题，让学生亲历解决问题的全过程，注意融知识传授，能力培养于一体，目的是使学生得到选择算法、编写程序、分析数值结果，培养使用计算机进行科学计算和解决实际问题的能力，为以后从事现代数学科研工作和实践打下良好的基础。为此，在课程的讲授过程中，要注意引入工程实例，启发学生思考问题，引导学生利用现有知识探索解决问题的方法。
4 结语
数值计算方法面向算法，是利用计算机快速解决问题的一门学科，这一特点决定了教学中的授课模式，在理论教学的同时要注重与实践的结合。基于MATLAB的数值计算方法辅助教学，不仅增强了课堂教学的直观性，使枯燥难懂的理论知识易于接受，而且优化了课堂教学内容，改变了师生对课程固有的传统认识，能真正实现教与学的良性互动，让学生在应用数学解决实际问题的过程中感受数学的魅力和作用。因此，不能光讲方法而不实践，那样只会过于理论，让学生摸不着，看不到，很难理解数值计算方法的精髓，只有通过边学习边实践才能更好地掌握数值计算方法，并将其应用于工程实践。
参考文献
[1] 张玉柱，艾立群.钢铁冶金过程的数学解析与模拟[M].冶金工业出版社，1997.
[2] 张光辉，任敏.MATLAB平台上《数值分析》课程教学的几点思考[J].甘肃联合大学学报（自然科学版），2012，26（5）.
数值计算方法 篇4
圆环电流中心轴线及圆平面上磁场解析计算推导过程较简单[3—7], 而任意空间点磁场解析计算推导过程较复杂, 数学基础要求较高[3—13], 借鉴“割圆法”求圆面积的思想, 从毕奥-萨伐尔定律出发, 推导了一种简单的圆环电流周围空间任意点磁感应强度数值计算方法, 以便于求解空间点磁场计算机程序实现。
1 公式推导
建立如图1所示右手直角坐标系, 圆环电流位于XOY平面, 圆心O点为坐标系原点, 圆环电流所在空间为真空。设图1直角坐标系中有一点P, 圆环电流大小为I, 方向为逆时针方向, 圆环半径为R。圆环微小电流段Idl几何中心到P点位移为r, 则由毕奥-萨伐尔定律可得微小电流段产生的磁感应强度d B为[5,7,12]
则圆环电流在P点产生的总磁感应强度为
“割圆法”求圆面积时, 使用圆内接正N边形面积逼近圆的面积。如图2所示, 借鉴“割圆法”思想将电流圆环分为N等份, 当N→+∞时, 每一段圆弧电流近似等于微小直线电流段, 使用N段直线段电流元在空间某点产生的磁场总和逼近圆环电流在此点产生的磁场。设N段直线段电流元在P点产生的总磁感应强度为Bp, 则有
将式 (4) 和式 (5) 两向量代入式 (3) 中, 设Bp三分量为 (Bx, By, Bz) , 则式 (3) 可变为
在图 (1) 坐标系中, 给定P点坐标, 则由式 (6) 可以计算圆环电流周围空间任意P点的磁感应强度Bp=Bxex+Byey+Bzez, 其中ex、ey及ez为XYZ三轴向单位向量。
2 数值算例
如图2所示, 设圆环电流大小I为1 A, 圆环半径R大小为1 cm, 圆电流离散度N为1 000, μ0为真空磁导率。根据圆环电流空间点位置分布特征对式 (6) 简化, 分析圆环电流中轴线、圆电流平面、圆电流平面平行面及垂直面上的磁场分布特征。
2.1 圆环电流轴线上的磁场
在圆环电流轴线上, P点坐标为 (0, 0, Pz) , 式 (6) 可简化为
计算圆环电流中轴线上磁感应强度大小如图3所示。由图3可知, 在圆环电流中轴线上, 磁感应强度为Z轴正方向, 其大小在圆心处最大。
2.2 圆环电流平面的磁场
图4中, 圆环电流磁感应强度Bz分量在环平面上成圆对称分布, 过圆心的任意一条轨迹上圆环电流磁感应强度分布相同, 这是由圆环电流对称分布决定的。计算圆环电流平面上Py=0, -0.004
2.3 环平面平行面上的磁场
求电流环平面平行面上磁场, 即令式 (6) 中Pz等于一恒定值, 求平面上P点磁感应强度。计算平行面Pz=0.002, -0.004
由于平行面上磁感应强度分布的对称性, 平行面上过 (0, 0) 的磁感应强度大小剖面图相同。分别令Pz为0.001 m、0.002 m、0.004 m、0.01 m、0.02m, 做平行面上过 (0, 0) 的磁感应强度大小的剖面如图7所示。由图7可知, 圆环电流平行面上磁感应强度大小剖面曲线在距电流环较近时, 由于受到局部电流段磁场过大影响, 在对应电流环处会出现两个峰值。随着Pz的增加, 局部电流段磁场对剖面曲线的影响逐渐减小, 对应电流环处的两个峰值逐渐消失, 而平行面 (0, 0) 点因距圆环电流中心最近而逐渐出现磁场峰值。随着Pz的增大, 平行面相同 (Px, Py) 处的磁感应强度逐渐减小。
(Pz=0.002 m)
2.4 环平面垂直面上的磁场
为了研究的便利, 不妨取垂直于X轴的一组平面作为研究对象, 当Px分别为0.005 m和0.012 m时, 由式 (6) 计算该平面-0.004
图8中, 当垂直面到圆环电流圆心距离小于圆环半径时, 垂直面上无限趋近电流环区域磁感应强度大小会趋近于无穷, 而使得图8 (a) 中等值线出现两个峰值;图8 (b) 中由于垂直面上没有无限趋近于电流环的区域而只存在一个离电流环最近的点, 使得图8 (b) 中等值线只有一个峰值;在图8中, 垂直面上圆环电流外磁感应强度大小等值线成椭圆形分布, 在圆环电流近场区域, 等值线Z轴方向为长轴, 随着垂直面上点到圆环电流圆心距离的增加, 等值线椭圆形长轴逐渐由Z轴方向变为Y轴方向。这是由于垂直面上圆环电流近场区域磁感应强度受局部电流段影响大于全部电流段的影响, 而使得椭圆形等值线长轴方向沿着电流环径向;在垂直面上圆环电流远场区域, 结合图7结果可知, 局部电流段对圆环电流空间磁场的影响越来越小, 而在垂直面上距圆环电流轴线最近的平行线方向, 全部电流段所致磁感应强度综合影响越来越大, 使得椭圆形等值线长轴由圆环电流径向变为圆环电流轴向。
3 准确性分析
由对称性分析可知, 过原点的电流环平面垂直平面上磁场分布相同, 研究XOZ平面上点的磁感应强度数值计算准确性, 即能够说明式 (6) 在整个圆环电流周围空间磁场数值计算的准确性[8]。在XOZ平面上, 圆环电流空间磁场椭圆积分解析计算式为[3, 5, 8—12]
图9中, 横坐标为轨迹上P到原点的距离L, 纵轴为式 (6) 数值计算结果与式 (7) 解析计算结果之间的相对误差。由图9可知, 与式 (7) 解析计算结果相比, 式 (6) 数值计算圆环电流磁感应强度分量相对误差, 随着点P到原点距离的增加而逐渐增大, 但总体趋势上相对误差不大于0.004%。
圆环电流轴线上, 式 (6) 数值计算磁感应强度在N→+∞时, 与解析计算结果完全一致, 说明在圆环电流轴线上随着N增大, 数值计算越准确。同时, 随着N增大, 数值计算量也越大, 合理的圆环电流离散度是其空间磁感应强度数值计算准确度与计算量优化的关键。以式 (7) 计算结果作为准确值, 选取XOZ平面上三个不同位置点分析圆环电流离散度N分别取6至100时, 式 (6) 数值计算结果准确性与圆环电流离散度N的关系, 结果如图10所示。所选三点到原点的距离分别为0.005 m、0.03 m、0.2 m, 到原点的连线与Z轴夹角分别为由图10可知, 圆环电流磁感应强度分量数值计算结果的相对误差因空间点位置的不同而不同, 在离散度N取值较小时其相对误差受电流环离散度影响较大。随着离散度N增大, 所有空间点磁感应强度分量数值计算结果相对误差逐渐减小并趋于0, 结合图9结果, 在N=1 000时, 点处数值计算结果相对误差不大于0.001%。
4 结论
借鉴“割圆法”推导的一种圆环电流空间任意点磁感应强度数值计算方法, 与传统解析计算方法相比, 推导过程简单, 数学理论要求较低。基于该数值计算方法计算圆环电流中轴线、圆电流平面、圆电流平面平行面及垂直面上的磁场分布, 数值算例结果与已知结果及理论分析一致。在圆环电流离散度为1 000时, 圆环电流周围空间磁感应强度分量数值计算结果相对误差小于0.004%, 说明本文推导的圆环电流空间任意点磁感应强度数值计算方法准确度高。
摘要：圆环电流是最基本的理论磁体单元, 其周围空间任意点磁场解析计算推导过程复杂。借鉴“割圆法”原理, 推导了一种圆环电流周围空间任意点磁感应强度数值计算方法, 应用该数值计算方法进行了圆环电流中轴线、圆电流平面、圆电流平面平行面及垂直面上磁场特征数值分析, 并通过对比圆环电流完全椭圆积分磁感应强度解析计算结果验证了数值计算方法的准确性。研究结果表明圆环电流空间磁感应强度数值计算方法推导过程简单, 能够方便计算空间任意点磁感应强度且准确度高。
数值计算方法 篇5
优选数值计算方法用于不对称街道峡谷的研究
分别将迎风、混合及调和QUICK格式用于对称街道峡谷的汽车排放污染物浓度预测,并与风洞实验结果[1]相比较,研究表明数值模拟与风洞实验吻合较好,调和QUICK格式的预测效果最佳.将调和QUICK格式用于不对称街道峡谷的研究,发现当迎风面建筑物高于背风面建筑物时,峡谷内的流场和浓度场与对称街道峡谷相似,即峡谷内有一个顺时针方向的强漩涡,使得背风面污染物浓度高于迎风面污染物浓度.同时,随着迎风面建筑物高度的增加,漩涡中心的`位置成垂直向上分布,且与迎风面建筑物的高度基本成线性关系;当迎风面建筑物低于背风面建筑物且差别较大时,街道峡谷内出现了两个反时针方向旋转的强漩涡,使得峡谷内的流场与浓度场都发生较大的变化,迎风面污染物浓度会高于背风面污染物浓度.同时表明,峡谷上方的顺时针方向漩涡的中心呈抛物线向右上方向发展,而峡谷下方的逆时针方向旋涡的中心呈抛物线向左上方向发展.
作 者：汪立敏 王嘉松 黄震 谢晓敏 WANG Li-min WANG Jia-song HUANG Zhen XIE Xiao-min 作者单位：上海交通大学机械与动力工程学院,上海,30 刊 名：水动力学研究与进展A辑 ISTIC PKU英文刊名：JOURNAL OF HYDRODYNAMICS 年，卷(期)： 20(1) 分类号：X169 关键词：数值模拟 街道峡谷 差分格式 风洞实验某滑坡稳定性计算与数值模拟研究 篇6
摘要:文章以某滑坡为研究对象,分析了该滑坡的稳定性,并借助于FLAC数值模拟软件,对其变形中的应力和应变进行了分析,研究结果表明该滑坡处于欠稳定状态,需水后变形加大,因此需要工程治理。
关键词:滑坡;稳定性;数值模拟
中图分类号:P642 文献标识码:A文章编号:1006-8937(2009)20-0140-01
1工程概况
该滑坡为一相对平缓的凹形坡地,总体倾向西,坡形为折线型,滑坡中部地形坡度较缓,坡角在10～2°之间,中前部及前缘和后部地形坡度较陡,前部地形坡度在25~30°之间,前缘近河边部位多形成5～10 m的陡坎,后部地形坡度在17~22°之间。滑坡区北侧山体倾向沿渡河的坡地呈直线型,坡度35～40°,倾向滑坡区一侧山坡坡地,多形成2~5 m的岩质陡坎,地形坡度较陡,坡度在35～40°之间,滑坡区南侧坡体,坡面形态不很规则,地形坡度在17～30°之间。
研究区地层在区域上划属扬子区巴东利川小区,区内出露地层主要为三叠系和第四系,三叠系为一套碳酸盐岩和细碎屑岩互层构成;第四系有残坡积、崩坡积、冲洪积和滑坡堆积层。
研究区在大地构造上位于新华夏系一级隆起带的第三隆起带,以及淮阳山字型西翼反射弧和长江中下游区域东西向构造带交汇部位,区域上经历了晚元古代末的晋宁、侏罗纪末的燕山、三叠纪的印支和第三纪以来的喜马拉雅四次大的构造运动,随着中三叠世以后构造变动的加剧,形成了北西向、东西向、北东向、北北东向及北北西向等一系列不同方位、不同性质和不同特征的构造形迹,其中的东西向构造带是极为重要的构造体系之一,主导着研究区所在区域的基本构造构架,该滑坡就位于该构造带奉节复向斜的黄家河向斜南翼。
2稳定性计算
根据野外地质调查和宏观分析,滑坡可能的变形破坏模式为沿岩土界面滑移。其他坡段基本上不存在整体滑移的几何边界条件,不会发生整体的变形破坏,但存在表层风化破裂岩体的塌落、掉块等。
考虑滑坡区域可能遇到的各类情况,特别是最危险的情况,由于区内基本地震烈度为6°,可不考虑地震的影响,故综合确定以下两种计算工况:
①天然状况(坡体自重);②虚水到490。本滑坡坡的安全等级为三级,根据《建筑边坡工程技术规范》规定,采用平面滑动法和折线滑动法计算安全系数为1.25,采用瑞典圆弧法计算安全系数为1.20。由现场勘察和稳定性计算结果分析,可得出以下结论:现状条件下,滑坡整体处于稳定或基本稳定状态,仅局部位置出现表层岩体的剥落、崩塌等小规模的破坏现象。
②在工况二的条件下,稳定性系数最低为0.86,安全储备不足,需要进行加固处理。
从分析结果来看,目前坡体处于欠稳定状态,预测评价滑坡易发展并将失稳。
3滑坡应力与变形数值模拟
为了了解天然状态和蓄水后区内斜坡岩土体的应力分布特征及变形破坏特征,并进行分析评价,特选用岩土工程数值分析的通用软件FLAC 3D进行模拟分析。可以得到各种不同工况下的 应力和应变分布情况,部分成果见图1。
从位移等值线图可知,天然状态下的滑坡变形主要为水平向,主要分布在滑坡前缘,而蓄水后,斜坡水平位移急剧增大,但分布规律不变,为前缘大、后缘小,表明该滑坡在蓄水后可能在前缘出现蠕滑,进而诱发后缘的滑移变形。由此可知,滑坡前缘是潜在失稳的关键部位,在治理和监测过程中应重点考虑。
4结 语
研究表明,目前滑坡处于欠稳定和不稳定状态,滑坡易发展,特别是蓄水后,其潜在危害性很大。建议采用治理措施。根据气象预报和水雨情,对水库运行做到科学调度,精心指挥,使库水位变化尽量平稳,不出现大幅骤降和迅速抬升,以减小对滑坡稳定性的影响。
参考文献:
[1] 刘春原,朱济祥,郭抗美.工程地质学[M].北京:中国建材工业出版社,2000.
数值计算方法 篇7
针对《数值计算方法》课程教学作了一些尝试和改革,主要包括优选教学内容,并做适当合理的补充,重点建设实验课程,熟练掌握使用MATLAB软件,强化数值方法与计算机技术的应用能力训练,彻底改革作业形式,养成学生动手又动脑的良好学习习惯,将数学建模思想贯穿整个教学过程中,激发学生学习兴趣等措施。
1 选择合适的教学定位,优选理论教学内容
针对《数值计算方法》课程的基本要求和我校信息与计算科学专业学生的实际情况,我们选用中国科学技术大学张韵华等编著的《数值计算方法与算法》[1]作为教材,J H.Mahews等著,周露等译的《数值方法(MATLAB版)》(第四版)[4]一书作为参考书。考虑到我校信息与计算科学专业主要培养应用技术型人才,我们在保证课程体系完整,内容科学的前提下对教材进行了适当的处理,例如,我们授课中不包含Hermite插值,Gauss型积分,计算实对称矩阵特征值和特征向量的Jacobi方法和QR方法以及常微分方程的稳定性这些理论内容。这部分内容理论性较强,概念抽象不容易理解,我们建议基础好的同学自己钻研,其他同学会借助计算机软件用相应算法来解决问题即可。我们理论授课的内容包括插值,数值微分和数值积分,曲线拟合的最小二乘法,非线性方程求根,解线性方程组的直接法,解线性方程组的迭代法,计算矩阵的特征值和特征向量,常微分方程数值解。
2 适当补充教学内容,铺平解决实际问题的道路
随着计算机技术的飞速发展,数值计算方法也得到了蓬勃的发展,有些技术已经很完善,已经成为解决实际问题的有力工具,但是没有教科书介绍,造成理论与应用的脱节。例如,2011年全国大学生数学建模竞赛A题,给了某个地区的一些地点的重金属污染浓度数据,存储在EXCEL里面,题目需要我们给出重金属污染浓度的空间分布。要解决这一问题,除了需要二维插值的知识外,还需要知道使用MATLAB编写程序时如何调用EXCEL的数据。这是一个很实用的技巧,但是教科书中没有介绍。为此在讲一维插值的时候我们补充了如何使用MATLAB调用外部数据,处理数据以及二维插值等内容,而对于二维插值的理论我们不做介绍。又例如,在讲授多项式插值时,我们知道采用等距节点,高次插值会出现龙格(Runge)现象,解决Runge现象的一个传统方法是采用分段等距插值方法。我们补充介绍另外一种解决Runge现象的方法—Chebyshev插值。如考察函数f(x)=11+25x2,x∈[-1,1],取等距节点,分别构造5次插值多项式和10次插值多项式,会出现Runge现象,如图1所示。若取10次Chebyshev多项式插值,效果比等距10次插值多项式要好很多,如果取20次Chebyshev多项式插值,效果会更好,如图2所示。
3 重点建设实验课程,熟练使用MATLAB软件
数值计算方法课程是一门实践性很强的课程,各种算法最终是为解决实际问题服务的,所以我们更看重的是算法在计算机上运行的效果,为此我们增设了16个课时的实验课程。另外,《数值计算方法》的许多内容在理论和实践中都非常成熟,很多算法都已经被开发并集成到专门的数学软件,这些软件具有强大的数值计算功能,易学且具有开放性,其中最具代表性的就是MATLAB软件。
在实验课程里,我们使用MATLAB软件实现理论课中所有的算法。包括插值,数值微分,数值积分,曲线拟合的最小二乘法,非线性方程求根,解线性方程组的直接法,解线性方程组的迭代法,计算矩阵的特征值和特征向量,常微分方程数值解法等。另外我们还补充介绍MATLAB强大的图形展示功能,曲线拟合工具箱丰富的GUI界面以及非线性方程组求零点。实验课不仅提高学生解决实际问题的能力还能帮助学生加深对理论知识的理解。例如,考虑估算山崖高度的问题,如图3所示。在考虑了空气阻力,反应时间,回声传播时间等因素之后,引导学生建立了如下数学模型,其中h,t1,t2是未知数。
这是一个看似简单的三元非线性方程组,Newton迭代法数值求解它需要初值,有些同学虽然会使用MATLAB求解方程组零点,但因为初值选取不好,一直找不到解。从这个实验,加深了学生对Newton迭代法严重依赖初值的理解。
最重要的是我们将MATLAB软件介绍给学生,引导他们入门,激发他们自己学习的兴趣,鼓励他们自学MATLAB其他功能,熟练使用MATLAB解决各种计算问题。
4 改革作业模式,动手又动脑
作业是教学改革的重要部分,作业布置得不好会让学生更加讨厌这门课程,相反,作业布置得好可以激发学生更大的学习热情。《数值计算方法》教材和参考书都有很多题目可供学生练习,但是这些题目无论从形式上,还是从内容上都很陈旧,题目的答案也很容易找到,学生大多彼此抄袭,敷衍了事,根本达不到预期的作业效果。针对这一情况,我们设计了形式和内容都很新颖的作业题。
例如,数值积分部分的作业题是发给每个人一个形状不规则的卡片如图4所示,让他们分别用梯形公式,Simpson公式,复化梯形公式,复化Simpson公式计算其面积。作业最后以小论文的形式上交,作业内容包括设计算法,编写代码,图像展示数值结果,估计误差。由于每个人的卡片不同,坚决杜绝了作业抄袭的现象。另外由于形式新颖,且需要动手测量,极大的调动了学生的学习热情。
还有的作业题目需要学生合作完成,这时我们将学生分成若干组。例如,插值部分的作业题是观测一天中某一时段的气温变化,采集数据,然后估计某一时刻的气温。我们将学生分成六个组,每组配置一个气温表。第一组记录上午7,8,9,10,11点的气温,第二组记录7:10,8:10,9:10,10:10,11:10的气温,依此类推。这样每一组同学得到的数据都不相同,每组同学需要根据自己测得的数据分别用四次多项式插值,分段线性插值,三次样条插值估计10:15时刻的气温,每组上交一个小论文。为了弄清楚龙格现象,我们让第一组和第四组,第二组和第五组,第三组和第六组共享数据,用9次多项式插值来估计11:35的气温。
最小二乘拟合部分的作业是六个小组共享数据,每个学生用所有数据拟合三次多项式估计10:05的气温。
通过这样形式新颖的作业,极大调动了学生的学习热情,学生反响良好,得到了很好的教学和学习效果。
5 将数模思想贯穿整个教学,鼓励学生参加数模竞赛
《数值计算方法》课程理论性较强,背景知识较少,在授课过程中我们着重加强背景知识的介绍,精选教学实例,将数学建模思想贯穿到整个教学过程中,从提出问题,分析问题,建立模型,数值求解,结果展示,误差分析,力求完整的解决实际问题。另外,我们鼓励学生积极参加校内数学建模竞赛,网络挑战赛,全国大学生数学建模竞赛,美国大学生数学建模竞赛,建议每个学生毕业前都要至少参加一次数学建模竞赛。通过参加数学建模竞赛活动,学生更加认可了《数值计算方法》课程的重要地位,激发了学生的学习热情,有效地提高了学生解决问题的能力。
6 改革教学方法,更新教学模式
《数值计算方法》课程理论性较强,在教学过程中,我们采用启发式、讨论式等多种教学方法,营造良好的课堂气氛,加强师生之间的交流。由于《数值计算方法》课程涉及较多的概念、公式和定理,传统的教学方法,在算法推导、理论分析等方面能更好地引导学生去感受和思考数学逻辑的过程以及创造性的思维过程,加深对数学理论的理解和认识,培养学生的逻辑和思维能力。而在讲述背景知识,算法的应用,算法的程序实现的时候最好用多媒体课件进行演示。所以,我们认为需要将传统的教学方法和现代的教学手段结合起来,充分发挥各自的优势,在传统教学中穿插使用多媒体课件,根据教学内容选择合适的教学手段。
7 结束语
我们在《数值计算方法》课程教学改革方面作了以上的探索和尝试,但课程教学改革是一项艰巨的,长期的工程,我们仍然任重而道远。
摘要：针对一般本科师范院校的实际情况,结合当前社会需求和就业形势,探讨《数值计算方法》课程的教学改革,以期提高学生的动手实践能力和解决实际问题的能力。改革的内容主要包括准确的教学定位,优选教学内容,重点建设实验课程,熟练使用MATLAB软件,彻底改革作业形式以及将数学建模思想贯穿到教学过程等方面。
关键词：数值计算方法,教学改革,MATLAB,数学建模,作业改革
参考文献
[1]张韵华,奚梅成,陈效群.数值计算方法与算法[M].北京:科学出版社,2006.
[2]杨韧,张志让《微分方程数值解》课程教学改革与实践[J].大学数学,2011,27(4):19-22.
[3]张韵华,陈效群.数值计算方法课程改革初步[J].大学数学,2003,19(3)23-26
[4]Mahews J H.等著,周露等译.数值方法:MATLAB版[M].北京:电子工业出版社,2005.
汽油机进气歧管数值计算方法研究 篇8
进气歧管是发动机的重要部件之一, 是影响发动机进气阻力和各缸进气均匀性的关键因素[1,2]。而汽油的进气阻力和进气均匀性对整机的经济性、动力性和工作过程的可靠性有着重要的影响。随着计算机技术的飞速发展, 数值模拟分析技术也被引入到发动机进气歧管的研究中, 并发挥着愈来愈重要的作用。通过数值模拟计算, 研究者能够更直观地观察管道内气体的流场组织情况, 更合理地指导进气歧管的研究及优化工作[3]。
目前, 四缸进气歧管稳态CFD仿真主要存在两种数值计算方法: (1) “四口全通” (下文简称方法A) , 即4个支管全部设置成出口边界; (2) “三闭一通” (下文简称方法B) , 即只有一个支管设置出口边界, 其余3个支管设置成壁面。文献[4-6]采用方法A对进气歧管进行数值模拟, 分析其流场特性、进气阻力及进气均匀性, 进而对结构进行优化设计;而文献[7-10]则采用方法B对进气歧管进行数值模拟, 分析其流场特性等。
因此, 本研究以某四缸汽油机的两款进气歧管为研究对象, 研究方法A、B对进气歧管的进气阻力和进气均匀性计算结果的影响。此外, 还建立相应的发动机试验台架, 通过分析两款进气歧管对发动机的动力性、经济性及运行平稳性的试验数据来验证数值计算方法的可靠性。
1 模型建立
1.1 数学模型
发动机进气歧管曲线形状复杂, 进气过程为可压、粘性、非定常的三维湍流流动, 但它仍然满足质量、动量和能量守恒方程。各方程具体形式如下。
质量守恒方程:
式中:ρ—流体密度, t—时间, —流体速度矢量, Sm—质量源项。
动量守恒方程:
本研究采用RNG k-ε湍流模型。在RNG k-ε湍流模型中, k和ε是两个基本未知量, 与之相对应的输运方程为:
式中:Gk—由层流速度梯度而产生的湍流动能;Gb—浮力而产生的湍流动能;YM—由于在可压缩湍流中, 过渡的扩散产生的波动;C1ε, C2ε, C3ε—常量;αk, αε—k方程和ε方程的湍流Prandtl数;Sk, Sε—自定义的值。
1.2 数值模型
本研究采用UG建立了两款进气歧管的几何模型, 并采用Hypermesh软件对进气歧管的三维模型进行四面体网格的划分, 并在出口和入口处进行网格加密。该发动机两款进气歧管 (下文简称I、II型) 的网格模型如图1、图2所示。经网格无关性研究后, 最终划分的I、II型进气歧管网格模型的网格数分别为576 000个和803 000个。进出口边界条件均采用压力边界条件, 压力入口条件P=101 325 Pa;压力出口条件P=998 25 Pa, 壁面温度为固定温度293 K。
2 计算结果及分析
2.1 进气阻力分析
本研究采用方法A、B计算获得的Ⅰ、Ⅱ型歧管XOZ平面速度矢量迹线图如图3、图4所示。比较图3和图4, Ⅱ型歧管各支管的流速都较Ⅰ型歧管大, 也都未出现明显的回流现象。
本研究采用方法A、B计算获的Ⅰ、Ⅱ型歧管XOZ平面压力云图如图5、图6所示。如图5所示, Ⅱ型歧管的支管2、3在XOZ平面的末端平均压力较I型歧管大150 Pa, 支管1、4则只大于Ⅰ型歧管50 Pa;而如图6所示, Ⅱ型歧管4支管都较Ⅰ型大170 Pa左右。
采用2种数值计算方法获得的两型歧管总各个分支管出口的总压损失如表1、表2所示。由表中数据可知, II型歧管各支管的总压损失都较I型少。其中, 方法B获得的两型歧管各支管总压损失变化较均匀, 而方法A则中间支管2、3压损变化较大, 支管1、4压损变化较小。可见, 方法A获得的两型歧管进气阻力变化和方法B获得的趋势一致, 即II型歧管进气阻力相较于I型歧管更小, 但采用方法B获得的各支管进气阻力变化较均匀。
2.2 进气均匀性计算结果比较
发动机进气的最大不均匀度为σmax:
式中:Qmax—分支管最大出口质量流量, kg/s;Qmin—分支管最小出口质量流量, kg/s;Qm—平均质量流量, kg/s。
Ⅰ、Ⅱ型歧管采用方法A、B获得的各个分支管出口的质量流量如表3、表4所示。由表3可得, I型歧管各支管的最大的不均匀度σmax为2.44%, 而II型歧管中间两根支管的气流流量与两端支管相差很大, 各支管的最大不均匀度σmax达到了28.68%, 说明II型歧管的进气均匀性很差, 将会导致燃烧恶化、各缸工作不均匀等后果, 从而出现游车及转速波动数值较大等现象。从表4可知, I型歧管各支管的最大的不均匀度σmax为1.71%、II型歧管中各支管的最大不均匀度σmax仅为1.23%。由此可知, 不同数值计算方法的进气均匀性计算结果有很大的差异且两型歧管采用方法B得到的均匀度结果都不同程度的好于采用方法A计算获得的。
3 发动机性能试验
3.1 试验设备
台架试验系统结构如图7所示, II型歧管试验测试台架如图8所示。
1—空气滤清器;2, 9, 13, 18, 21, 22, 27—温度传感器;3, 10, 12, 19, 20, 24, 26—压力传感器;4—空气质量流量计;5—发动机;6—测功器;7—燃油箱;8—油耗仪;11, 25—冷却液体积流量计;14—冷却水箱;15—电子风扇;16—膨胀水箱;17—变频调速风机;23—原机水泵
3.2 试验结果及分析
本研究通过对两款进气歧管进行台架试验, 测试两款歧管的运行平稳性并获得两款进气歧管外特性下的扭矩、油耗数据。
I、II型进气歧管动力性 (功率) 对比图如图9所示。试验结果表明, 在1 800 r/min~3 600 r/min区间内, II型歧管功率平均提升4.21%。I、II型进气歧管经济性 (油耗率) 对比图如图10所示。从图中可以得到, 在1 800 r/min~3 600 r/min区间内, II型进气歧管相较I型, 油耗平均下降2.99%;除转速在1 800 r/min左右外, II型歧管其他转速下的外特性油耗均有下降。表明II型歧管相较于I型对发动机性能有较大的提
I、II型进气歧管在转速2 000 r/s下30 s内的速度波动对比图如图11所示。研究结果表明, 发动机工作时转速波动幅值都只接近6 r/s, 没有出现由于进气不均匀导致的转速波动数值较大现象, 说明两型歧管的进气均匀性都较好。由于发动机实际工作过程中是按进气行程、压缩行程、作功行程和排气行程的顺序不断循环反复的, 在同一时刻中只有一缸是处于进气行程。因此, 采用方法B得到的进气均匀性结果更加符合实际。
4 结束语
(1) 本研究采用“四口全通”和“三闭一通”这两种数值计算方法获得的进气阻力变化趋势一致且与发动机经济性、动力性试验结果吻合。因此, 两者均适用进气歧管流场特性分析。
(2) 不同数值计算方法的进气均匀性计算结果有很大的差异且采用“三闭一通”数值计算方法获得的均匀度不同程度的好于“四口全通”计算结果。
(3) 发动机运行平稳性试验结果表明, “三闭一通”高。由此可见, 采用方法A、B这两种数值计算方法分析获得的进气阻力变化趋势与实验结果一致。因此, 两者均适用于进气歧管流场特性的分析。数值计算方法更加适用进气歧管进气均匀性分析。
参考文献
[1]BRENNAN S L, KEE R J, KENNY R G, et al.A Theoretical and Experimental Study of Resonance in a High Performance Engine Intake System:Part2[N].SAE Paper, 2007-01-1399.
[2]HAMILTON L J, ROZICH J, COWART J.The Effects of Intake Geometry on SI Engine Performance[N].SAE Paper, 2009-01-0302.
[3]胡景彦, 苏圣, 洪进.某缸内直喷发动机进气歧管CFD模拟分析[J].液压气动与密封, 2009 (9) :25-28.
[4]张强, 李娜.沼气发动机进气均匀性数值分析[J].农业工程学报, 2010, 26 (8) :145-149.
[5]王晗, 蔡忆昔, 毛笑平.发动机进气系统不均匀性的三维数值模拟[J].小型内燃机与摩托车, 2007, 36 (3) :41-44.
[6]韩同群, 马祥宁.应用CFD/CAD技术对柴油机进气歧管进行优化设计[J].内燃机, 2007 (1) :13-16.
[7]许元默, 帅金石, 王建昕.电喷汽油机进气歧管的CAD/CFD设计[J].汽车工程, 2002, 24 (4) :314-321.
[8]张继春, 李兴虎, 杨建国, 等.壁面函数对进气歧管CFD计算结果的影响[J].农业机械学报, 2008, 39 (7) :47-50.
[9]江国华, 温苗苗.EQD180N-20型发动机进气不均匀性分析[J].武汉理工大学学报:交通科学与工程版, 2006, 30 (6) :1008-1011.
数值计算方法 篇9
1 曲面在数值模拟中的表示方法
三维空间中的一张二维曲面S, 我们可以用曲面上uv两个方向上的离散格点X (u, v) 来表示, 其中X (u, v) = (x (u, v) , y (u, v) , z (u, v) ) 。若uv两个方向上的阶数分别为m和n, 则该曲面由m×n个格点形成的网络来表示。如果将Xi (u, v) 和Xi (u, v) 分别定义为格点坐标在i (i=u, v) 方向上的一阶和二阶差分, 将Xuv (u, v) 定义为混合差分, 则曲面S的法向单位向量为n= (Xu×Xv) /||Xu×Xv||。曲面的平均曲率为曲面的两个主曲率之和, 对于使用上述方法表示的曲面, 其平均曲率可以表示为
对于图1 (a) 中所示的具有三维轴对称性的曲面, 若与对称轴z轴垂直的为曲面的v向, 曲面u向与v向正交, 则曲面与yz平面的交线是曲面内的一条u向曲线 (表示为y (z) 或z (y) ) , 由于对称性, 数值计算中只需记录这条曲线上格点的位置, 就能反映整张曲面的形状图1 (b) 。
2 三维轴对称情况下的曲率算法改进
进一步, 由于〈Xu, Xu〉=yu2+zu2, 〈Xv, Xv〉=xv2, 〈Xu, Xv〉=0, 则根据公式1) , 曲率可以表示为
接下来对计算精度进行分析。对于底面半径为1的圆柱曲面, 在数值模拟中使用我们改进得到的公式3) 进行计算, 得到的曲面曲率为1, 是完全准确的。而如果按照传统方法使用边型棱柱来拟合该圆柱, 并完全通过数值差分计算曲率, 得到的曲率如表1所示。可以看出才能使计算得到的曲率的误差小于10%, n≥32才能使误差小于1%。对于非圆柱曲面, 使用公式3) 计算得到的曲率的精度同样比传统方法有显著提高。使用单条u向曲线表示具有三维轴对称性的曲面, 还节省了数值计算所使用的存储空间和计算量。
3 结论
本文在uv方向离散格点表示的曲面的曲率计算方法基础上, 针对三维轴对称性曲面的特性, 使用特定的u向曲线来表示曲面, 并改进了曲面平均曲率的计算方法。改进的计算方法比完全通过数值差分计算的传统方法具有更高的精度, 同时还减小了计算量, 这对涉及三维轴对称性曲面曲率计算的数值模拟具有重要意义。本文使用的推导方法还可以进一步衍生出曲率梯度、曲率梯度的散度的改进算法。
摘要：三维轴对称模型是材料微观结构研究中的常见模型。本文以数值模拟中界面的曲率计算为基础, 对三维轴对称情况下的计算方法进行了改进, 提升了效率和精度。
关键词：三维轴对称,界面,曲率
参考文献
[1]沈志强, 沈耀, 蔡珣.纳米金属多层膜的微结构热稳定性研究[J].纳米科技, 2009.
[2]王奎武, 陈发来, 陈意云.基于点表示的曲面曲率计算方法[J].小型微型计算机系统, 2005.
数值计算方法 篇10
1 接地网雷电流冲击特性模拟模型计算方法
该计算方法:用傅里叶级数把雷电流分解成不同频率下幅值;利用电磁场理论计算出自互阻和自互感;利用节点电压法建立接地网模型;根据地网模型建立适用于编程的矩阵方程, 从而求解出频率分量的电流响应。
1.1 傅里叶分解雷电流
图1是中国电工部门推行的标准雷电流波形, 波头为tf, 波长时间为tt, 在波形上取两点, 分别是10%倍和90%倍的幅值, 以这两点做直线交与时间轴所得的值即是t0。实际工程中, 大多都采用直角波形代替标准波形, 基本可以满足实际需求。
图1所示的雷电流可采用双指数函数来表达:
式中:I0是雷电流幅值;t1和t2是待求的参数。文献[6]中利用闪电的三个特性得出t1和t2的值。求出t1、t2后, 通过傅里叶级数对式 (1) 进行转化, 使雷电流分解为各个频率下随时间变化的函数。
傅氏级数的转化公式如下:
式中:f (x) 是具体的函数;l是函数的周期, 根据实际情况选取。
图2显示的为30 kHz的频谱。幅值为1.4 kA, 波形为2.6/50μs的雷电流通过上述的傅里叶分解得到频域中的各个频率分量。
1.2 节点电压法在接地网中的应用
接地网假设导电材料形状为圆柱体, 且导体长度远远大于导体半径。计算时各种接地极形状均可等效成圆柱体。本文分析接地网模型是在假设埋设处于均匀土壤的情况下, 土壤的电阻率为ρ, 介电常数为ε=εr·ε0, 磁导率为μ0。
文献[7]中将接地网的接地体分成N段, N值越大计算结果越接近真实值, 但N取得过大, 不仅需要较大的计算机内存而且计算时间过长。文中提出多个假设:
1) 散流到土壤的电流都是从该导体的中点入地。
2) 用端点电压来求得支路电压。
3) 把接地网分成r段及n个节点。
图3是某个节点的等效电路, 每段导体k上有一个沿导体方向的轴向电流和一个散流到土壤的散流电流。
接地网中任意点注入电流F后, 地网中每个节点就会产生电位。等电位计算时, 假设整个地网节点电位是均等的。令节点j的电压为Vj (以无穷远点为参考点) 。定义第k段支路电压为[8]
式中:l和m为k支路的两个端点。所有支路可以写为矩阵形式:
式中:[U]为支路电压列向量;[V]为节点电压列向量;[K]是一个二维矩阵, 当支路i和支路j相连时, 矩阵中的元素kij=0.5, 否则为0。
每段导体上除了轴向电流分量外, 同时有泄露电流I散流到土壤中, 因此所有支路有
式中[G]为自互阻阵, 文献[8]中已经给出了计算方法。
将支路泄露电流I等分到与之相连的每个支路上有
式中:[J]为等效节点泄露电流列向量;[K]T为[K]的转置。
对于整个接地网, 可得到
式中[Y]为节点导纳矩阵, 文献[9]已经给出计算方法。
综合式 (2) —式 (6) 有
由式 (7) 中可以得出, 输入注入电流[F]就可以求出[I]、[V]和[U]。
2 计算实例与分析
程序运行过程中, 在工频情况下输入注入点电流, 得出的注入点电压即为接地电阻。计算冲击接地电阻时, 按照上述方法, 将不同频率分量下的电位经过叠加, 算出入地点的最大电位, 从而求出冲击接地阻抗。上述过程, 已经用MATLAB编程中得到实现。在编程过程中, 基本频率采用1 kHz, 最高频率为300 kHz。
将本文计算所得到的结果与运用CDEGS计算得到的结果进行比较。一接地网100×100 m2钢材料组成, 接地极直径为0.01 m, 埋深为0.6 m, 土壤电阻率为1000Ω·m, , 导体电阻率为1×10-7Ω·m, 真空磁导率4×π×10-7, 雷电流幅值为1 kA, 电流由角网孔中间入地。本程序运行结果为7.2592Ω, CDEGS为7.4571Ω。本文所编程序运行和CDEGS运行的电压震荡图如图4所示。
从图4可以得出, 由两者数值计算得到的结果非常接近。
3 结语
本文利用节点电压法, 建立接地网在雷电流下的接地模型, 并且给出了相应的解析表达式。利用傅里叶级数分解雷电流, 分析出频域特性, 从而求出需要的接地参数。将所得数据与CDEGS计算比较, 可以得出本文计算方法的正确性。
参考文献
[1]GUPTA B R, THAPAR B.Impulse impedance of grounding grids[J].IEEE Transaction on Power Apparatus and Systems, 1980, 99 (6) :2357-2362.
[2]MCLIOPOILOS A P, MOHARAM M G.Transients analysis of groun-ding systems[J].IEEE Transaction on Power Apparatus and systems, 1983, 102 (2) :389-399.
[3]MAZZETTIE C, VECA G M.Impulse behaviour of grounding elect-rodes[J].IEEE Transaction on Power Apparatus and systems, 1983, 102 (9) :3148-3154.
[4]余绍峰.考虑时频特性的地网导体的暂态分析[D].北京:清华大学, 2003:35-43.
[5]NEKHOL B, C.GUERIN C, LABIE P, et al.A finite element methed for calculating the electromagnetic fields generated by substation grounding systems[J].IEEE Transaction on Magnetics, 1995, 31 (3) :2150-2153.
[6]GOLDE R H.雷电[M].北京:电力工业出版社, 1982:24-25.
[7]李中新, 袁建生, 张丽萍.变电站接地网模拟计算[J].中国电机工程学报, 1999, 19 (5) :76-79.
[8]李中新.变电站接地网模拟计算[D].北京:清华大学, 1999:76-77.
岩土工程数值模拟方法的发展 篇11
关键词：岩土工程 数值模拟有限差分有限元边界 元离散 元无界元
1.引言
近几十年来，随着计算机应用的发展，数值计算方法在岩土工程问题分析中迅速得到了广泛应用，大大推动了岩土（体）力学的发展。在岩土（体）力学中所用的数值方法主要有以下几种：有限差分法、有限元法、边界元法、加权余量法、半解析元法、刚体元法、非连续变形分析法、离散元法、无界元法和流形元法等。下面就对这些方法进行简要的介绍和分析。
2.有限差分法
有限差分法是一种比较古老且应用较广的一种数值方法。它的基本思想是将待解决问题的基本方程和边界条件近似地用差分方程来表示，这样就把求解微分方程的问题转化为求解代数方程的问题。亦即它将实际的物理过程在时间和空间上离散，分解成有限数量的有限差分量，近似假设这些差分量足够小，以致在差分量的变化范围内物体的性能和物理过程都是均匀的，并且可以用来描述物理现象的定律，只是在差分量之间发生阶跃式变化。有限差分法的原理是将实际连续的物理过程离散化，近似地置换成一连串的阶跃过程，用函数在一些特定点的有限差商代替微商，建立与原微分方程相应的差分方程，从而将微分方程转化为一组代数方程，通常采用“显式”时间步进方法来求解代数方程组。
3.有限单元法
有限元法将连续的求解域离散为有限数量单元的组合体，解析地模拟或逼近求解区域。由于单元能按各种不同的联结方式组合在一起，且单元本身又可有不同的几何形状，所以可以适应各种复杂几何形状的求解域。它的原理是利用每个单元内假设的近似函数来表示求解区域上待求的未知场函数，单元内的近似函数由未知场函数在各个单元节点上的数值以及插值函数表达。这就使未知场函数的节点值成为新未知量，把一个连续的无限自由度问题变成离散的有限自由度问题。只要解出节点未知量，便可以确定单元组合体上的场函数，随着单元数目的增加，近似解收敛于精确解。按所选未知量的类型，有限元法可分为位移型、平衡型和混合型有限元法。位移型有限元法在计算机上更易实现，且易推广到非线性和动力效应等方面，故比其他类型的有限元法应用广泛。
4.边界元法
边界元法出现在20 世纪60 年代，是一种求解边值问题的数值方法。它是以Betti 互等定理为基础，有直接法与间接法两种。直接边界元法是以互等定理为基础建立起来的，而间接边界元法是以叠加原理为基础建立起来的。边界元法原理是把边值问题归结为求解边界积分方程的问题，在边界上划分单元，求边界积分方程的数值解，进而求出区域内任意点的场变量，故又称为边界积分方程法。边界元法只需对边界进行离散和积分，与有限元法相比，具有降低维数、输入数据较简单、计算工作量少、精度高等优点。比较适合于在无限域或半无限域问题的求解，尤其是等效均质围岩地下工程问题。边界元法的基本解本身就有奇异性，可以比较方便地处理所谓奇异性问题，故目前边界元法得到研究人员的青睐。
5.加权余量法
加权余量法也是一种求解微分方程的数值法，它在流体力学、热传导以及化学工程等方面应用较广。它具有两个方面的优点：①由于加权余量法是直接从控制方程出发去求解问题，理论简单，不需要复杂的数学处理，且它的应用与问题的能量泛函是否存在无关，因而它的应用范围较广，利用加权余量法这一优点去建立有限单元的刚度矩阵，可以大大扩展有限元法的应用范围；②加权余量法的计算程序简单，要求解的代数方程组阶数较低，对计算机内存容量要求不高，计算所需要的原始数据较少，这样就大大减轻了准备工作量。
6.半解析元法
半解析元法是Y. K. Cheung 于1968 年提出来的，同有限元法一样，它也是基于变分原理的。不同点是半解析元法根据结构的类型和特点，利用部分已有的解析结果，选择一定的位移函数，使解沿某些方向直接引入已知解析函数系列，而不再离散为数值计算点，因此自由度和计算量大大降低。这几年半解析法发展很快，种类很多，主要包括有限条法、有限层法、有限厚条法、有限壳条法、样条有限元法以及无限元法等。这类方法适用于求解高维、无限域及动力场等较复杂的问题。
7.无界元法
无界元法是P. Bettess 于1977 年提出来的，用于解决用有限元法求解无限域问题时，人们常会遇到的“计算范围和边界条件不易确定”的问题，是有限元法的推广。其基本思想是适当地选取形函数和位移函数，使得当局部坐标趋近于1 时，整体坐标趋于无穷大而位移为零，从而满足计算范围无限大和无限远处位移为零的条件。它与有限元法等数值方法耦合对于解决岩土（体）力学问题也是一种有效方法。上述介绍的几种数值法都是针对连续介质的，只能获得某一荷载或边界条件下的稳定解。
8.离散单元法
离散单元法随着非连续岩石力学的发展而不断进步，与现有的连续介质力学方法相比，还有以下问题需要研究：
（1）刚体离散单元法是基于非连续岩石力学的，更适合于低应力状态下具有明显发育构造面的坚硬岩体的变形失稳分析。对于软弱破碎、节理裂隙非常发育和高应力状态下的岩体变形失稳分析，则不适合。
（2）岩体介质种类繁多，性质非常复杂。在通常情况下，节理岩体或颗粒体表现为非均质和各向异性，并且常表现有很强的非线性，所处的地质环境不尽相同，这就使得岩土工程计算有很多不确定性因素。离散元的主要计算参数（如阻尼参数、刚度系数），影响到岩土工程稳定过程的正确模拟以及最终结果的可靠性，尤其是离散元计算中的参数选取，没有统一和完善的确定方法。
（3）计算时步的确定。现在的选取原则是出于满足数学方程趋于收敛的条件，与实际工程问题中的“时间”概念如何联系起来，合理地考虑时间效应，是今后需要研究的问题。
（4）迭代运算的时间较长。用计算机进行离散元计算时，CPU 占用时间较多，特别是在考虑岩块变形的情况下，模型划分单元数受到限制，对迭代方法需做进一步的改进。
9.刚体节理元法
刚体节理元法是Asai 在1981 年提出的，它是在Cundall 刚体离散元间夹有Goodman 节理单元的组合单元，但此节理单元有一定厚度而使离散元间不能“叠合”。刚体节理元法也可考虑不含节理单元的情况，即所谓的单一三角形刚体元非连续变形分析法，是石根华博士和古德曼教授于1984 年首次提出的一种新型数值分析方法，至1988 年该方法已形成了一种较为完整的数值计算方法体系。非连续变形分析方法以严格遵循经典力学规则为基础，是一种平行于有限元法的数值计算方法。
10.流形元方法
数值计算方法 篇12
结构拓扑优化是指在给定的设计空间、支撑条件、载荷条件和某些工艺设计等要求下, 确定结构构件的相互连接方式, 结构内有无孔洞、孔洞的数量、位置等拓扑形式, 使结构能将外载荷传递到支座, 同时使结构的某种性能指标达到最优。
结构拓扑优化研究最早从离散结构中最具代表性的桁架结构开始, 可以追溯到1904年, 由Michell提出的Michell桁架理论[1], 但其只适用于单工况, 并且依赖于选择适当的应变场, 不能适应工程实际问题。Bendse和Kikuchi[2]首次将均匀化方法成功用于连续体结构的拓扑优化设计中, 这标志着拓扑优化设计技术在理论和应用上的发展进入一个新的阶段。
拓扑优化研究的内容:材料插值方法、优化算法研究、数值计算不稳定性的消除方法、拓扑优化在工程中的应用等。其中数值计算不稳定性的消除非常重要, 直接关系到数值计算的准确性、收敛性和计算结果的可制造性等问题。
1 数值计算不稳定现象
拓扑优化中容易出现的数值计算不稳定现象包括:多孔材料 (Porous) 、棋盘格 (Checkerboard) 、网格依赖性 (Mesh dependency) 和局部极值等问题 (Local minima) 。多孔材料和棋盘格导致优化结构的提取和制造变得困难, 大大降低了结果的可制造性;网格依赖性使计算结果中杆状单元数量增加、尺寸减小, 不便于制造, 容易产生可靠性下降, 造成失稳;局部极值问题导致计算得不到全局最优解或得不到工程可行解。
所谓多孔材料, 就是在计算结果中出现许多中间密度单元, 使计算结果的可制造性差。多孔材料的计算结果如图1所示, 图中灰色区域就代表多孔材料。
所谓棋盘格, 就是在计算区域中材料密度为1和0的单元呈现出周期性的分布状态, 优化结构中会出现一种类似棋盘的网格现象。具有棋盘格的优化结果在工程上可制造性很差, 没有实际意义。棋盘格计算结果如图2所示。
注:图片来源于左孔天.连续体结构拓扑优化理论与应用研究:[博士学位论文].武汉:华中科技大学, 2004。
所谓网格依赖, 就是指用不同数目的单元离散化初始设计域, 最终会产生不同的优化结构, 随着网格剖分密度的增加, 优化结构的几何复杂性增加, 优化结构中杆状单元的数量逐渐增加, 而杆状单元的尺寸逐渐减小、数量变多且变得细长, 不便于制造, 在受压时容易失稳。网格依赖现象在不同计算网格下的计算结果如图3a、图3b、图3c所示。
所谓局部极值, 是指在相同离散格式下, 选择不同的算法参数以及不同的计算初始点时, 得到不同的优化结果。初始参数的一些小变化 (如几何设计域、单元数、周长约束值、过滤参数等) , 就可能导致优化设计结果的较大变化。如果是求解凸问题, 优化算法通常能够保证计算结果收敛, 但如果是非凸问题, 则往往只能收敛到附近的局部极值点而非全局极值点。
2 消除数值计算不稳定的方法
拓扑优化过程中出现的一些数值收敛问题实际上是数值计算中场函数的一种数值不稳定现象或数值奇异解造成的。国内外很多学者通过研究, 提出了一些解决问题的方法, 如周长约束法、局部梯度约束法、网格过滤法等。这些方法对于解决棋盘格和网格依赖现象有一定作用, 然而, 对于解决拓扑优化中的局部极值和数值奇异解问题, 仍然缺乏一种有效的方法和手段。
多孔材料的去除, 一般通过在材料模型中引入惩罚因子来达到, 这种方法比较容易实现, 主要是惩罚因子的选取要合理。例如密度法中常见的密度-刚度插值模型SIMP和RAMP, 就是通过引入惩罚因子的方法实现去除多孔材料。
局部极值问题导致在计算上得不到最优解, 或在工程上得不到可行的局部极小值解。目前还没有一种有效的克服局部极值问题的方法, 一般采用两种措施来减少局部极值问题的影响:一是从优化算法上考虑, 寻找一种合适的全局优化算法。二是从迭代初值上考虑, 采取选用不同的初始计算值来进行试算, 选取较好的优化结果。
对棋盘格现象, 常常采取增加几何约束的办法来限制数值解空间, 保证有限元求解的光滑性和稳定性;对网格依赖性, 常采取一些对密度变量的全局和局部约束, 以防止精细尺度的形成, 具体方法包括对优化问题增加约束以减少参数空间, 或采取过滤措施。
2.1 几种常见克服拓扑优化中数值计算不稳定性的方法
2.1.1 高阶单元法和非协调元法
Diaz和Sigmund[3]用高阶单元, 但对于有大量设计变量的优化问题会大量增加有限元求解时间, 造成求解的困难;吴长春等[4]认为低阶等参元的应力场难以满足局部坐标下的一次完备性, 而将非协调元和杂交元引入拓扑优化, 但该方法更多的是用在强度拓扑优化中。
2.1.2 周长控制法
Haber和Petersson提出周长控制法[5], 通过在优化问题中加入额外的周长约束来保证优化问题解的存在。离散结构的边长可以定义为:
这里lk是相邻单元i和j间的接触面长度, xi, xj是单元和的单元密度变量, ε是一个小的正数, 引入的目的是确保周长具有可微性。通过在优化模型中引入周长的上限约束P≤P*, 来对棋盘格等数值不稳定性现象进行控制。但周长控制法是一种全局约束方法, 并不能保证消除局部效果, 且周长约束值由设计人员的经验确定, 取值的合理与否直接影响到解的存在性。
2.1.3 局部梯度约束法
该方法由Petersson和Sigmund提出[6], 目的是对单元密度局部梯度进行约束, 以确保有限元离散后的精度。在二维网格分划中定义局部单元梯度约束为:
这里指标j和k表示单元在坐标系中的位置, s为密度斜率的上边界, h为单元网格尺寸。这种方法虽能较好地控制优化结构的几何复杂性, 消除棋盘格式和网格依赖性, 但在优化问题中引入2N或3N个额外的约束, 计算的时间大大增加。
2.1.4 网格过滤法
该方法由Sigmund提出[7], 用于修改目标函数的敏度信息:
其中H=rmin-dist (i, k) , {i∈N|dist (i, k) ≤rmin}, rmin是预先定义的最小单元直径dmin的一半, dist (i, k) 为邻近单元i和k的距离。Hi确定的一个圆形区域称为过滤区域, Hi的数值在过滤区域之外为零, 在过滤区域之内非零。式 (3) 表明:在单元xe影响区域rmin之内的单元对响应函数的敏度计算有影响, 在单元xe影响区域rmin之外的单元对响应函数的敏度计算没有影响。并且当rmin趋近于0时, 修改的敏度值收敛于原始敏度, 当rmin趋近于无限大时, 各敏度相等。通过对目标函数敏度的修改, 可以去除数值计算中的不稳定现象。
2.1.5 其他方法
Bendsoe[8]提出了全局梯度约束法, Kikuchi[9]提出了消除均匀化方法中数值计算不稳定性现象的引力函数法。Petersson[10]提出了一种显式约束方法以控制中间密度单元的生成。Bourdin[11]基于SIMP插值模型, 提出一种改进的柔度最小化的拓扑优化模型。在该模型中, 用卷积函数同点态密度 (Pointwise Density) 的方法来提高密度函数的光滑性和可微性。左孔天等[12]提出了一种基于卷积因子前处理和小波分析后处理的综合过滤算法。在Sigmund提出的一重敏度过滤方法的基础上, 罗震[13]提出了一种二重敏度过滤方法, 进一步对过滤范围内各单元的设计敏度重新进行加权均化, 进一步提高敏度过滤技术的初始网格无关性。郭旭等[14]以水平集函数作为拓扑设计变量求解连续体问题。该方法具有非局部效应特点, 可减缓棋盘格式等。
2.2 克服数值不稳定性各种方法小结
对于结构拓扑优化过程中的数值不稳定性现象, 目前还没有一种能够适用于所有情况的方法来消除该问题。以上描述的一些方法和技术, 都有各自的本身优缺点和适用范围, 特别是在实际工程应用中, 不同的数值计算方法对计算收敛性和计算速度影响很大, 根据实际的情况选择一种合适的克服数值不稳定性的方法显得尤为重要。
3 结语
结构拓扑优化是结构优化领域相对比较难的问题, 技术方法也相对不够成熟, 目前虽然已经取得了长足的进步, 但是还是有很多问题亟待解决。数值计算的稳定性是整个结构拓扑优化过程中最重要的部分, 直接关系到数值计算的准确性、收敛性和计算结果的可制造性, 对于实际工程应用方面, 显得尤为重要。
目前所出现的各种去除数值计算不稳定性的方法, 或多或少都存在各自的局限性, 对于结构优化过程中的去除数值不稳定性的研究, 还有很长的路要走。随着计算硬件和各种算法的不断发展, 也必将出现越来越多先进的方法和技术, 用于去除拓扑优化过程中的数值计算不稳定性。
摘要：介绍了结构拓扑优化的研究内容和方法, 分析了拓扑优化计算中常出现的多孔材料、棋盘格、网格依赖性、局部极值等数值不稳定现象的表现形式以及产生原因。提出了去除结构拓扑优化中数值不稳定性现象的各种方法, 及其各自的特点。
【数值计算方法】推荐阅读：
数值计算方法结课心得07-06
耦合数值计算06-09
电场数值计算07-28
仿真数值计算08-01
数值计算模型08-13
数值计算论文08-29
MATLAB数值计算06-02
CFD数值计算07-22
数控加工中的数值计算09-09