数值仿真十篇

数值仿真篇1

关键词：数值仿真；工科课程；应用研究

数值仿真是一类基于计算机求解，针对实际工程问题采用仿真模拟，从而进行数值实验的方法，已成为与实验研究、理论研究相并列的研究科学问题的必备方法。基于数值仿真，可通过图像显示所研究问题发生的物理、化学过程，如天气预报中，即可通过计算机仿真的方式预测某地区的温度变化过程，具有直观性强、成本低、可研究问题范围广的优点。如今，数值仿真方法已被应用于力学、机械工程、能源动力工程等工科类专业课程教学中，为复杂或不易开展的实物实验提供数值实验途径，加深学生对课程中复杂的物理、化学原理的理解。

一、数值仿真应用于工科类课程理论教学的具体方法与特点

数值仿真是利用数值计算原理与计算机求解技术开发出的现代科研方法。基于所求解问题的物理、化学原理，采用有限元等数值分析方法，将连续的物理、化学问题离散求解，获取所求解对象的物理、化学量分布特性。目前，数值仿真方法广泛应用于科学研究，如高温燃烧热场的温度分布、高超音速流动的流场压力分布、天气温度预测等。在工科类课程教学中，由于理论、概念多，物理、化学过程复杂，常制约着学生的理解，而数值仿真方法可通过图像显示的方式展示所研究问题发生的物理、化学过程，具有可操作性强、结果显示直观的特点，因此，可将数值仿真方法引入工科类课程教学，提升学生对课程知识的理解。一般的数值仿真过程分为预处理、求解、后处理三部分，学生在掌握基本数理知识的基础上，即可理解数值仿真的基本流程。在预处理中，需建立所分析问题的几何模型，并给定求解条件。例如，分析钢结构塔在承受一定重量时是否具有足够的强度，需首先基于分析软件建立钢结构塔的外形几何模型。在求解环节，基于数值计算方法，利用分析软件中求解仿真模型以获取需掌握的物理、化学量。在后处理环节，完成分析结果的可视化，以图像形式直观展示分析结果，如天气预报中展示的温度云图，即是后处理结果。在工科类课程中，涉及大量有关物理、化学现象的理解，而某些环节无法通过直观实验复现，也无法通过理论分析获得结果，对学生的理解带来困难。例如，在锅炉设计中通过改进折焰角的尺寸，即可优化锅炉炉膛内部的燃烧热场，而炉膛内部燃烧的温度无法通过理论计算获取，也无法直接进行实验测试，但可基于数值仿真获取其燃烧温度分布图，直观展示折焰角尺寸调整对热场分布的影响，有利于学生通过实践操作，获取分析结果，加深理解。

二、数值仿真应用于理论教学的具体实例———材料力学弯曲应力求解教学

材料力学是机械、土木、能源等工科类专业的必修专业基础课程，是结构设计的必要基础，基于材料力学基本原理，可分析部件的强度、刚度及稳定性，为工程结构设计提供指导。材料力学中的杆件弯曲问题是重要的教学知识点，当矩形截面杆件受弯时，杆件横截面上的正应力呈线性分布，在横截面上、下边缘正应力最大，在中部中性轴处最小。在理论求解中，需首先根据给定约束条件与外部受力特点，分析杆件受到的弯矩，然后利用理论公式求解正应力。弯曲应力求解为材料力学中的重点与难点内容，区别于杆件的轴向拉压与扭转问题，弯曲问题的应力分布更复杂，学生在学习时不易获取直观理解，故可借助数值仿真方法。利用受力求解的数值仿真软件，在预处理环节，首先建立杆件的几何模型，根据求解要求设置约束，其后在给定位置施加受力，最后，指定杆件所使用的材料并由数值仿真软件自动完成求解。该问题为静力学求解问题，求解精度较高。在完成求解后，进入后处理环节，在该环节可通过软件直接显示梁在弯曲作用下的变形，也可以通过图像显示梁的应力分布。藉此，直观清晰地向学生描述了相关力学概念，展示了复杂的受力情况，使学生对课程理论学习产生更直观的认识，既提升了学习兴趣，又进一步加深了对课程内容的理解。

三、数值仿真应用于理论教学的结论

数值仿真篇2

关键词：石油工程；岩石力学；仿真；应力；应变

引言

为使学生更好地认识和理解石油钻井、开发过程中普遍存在的岩石力学问题，在石油工程岩石力学实验教学中，岩石的拉、压、剪基本实验及岩石的破裂与失稳过程是一个复杂而重要的基本实验教学内容[1]。由于井下高温高压并存，井下岩石材料的非均匀性、非连续性，以及外荷载作用下岩石组织结构之间相互作用的复杂性，当前的实验教学方法难以直接演示岩石变形及破坏的复杂现象，已有的理论解析方法仍然缺少对此过程有效的研究手段，且目前的物理实验不易对岩石的破裂与失稳过程等现象做准确的描述。石油工程岩石力学物理实验教学尚存些许不足之处，而数值实验方法可有效地克服此类问题。

1石油工程岩石力学仿真实验的必要性

（1）石油工程岩石力学物理实验教学内容时效性较差，创新性不足。随着大数据时代的来临，油气钻井、压裂技术已数据化、信息化，并朝向智能化发展，高校石油工程岩石力学实验也应与时俱进，及时更新实验内容、方法及手段[2]。而石油工程岩石力学出于成本及资源利用率的考量，高校不可能将最新的岩石力学试验设备技术实时引进实验室，故如何使学生及时接触最前沿的石油工程岩石力学实验技术，有效拓宽思维空间，就成为一个亟待解决的问题。（2）石油工程岩石力学物理实验成本高，可重复性差。受到现场条件、人力、物力和财力的限制，通常的岩石力学实验教学很难通过大量的物理实验向学生直观演示各种岩石变形、破坏的复杂现象。（3）石油工程岩石力学物理实验过程演示直观性不足。由于岩石介质的复杂性，传统的物理实验测试方法很难全面地反映岩石在变形损伤演化和宏观破坏过程中的应力场、变形场等重要信息，而数值实验法可以取得较好的效果。（4）石油工程岩石力学物理实验缺乏必要的井下工况，可理解性较差。尽管岩石破坏试验可得到拉、压、剪等更基本数据，但与井下岩石结构破坏密切相关的工程力学现象，如井壁失稳、井下工具破岩、水力压裂等，通过实验室小型试样的加载试验是难以理解的，因为这些破坏不仅与构成岩石结构的材料性质有关，而且还涉及复杂的井下高温高压环境。

2石油工程岩石力学仿真实验实施方法

2.1融合仿真技术，拓宽石油工程岩石力学实验对象

了解当前国内外岩石力学发展的现状和进展，考虑当前石油钻井、开发过程中普遍存在的岩石力学问题，掌握解决这些问题的岩石力学基本理论和实验方法，基于大数据时代的计算机仿真技术，将实验新技术与仿真技术相融合，升级实验项目，培养学生的创新能力。

2.2考虑井下工况，构造岩石仿真模型

根据所要求的井下工况，虚拟制造岩石材料，加工其几何形状。将所要建立的岩石模型简化为复杂多面体。首先，初始化结构体系统；其次，切割岩体模型，实现单元的离散；然后，将结构面进行统一编号，对现存的单元进行判断；最后，去除虚拟结构面，形成岩石仿真模型。此过程可重复进行。

2.3岩石模型应力应变分析，渐进破坏演示

（1）在上述岩石几何模型基础上，考虑井下岩石的矿物成分、结构构造、结晶情况、试样尺寸、围压、加载速率、应力路径、孔隙水压力、温度及湿度（含水率）等因素影响，对不同围压下的岩石变形与强度特性进行研究，模拟弹性、塑性、硬化、软化和摩擦等阶段应力状态下岩石的应力—应变关系。（2）模拟岩石宏观力学破坏行为，实现相应的力学及多场（应力场、应变场）耦合分析演示（如图1所示），再现细观结构上岩石的损伤演化、裂纹的扩展及应力的重分配等特征，使岩石破坏过程及蠕变现象能以3D动画的形式展示。

2.4仿真结果分析，实现岩石力学实验虚拟仿真

（1）对照岩石模型结构（微观、宏观）和组织特点，结合仿真结果，分析岩石的物理性质，弄清影响岩石力学性质的影响因素，理解岩石的流变性质（已知现象）。（2）探索井眼围岩的应力状态与岩石的破坏准则（未知现象）（如图2所示），实现岩石力学实验抗拉、抗剪、抗压虚拟仿真。促使学生能够独立思考，对岩石力学实验相关机理充分吸收，达到较为理想的效果。

数值仿真篇3

关键词：三维刚体弹簧元 拱坝地基破坏 拱坝 仿真

近年来，用非连续介质力学模型研究高坝地基、高边坡与地下工程的变形稳定有长足进展，离散元[1～4]、DDA[5]、刚体弹簧元[6]等的发展提供了求解非连续介质力学问题的有效工具.关于高坝-地基系统的破坏机理与仿真分析，由于问题的极端复杂性，迄今仍未见诸报导，本文利用改进的刚体弹簧元研究玛尔帕塞拱坝溃坏的机理与过程，是试图用连续-非连续介质统一的数值模型研究高坝-地基系统破坏过程的尝试.

法国玛尔帕塞拱坝1959年12月2日的溃坝事故是坝工史上的重大事件[7，8]，引起了坝工界的极大重视.世界各国专家和学者从未停止过对玛尔帕塞拱坝失事原因的探究，并提出了各自见解[9]，归纳起来有以下3种看法：(1)坝基变形，认为左岸坝基修建在弱风化顶层上，岩体弹性模量低，仅1000MPa左右，地基的过大变形使拱的推力向上向下转移，其中向上转移的推力使左岸重力墩超载而产生位移，导致拱坝失去支撑而破坏.(2)坝体上滑失稳，认为岩体浅层的裂隙在不利的条件下不能提供足够的抗滑力，致使坝体沿建基面向上和向下游滑动，拱圈拱弦被拉长，使坝沿一岸或两岸拱端转动，导致拱冠断裂，坝体崩毁.(3)坝肩岩体滑动，认为左岸坝基存在一个潜在的由下游断层、发育的节理面，以及下游临空面组成的滑动楔形体.由于地质条件在楔形岩体上游侧形成很高的静水扬压力使岩体滑移失稳.本文利用改进的三维刚体弹簧元对玛尔帕塞拱坝的失稳破坏过程进行仿真，对该坝的失稳机理作进一步的探讨.

1 拱坝破坏仿真模型

1.1 三维刚体弹簧元模型 假定三维刚体弹簧元的单元块体为不可变形的刚性体，在引入了刚度、强度等效准则后，将块体变形与强度指标统一集中于接触弹簧上.块体的运动符合刚体的平动和转动定律，块体与块体依靠相互之间连接的切向和法向弹簧和阻尼器来传递相互之间的作用，如图1.弹簧的作用力由块体间的相对位移确定.如图2，块体P对块体Q的相对位移为Δu，则他们之间的弹簧力变化量为

式中：Ks、Kn分别为切向和法向弹簧的刚度；Δs、Δn分别为相对位移增量的切向和法向分量；Δs、Δn分别为切向和法向弹簧力的变化量.

图1 接触模型

图2 相对位移增量示意

1.2 使用三维刚体弹簧元模拟非连续介质 三维刚体弹簧元模型基于离散元方法提出，是非连续介质模型的1种，由于块体之间可以错动、分离，因此可以方便地模拟非连续体如裂隙发育岩体的大变形、大位移问题.离散元方法中接触关系一般可分为角面、角边、角角和边边四类，关于各类接触关系的检索和其它离散元的详细内容已在Cundall和鲁军等人的文章中有详细叙述，参见文献[1～4].

1.3 使用三维刚体弹簧元模拟连续介质

图3 刚体弹簧元模拟连续介质原理示意

图4 中心受拉杆件

1.3.1基本原理 如图3，通过将连续介质的变形和应力凝聚到单元之间连接的法向和切向弹簧上，可以用刚体弹簧元模型来模拟连续介质.即由弹簧的变形刚度(Kn、Ks)来等效连续介质的变形模量(E、G)，由弹簧力系统(Fn、Fs)来推求与之静力等效的连续介质的应力(σ、τ)分布，从而将刚体弹簧元系统和连续介质联系起来.(1)刚度等效.如图4所示，一连续的l×l1×l2方杆，弹性模量为E，泊松比为υ.将杆沿长度方向等分成n个l3×l1×l2的块体，其中l3=l／n，每个块体间的面面接触被视为四个角点处的点面接触，当n足够大，取弹簧的法向刚度Kn=l1×l2×E/4l3时，使用三维刚体弹簧元计算它在中心作用力P下的变形与使用材料力学方法的结果是一致的.弹簧的切向刚度借用弹性模量和剪切模量的关系可取Ks=Kn/[2(1+υ)].单元应尽可能采用长方体，对于不规则的块体可以将其均化成长方体使用上述估算式，或另外推求相应的计算式.

(2)强度等效.如果刚体弹

簧元网格划分得比较密，则可以取单元面上每个弹簧力合力除以该面的面积求得其平均应力作为该面上的应力值.因为单元网格较密，可以认为当面上的平均应力达到材料的强度时该面上的弹簧全部断开，则单元在该面发生破坏.由于刚体弹簧元的计算工作量比较大，当为了节省计算工作量，网格划分得较稀疏时，可以按照材料力学方法假定单元面上的应力分布方式，由弹簧力根据静力等效原则求出该面上的应力分布.

(3)荷载等效.对于集中荷载，将荷载作用在相应的块体形心上，并将荷载对中心的力矩加入到该块体的合力矩中即与原荷载的作用等效.对于分布荷载也要等效成集中荷载加入到相应的块体中.

图5简单算例

1.3.2简单算例 如图5，一长40m，高4m，宽1m的悬臂梁受均匀荷载P=105N/m2作用，支座A固定，不考虑自重，坐标见图中所示.悬臂梁弹性模量E=1010N/m2，泊松比υ=03，密度ρ=2.4×103kg/m3，离散网格如图5(在厚度方向为一层，厚1m).使用上述的刚体弹簧元方法计算悬臂梁的变形和应力并与有限元的结果进行对比，图6为悬臂梁沿y轴的位移对比结果，图7和图8分别为x方向的正应力σx和剪应力τxy的对比结果.从图中可以看出两种方法计算的位移和应力结果除有限元法在固定端能反应局部应力集中外，整体应力具有很好的近似性.

图6 y方向位移等值线(单位：m)

图7 x方向正应力σx等值线(单位：MPa)

图8 剪应力τxy等值线(单位：MPa)

刚体弹簧元通过引入单元表面应力分布假设和变形近似等效来模拟连续介质，利用这些假设，刚体弹簧元通过上述的刚度等效、强度等效和荷载等效，能够把连续介质和非连续介质统一在一个模型之中，为研究坝肩失稳导致拱坝断裂破坏需要仿真连续与非连续耦合介质的大变形与破坏断裂的复杂过程提供了有力的工具.

2 玛尔帕塞拱坝失稳机理的研究

2.1 模型的建立

2.1.1 坝肩和坝体的离散和计算条件的设定 法国玛尔帕塞拱坝坝高66m，坝顶高程102.55m，坝顶长222.7m，长高比3.3，坝顶厚1.5m，底厚6.78m，厚高比为0.1.坝址岩体由带状片麻岩组成，含千节理倾向下游右岸，倾角30°～50°.该坝于1954年建成，蓄水历时5年至1959年12月蓄至最高水位约100m时，坝在短瞬间溃决，左岸与中部坝体全部冲走，右岸底部坝体残存.根据玛尔帕塞拱坝的地形和地质资料以及拱坝的体形和横缝设置将坝址处的岸坡简化成如图9所示，坐标系设置和拱坝离散网格如图10.在分析中边界条件设定为除“二面沟”内的岩体外，其他岩体均被固定.坝体、右岸坝肩和中间岩体依照连续介质进行模拟，坝体弹性模量取E=1010Pa，坝体材料的抗拉强度取2.0MPa，抗压强度为20MPa，摩擦系数f=1.0；考虑到右岸坝基事后调查仍有坝体残存，表明该区岩体强度较高，故将坝体右岸及中间部分与基础的交界面强度较坝体提高1倍.左岸坝肩依照非连续体计算，重力墩和地基之间的摩擦系数设定为1.0，坝体和左岸坝肩之间以及“二面沟”内的岩体两组节理面与水平层面取相同的摩擦系数f1，根据不同的工况选取，上、下游断层摩擦系数为f2，亦根据不同的工况进行选择.计算中未考虑粘结强度c的作用.主要荷载为库水作用，稳定分析中库水水位均取100m，考虑左坝肩扬压力时，按静水压力作用施加在沿上游面坝基线的裂隙处，作用面见图9中虚线段所示.为了简化计算，所有荷载均一次施加.计算中阻尼参照Cundall提供的阻尼模式选取刚度与质量成比例的Rayleigh阻尼型式，阻尼系数α、β由试算确定，试算表明：如不出现过阻尼或阻尼过小情况计算结果均收敛于一稳态解的近似值，详见参考文献[3].

2.1.2位移分析 计算拱坝位移时左坝肩暂按连续介质考虑，并且不考虑扬压力作用，得到拱坝在1955年9月(水位79.75m)蓄水至1959年7月(水位94.19m)拱坝产生的位移，并与相应实测位移值进行比较，图11列出了顶拱和拱冠梁处的比较结果.可以看出二者最大值基本接近，但在拱冠梁坝基附近，实测的位移值较计算大1cm左右，顶拱计算结果较实测结果略偏大.这一比较结果表明选择的弹簧刚度参数Kn、Ks以及单元等效模型基本符合实际.

2.2 玛尔帕塞拱坝稳定性分析及失稳过程的仿真 通过改变坝体和左岸坝肩之间以及左岸坝肩“二面沟”内岩体的摩擦系数f1和上下游断层面的摩擦系数f2可以模拟和分析拱坝在各类参数条件下失稳破坏的情形.

图9岸坡离散

图10刚体弹簧元离散网格

图11拱坝位移计算结果与实测结果比较

2.2.1工况1——不考虑扬压力的浅层滑动 不考虑左岸坝肩处的扬压力，将坝体和左岸坝肩之间以及左岸坝肩“二面沟”岩体的摩擦系数f1和上、下游断层面的摩擦系数f2设为相同，从1.0开始折减，逐步计算，直到拱坝发生破坏.计算结果表明当摩擦系数f1和f2折减到0.6时坝肩失稳，拱坝破坏.图12列出了左坝肩附近的块体1(见图10)的位移过程.f=0.7时，块体位移量陡增，接近30cm，此时拱坝处于失稳的临界状态，f=0.6时，块体位移不再收敛，开始沿着坝肩向下游、向上滑动，整个拱坝随之破坏.

图12 块体1的位移过程(工况1)

图13 f=0.6时左坝肩重力墩位移过程(工况1)

分析图13～图15，拱坝的整个破坏过程可以分为两个阶段，第一阶段是大约3s以前拱坝整体位移，第二阶段是3s以后由于过大的变位使拱坝的工作方式发生改变，应力恶化，坝体溃毁.描述如下：在第一阶段，拱坝利用拱的作用将水压力转化为水平作用在坝肩上的推力，沿坝肩向上的分量使拱圈产生向上滑动的趋势，由于左坝肩岩体的摩擦系数较小，不能提供足够的抗力，巨大的拱推力转移到了左岸的重力墩上，使重力墩沿着墩轴线偏向下游移动了近2m(见图13，3s前重力墩主要以平动为主)，同时整个坝体沿建基面向下游、向上滑动，带动坝体以右岸为轴向下游转动(由图14中y方向的位移分布，坝体向下游位移由右到左依次增大，显示以右坝肩为轴向下游转动)和向上旋转(见图14中z方向的位移分布，坝体上抬量由右到左依次增大，显示以右坝肩为轴向上转动).此时，整个拱坝基本上保持整体发生变位，同时有沿右岸转动趋势.

(a)y方向位移等值线

(b)z方向位移等值线

图14 f=0.6，3s时拱坝坝体位移分布(工况1，单位：m)

拱坝在经过了第一阶段后，拱圈拱弦拉长，逐渐被压扁，承载能力极度降低，应力发生恶化，坝体破裂，拱坝的破坏进入到第二阶段——破坏阶段.由于中、下部块体受水压力较大，拱坝坝体从该部位偏左岸发生断裂，向下游弯折，首先被冲出，随后其周围坝块跟着折向下游，上部坝块受水压力较小，发生翻转后随在中下部坝块之后冲出，拱坝彻底破坏.与重力墩连接的坝段在被向下游冲出的过程中又带动重力墩发生了向下游的位移和旋转(见图13中3s以后的位移，向下游移动1m左右，在平面上旋转约3°).图15给出了上述过程

一些典型时刻拱坝的状态，拱坝破坏后最后残留状态见最后一幅.

图15 f=0.6时拱坝破坏过程(工况1)

上述拱坝破坏过程的计算结果和拱坝失事后的调查结果基本相符：(1)f值的大小，“实有摩擦角为30°～35°”，本文计算表明f值取0.7时拱坝处于临界状态，对应的摩擦角为35°，当拱坝左坝肩的实有平均摩擦角小于该值时，拱坝会发生上滑失稳，与拱坝破坏事实相符；(2)坝体的残留部分，拱坝失事后“右岸部分的坝体和中间部分的坝体底部基础留存了下来”，坝基由于嵌固作用，一般不易破坏，当拱坝上滑失稳时，应力发生恶化，拱坝沿横缝和工作缝等薄弱处断裂，同时底部的坝体得到卸载，从而保留了下来，计算中将坝体右岸及中间部分与基础的交界面强度较坝体提高了1倍，坝体的残留部分如图15，和拱坝破坏后的结果相符；(3)重力墩的位移，“左按重力墩在沿墩轴线方向上已移动约2m，沿垂直于墩轴线方向，向下游移动了约80cm，并略有倾侧”，本文计算的过程解释了上述重力墩先沿轴线方向移动(拱坝整体滑动时的推力所致)，又向下游移动并转动(拱坝破坏后，重力墩附近的坝段向下游破坏时的带动)的现象，位移量也基本接近；

(4)残留部分的位移规律，“AB坝段的位移很小，各坝段的位移值自右岸向左岸逐渐增大，至JK坝段处增达80cm；位移方向与坝轴线略呈斜交，并沿上游至下游方向略偏向于左岸.”，“左右岸，建筑物上游的岩石从混凝土与基岩接触点的几十cm处发生破裂，并有一条张开裂隙从底部沿岸坡脚向上延伸至接缝B处；这条裂隙在底部宽40～50cm，向上续渐减小.在右岸坝的下游可以看到若干在建筑物发生位移时所压碎的岩石区”，由于本文计算原来主要希望研究左岸坝肩滑动安全，未对中间坝基和右岸坝肩做离散和参数研究，因此未能计算出上游的裂缝和残留坝段80cm的位移量，但3s时的位移分布可以看出整个拱坝沿右岸的旋转和向下游的倾侧现象，这和上述调查结果相吻合，是产生上述裂隙和位移的原因，如果坝基岩体强度不够，整个拱坝沿右岸的旋转以及沿坝基向下游的倾侧必然产生很大的对上游面岩体的拉力和对下游面岩体的压力，使其开裂或破坏.

2.2.2 工况2——考虑扬压力时的浅层滑动在左岸坝肩处施加扬压力，重复工况1的计算.其结果与工况1基本一致，f=0.7时，拱坝处于失稳的临界状态，f=0.6时，拱坝破坏.表明扬压力对于浅层滑动的情况影响较小.

2.2.3 工况3——不考虑扬压力的深层滑动取坝体和左岸坝肩之间以及左岸坝肩“二面沟”岩体的摩擦系数f1为1.0并保持不变，在不考虑左岸坝肩处的扬压力情况下，将上、下游断层面的摩擦系数f2从1.0开始折减，逐步计算，直到拱坝发生破坏.计算结果表明：直到f2折减到0.1时坝肩仍然保持稳定，拱坝没有溃毁.因此，如果没有扬压力的作用，拱坝发生深层滑动的可能性极小，依靠岩体自身的重量基本上就可以保持稳定.

2.2.4 工况4——考虑扬压力时的深层滑动在左岸坝肩处施加扬压力，重复工况3的计算.计算结果表明当上、下游断层的摩擦系数f2折减到0.3时坝肩发生滑动失稳，拱坝破坏，见图16.分析图17、图18，与拱坝发生浅层滑动相对应，整个破坏过程也可以分为拱坝整体位移和坝体溃毁两个阶段.其过程如下：在第一阶段，作用在坝肩上的拱端推力向坝肩岩体深层传递，在最不利的滑动面(下游断层面)上，坝肩推力沿该面向上的分量和作用在坝基岩体上的扬压力有使坝肩岩体沿滑动面产生滑动的趋势，当岩体自重和滑动面上的抗滑力不足以阻止岩体的滑动时，整个坝肩岩体将沿着断层面产生滑动，巨大岩体的滑动带动整个左坝肩和重力墩运动(见图18中6.25S时的状态，坝肩岩体向上和下游滑动，重力墩折断破坏)，使重力墩发生平动和旋转(见图17，625S前重力墩的旋转)，拱坝坝体左边坝肩部位被抬起并推向下游，带动拱坝整体沿着右岸发生向下游转动.以后的拱坝运动、破坏的过程和工况1十分类似，不再赘述.

图16 块体1的位移过程(工况4)

图17 f2=0.3时坝肩重力墩位移过程(工况4)

深层滑动的计算结果和拱坝失事后的调查对照如下：(1)f2值的大小，“下游断层面中充填的粘土(摩擦角)的试验值只有30°”，本文计算的结果表明f2值取0.4时拱坝处于临界状态，对应的摩擦角为22°，当拱坝左坝肩上下游断层的平均摩擦角小于该值时，拱坝会发生深层滑动；(2)坝体的残留部分，深层滑动计算结果亦与其相似.(3)重力墩的位移，深层滑动计算的过程是重力墩一直在发生转动，稳定后在平面内的位移量和转动与实际调查结果接近；(4)残留部分的位移规律，这一现象和浅层滑动相同，不再重复.

图18 f2=0.3拱坝破坏过程(工况4)

3 结论

(1)计算表明玛尔帕塞拱坝在坝肩岩体抗滑强度较低时两种滑动失稳模式均有可能发生.当左坝肩岩体节理、裂隙的摩擦系数小于0.7时拱坝将发生浅层滑动失稳，而左坝肩上、下游断层的摩擦系数小于0.4，且坝肩作用有扬压力时拱坝会发生深层滑动失稳，因此，似乎拱坝发生浅层滑动的可能性更大一些.

(2)对拱坝沿坝肩岩体的浅层滑动和深层滑动失稳破坏过程的仿真计算均较好的解释了拱坝失事原因，再现了拱坝溃坏的过程和坝址处的残留情况.浅层滑动失稳过程似乎与调查结果符合得更好，然而由于计算资料并不十分充分，考虑的因素也不可能很全面，另外在计算中有所简化，深层滑动过程与浅层滑动过程亦有很多相似，因此很难排除或肯定拱坝失稳确切是那一种滑动模式.

(3)扬压力对浅层滑动影响很小，它主要是由于库水荷载的作用产生沿坝肩的滑动，但扬压力对深层滑动起着重要的影响，如果没有扬压力，对于玛尔帕塞拱坝发生深层滑动失稳的可能性很小.

数值仿真篇4

关键词：微型涡喷发动机 离心式压气机 数值仿真

中图分类号：TH452 文献标识码：A 文章编号：1674-098X（2017）01（b）-0022-04

Abstract：In this paper， by using BladeGen， Turbogrid and CFX commercial software for 3D modeling divided hexahedron grid and numerical simulation to the configuration 66mm micro centrifugal compressor. Finally get the compressor flow characteristic curve of 20 000 rpm to 140 000 rpm. It found that the efficiency of 94.7%. type compressor has a stable working area is the width is between 60 000 rpm and 100 000 rpm， the results coincide with the original engine performance parameters. At the same time， this paper selects the optimum condition of 100 000 rpm the speed of the compressor at near stall condition， the blocking condition three cases， on the flow field inside the compressor of the three cases were analyzed， to provide a reference for improvement the future of the compressor.

Key Words：Micro Turbojet Engine； Centrifugal Compressor； Numerical Simulation

浩机是涡轮喷气式发动机的重要部件。作为微型涡喷发动机的核心，压气机的性能决定了发动机的整体性能。该课题以一款微型涡轮喷气式发动机所使用的直径为66 mm的微型离心压气机为研究对象，分析其各工况流动特性与性能。

西安交通大学的赵会晶等人运用CFD方法模拟了离心压气机叶轮前缘倾掠对压气机性能的影响[1]，韩国仁荷大学与三星泰科能源设备研发中心通过进化算法对离心式压气机叶轮进行多目标优化[2]，由此可见对于离心式压气机的研究一直是热门话题。

该文主要采用数值模拟的方法对某一构型微型离心压气机进行流场分析以及性能测试。数值模拟的方法可以在短时间内获得较为精准的流场以及实验结果，可以模拟各种复杂的流动。通过改变初始条件与边界条件获得不同工况下该构型离心式压气机的流量特性曲线，同时获得不同转速、不同流量下叶轮流道内气体流动情况、叶片表面的压力与速度矢量分布和不同叶高的压力分布与速度矢量分布。

1 压气机三维构型建模

根据扩压器尺寸以及机匣尺寸确定子午面几何形状，确定子午面尺寸参数后使用CFX-BladeGen组件进行叶轮建模。经过与创新创业项目的压气机实物进行比对与实验，确定了叶片构型。

叶片数量：5组叶片，每组包含一片主叶片与一片小叶片；叶片角度：主叶片进气角49°，出气角12°；小叶片进气角69°，出气角5°。通过测量各截面叶厚及叶片安装角可绘制出叶轮模型。

确定各项参数后即可生成叶轮的三维模型并确定最终构型，最终叶轮构型如图1所示。

2 压气机流场数值计算

此次课题的CAD模型已经使用ANSYS BladeGen建立完成，选择使用ANSYS TurbiGrid进行网格划分，使用ANSYS CFX进行求解并使用ANSYS CFD-Post进行后处理分析。

2.1 导入模型

将之前在BladeGen组件中生成的叶片模型导入CFX中的TurboGrid组件，准备进行网格划分。

2.2 网格划分

将模型导入TurboGrid后，软件会根据BladeGen的模型信息自动划分出计算域，包括进气口、出气口、主页片、小叶片、机匣面、轮毂面以及两个切面共8个部分，如图2。

导入模型后输入叶尖间隙为固定距离0.2 mm。

为获得较好的网格质量，在Topology set中Placement选择Traditional with Control Point，通过控制点控制网格形状。在Mesh Data菜单中，Node Count选项选择Fine（250 000），使生成的网格包含大约250 000个节点，提高网格密度。

生成网格后通过添加与移动控制点来获得一个较好的拓扑学模型，如图3，从而得到质量更好的网格。最终生成的网格如图4，共包含251 137个节点与232 320个网格单元，全部为六面体网格。

2.3 初始条件设定

（1）初始计算域。

基本设置：Material设置为“Air Idea Gas”；Reference Pressure设置为“1 atm”；Domain Motion中Option设置为“Rotating”；Angular Velocity设置为计算所需求转速。

（2）湍流模型的选择。

对于离心式压气机的计算，比较适合的有模型与SST模型。比较两种模型的优势，该课题选择使用模型作为计算的湍流模型。

（3）热传导模型：在此次计算中将压缩过程理想化为绝热过程，因此热传导模型选择Total Energy。

（4）确定边界条件。

在该次计算中，边界条件如下所示。

Inlet：固定总温298 K，总压1atm，进气方向与转轴平行；Outlet：平均静压，该文通过改变出口平均静压值以模拟不同工况；Wall：机匣，轮毂；叶片均设置为无滑移、不导热壁面；Interface：设置叶间隙平面为接合平面；Rotational Periodicity：为减少计算量，取叶片的1/5计算，其两侧壁面设置为旋转周期边界。

2.4 求解设置

完成初始设置后，在模型的出口面上添加一个监测点，以便在计算时观察收敛情况。在求解设置时需要确定计算收敛残差与迭代步数。此次计算中将收敛残差设置为10-5，迭代步数为最大500步，无论哪个条件先到达即认为结果收敛。前处理完成后生成相应的def文件，在求解器中进行计算。

2.5 求解器计算

在求解过程中，监视器会显示计算残差的折线图，同时可以监控在上文中加入的检测点的总压数值，以辅助判断计算是否收敛，若总压趋于平稳可判断计算收敛。

3 数值计算结果与分析

在该章中对叶片的计算结果进行分析，绘制出该构型叶片的流量特性图，并对叶片在3个不同工作状态下的流场进行分析。

3.1 压气机流量特性

表征压气机性能最主要的描述是压气机的流量特性图[3]。压气机的增压比随着压气机流量的减小而增大，但当流量减小至某一值时，压气机效率急剧下降。在流量特性曲线中，将所有转速下的喘振边界连起来即可得到整_压气机的喘振边界。喘振边界将流量特性曲线分为了两个区域，左侧为不稳定工作区，右侧为稳定工作区，当压气机工作状态落入不稳定工作区时，压气机进入喘振状态，此时增压比下降，压气机效率下降，在压气机工作中要极力避免[4]。

该文中使用CFD方法计算出压气机的流量特性，可以节省大量时间与经费，在压气机初期设计时有较好的应用，可以快速地对压气机效率进行评估。在该文中计算了7个转速条件下（20 000～140 000 rpm）通过改变出口静压改变压气机工况，通过观察压气机是否出现分离涡以及是否出现堵塞以分辨压气机是否在稳定工作状态，以此绘制出压气机的流量特性曲线，如图5与图6。

由图5和图6中可以得出以下几点。

（1）在20 000～40 000 rpm的转速区间内，压气机的增压比较低，同时稳定工作区也很小，在压气机流量变化较大时，压气机效率下降明显。

（2）在60 000-100 000 rpm的转速区间内，压气机增压比较高，稳定工作区也比较大，压气机效率随流量变化较缓慢。

（3）在120 000-140 000 rpm的转速区间内，压气机增压比很高，但由于此时叶片吸力面气流已超音速，在叶片吸力面产生斜激波，同时由于叶片攻角为负，在叶盆部分出现气流分离，这两种现象堵塞叶片流道，使叶片几乎无法稳定工作。表现为即使出口平均静压持续降低，压气机流量也不上升，此时压气机出现堵塞，增压比与效率均直线下降。

由此可以判断出压气机在各个转速下的稳定工作区。此结果与原发动机慢车转速30 000 rpm，最大转速110 000 rpm的性能参数相吻合，当飞机在巡航状态时，压气机正好处于60 000～100 000 rpm这一稳定工作转速之间。

3.2 不同流量下压气机内部流场与性能

取100 000 rpm这一特定转速分析在不同流量条件下压气机的流场变化以及性能变化。

从图7可看出，压气机的增压比在进入喘振状态以前随着流量的下降几乎直线上升，但是从效率曲线可以看出，压气机的效率存在一个最大值，在该点左边时，压气机向喘振边界靠近，效率降低；在该点右边时，压力面开始出现气流分离堵塞通道，导致压气机效率下降。当流量超过某一值时压气机通道堵塞严重，此时改变出口静压压气机流量变化不明显，但效率与增压比却直线下降。

（1）最佳工况点压气机内部的流场。

100 000 rpm下效率最高的工况代表在此转速下，该工况点时压气机工作最稳定，流场最通畅，此时压气机的各项损失最小，是压气机的最佳工作状态。

（2）喘振边界点压气机内部的流场。

当流量减小到喘振边界时，叶片吸力面开始出现气流分离，气流分离产生的分离涡堵塞气流通道使压气机效率下降，此时若继续降低流量，压气机将进入喘振状态，吸力面出现大面积气流分离，流道被堵塞，增压比与效率均急剧下降，压气机进入不稳定工作状态。

（3）大流量工况下压气机内部流场。

流过压气机的流量过大同样会导致压气机效率下降。在流量过大时，叶片出现负攻角，在压力面出现气流分离。由于压力面空气压力较大不容易产成类似吸力面气流分离时那样的分离涡，但同样会堵塞气流通道。

4 结论

使用BladeGen软件进行建模，并使用CFX TurboGrid软件对模型进行网格划分。并使用ANSYS CFX软件对划分完毕的网格进行计算。在CFX中通过改变压气机出口平均压力以及压气机转速得到一系列工况点绘制出流量特性曲线并对其进行分析。并选择典型转速分析不同流量条件下压气机内部流场以及压气机性能变化。经过以上工作可以得出以下结论。

（1）压气机的增压比与压气机转速与流过压气机的流量密切相关，在一定范围内转速越大，压气机的增压比越高，流过压气机的气体流量越大，压气机增压比越低。

（2）压气机在每一个特定转速下其效率均随着气体流量的改变而改变，在每一个特定转速下均有一个最佳工况点，此时压气机效率达到最高，无论工况点移动到该点左边或是右边，压气机效率都将下降。

（3）压气机工况点向左移动至喘振边界附近时，最先在叶尖附近出现分离涡，同时尾迹涡流进一步扩大。若工况点进一步向左移动，此时叶尖泄漏涡、叶片吸力面出现的分离涡、尾迹涡可能混合在一起，产生一个较大的湍流区阻塞流道，使压气机效率与增压比急剧下降。同时由于湍流区的阻碍，其前方气流发生偏转，可能出现旋转失速现象。压气机工况点向右移动时流量过大，压气机可能出现堵塞现象，叶片进入负攻角状态，气流在压力面出现分离，在压力面位置产生湍流区。当流量进一步增大，湍流区也不断扩大，最终与尾迹涡流混合，完全堵塞压气机流道。

使用CFD方法的优点有成本低、准确度高、实验周期短等，但其与实际情况依旧存在差别。现今3D打印技术飞速发展，可以利用3D打印快速成型技术在CFD计算之后将叶轮模型制作出实体，在试验台上进行测试并计算结果。

参考文献

[1] 赵会晶，王志恒，席光，等.叶片前缘倾掠对离心叶轮气动性能的影响[J].西安交通大学学报，2015，49（11）：1-7.

[2] JH Kim，JH Choi，et al.Multi-objective optimization of a centrifugal compressor Impeller through evolutionary algorithms[J].IMechE Part A：J. Power and Energy，2010（224）：711-721.

数值仿真篇5

关键词：简支工字型钢梁；内力分析；SAP2000

0 前言

某简支工字型钢梁，利用有限元软件分析一型钢梁的应力，钢梁的材料特性如下：弹性模量E=2.1e11N/m2，泊松比∪=0.3,屈服强度fy=3.45e8N/m2；长度为6米，翼板宽为0.2米，腹板厚度为0.08米，高度为0.3米。由于工字型钢梁变形复杂受力多变，为较好把握该结构的受力特点，现基于通用有限元程序对该结构进行内力分析，现把该结构计算模型的简化、有限元模型的建立以及求解做一介绍[1]。

1 创建三维模型

三维模型的建立是通过几何模型的建立，在几何模型的基础上根据结构的特点进行几何模型的单元划分，从而建立与几何模型相应的有限元模型。根据该结构的特点，采用杆系结构对该结构进行几何模型简化得到相应的计算几何模型， 几何模型建立以后就可以根据几何模型的尺寸进行单元的划分。

2 分析条件的引入

2.1 材料和截面属性

在Module列表中选择Property功能模块，此模块中可以定义钢梁的材料本构模型和截面属性，并将截面属性赋予在相应的区域上填入弹性模量，泊松比、屈服强度等点击确定即可创建材料，之后就可以在Create Section对话框中完成截面属性的定义，并给部件赋予截面属性，此时模型由白色变成绿色。

2.2 定义荷载及约束条件

在环境栏的Module列表中选择Load功能模块进行荷载及边界条件的定义。第一步是施加荷载，分析步的载荷类型设为Pressure，点击continue，在Magnitude后面输入3.5e5，保持其他参数不变，退出Edit load对话框，完成型钢梁荷载的定义，加载后便是定义边界条件，将step设为initial，将施加边界条件的方式选为位移/转角形式，保持剩余参数默认值不变，点击continue，在Edit Boundary condition对话框中，选中U1=U2=UR2=UR3=O,对选中面施加铰接约束，点击OK，完成对所选截面边界条件的定义。该结构约束条件是一个固定铰支座其余采用活动铰支座，固定铰支座设置在梁的左端。结构约束的施加在梁的节点上，是通过节点约束选项进行施加完成[2]。

2.3 定义装配件及分析步骤

在环境栏的Module列表中选择Assembly功能模块，，在Create Instance对话框中完成装配件的定义。 在环境栏的Module列表中选择Step功能模块，出现如图对话框，选择step，在下拉菜单中选择static general项，然后保持其他参数不变，点击continue，进入下图第二个界面，保持参数不变，点击确认即可。

3 有限元模型的建立

在环境栏的Module列表中选择Mesh功能模块，对模型进行网格划分。将环境栏中的Object项设为Part：beam，即为部件beam划分网格，首先是要对型钢梁进行分割，将上翼缘和腹板，下翼缘和腹板分割开来，之后在边界上布置种子，在Approximate glolbal size后面输入0.05，保持其余参数默认值不变，点击Apply，模型即按要求布满种子，点击OK后，点击划分网格按钮，模型按照前面定义的种子自动划分网格，模型由绿色变成青色。

在环境栏的Module列表中选择Job功能模块进行作业提交。首先是创建分析作业，然后是提交分析，点击Submit，对话框中的状态提示由Submitted变为Running最终显示为Completed，点击对话框中的Results，自动进入Visualization模块。分析后需要对分析的结果进行处理，对型钢梁已经进行过运算分析，节点变形后的网格模型如上图5所示。在显示出应力云纹图的同时，还可以显示节点的Mises应力随时间变化的的关系曲线，选择型钢梁跨中底部中间节点，即最先屈服节点。在XY Data from ODB Field output 对话框中的Position标签页下，点击U：Spatial displacement旁边的三角形，选中U2，点击Plot，绘图区显示此节点也即型钢梁跨中底部中间节点Y方向挠度随时间变化的关系曲线：

5 结语

通过对该简支工字钢梁的数值模拟分析可以得出以下结论：

1）有限元软件可以方便地对桥梁结构进行仿真数值分析获得较理想的效果，与传统的手算比较获得较大的经济效应。

2）在钢结构的分析中数值分析具有更广泛的应用，比如静动力分析、非线性分析以及疲劳分析等都能进行全面的分析与获得满意的结果，建议提出更多的合适的单元来发展有限元程序在结构分析中的应用。

参考文献：

数值仿真篇6

关键词： 物质点法； 无网格法； 超高速碰撞； 侵彻； 爆炸； MPM3D

中图分类号： O389；TB115.1文献标志码： B

3D simulation based on material point method for

impact and explosion problems

ZHANG Xiong, LIAN Yanping, YANG Pengfei, LI Jinguang,

ZHANG Yantao, WANG Hankui, LIU Yan

(School of Aerospace, Tsinghua University, Beijing 100084, China)

Abstract： Based on the characteristics of Material Point Method(MPM) in the simulation of hypervelocity impact and explosion problems, the extension on MPM and its application are introduced, including the application hypervelocity impact problems by using MPM, material point finite element method(MPFEM), hyprid MPFEM, an adaptive particle splitting scheme for MPM, contact algorithm based on local multiple background mesh, and parallel MPM algorithm. Based on the improvement, a 3D explicit parallel simulation software MPM3D is developed for impact and explosion problems. C++ is used to develop MPM3D and the graphical user interface PeneBlast is developed by Qt and VTK. MPM3D can run on different platforms such as Windows, Linux, Mac OS, and so on. Many examples about hypervelocity impact, penetration, explosion, slope failure and metal cutting verify the reliability and accuracy of MPM3D. MPM3D can be an effective design tool for spacecraft protection on space debris, conventional weapon development and protection, and so on.

Key words： material point method; mesh-free method; hypervelocity impact; penetration; explosion; MPM3D

0引言

材料与结构在冲击爆炸载荷作用下的动态响应问题是几何、材料和边界条件均为非线性的多物理场强耦合问题，涉及高应变率、高压、高温、相变乃至化学反应，气体、液体和固体等多种物质间相互耦合甚至混合，材料不但会严重扭曲和破碎，还会熔化甚至气化.在这类问题的研究中，数值模拟比其他研究手段更经济，更便于观察，能突破试验和理论研究的局限性，正发挥越来越重要的作用.

数值模拟方法可分为拉格朗日法和欧拉法两大类.拉格朗日法中计算网格随物质一起变形，可方便地跟踪材料界面和引入与变形历史相关的材料模型，但对于涉及特大变形的问题会因网格严重畸变而产生数值求解困难，且难以有效模拟材料的破碎、熔化和汽化等行为；此类方法的代表性程序为DYAN.欧拉法中计算网格固定在空间中，不存在网格畸变问题，但不易跟踪材料界面，且非线性对流项也导致数值求解困难.国内学者在此方面做了大量研究工作，宁建国等［1］开发出欧拉型三维爆炸与冲击问题数值模拟软件EXPLOSION-3D；北京计算数学与应用物理研究所开发出流体动力学程序MEPH2Y和MEPH3D ［2-3］；王景焘等［4-5］基于高精度时空守恒元解算法开发出SUPER CE/SE程序.

近年来，无网格算法在国内引起广泛关注并被应用于冲击爆炸领域.LIU等［6］用光滑质点流体动力学（Smoothed Particle Hydrodynamics，SPH）方法对爆炸问题进行大量工作；杨秀敏［7］用SPH方法并结合有限元法开发出爆炸冲击问题数值模拟软件EF3D和EP3D；王宇新等［8-9］用物质点法［10］（Material Point Method，MPM）模拟冲击爆轰等问题.

本文简要介绍MPM和本课题组针对超高速碰撞、侵彻和爆炸问题在MPM方面的研究，重点介绍所研发的三维显式并行MPM数值仿真软件MPM3D及其在超高速碰撞、侵彻和爆炸等方面的应用.

1MPM

基于FLIP(Fluid Implicit PIC)方法［11］发展起来的MPM采用拉格朗日法和欧拉法双重描述，将物体离散为一组在空间网格中运动的质点.这些质点携带所有物质信息，如质量、速度、应变和应力等，其运动表示物体的变形和运动.空间网格即背景网格用于动量方程的求解和空间导数的计算，不携带任何物质信息.背景网格可在空间中固定，也可根据问题的特点按照某种方式自由布置.在每个时间步中，将物质点和背景网格完全固连，将物质点的物理量映射到背景网格节点上建立背景网格节点的运动方程，求解运动方程后再将背景网格节点的物理量映射回物质点，得到下一时刻物质点所携带的物质信息.这一步物质点和网格节点一同变形，没有相对运动，是拉格朗日求解，避免欧拉法中处理对流项的困难，且很容易跟踪物体的界面.物质点携带所有物质信息，故在下一个时间步中可以丢弃变形后的背景网格，仍采用规则的背景网格，从而避免拉格朗日法中因网格畸变而产生的数值困难.可见，MPM结合拉格朗日法和欧拉法各自的优势，适于分析特大变形及流动问题.

MPM采用线性形函数，其导数在格子(cell)之间不连续，因此当质点跨越格子时会对系统产生扰动，称为跨格子噪声（cell-crossing noise），导致精度降低，甚至会使部分问题结果完全失真.BARDENHAGEN等［12］针对此问题发展广义插值MPM（Generalized Interpolation MPM，GIMPM），有效抑制由于质点跨越格子而引起的数值噪声.

由于速度场的单值性，MPM自动满足无滑移接触条件，便于处理碰撞接触问题，特别是在超高速碰撞问题的模拟中，弹靶间的滑移可以忽略，采用MPM无须采用特别的接触算法即可达到较好的模拟效果.为处理滑移接触，HU等［13］引入多重背景网格概念，各物体在接触节点处的法向运动量在公共背景网格上求解，而切向运动量则在各物体的背景网格上求解，既保证各物体在接触点处的法向速度和法向加速度连续，又允许它们的切向运动相互独立；BARDENHAGEN等［14］基于多重速度场的思想给出MPM的接触/摩擦/分离算法，用于颗粒材料的数值模拟；NAIRN［15］将接触算法扩展到物体内部表面之间的自接触中，用于分析裂纹扩展问题.目前，MPM已应用于冲击、爆炸、裂纹扩展、材料破坏(脆性破坏、分层复合材料脱层破坏、层裂破坏)、多孔材料与颗粒材料和生物力学等的数值模拟中.

2对MPM的扩展

2.1超高速碰撞

在空间碎片防护设计中，研究碎片云对航天器的作用非常重要.通过侵蚀算法，有限元法可模拟弹丸对靶板的侵彻过程，但无法得到超高速碰撞所形成的碎片云形貌.通过引入Johnson-Cook材料模型和Mie-Gruneisen状态方程，将MPM应用于弹丸超高速碰撞薄板问题，能较好地模拟薄板背后形成的碎片云形貌，展示MPM在该领域的应用前景.［16］

2.2物质点有限元法

MPM虽然能避免拉格朗日法的网格畸变问题，但对于小变形问题，其效率和精度均低于有限元法，因此有必要将两种方法结合.通过在大变形区域中引入背景网格，提出物质点有限元法(Material Point Finite Element Method, MPFEM).［17］物体初始用有限元离散，但在可能发生大变形的区域预先布置背景网格；将进入背景网格区域的有限元节点视为物质点，动量方程在背景网格上求解，而在其他区域，动量方程仍在有限元网格上求解.MPFEM充分发挥MPM处理特大变形的能力和有限元法的高效性.

2.3杂交MPFEM

在将MPM扩展应用于钢筋混凝土的侵彻问题时，钢筋尺寸相对于混凝土尺寸较小，给离散带来困难.为此，将有限元（如杆单元）引入到MPM中，提出杂交MPFEM.［18］单元的节点与质点相同，将其动量方程在背景网格上积分，但其本构则在单元的高斯点上计算而不在质点上.通过引入杆单元，可将MPM成功应用于钢筋混凝土侵彻的数值模拟中.

2.4质点自适应分裂

在MPM中，当2个物质点间的距离大于背景网格的节点间距时，它们之间将没有相互作用而发生数值断裂.针对该问题提出MPM的自适应分裂质点方案［19］，使质点可根据累计应变在变形最大的方向上自适应分裂，从而避免数值断裂，并将其成功应用于激波管和聚能射流等爆炸问题的数值模拟.

2.5接触算法

对于中低速侵彻问题，标准MPM中的无滑移接触约束将产生较大的侵彻阻力，因此在分析侵彻问题时需引入接触条件.［20］建立基于局部多重背景网格的接触算法［21］，即仅在接触区域建立局部多重背景网格并引入多重速度场的思想，在背景网格上调整速度场并求解接触力.同时，为避免MPM中提前接触的现象，给出计算即将接触的物体间物理距离的算法.该方法减少对存储空间的需求并提高接触算法的效率，改善接触判断的精度.

2.6并行MPM

使用并行计算技术可实现更大规模问题的模拟.MPI为分布式内存并行技术，适于在Cluster机器上运行.基于MPI给出针对Cluster机器的MPM并行算法，该算法涉及局部背景网格分解，每个时间步均需要进行子区域重复节点之间的通信，而且还要建立复杂的并行I/O.OpenMP为共享内存并行技术，适于在多核计算机或者SMP机器上运行.基于OpenMP提出针对共享内存计算机的MPM并行算法.［22-23］OpenMP支持增量并行，且不涉及复杂的区域分解，因此已逐步在显式有限元计算、分子动力学计算和流体力学计算并行化中得到应用.

3MPM数值仿真软件MPM3D

在算法研究的基础上，采用C++开发出冲击爆炸问题的三维显式并行MPM数值仿真软件MPM3D，并为其开发GUI系统PeneBlast，见图1.二者均采用跨平台的开发工具和程序库，可运行于Windows，Linux和Mac OS等多种平台上.

3.1PeneBlast

PeneBlast基于Qt和VTK开发.Qt是跨平台的C++ GUI应用程序框架，跨平台特性优良、模块化程度高、可重用性好.VTK是开源的跨平台可视化工具包，可用于处理三维计算机图形、图像和可视化.

PeneBlast主要用于前处理和模型求解.通过PeneBlast，用户可完成从建立三维实体模型、物质点离散、从材料库中选择材料、施加边界条件、设置求解参数到生成供求解器读取的输入文件等一系列前处理过程.其中，材料库中包含常用的材料参数，用户可直接在材料库中选择使用某种材料用于仿真，也可拷贝、修改或添加新的材料到材料库中.

生成输入文件后，可以在界面上调用求解器MPM3D对所建立的物质点模型进行求解，并实时监控整个求解过程；可以暂停、重启动或中断求解器的运行.通过PeneBlast生成的供MPM3D读取的输入文件采用XML格式，易于扩充和读取.同时，PeneBlast还可导入已有的输入文件并在此基础上继续建模，PeneBlast也可对输入文件进行格式检查，并直接编辑输入文件.

后处理采用开源的多平台数据分析和可视化应用程序ParaView，它可定性或定量分析数据并实现可视化.MPM3D将计算结果保存为h5格式和VTK PolyData格式，供ParaView进行数据分析与可视化处理，绘制和导出变形动画、云图，绘制时间历程曲线以及提取某些点的变量值等.

3.2MPM3D求解器的主要功能

MPM3D采用C++编写，使程序具有更好的可维护性和良好的可扩充性，且易于添加新类以增加新的功能.目前，MPM3D中具有多种材料模型，如弹性及弹塑性模型、Moony-Rivlin超弹性模型、Johnson-Cook模型、Drucker-Prager岩土模型、HJC混凝土模型、高能炸药模型、Gurson模型以及可用于流体模拟的空材料模型等.冲击爆炸问题涉及高温、高压和高应变率，相对于低速下的材料响应，热力学效应更明显，故需用状态方程描述压强、密度和内能之间的关系.MPM3D中包含用于高压固体模拟的Gruneisen状态方程、理想气体状态方程、多项式状态方程，用于爆轰产物模拟的JWL状态方程和用于模拟混凝土的HJC状态方程及Palpha状态方程.为处理材料断裂破坏问题，MPM3D具有最大静水拉力失效、最大等效塑性应变失效、最大主应力/剪应力失效，最大主应变/剪应变失效以及Johnson-Cook损伤失效模型.为处理材料的随机特性导致其损伤破坏的随机性问题，MPM3D支持随机失效设置，目前支持Gauss和Weibull 2种随机分布.

MPM3D中爆炸采用基于CJ理论的理想爆轰，起爆方式为程序起爆.用户可定义多个起爆点，也可定义起爆面.为处理冲击波，添加人工体积黏性.

MPM3D可对某物体设置初速度、常速度或者施加随时间变化的载荷及加速度场（如重力场等）.

MPM3D的背景网格可固定在空间不动，也可根据物体所占的区域动态变化.用户可设置是否使用GIMP算法和是否启用接触算法.时间积分为可变时间步长的显式积分，支持USF，USL或MUSL格式.［15］用户既可采用OpenMP并行设置实现在共享内存计算机上的并行计算，又可采用MPI并行设置实现在计算机集群上的并行计算.

4应用实例

4.1超高速碰撞

4.1.1弹丸超高速碰撞薄板――碎片云

一个半径为7.5 mm，质量为20 g的铅弹以6.56 km/s的速度撞击厚度为6.35mm的铅靶.30.6 μs时碎片云的数值计算结果和试验结果见图2.［24］由MPM3D得到的30.6 μs时碎片云前端与靶板的距离和碎片云宽度分别为198 mm和142 mm，其试验结果分别为200 mm和145 mm.

4.1.2弹丸超高速碰撞厚板――成坑

一质量为0.5 g球形铜弹丸以6.0 km/s的速度撞击厚度为40 mm，长和宽均为80 mm的铜靶.弹坑参数的数值计算结果与试验结果［25］见表1，可知二者吻合较好.32 μs时碰撞变形结果切片见图3.

4.2侵彻问题

4.2.1钢珠侵彻薄板

一个直径为10 mm的钢珠正侵彻一个厚度为1 mm，直径为178 mm的薄圆板.靶体的最终变形数值计算结果与试验结果比较见图4.［26］试验所得靶体的变形高度与靶体的穿孔直径比值h/D为0.84，MPM3D计算得到的值为0.85.

(a)试验结果

(b)MPM3D 计算结果

4.2.2金属靶斜侵彻

一长为88.9 mm，直径为12.9 mm的尖拱钢弹以575 m/s的速度撞击倾角为30°，厚度为26.3 mm的铝合金靶.194.4 μs时的数值计算结果和试验结果见图5.［27］由MPM3D得到的剩余弹速为453.4 m/s，与试验值455 m/s吻合很好.

(a)试验结果(b)MPM3D 计算结果图 5194.4 μs时尖拱弹侵彻铝靶变形结果比较

Fig.5Deformation comparison for ogive-nose projectile penetrating aluminum target at 194.4 μs

4.2.3钢筋混凝土靶正侵彻

一个长度为144 mm，直径为25.4 mm的尖拱钢弹正侵彻厚度为178 mm的钢筋混凝土靶板.弹体以不同速度入射时，穿过靶体后剩余速度的数值计算结果与试验结果见图6［28］；当弹体以初速度749 m/s打中钢筋时，钢筋混凝土的损伤结果见图7.

4.3爆炸问题

4.3.1激波管

激波管由薄膜分为2个初始静止的、具有不同密度和压力的区域（该算例中各参数和结果值采用无量纲描述），见图8.在时间t>0时，薄膜破裂，激波和接触面会以不同的速度从左向右运动.通常观察t=0.143时刻的结果，此时激波前进大约0.25.

4.3.2一维板条TNT爆轰

一维爆轰过程经常用于测试爆轰模拟的数值方法的基准检验，长为0.1 m的板条TNT炸药在一端起爆，并在起爆端固定.0~14 μs时间内沿板条方向的压力曲线见图10，相邻曲线时间间隔为1 μs.数值模拟所得各曲线的峰值与解析解吻合较好.

4.3.3爆轰驱动飞片

爆轰驱动飞片在多个领域应用广泛，该过程涉及爆轰波的传播、爆轰产物的急剧膨胀及其与周围材料相互作用产生的大变形.非对称构型装药一维和二维模型的飞片终态速度与炸药爆轰速度比值随质量比变化的比较见图11，数值计算结果与Gurney公式及其修正的经验公式［29-30］预测值吻合很好.

4.3.4爆炸碎片模拟

金属壳在爆炸作用下生成的碎片极具破坏性，研究其机理对公共安全有重要意义.采用Gurson模型及Weibull随机失效机制模拟材料断裂破碎过程的随机性.采用Gurson模型的柱壳破碎过程见图12，三维球壳在爆炸作用下产生的裂纹和碎片形状见图13.

4.3.5聚能射流模拟

聚能射流在穿甲、石油开采等领域有着重要的应用，为此进行聚能射流模拟.聚能装药的初始构型见图14，半锥顶角为38°，炸药为TNT，外壳和药型罩材料为铜.射流尖端的速度和杵的速度可由BIRKHOFF等提出的理论模型进行估计.根据理论估计，本例中射流尖端速度应为3.7 km/s，杵的速度应为0.6 km/s.［29］模拟得到的21 μs时x方向速度分布见图15(a)，与理论估计基本吻合；射流中温度的分布情况见图15(b).

4.4其他问题

4.4.1边坡失效问题

在岩土工程中，常采用铝棒堆积物的流动模拟沙和土的颗粒流动.用Drucker-Prager模型对BUI［31］进行的干燥铝棒堆积物的坍塌流动试验进行模拟.数值计算与试验的比较见图16（a）和16（b）；试验得到的失效区和计算结果的定量比较见图16（c），表明坡面及坡体失效区的数值计算结果与试验结果吻合.

4.4.2金属切削问题

切削加工技术是工业装备过程中最重要和普遍的技术之一.切屑过程伴随着高速、高应变率、断裂和高温，极易产生绝热剪切带进而产生锯齿状切屑.模拟一个简单的切屑过程，见图17.在0.01 ms时第1条剪切带触发，随着更多的剪切带触发，锯齿状切屑形成.

4.4.3碳纳米管复合材料

碳纳米管有着突出的力学性能,将其掺杂于基体中可形成高性能的复合材料.但碳纳米管细长比大，在基体中均以弯曲形式存在，给基于网格的建模带来极大的困难.在MPM3D中，利用塌落生成的碳管几何模型布点，可构建出碳管的计算模型；管间正交布点，可构建出基体材料模型.

碳纳米管复合材料的计算模型见图18，基体材料点未画出，碳管间黑点用于代表碳管间近距离的相互作用.管间不同连接强度的力-位移曲线见图19，可得出碳管间相互作用对复合材料整体模量的影响.

5结论

介绍冲击爆炸问题的显式三维MPM数值仿真软件MPM3D及其图形用户界面系统PeneBlast，可运行于Windows，Linux和Mac OS等多种平台上.大量的超高速碰撞、侵彻爆炸数值算例说明该软件的可靠性和准确性.该软件可用于超高速碰撞、侵彻和爆炸等问题的数值仿真，为航天器空间碎片防护设计、常规武器研发与防护工程设计等领域提供有效的数字化设计工具.

致谢：感谢马上博士、黄鹏博士及出站博士后马志涛对本课题所做的出色工作.

参考文献：

［1］宁建国, 王成, 马天宝. 爆炸与冲击动力学［M］. 北京: 国防工业出版社, 2010: 348.

［2］刘军, 何长江, 梁仙红. 三维弹塑性流体力学自适应欧拉方法研究［J］. 高压物理学报, 2008, 22(1): 155-178.

LIU Jun, HE Changjiang, LIANG Xianhong. An Eulerian adaptive mesh refinement method for three dimensional elastic-plastic hydrodynamic simulations［J］. J High Pressure Phys, 2008, 22(1): 155-178.

［3］何长江, 于志鲁, 冯其京. 高速碰撞的三维欧拉数值模拟方法［J］. 爆炸与冲击, 1999, 19(3): 216-221.

HE Changjiang, YU Zhilu, FENG Qijing. 3D Eulerian numerical simulation method of high speed impact［J］. Explosion & Shock Waves, 1999, 19(3): 216-221.

［4］王景焘, 张德良, 刘凯欣. 基于CE/SE方法的二维Euler型多物质流体弹塑性问题计算［J］. 计算物理, 2007, 24(4): 395-401.

WANG Jingtao, ZHANG Deliang, LIU Kaixin. A Eulerian approach based on CE/SE method for 2D multimaterial elastic-plastic flows［J］. Chin J Comput Phys, 2007, 24(4): 395-401.

［5］WANG Jingtao, LIU Kaixin, ZHANG Deliang. An improved CE/SE scheme for multi-material elastic-plastic flows and its applications［J］. Computers & Fluids, 2009, 38(3): 544-551.

［6］LIU Guirong, LIU Moubin. 光滑粒子流体动力学――一种无网格粒子法［M］. 韩旭, 杨刚, 强洪夫, 译. 长沙: 湖南大学出版社, 2005: 190-292.

［7］杨秀敏. 爆炸冲击现象数值模拟［M］. 合肥: 中国科学技术大学, 2010: 308-309.

［8］王宇新, 陈震, 孙明. 滑移爆轰问题无网格MPM法模拟［J］. 力学与实践, 2006, 29(3): 20-25.

WANG Yuxin, CHEN Zhen, SUN Ming. Numerical simulation of slippage detonation by material point method-MPM［J］. Mech Eng, 2006, 29(3): 20-25.

［9］王宇新, 顾元宪, 孙明. 无网格MPM法在冲击载荷问题中的应用［J］. 工程力学, 2006, 23(5): 46-51.

WANG Yuxin, GU Yuanxian, SUN Ming. Application of material point method to shock load problems［J］. Eng Mech, 2006, 23(5): 46-51.

［10］SULSKY D, CHEN Z, SCHREYER H L. A particle method for history-dependent materials［J］. Comput Methods Appl Mech & Eng, 1994, 118(1-2): 179-196.

［11］BRACKBILL J U, KOTHE D B, RUPPEL H M. FLIP: a low-dissipation, particle-in-cell method for fluid flow［J］. Comput Phys Communications, 1988, 48(1): 25-38.

［12］BARDENHAGEN S G, KOBER E M. The generalized interpolation material point method［J］. Comput Modeling Eng & Sci, 2004, 5(6): 477-496.

［13］HU W, CHEN Z. A multi-mesh MPM for simulating the meshing process of spur gears［J］. Computers & Structures, 2003, 81(20): 1991-2002.

［14］BARDENHAGEN S G, BRACKBILL J U, SULSKY D. The material point method for granular materials［J］. Comput Methods Appl Mech & Eng, 2000, 187(3-4): 529-541.

［15］NAIRN J A. Material point method calculations with explicit cracks［J］. Comput Modeling Eng & Sci, 2003, 4(6): 649-663.

［16］马上, 张雄, 邱信明. 超高速碰撞问题的三维物质点法［J］. 爆炸与冲击, 2006, 26(3): 273-278.

MA Shang, ZHANG Xiong, QIU Xinming. Three-dimensional material point method for hypervelocity impact［J］. Explosion & Shock Waves, 2006, 26(3): 273-278.

［17］ZHANG Xiong, SZE K Y, MA Shang. An explicit material point finite element method for hypervelocity impact［J］. Int J Numer Methods Eng, 2006, 66(4): 689-706.

［18］LIAN Yanping, ZHANG Xiong, ZHOU Xu, et al. A FEMP method and its application in modeling dynamic response of reinforced concrete subjected to impact loading［J］. Comput Methods Appl Mech & Eng, 2011, 200(17-20): 1659-1670.

［19］MA Shang, ZHANG Xiong, LIAN Yanping, et al. Simulation of high explosive explosion using adaptive material point method［J］. Comput Modeling Eng & Sci, 2009, 39(2): 101-123.

［20］HUANG Peng, ZHANG Xiong, MA Shang, et al. Contact algorithms for the material point method in impact and penetration simulation［J］. Int J Numer Methods Eng, 2011, 85(4): 498-517.

［21］MA Z T, ZHANG Xiong, HUANG Peng. An object-oriented MPM framework for simulation of large deformation and contact of numerous grains［J］. Comput Modeling Eng & Sci, 2010, 55(1): 61-88.

［22］HUANG Peng, ZHANG Xiong, MA Shang, et al. Shared memory OpenMP parallelization of explicit MPM and its application to hypervelocity impact［J］. Comput Modeling Eng & Sci, 2008, 38(2): 119-147.

［23］ZHANG Yantao, ZHANG Xiong, LIU Yan. An alternated grid updating parallel algorithm for material point method using OpenMP［J］. CMES: Computer Modeling Eng & Sci, 2010, 69(2): 143-165.

［24］Jr ANDERSON C E, TRUCANO T G, MULLIN S A. Debris cloud dynamics［J］. Int J Impact Eng, 1990, 9(1): 89-113.

［25］WINGATE C A, STELLINGWERF R F, DAVIDSON R F, et al. Models of high velocity impact phenomena［J］. Int J Impact Eng, 1993, 14(1-4): 819-830.

［26］SEO S W, MIN O K, LEE J H. Application of an improved contact algorithm for penetration analysis in SPH［J］. Int J Impact Eng, 2008, 35(6): 578-588.

［27］PIEKUTOWSKI A J, FORRESTAL M J, POORMON K L, et al. Perforation of aluminium plates with ogive-nose steel rods at normal and oblique impacts［J］. Int J Impact Eng, 1996, 18(7-8): 877-887.

［28］HANCHAK S J, FORRESTAL M J, YOUNG E R, et al. Perforation of concrete slabs with 48 MPa and 140 MPa unconfined compressive strengths［J］. Int J Impact Eng, 1992, 12(1): 1-7.

［29］MAYERS M A. Dynamic behavior of materials［M］. 张庆明, 刘彦, 黄风雷, 等, 译. 北京: 国防工业出版社, 2006: 165-170.

数值仿真篇7

关键词：有限元模型；仿真结果与讨论；光电系统

中图分类号：TP211 文献标识码：A

光电技术是以激光、红外和微电子技术为基础，由光学、电子、精密机械和计算机技术相结合而形成的一个高新技术领域。光电经纬仪或其它地基光电跟踪设备，在结构系统、伺服控制系统和光学系统设计等方面，均形成了较为完整的理论。而机载光电系统不仅包括精密仪器的设计，还要考虑精密仪器的载机环境，尤其是高低温环境的适应性问题。光机系统中的光学基座设计，作为光电元器件的承载平台，不仅要考虑安装在其上的光电元器件的安全性，保证光机系统的高精度与可靠性，还要考虑自身在机载环境下的热变形和温度场分布是否满足设计要求，实现系统在机载高低温环境中，清晰稳定成像跟踪需求。

目前，机载设备的高低温环境适应性问题主要通过两种途径解决，一种是采取改善环境或减缓环境影响的措施，通过添加辅助设备，如环控系统来改善光机所处的机载环境；另一种是选用耐环境能力强的结构、材料、元器件或工艺等，即通过改善系统的结构设计，提高光机系统的耐环境能力。本文主要通过第二种途径的研究，来提升光机系统高低温环境适应性。计算机数值仿真，可在设计初期就获取丰富的分析数据，对设计进行优化，且周期短、成本低，便于实现。

为此，基于ABAQUS结构分析模块对某机载光电系统光学基座高低温环境适应性进行了数值仿真研究。通过光学基座的数值仿真分析，在综合考虑结构热形变、强度、重量、及加工工艺等因素下，经过设计改进，获取一个满足系统要求的光学基座设计方案。

1 结构有限元模型建立

建立有限元分析模型包括三维结构模型的简化、网格划分、材料属性设置。之后依据具体分析类型，进行约束与载荷施加，选择计算类型进行仿真，分析结果。

1.1 模型简化

根据结构热传导、热对流，以及热辐射等方面的理论依据，结合具体的分析对象和分析内容简化有限元分析模型，确定被分析部件的简化方式与程度。本文研究光学基座的高低温环境适应性，所以文中光学基座为详细模型，其上光电元件均简化为等效质量块，系统有限元模型如图1所示。

1.2 网格划分

该系统三维模型简化后，以iges格式导入ABAQUS中，采用四面体网格和六面体网格进行划分。三个光学基座方案的网格划分结果如图2（a）～（c）所示。

1.3 材料属性

本文中设计方案的材料主要有铟36，硬铝合金7075，以及各简化部件的等效材料，见表1。由于各光学支路的镜座和光电器件的安装支架均采用硬铝合金材质，所以等效材料时，使用重量等效原理，与铝合金的材料性能相近。

1.4 边界与载荷施加

机载高低温环境有限元仿真的工况边界与载荷见表2。系统通过6个安装点处的不锈钢螺栓与机体固定约束，所以在结构热分析中，将根据不同材料的热膨胀系数不同，将该约束处理为螺纹孔内壁的平面内约束。光学基座的下表面安装光电元器件，上表面有载机提供的5℃气源吹过进行强迫对流，所以光学基座上表面的结构设计会影响整个系统的散热效率。

2 仿真结果与讨论

该光电系统正常工作时受载机高低温环境影响，可能会出现成像质量差，不能稳定跟踪与瞄准等，甚至结构的破坏。对于三种方案的结构热分析结果见表3。

该光电系统机载环境温度范围-55℃～ 70℃：系统设计之初光学系统对光机座热变形要求严苛，因此方案一选用了线膨胀系数较小的铟36钢材料。但该材料密度较大，经详细减重设计后，即使在结构局部加强筋仅2mm厚，且局部热应力已超过材料许用应力的情况下，其重量仍不满足要求，具体分析结果见表3。

方案一中的铟钢光机座虽热膨胀系数小，但其热传导性也差，温度场分布不均匀，温度梯度大，产生较大的热应力集中。且其重量超过了20kg的设计限制，固不满足设计要求。因此在方案一的基础上，综合材料力学和热力学特性，选用硬铝7075。该种材料是一种高强度的铝合金，尤其在T6状态下，可加工性好，且性价比高。方案二着重减小温度梯度和重量，使温度场更为均匀，重量满足设计要求。但是由于铝合金热膨胀系数较大，导致其最大热形变值达到0.42mm，超过0.2mm的设计要求，固亦不满足设计要求。

在方案二的基础上，针对热形变较大的问题，对光学基座的上表面加强筋进行了散热优化设计。并进一步进行减重设计和提高加工工艺性，减轻加工应力对光机座精度的影响，最终完成了方案三的设计。有限元结构70℃环境下的热分析温度场和热形变结果如图3所示。

综合以上三种方案，最终选择各方面满足设计要求的方案三，开展该光电系统的光学基座设计。该基座在满足最大热形变小于0.2mm的情况下，比原设计指标减重2.3kg，说明有限元分析手段是结构在严酷环境下设计和优化的有效途径。

结语

综上所述，本文使用有限元分析手段，对光电系统中的光机座进行有限元高低温环境下的有限元分析，根据分析结果，得出在有强迫对流的条件下，热膨胀系数小不是降低热形变的唯一途径，根据流场情况，进行合理的导热部件结构设计是降低设备重量和热变形的良好选择。且确定最终方案进行指导设计，样机在高低温试验中，可靠稳定工作。

参考文献

[1]王海晏.光电技术原理及应用[M].北京：国防工业出版社，2009.

[2]殷纯永，方仲彦.光电精密仪器设计[M].北京：机械工业出版社，1996.

[3]毛英泰.误差理论与精度分析[M].北京： 国防出版社，1982.

[4]吴晗平.光电系统环境与可靠性[M].北京：科学出版社，2009.

[5]李薇.某天线座砖台结构的计算机辅助设计[J].电子机械工程，1996（06）：33～40.

数值仿真篇8

关键词 声场；耦合；有限元；封闭声腔

中图分类号X5 文献标识码A 文章编号 1674-6708（2011）50-0170-02

1 概述

1.1研究背景与研究目的

随着近代工业的发展，环境污染也随着产生，噪声污染就是环境污染的一种，已经成为对人类的一大危害。噪声污染与水污染、大气污染、固体废弃物污染被看成是世界范围内3个主要环境问题，所以控制噪声污染已成为环境保护的重要内容。噪声污染按声源的机械特点可分为：气体扰动产生的噪声、固体振动产生的噪声、液体撞击产生的噪声以及电磁作用产生的电磁噪声。按照来源分，则可分为交通噪音、工业噪音、建筑噪音、社会噪音、家庭生活噪音污染。

1.2 封闭声腔结构-声耦合分析的一些求解过程

有限元分析法当今最常见的是基于格林函数法的分析与研究，以封闭声腔为模型,在考虑流固耦合作用的基础上,结合流体格林函数和Helmholtz方程及其边界条件,导出了各阶声压模态对应的声压振幅响应公式;结合结构格林函数和板的振动方程及其边界条件,导出了各阶板模态对应的速度振幅响应公式。Dowell[1]等建立了弹性薄板声腔系统的结构-声耦合理论模型，分析了耦合系统的固有特性，并进行了实验验证。Kim[2]等人发展了Dowell的理论，在前人的基础上，用阻抗和导纳方法分析了结构-声耦合问题，但其却没有对系统耦合特性以及影响系统耦合程度的因素作具体研究。1984年,美国通用汽车的Sung和Nefske[3]应用有限元方法对完整车身内部结构噪声进行了分析,并首次考虑了车身结构和声场的耦合作用。Kompella[4]从结构-声腔耦合的角度建立了车内声辐射数学模型，很明显，在这个问题的研究上，国外科学工作者确实领先了我国科学家一步。

2 建模和分析软件

2.1声固耦合问题概述

所谓声固耦合问题，简单地说，就是在外加载荷的作用下，使弹性结构振动，并通过振动辐射产生周围的声场，而辐射出的声场再反过来对结构产生作用力，这就是所谓的声固耦合。

3具有弹性板的矩形封闭声腔的声模态分析

3.1具有弹性板的矩形封闭声腔结构模型建立

本课题直接运用ANSYS建立了具有弹性板的矩形封闭声腔结构模型，其长为2m，宽为3m，高为2m，弹性板的厚度为1cm。

3.2模型的薄板定义和网格划分

首先，定义声腔薄板材料属性，设置为密度dens为“7 800kg/m3”，弹性模量ex为“200GPa”,泊松比nuxy为“0.3”，设置成封闭的矩形声腔结构。矩形声腔结构以空气为介质，在定义材料性能参数时设置单元类型为“fluid 30 3D”,材料属性设置为密度为“1.21kg/m3”，声速为“340m/s”。

然后进行“网格划分”，在对矩形封闭声腔内声场进行网格划分时，最大声场流体单元的尺寸应小于波长的1/12，每个声波波长内的声场单元数不应小于8。如果网格划分越密，用有限元方法得到的求解精度越高，但对计算机的性能要求也越高，计算时间长。综合考虑以上因素，设置“网格单元尺寸”为0.1m，用“mesh volumesfreepick all”命令对体自由划分网格，共有121个节点，1 321个单元。

下面一步则是声固耦合设定，具体操作为在命令栏输入“asel,u,loc,y,width,sfa,all,,fsi alls”，其含义为在腔内介质与弹性板之间设置耦合界面。

3.3矩形封闭声腔模态分析

3.3.1声模态分析步骤

第一步是施加约束，在命令栏输入“d,all,,,,,,ux,uy,uz,”即可，其具体含义为固定X轴Y轴Z轴。

接着进入分析计算模块对其求解。分析类型设置为“modal（即模态分析）”，在选项里选择“Unsymmetic（即非对称分析）”。

1）New Analysis[ANTYPE]――Analysis Type――Modal

2）Analysis Options――选择Unsymmetric

3.2.3模态计算与结果分析

现列举前十阶模态频率和模态振型，可从中看到具有弹性板的矩形封闭结构声腔在不同模态频率下的声压分布情况。

软件分析得到具有弹性板的矩形声腔结构有限元声模态振型，做出以下分析讨论：

1）从各阶模态振型中可以看出,具有弹性板的矩形封闭结构声腔呈横向对称。

2）第1阶模态，声腔中间偏上部位Y向声压最大，X向声压最小，整个声腔的声压较小；第2阶模态，声压中间偏上部位Y向声压较大，其余部位声压较小，声腔的声压也较小；第3阶声模态时，声腔中间偏上部位X向声压很小，总体声压也很小；第4阶声模态，与第一阶模态相似，只是由靠上部位移到了靠下的部位；第5阶声模态的声压呈纵向向上递减，整体声压较大;之后略同。

4使用谐波分析法对声压进行分析

4.1分析不同频率下声压分布的具体步骤

1）在命令栏输入antype,harmichropt,fullf,121,fY,100。点击回车确认，即在编号为121点处，施加一个方向为Y正方向，大小为100N的力。并且求出空气介质在此载荷下，20Hz~300Hz之间的声压分布；

2）在命令栏输入alls nsubst,10kbc,1HARF,20,300SOLVE。点击回车确认。即在20Hz~300Hz之间选取十个频率作为分析点，查看每一个点频率的声压分布，从20Hz起每增加28Hz进行一次仿真分析，计算十次，得出结果；

3）在ANSYS软件中点击General PostprocRead ResultsFirst Set，然后点击Load CasePlot ResultsContour PlotNodal Solu，在弹出的菜单中选取DOF Solution菜单下的Pressure即可，声腔内空气介质在20Hz~300Hz下时的声压分布。

下面一步则是使用谐波分析法对20Hz~300Hz频率之间进行频率扫描计算。

结果与数据：

4.2结论

由上述曲线可以看出，每隔约为140HZ则产生一次最大值，而只有简谐波频率和固有频率产生叠加时才会产生峰值，所以可以求得其固有频率约为140Hz，而总体上，曲线随着频率的增加而增大，所以，可以得出结论，在与固有频率叠加的点，声压会间歇性的达到峰值，而总体上，声压随频率的增加而增大。

4.3不固定X轴的边界条件下腔内声压分布

将之前计算声模态的步骤“d,all,,,,,,ux,uy,uz，”改变为“d,all,,,,,, uy,uz”，其含义是不限制X轴的边界条件，然后在命令栏输入antype,harmichropt,fullf,121,fY,100。点击回车确认，即在编号为121点处，施加一个方向为Y正方向，大小为100N的力。

4.4结论

在不限定X轴的情况下，在250HZ左右时，声压达到最大值，整体声压值与限定X轴比较，总体幅度下降。

4.5不同边界条件下的结果比较

限定边界条件时，弹性板所吸收的能量最小，所以辐射出的内声场声压最大，而限定的条件越少，弹性板因为其自身旋转或呈波浪移动，故吸收了较多的能量，声场内声压明显降低。故如若将矩形腔假象为工作室，则将弹性板的固定方式越牢固，内声场声压越高，而适当的减少连接强度，即可以达到吸收震动能量，减少噪音辐射的效果。

参考文献

[1]Dowell E H，Gorman G，Smith D A.Acoustoelasti-city General theory ,acoustic natural modes and forced response to sinusoidal excitation including comparison with experiment [J].sound Vib.1977：52(4)：519-542.

[2]Kim S M Brennan M J.A compact matrix formulation using the impedance and mobility approach for the analysis of structural-acoustic systems[J].SoundVib,1999;223(1):97-113

数值仿真篇9

关键词：仿真系统 数据融合 性能评估 指标体系

中图分类号：TP391.9 文献标识码：A 文章编号：1007-9416（2013）04-0189-02

数据融合技术在现代C4ISR系统中得到广泛的应用，共享多个平台上多种传感器的量测数据，可获得任何单一平台所无法获得的更高质量的信息[1]。多传感器数据融合可实现多目标复合跟踪，生成编队统一态势图，提高了武器打击的效能。共享信息的潜在能力被充分挖掘出来，信息优势可有效地转化为作战优势[2]。

随着目标和信源数量的增加，数据融合系统复杂度急剧上升[3]。需要对多传感器数据融合系统进行反复、全面的测试，验证其正确性，评估其性能和效能。依据此研究目标设计并实现多传感器数据融合仿真系统，目的是经济且高效的对数据融合算法、系统进行调试、测试和评估。

1 仿真系统的需求

一个完整的仿真系统，可以模拟战场环境，为数据融合系统提供输入，记录数据融合结果并评估，并把仿真过程和结果以可视化方式呈现，辅助对系统各个环节和存在问题进行分析。综上，数据融合仿真系统应具备以下主要功能：（1）想定元素的编辑，剧情和真值数据生成，仿真系统控制等；（2）各种作战元素的仿真，包括传感器仿真、平台仿真等；（3）战场信息的采集和回放；（4）战场态势显示，仿真系统状态等信息的可视化；（5）数据融合性能、效能的评估。

2 仿真系统分析与设计

该仿真系统应遵循可扩展、易维护和界面友好等设计原则[4]，具备灵活的配置方式，各类模型和算法可以通过更改配置实现“无缝”选择和变更；具备良好的通用性和开放的系统结构，支持多平台、多传感器的扩展和即插即用，支持多个数据融合系统并行测试和评估。

2.1 系统组成

根据需求和功能，仿真系统可划分为八个子系统，其组成如图1所示。

剧情生成子系统是整个仿真系统的基础，其功能主要包括：想定元素编辑，剧情生成，为平台模拟子系统、传感器模拟子系统提供初始化参数，并发送目标真值数据；驱动网络工作， 控制整个仿真平台的启动与终止；反馈仿真系统状态，生成仿真系统的时序控制信息。

数据采集与管理子系统可实时接收各子系统发出的数据，包括真值数据、模拟数据、仿真系统控制命令和状态信息，并存储到数据库相应的表中；可分段选择并变速回放各类数据；可对数据库进行维护与管理。

传感器模拟子系统可由若干个传感器模拟器组成，其数量由剧情配置。每个传感器模拟器接收剧情真值数据，变换到该传感器的测量坐标系下，并按照一定扫描周期发出仿真数据。

平台模拟子系统可由若干个平台模拟器组成，其数量由剧情配置。每个平台模拟器接收剧情真值数据，根据惯导、卫星导航或罗兰C导航的模型，模拟平台定位和姿态数据并发送。

网络管理控制子系统可筛选真值数据、测量数据和控制数据，实现网络数据统一管理。

数据融合评估子系统读取仿真数据库中存储的数据，包括目标真值、传感器测量值、平台定位值和融合输出值，定量评测融合系统的性能。

态势可视化子系统可动态显示战场二维统一态势，包括目标、平台的运动轨迹和数据融合子系统输出航迹，实时显示仿真系统状态和仿真中间结果等。

数据融合子系统由若干个可独立工作的数据融合模块组成，其接收传感器模拟器输出点迹、航迹数据和平台模拟器输出定位和姿态数据，对多目标进行跟踪形成航迹，支持从融合算法库动态配置多种融合算法，支持点迹、航迹、点航混合融合等多种处理模式。

2.2 系统架构

由于剧情真值和仿真模拟数据输出频率高、数据量大，如果和被测多个数据融合模块放在一个网中，会产生不可预料的延迟。为了保证系统的实时性，物理上分成A、B两个局域网，由网络管理控制子系统连接。可以过滤A网的真值、控制和状态信息，只发送仿真模拟数据到B网，减轻了网络负载，提高了系统的灵活性。系统数据流图如下：

3 仿真系统的关键技术

3.1 传感器模拟器的构建

仿真系统是否能够对实践起指导作用，关键在于是否真实还原了作战环境。因此，对各类传感器的模拟成为重点和难点。在设计层面，传感器模拟子系统需要尽量模拟实际的外场工作环境，如杂波、电磁干扰、大气干扰、多径等。实现层面，需要采用面向对象思想，把每种传感器作为一个类，根据剧情生成的配置信息，创建相应的类的实例，并初始化参数，使该子系统易于继承和扩展。下图为传感器模拟器处理流程。（如图4）

3.2 融合评估指标体系的构建

要全面、合理的评估多传感器数据融合系统的性能和效能，必须建立相对完备的评估指标体系，建立数学模型，并运用计算机算法实现。一般情况下，评估指标的选取与系统应用背景、融合体系结构和融合算法密切相关[5]。在建立指标体系过程中，必须考虑指标体系的完整性、层次性和针对性等问题。下图为评估指标体系。

4 仿真系统的实现

该仿真系统由多台计算机组成独立的A、B两个以太网，并通过网络管理控制计算机连接。A网中真值信息采用广播方式发送，控制信息采用TCP方式发送。B网中仿真模拟数据和数据融合输出采用组播方式发送。软件运行在WindowsXP-SP3操作系统下，在Visual Studio 2010环境下，采用面向对象思想使用C++和C#语言开发，数据库采用Oracle 10i。

运行一个仿真系统的实例，剧情生成呈倒三角形分布的3个平台，分别搭载3部传感器，两批目标在分别在西北和东北方向做机动运动。

通过此实例的运行得出：该仿真系统配置管理方便，符合实际，为多传感器数据融合调试、测试和评估提供了动态、实时、直观的仿真环境。

5 结语

本文设计和实现了一种开放式的数据融合仿真环境，该仿真系统架构合理，操作方便，通用性良好。一方面，可用来测评不同的融合算法的性能，给出定量评价；另一方面，可验证仿真中机动模型、杂波模型和系统误差模型等的正确性，进行参数寻优。今后，可在此基础上将通信链路的仿真加入网络管控子系统，将可视化子系统扩展到三维方式等，以期实现一个实用性更强的仿真系统，从而指导工程实际应用。

参考文献

[1]韩崇昭，朱洪艳，段战胜.多源信息融合[M].第2版.北京：清华大学出版社，2010.

[2]吴泽民，张磊，刘熹.多源多目标信息融合仿真平台的分析与实现[J].军事通信技术，2008，29（3）：12-16.

[3]邹伟，刘兵，孙倩.多源信息融合能力评估关键技术综述[J].计算机与数字工程，2010，38（3）：1-5.

[4]程云鹏，肖兵.一种多传感器数据融合算法评估平台的设计[J].空军雷达学院学报，2006，20（1）：11-13.

数值仿真篇10

关键词：电力系统；数值仿真；参数调整；动态负荷中图分类号：TM32 文献标识码:A

1 引言

电力系统暂态稳定的研究是电力系统安全分析、规划设计和运行等工作的重要内容。目前在电力系统稳定性和安全性分析中时域数值仿真法由于直观性，是目前应用最为广泛的方法[1]。但是国内外几次事故和试验均表明实测数据和仿真数据有较大偏差[2]。因此电力系统模型和参数的校正研究就显得越来越重要。

为验证仿真模型的准确性，中国东北电网公司分别进行了两次500kV的大扰动试验[3]，利用实验数据开展了不同模型对电网输电能力的影响[4]和不同模型参数对时域仿真造成的影响[5]等等大量的研究工作。与负荷模型相比，发电机、励磁系统、PSS、发电机调速器、直流控制系统已经建立了详细的数学模型，其参数的不确定性相对较小，因此在仿真计算时首先应对这些元件参数进行灵敏度分析及参数调整，在此基础上进行负荷模型的校正。负荷模型通常是根据经验确定，并定性的选定模型中的参数，然后通过对本网发生典型事故的仿真计算对负荷模型进行校正。

现代电力系统广域测量系统（wide area measurement systems，WAMS）的建立，能够直接量测到发电机相角轨迹。WAMS系统为分析电力系统的内部特性，提供了海量的数据。就如何充分利用WAMS系统及其量测数据的研究，已经有了一定的成果。如运用实测功角轨迹来判定和衡量系统的稳定程度[6，7]；将WAMS系统运用在实时监测中[8]；将其用于参数量测中[9]；基于量测轨迹求取轨迹灵敏度的方法[10]等等。

本文提出了一种基于量测发电机相角轨迹调整感应电动机负荷模型参数的方法。由于PSASP（Power System Analysis Software Package）是中国电力科学研究院开发的一款优秀的电力系统综合分析软件，是目前国内电力设计和运行部门应用最广泛的程序，因此本文选用它作为仿真工具进行分析和应用。本文首先介绍了量测和仿真相角轨迹的误差指标；然后介绍了调整感应电动机参数的方法；最后通过电力系统动态模拟实验验证了该方法的有效性。

2 相角轨迹误差指标体系

量化电力系统数值仿真与实测轨迹的误差是电力系统仿真模型参数调整的基础。本节根据仿真和实测相角轨迹的振荡特点，结合物理本质，从整体和局部两部分对误差进行介绍。

2.1整体指标

对相角全部变化过程进行比较和分析，定义为误差能量指标（Error Energy）[11]

（1）

式（1）中是仿真变量序列，是实测变量序列，是变量的稳态值，N是实测变量与仿真变量个数。

2.2局部指标

2.2.1第一摆幅值误差（First Swing Magnitude Error）及第二摆幅值误差（Second Swing Magnitude Error）[12]

（2） （3）

式(2)（3)中FMagsimu是仿真值第一摆的摆动幅值,FMagmeasur是实测值第一摆的摆动幅值, SMagsimu是仿真第二摆的摆动幅值, SMagmeasur是实测第二摆的摆动幅值。

2.2.2第一摆摆动的周期误差（First Swing Period Error）及第二摆的周期误差（Second Swing Period Error）[12]

（4） （5）

式(4)（5）中FPersimu是仿真的第一摆的摆动周期,FPermeasur是实测第一摆的摆动周期, SPersimu是仿真的第一摆的摆动周期 , SPermeasur 是实测第一摆的摆动周期。

3 感应电动机负荷模型参数调整方法

与负荷模型参数相比，发电机及其控制系统模型参数容易测量和辨识。当除负荷模型以外的其它元件的模型和参数经过量测和辨识后，仿真与实测结果仍不符合，则考虑对负荷模型和参数进行调整。若动态部分采用感应电动机模型，则应对电动机模型参数进行调整。

3.1灵敏度分析

在调整感应电动机模型参数时，分析发电机相角轨迹对感应电动机参数的灵敏度，并按其数值大小对参数排序。该排序反应了参数对仿真轨迹影响大小。在调整幅度相同的情况下，优先调整灵敏度大的参数会对相角轨迹产生较大影响。

3.2 摄动分析

为分析感应电动机模型参数对仿真相角轨迹的影响，选感应电动机负荷模型任意一个参数在已知数值的基础上下浮动30％，在其它参数不变的前提下，做PSASP仿真并记录所得的发电机相角轨迹。以所得的轨迹族，作为分析此参数在这组参数数值下的影响数据。然后将被变化参数恢复原始值，选下一个参数上下浮动30％做同样的分析，直到将所有参数分析到。

为将参数摄动对误差指标的影响程度进行量化，将影响程度量化为：

（6）

y为参数影响程度，Ymax为参数摄动到上限时的相应的物理量，Ymin为参数摄动到下限时相应的物理量，Yinit为参数为初始值时相应的物理量。

3.3 确定参数调整方案

在确定感应电动机负荷模型参数的调整顺序后，对其参数做摄动分析，分析参数变化对相角轨迹的各个误差指标的影响方式。通过寻找使误差指标最小的变化方向和变化量，来确定参数的调节方式和调节量。

在上述分析中，选择参数的摄动域越宽摄动步长越小，得到的信息量越多，对该参数的调整越精确；但是计算量也会越大，计算时间也会越长。相反，分析中选择参数的摄动域越窄摄动步长越大，得到的信息量越少，对该参数的调整精度越差；然而计算量也会越小，计算时间也会越短。因此在参数摄动域和摄动步长的选择上，要考虑实际的调节需要和调节精度以便选择合适数值。

3.4 适应性分析

若基于一次故障轨迹对感应电动机负荷模型参数的调整也能减小其他故障时的仿真误差，则证明对其参数所做调整具有适应性。通过电力系统动态模拟实验发现，当故障发生在距动态负荷近的位置时，基于本文方法对电动机参数所做调整同样适应于当故障发生在距动态负荷远位置的情况。当故障位置在动态负荷附近时，动态负荷仿真模型所受激励剧烈，发电机相角轨迹对电动机参数的灵敏度比较大，易对电动机参数做出调整。由此认为在这种情况下对感应电动机负荷模型参数的调整应用于全网具有一定的适应性。

4感应电动机负荷模型参数调整方法验证

4.1 应用PSASP对图1系统进行仿真

图1 实验电气主接线图

按照图1系统在电力系统动态模拟实验室进行实验，4#为隐极发电机、3#为凸极发电机，静负荷为一组总额定功率4.5kw的电炉，动态负荷为两台大小不同总额定功率3.7kw的异步电动机，电动机的负荷为直流发电机带电炉。在整个实验过程中利用WAMS系统对系统各信息量进行采集和记录。

以母线6处发生200ms三相短路故障为例，分析发电机出口母线电压相角仿真与实测轨迹偏差。4#发电机模型采用PSASP自带的六阶模型，3#发电机模型采用PSASP自带的5阶模型，电动机采用PSASP自带的三阶感应电动机模型。该实验仿真所用发电机及其励磁系统的模型参数是实验室以往通过专门实验量测和辨识所得，具有较高准确性；三阶感应电动机模型则直接采用PSASP自带的经验参数。本文中所提到的发电机出口电压相角均以无穷大母线1的电压相角为基准值。

图2 4#发电机相角图3 3#发电机相角

仿真和实测相角轨迹对比如图2、3所示。通过比较仿真和实测相角轨迹，实测轨迹和仿真轨迹存在偏差，说明感应电动机模型参数存在一定的误差。

4.2 基于一次故障的感应电动机负荷模型参数调整方案

现将感应电动机负荷模型所有参数进行灵敏度分析和摄动分析，确定的参数调节方案如表1 所示。图4、5为参数调整前后相角仿真轨迹的对比。表2为参数调整前后的误差指标。通过比较认为调整后的仿真相角轨迹和实测相角轨迹的误差明显变小，说明调整后的参数比原参数对数值仿真更加有效。

感应电动机负荷模型参数调整方案

表1 原始参数与调节后参数对比

图4 4#发电机相角图5 3#发电机相角

表2 参数调整前后各种误差指标的对比

4.3 调整方案在其他故障时的适应性验证

按照如图1所示系统，在5母线发生250ms三相短路并记录动态轨迹。分别用表1中原始参数和调整后的参数仿真发电机相角轨迹，其结果如图6、7和表3所示

图6 4#发电机相角图7 3#发电机相角

表3 参数调整前后各种误差指标的对比

通过对5母线三相短路的数值仿真，发现感应电动机负荷模型参数调整前后发电机功角轨迹的变化小于6母线故障时发电机功角轨迹的变化。但是，基于6母线短路故障时对感应电动机负荷模型参数的调整方案对于5母线故障仿真也适用，增强了5母线故障时仿真的有效性。当故障位置越靠近需调整的动态负荷模型时，感应电动机负荷模型发生的动态响应越大，其模型参数对功角轨迹的影响量越大，越有利于调整其模型参数。

5结论

在发电机及其控制系统模型参数准确的情况下，本文提出了一种基于实测相角轨迹，运用误差指标调整感应电动机负荷模型参数的方法。本方法仅需要根据电力系统某一次故障时的实测相角数据，然后据此调节感应电动机仿真模型参数，只要保证在这组参数下仿真轨迹和实测轨迹误差指标充分小；那么就可以保证在调整后的参数基础上，所做仿真能充分反映实际系统的行为特性。最后通过电力系统动态模拟实验的实例分析验证了该方法的有效性。

参考文献

1 倪以信等.动态电力系统的理论和分析[M].北京：清华大学出版社.2002年.

2 Dmitry N. Kosterev, Carson W.Taylor.Model validation for august 10, 1996 WSCC System Outage.IEEE Trans on Power System[J],14(3). 1996.

3 王钢，陶家琪，徐兴伟等.东北电网500kV人工三相接地短路试验总结[J].电网技术,31（4）,2007.

4 张红斌，汤涌，张东霞等.不同负荷模型对东北电网送电能量的影响分析[J].电网技术, 31（4）,2007.

5 王守相，郑志杰.模型参数不确定性对电力系统时域仿真的影响[J].继电器,35（4）,2007.

6 宋方，毕天，杨奇逊.基于广域量测系统的电力系统多摆稳定性评估方法[J].中国电机工程学报,26（16）,2006.

7毛家安，郭志忠，张学松.一种基于广域量测系统过程量测数据的快速暂态稳定预估方法[J].中国电机工程学报,26（17）,2006.

8 罗建裕，王小英等.基于广域量测技术的电网实时动态监测系统应用[J].电力系统自动化,27（24）,2003.

9 范琦，穆钢，王克英.基于同步相量量测的线路参数在线量测的实验研究[J].东北电力学院校报,22（4）,2002.

10 刘洪波，穆钢等.根据量测轨迹计算轨迹灵敏度的卷积法[J].电力系统自动化,31（5）,2007.