多场耦合系统动力学仿真方法研究进展.pdf

2023力学进展版权所有ling systems.Advances in Mechanic3,53(2):468-495Wang G P,Tao L,Rong B,Rui X T.Review of dynamics simulation methods for multi-fieldcoup引用方式：王国平陶玲，戎保，芮筱亭，多场耦台力学进展，2 0 2 3,53(2）：46 8-49 5研究综述力展进学2023年6 月第53 卷第2 期多场耦合系统动力学仿真方法研究进展王国平1陶 玲 2,*戎 保1芮筱亭11南京理工大学发射动力学研究所，南京2 10 0 942南京理工大学能源与动力工程学院，南京2 10 0 94摘要现代工程系统往往是以复杂结构/机械系统为主体，融合热、流、电磁等若干子系统的多场耦合系统.此类系统动力学建模复杂、计算难度大，给系统动态特性高效精确评估与设计优化带来前所未有的挑战，有关其高效精确动力学仿真方法的研究愈发受到关注.本文详细回顾了复杂工程系统多场耦合动力学仿真方法研究成果和进展，包括：多场耦合动力学建模与数值求解基本策略、网格变形处理方法、耦合数据交换技术、数值计算效率等问题，在此基础上详细讨论了单一和混合不确定性条件下多场耦合系统不确定性分析及可靠性评估方法，以期为相关研究提供有益的借鉴和参考.关键词多场耦合，动网格，分区求解，不确定性，可靠性，耦合数据传递中图分类号：0 30 2文献标识码：AD0I:10.6052/1000-0992-22-043收稿日期：2 0 2 2-0 9-2 8；录用日期：2 0 2 2-10-2 1；在线出版日期：2 0 2 2-10-2 2*E-mail:469王国平，仿真方法研究进展陶玲第2 期1引言兵器、航空、航天、船舶、核能、交通、通用机械工业等国防和国民经济建设领域诸多现代工程系统往往是以复杂结构/机械系统为主体，融合热、流、电磁等若干子系统的多场耦合系统.此类系统动力响应和动态特性分析，一般涉及热、流、电磁、结构/多体等至少两个以上相互耦合的学科领域，动力学建模复杂，数值计算难度大、效率低；特别是系统存在不确定性因素时，不确定性和场间耦合效应将进一步导致建模求解难度与计算量急剧增加，给系统动态特性高效精确评估与多学科设计优化带来极大挑战（Felippa&Park2005，Ne r o n&La d e v e z e 2 0 10,Rongetal.2019，A r n o l d&H e c k m a n n 2 0 0 7）如何实现复杂多场耦合系统高效精确动力学建模与仿真愈发受到关注，已成为相关领域研究热点和难点.本文将从多场耦合系统动力学建模与数值求解基本策略、网格变形处理与场间耦合数据交换、数值计算效率以及不确定性分析与可靠性评估等方面回顾多场耦合系统动力学建模与仿真研究成果和最新进展，以期为相关研究提供借鉴和参考.2多场耦合系统动力学建模与数值求解基本策略多场耦合系统动力学建模和求解，首先根据工程需求简化实际物理模型，分析研究对象涉及的物理场及相互耦合关系、耦合变量等，而后基于各物理场基本定律和耦合关系建立相应的动力学方程，如：线弹性结构场的运动微分方程、物理方程和几何方程，瞬态温度场的热传导方程，流场的纳维-斯托克斯方程等。上述方程一般为偏微分方程，可基于有网格（如：有限差分、有限元、有限体积法等）（张文生2 0 0 6，Zienkiewiczetal.2013，高效伟等2 0 15，刘君等2 0 16，张雄和刘岩2 0 0 4,Zhang etal.2017）或无网格方法（如：光滑粒子法SPH、运动粒子半隐式法、无格子法、无网格Galerkin法等）（孙鹏等2 0 17，葛杨俊2 0 19）离散求解区域将其转化为常微分（代数）方程后选择适宜的数值方法求解.以各种方式铰接的多个大运动部件所组成的多体系统与热、流、电磁等物理场的耦合也是现代工程系统常见的耦合形式.耦合场作用下多体系统运动学和动力学描述复杂，需在适当参考框架（如：浮动标架法、旋转坐标法、绝对坐标框架等）和柔体变形描述（如：有限段、有限元、模态综合法、无网格法等）基础上，应用Witenburg方法、Kane方法、Lagrange方法、绝对节点坐标法、多体系统传递矩阵法等多体动力学方法建立既能描述部件运动/变形又可刻画部件间相互约束关系的多体系统总体动力学方程或总传递方程（芮筱亭等2 0 0 8，Rongetal.2019，A r n o l d&H e c k ma n n 2 0 0 7）.除上述基于明确物理规律及本构关系构建的多场耦合系统动力学模型外，对复杂问题还可结合试验、经验等建立经验关系式或近似数学模型.由于场间耦合作用，上述方程通常不仅包含学科自身的物理量，还包含其他学科的物理量，如图1所示，需结合初始和边界条件耦合求解，求解难度远大于单一物理场数值求解难度（Rong etal.2019，王思莹等2 0 14,Hron&Turek2006,Degrooteetal.2009，宋少云2 0 0 7,Dietz et al.2001,Fernandez 2011,Erbts et al.2015,Degr0ote&Vierendeels 2011).目前，集成解法和分区解法是多场耦合数学模型求解的两种主流方法（HronT u r e k 2 0 0 6,Degroote et al.2009,宋少云 2 0 0 7,Dietz et al.2001,Fernandez 2011,Erbts et al.2015,Degroote力470展学2023年第53 卷进系统总体动力学方程或总传递方程等纳维一斯托克顿方程等压力多体流场系统位移、转动、变形等温度流速等热传导方程位移、压力（若双向耦合）转变形等温度场温度变形电磁场几何、平结构场(若双向耦合）电场力、磁场力衡、物理方程等变形等(若双向耦合）麦克斯韦方程等图1多物理场典型耦合关系和耦合变量示意图（Rongetal.2019，王思莹等2 0 14，Hron&Turek2006,Degroote et al.2009,宋少云 2 0 0 7,Dietz et al.2001,Fernandez 2011,Erbts et al.2015,Deg-roote&Vierendeels 2011)&Vierendeels2011）.集成解法在同一时间步内将各场动力学方程集总同步求解，不存在时间滞后，适于强耦合、高非线性多场耦合动力学问题.但多场耦合系统通常具有不同的数值特性和时间尺度，集成解法往往需面对病态矩阵难题，大规模问题系统总体方程涉及矩阵阶次高、计算成本大（Dietzetal.2001，W a l l 2 0 12），且其一般需针对具体问题建立相应的计算模型、方程和仿真方法，难以充分利用各学科现有专用求解器和最新研究成果.分区解法在各时间步独立计算各物理场（域），各场（域）可使用独立的空间离散方法和数值积分策略，场（域）间交互通过数据传递进行，求解过程中各时间步内一般需迭代计算以满足指定精度要求，其便于利用已有代码并充分使用各学科已有研究成果.分区求解可分为显式和隐式两种方案（Felippa&Park2005,Fernandez2011），前者按序依次求解各物理场控制方程，并且在各时间步中只求解一次或固定次数，不强制场间耦合参量平衡，出于计算精度和稳定性考虑对时间步长有严格限制，为保证计算收敛一般选择较小步长，但过小步长会产生由附加质量效应导致的计算稳定性问题(Brummelen 2009）；后者在各求解时间步内利用Gauss-Seidel迭代（Joosten et al.2009）、Ne w-ton/Quasi-Newton迭代（Michleretal.2005）、加速松池法（高行山和林胜勇2 0 0 1）、定点迭代(Kittler&Wall 2008）等方法多次计算，以强制场间耦合参量平衡（或达到指定收敛值）、提高耦合分析精度.分区求解削弱了场间耦合作用，求解域划分方式和各场（域）的时间积分顺序都会对数值计算产生较大影响，且单场求解异步性可导致数值误差、收敛困难或计算失稳等问题宋少云2 0 0 7），与整体方法相比无法保证一定具有更高的计算效率.以流场-结构场耦合问题为471王国平陶玲戎保，芮筱亭第2 期多场耦合系统动力学仿真方法研究进展例，其分区解法典型计算流程如图2 所示.同样对多场耦合多体系统分区解法，针对所研究问题的特点也出现了各具特色的计算方法，如：Rumold（2 0 0 1）针对充液多体系统动力学问题，将系统分割成多刚体子系统和流体子系统，多刚体子系统采用多体动力学方法求解，流体子系统采用基于有限体积法的计算流体力学求解，在每一个时间步将求解所需的多刚体子系统和流体子系统间耦合变量值近似为前一时刻值以实现耦合变量数据传递和两个子系统的独立求解，如图3所示（Rumold2001）与流场-结构场耦合不同，多刚体子系统和流体子系统间耦合变量包括充液部件连体坐标系与惯性系的坐标转换矩阵SF、角速度WF、角加速度WF、坐标原点相对惯性系的位置矢量rF及其加速度F、作用在充液部件上的液体力FMF和扭矩MMF等.若考虑充液部件柔性变形，则变形甚至流场网格变形/重构信息也需视情在各时间步子系统间传递.3网格变形处理与场间耦合信息交换技术3.1网格变形处理方法对大运动、大变形或时变不连续多场耦合问题采用有限元等有网格方法进行计算区域离散时，边界运动或计算域变形引起的网格变形/畸变严重影响耦合计算精度和数值稳定性（周璇等2011，张来平等2 0 10）.网格变形（如：虚拟结构法、偏微分方程法和代数法）（周璇等2 0 11）、网格重构（孙鹏等2 0 17，葛杨俊2 0 19,Yangetal.2011）及二者结合方法等动网格技术是解决上述问题的有效手段.虚拟结构法包括弹簧近似法（Degand&Farhat2002，M a r k o u e t a l.2 0 0 7）和弹性体法（Abgralletal.2014，仲继泽等2 0 16，张斌等2 0 18），其假设网格节点由弹簧或弹性介质连接，求解弹性问题获得网格变形，前者简单、占用存储空间小、适于小变形问题；后者变形能力优于前者，适于大运动/变形问题，但效率低、实现复杂，存在运动边界网格畸变甚至“负体积”非法单元问题偏微分方程法由位势方程求解网格位移，包括Laplace方程（Ogiso etal.2010）、双调和算子（Helenbrook2003）、四元代数（Samareh2002）等方法，适于小规模网格问题.代数法包括径向基函数插值（Estruchetal2013，孙岩等2 0 18）、Delaunay背景网格插值（范锐军和周洲2 0 10，Wangetal.2018）、距离函数插值（Zhao&Tai2002）和无限插值（石磊等2013）等，高效直接，仅根据节点位置插值得到其运动规律，无需考虑其邻接关系及网格类型，但存在径向基函数选取对变形后网格质量影响大（周璇等2 0 11,Estruchetal2013，孙岩等2 0 18）、Delaunay背景网格插值仅对凸域网格变形有效、可能网格堆积（范锐军和周洲2 0 10，Wang etal.2018）等问题.与网格变形法不改变原网格节点数和拓扑不同，网格重构法重新生成网格并将重构前相关变量插值映射到当前网格，存在网格重新生成计算量大和额外插值误差问题.提高网格变形能力及质量、减少网格重构并提高插值精度是动网格关注的焦点（周璇等2 0 11，张来平等2 0 10，Yangetal.2011）.周璇等（2 0 11）、张来平等（2 0 10）对主要动网格方法及其特点给予了详细比较，感兴趣的读者可参阅相关文献.为提高动网格生成效率和质量，上述基本方法结合或改进研究也层出不穷.孙岩等（2 0 13）、蒋波（2 0 17）研究了径向基函数和Delaunay背景图相结合的网格变形方法，提高了变形后网格的光顺性和质量（相关算法流程如图4所示），徐琳和宋万强（2 0 19）研究了基于Laplace网格变形算法与局部重构相结合的动网格技术，在单步小位力472展进学20233年第53 卷流场结构场FO流体域初始开始SO固体域初始条件和边界条件条件和边界条件F1流体域离散S1固体域离散;C.4)网格C.4)P网格变形i=0变形/重构Coupling/重构网格变F2离散格式的S2离散格式的形处理动力学方程动力学方程;C.5)接收与场间C.2)接收耦节点位移等耦合信合面压力F3流体域动力学息交换S3固体域动力方程送代求解响应分析C.3)发送结构变形等否判断耦合选代收敛C.1)发送耦合面压力i=i+1是否判断求解时是结束间是否结束图2典型流场-结构场耦合问题分区解法示意图F1N-S方程方程离散与迭代求解FO流体aFur+o(Fur:FV)Fur+2FWF X Fur-Fu r+FV p域初始at和边界条件=QSTIf-QFaF,(to)FVFur=O,流场couplingF（t o)couplingMcouplingrcouplingMt1时间积分=0,1,2,N0=22耦合信息耦合信息11(to)couplingMocouplingM(t1)111(to)couplingrcouplingF(t1)11couplingF(to)IcouplingF(t)11(to)couplingMcouplingl时间积分i=0,1,2，.NMO多M1多刚体子系统动力学方程求解体子系M(y,t)j+k(y,y,t)=q(y,y,t)统初始和边界条件多刚体子系统coupling.l=(rf),coupling(Fle,Mlih)(t)(t)(t)tirFWFWF(t)S(ti)SF图3典型充液多刚体系统动力学问题分区解法示意图(Rumold 2001)移下基于Laplace算法网格变形，当总位移量增大至网格质量低于预期时局部网格重构（仿真过程网格变形/重构示意图如图5所示）.动网格结合任意拉格朗日-欧拉（ArbitraryLagrangian-473王国平，陶玲仿真方法研究进展第2 期Bowyer-RBFWatson算法插值计算物面点初始变形后DClaunay加密后网格背景网格的远场边界点定位计算网格图4径向基函数和Delaunay背景图相结合的网格变形方法（蒋波2 0 17)上翼型（飞机）背景网格下翼型（导弹）outerbpx运动物体子网格图5双翼型轨迹模拟过程网格变形/重构示意图.(a)背景网格与运动物体子网格，(b)初始位置网格，(c)网格变形，(d)大变形的网格重构（徐琳和宋万强2 0 19)Eulerian，A LE）方法求解流-结构/多体等多场耦合问题愈发受到关注，如：充液刚体动力学（周宏等2 0 0 8）、飞机油箱液体晃动模拟（杨瑞2 0 15）、伞衣下落过程模拟（高兴龙2 0 16）等.ALE网格可承受更大畸变且解析度更精细，但ALE方程含网格与物质点间迁移速度引起的迁移项，计算复杂.除上述方法外，基于边界元、浸入边界法、无网格法等方法（Erbtsetal.2015,Doneaetal.2004，孙鹏等2 0 17，葛杨俊2 0 19）解决边界变形问题的研究也不断涌现.3.2场间耦合信息交换技术场间耦合信息高效精确传递是多场耦合分区求解必须面对的问题，其对计算稳定性、精度和效率影响大，确保能量守恒、荷载信息等效、基本未知量及其导数连续是保证场间信息准确传递的关键（刘深深2 0 14,Boeretal.2007）.根据样本空间范围，场间信息传递算法可分为全局、局部插值方法及上述两类方法的混合方法等.全局方法基于全部已知点信息拟合待插值点对应值（刘深深2 0 14,Boeretal.2007），如：样条函数法（HarderD e s ma r a i s 197 2,D u c h o n 197 7）、多面函数双调和法（Hardy1971，K a n s a 1990）、逆等参映射（Pidaparti1992）、Sh e p a r d 方法(Shepard 1968）、径向基函数法（BeckertW e n d l a n d 2 0 0 1，Smi t h e t a l.2 0 0 0）等.Harder和Desmarais（197 2）基于无限平板样条函数法（infinite plate spline，IPS）计算耦合场间平面法向位移和力的传递问题.Duchon（197 6）基于薄板样条函数法（thinplate spline，T PS）处理不规则曲面插值问题.IPS和TPS需满足薄板变形理论及板内位移连续且四阶可导假设，对复杂外形适力474展学2023年第53 卷进用性差.Hardy（197 1）、K a n s a（1990）研究了对各种曲面外形适用性较好的多面函数双调和法（M u l t i q u a d r i c-Bi h a r mo n i c，M Q）和逆多面函数双调和法.Pidaparti（1992）提出气动弹性耦合界面四边形单元插值的等参逆变换法.崔鹏和韩景龙（2 0 0 9）提出局部TPS/IPS插值方法，以提高计算效率和稳定性.Shepard方法按待插值点与各数据点距离倒数加权平均，对缓慢变化量插值精度低，引入各数据点高阶导数可提高精度.径向基函数插值包括全域和紧支径向基函数两类(Beckert&Wendland2001，Smi t h e t a l.2 0 0 0），可用于结构化、非结构化及混合网格，形式简单，前者采用全部已知数据点，计算矩阵通常满阵，精度高，大规模插值计算量大；后者仅采用部分数据点，计算矩阵稀疏、带状分布，效率高，但插值误差通常较前者大.吴宗敏（2 0 0 7）推演TPS,IPS和MQ方法获得统一的径向基函数插值格式.刘深深（2 0 14）提出基于几何外形缩比的改进径向基函数法，大幅减小了紧支C基函数的插值误差.尹果（2 0 14)研究了流固耦合数据传递问题径向基函数插值误差.局部方法基于待插值点周围已知点信息计算其对应值（刘深深2014，Bo e r e t a l.2 0 0 7），如：最邻近点插值法及其改进（Thvenzaetal.2009，汪学锋等2 0 0 9，杨敏等2 0 11）、反距离权重法（张福才2 0 15）、映射点插值法（DettmerPe r i c 2 0 0 6）、加权余量法（Steinetal.2000）、常体积法（Gouraetal.2001，程兴华2 0 12，徐敏和陈士橹2 0 0 4，霍霖2017）及局部形式表面拟合法（Cavagna etal.2007）等.最邻近点插值法（Thvenza et al.2009,汪学锋等2 0 0 9，杨敏等2 0 11）从已知网格的最近点获得相关插值信息，信息传递简单快捷，但两网格匹配较差时插值精度低.为降低最近邻点法网格匹配度限制，汪学锋等（2 0 0 9）、杨敏等(2011）研究了基于待插值点邻近若干已知点信息的改进邻近点加权平均法.反距离权重法（张福才2 0 15）按照数据点与待插值点距离大小与贡献直接插值计算，形式简单，不涉及矩阵求逆，适用性好.映射点插值法在耦合边界将插值域（发送方）中节点向被插值域（接收方）界面单元所在平面正交投影获得对应的映射点和主单元，基于映射点局部坐标和形函数进行插值，但外形复杂、网格质量差的局部区域插值精度低（DetmerPe r i c 2 0 0 6，蒋波2 0 17）.Stein等（2 0 0 0)发展了流固耦合问题的加权余量法，流体区域和固体区域形函数相乘的积分求解是其面临的难题.常体积转换法（constantvolumetransformation，CVT）（G o u r a e t a l.2 0 0 1）利用气动点及其正交投影所在结构三角形单元组成的四面体体积守恒，计算气动点和结构单元间的插值关系，广泛用于气-固耦合界面信息传递.为提高CVT算法稳定性，程兴华（2 0 12）提出内投影常体积转换法.Xu和Chen（2 0 0 4）提出引入面积限制的改进CVT方法，以提高插值精度.霍霖（2 0 17）提出基于非结构网格面元的内投影CVT法，用于高超声速飞行器体襟翼气动热弹性耦合分析.Cavagna等（2 0 0 7）研究了基于移动最小二乘近似的局部形式表面拟合方法，其插值形函数不具备狄拉克函数性质、所得曲面难以精确包含结构节点.全局方法本质是对散乱数据的无网格插值（拟合），更易处理耦合面网格不匹配问题，但相比于局部方法计算时间较长，且当网格规模大、形状复杂时，其系数矩阵条件数大影响插值精度；局部方法应用方便，但精度较低，整体位移协调性较差（刘深深2 0 14,Boeretal.2007）.综合各种耦合信息交换算法的优点，众多学者提出各种混合算法以改进数据交换的精度和效率.Kim和Kim（2 0 0 5）提出基于部分表面跟踪和部分表面拟合的气动弹性问题非结构网格混合插值方475王国平陶玲戎保，芮筱亭多场耦合系统动力学仿真方法研究进展第2 期法.宋少云和李世其（2 0 0 6）提出了自然邻居插值、滑动最小二乘及反距离移动平均法相结合的节点载荷混合插值算法.韩旭里等（2 0 0 7）结合B样条与径向基函数插值，研究了计算稳定、适于非矩形网格的散乱数据曲面拟合方法.An等（2 0 0 9）提出基于边界元法的网格插值方法.刘刚等（2 0 18）结合快速映射法和径向基函数法，提出了二维电磁场-流场-温度场耦合问题异构网格节点数据映射的混合算法，异常区域节点使用快速映射法，其他节点使用径向基函数法，提高了异构网格间数据映射的精度.Huang等（2 0 18）基于局部近邻搜索算法匹配流体节点和结构单元，基于薄板样条处理非匹配网格间数据传递，实现了高速列车单向流固耦合信息传递.目前尚无任一种方法能适用于所有耦合场间信息交换，如何根据具体问题选择适应的插值方法是值得研究的问题.4多场耦合系统动力学计算效率问题多场耦合问题计算效率是计算力学领域长期关注的焦点问题，众多学者分别从耦合问题建模方法、数值积分与迭代算法改进、并行计算等不同角度提出了一系列改进策略，极大促进了多场耦合数值计算的发展.在耦合问题建模方法方面，各学科自身高效建模理论和数值方法的发展是提高耦合问题计算效率的基础，如：多体系统传递矩阵法、递归建模方法等一系列新型建模方法（Rongetal.2019，Ru i e t a l.2 0 2 2），大幅度提高了多场耦合多体动力学的计算效率.芮雪等(2020）基于多体系统传递矩阵法和VanderPol尾流振子模型建立了海洋热塑性增强管的流固耦合动力学方程，实现了海洋热塑性增强管振动特性和涡激振动响应的快速预测（如图6 所示）：李书田（2 0 2 0）基于多体系统传递矩阵法建立了多管火箭燃气射流-多体动力学单向耦合模型，提高了计算速度快，为分析多管火箭发射动力学和发射装置设计提供理论支撑（如图7 所示），模型降阶和模型自适应的应用也有效平衡了多场耦合问题计算效率和精度.模型降阶在工程精度允许范围内寻找近似原高阶模型行为的低阶模型，以实现对复杂系统的快速计算.以多场耦合多体动力学模型降阶为例，诸多学者发展了诸如模态综合法、本征正交分解、Krylov子空间模型降阶法等一系列方法，有效提高了该类问题的计算效率（Heresetal.2007，I h r l e e t a l.2 0 13,Schilders&Lutowska2014，Su n e t a l.2 0 15）必须指出的是，多场耦合系统模型降阶须有效处理降阶过程中学科间耦合的问题，难度较单学科模型降阶大得多.Heres等（2 0 0 7）提出Krylov子空间方法和正交向量拟合相结合的模型降阶策略，用于平面双LC滤波器的准静态电磁分析.Ihrle等（2 0 13）将鼓膜和耳道以及鼓室中的空气视为弹性体，建立了中耳结构非线性弹性多体系统动力学模型，并应用基于Petrov-Galerkin投影的约简技术实现了声-结构耦合系统的模型降解.Schilders&Lutowska（2 0 14）在子系统降阶和耦合块低秩近似基础上提出耦合时不变线性模型的降阶技术.Erbts等（2 0 15）提出可用于强耦合系统、适用性广、计算成本低的分离基降阶算法.Sun等（2 0 15）基于多体动力学理论和模态降阶方法，建立了地铁轮-轨耦合振动分析模型.Yamashita等（2 0 19）采用Craig-Bampton部件模态综合方法降阶热力耦合模型，用于约束多体系统的瞬态热力分析，降低了计算时间.Tian等（2 0 2 2）将热力耦合系统划分为多个子结构，建立了基于自由界面部件模态综合法的热力耦合柔性多体动力学模型降阶方法，提高了仿真效率Sullivan等（2 0 2 2）利用本征正交分解（properorthogonaldecomposition，PO D）开发了轮胎-土力476展进学20233年第53 卷=0.1m/s,均匀来流ab0.2 m/s,均匀来流顶张力1.00.3m/s,均匀来流0.1m/s,剪切流剪切流0.80.2 m/s,剪切流20.3m/s,剪切流工刚性接头0.6UL0.40.2均匀来流00.10.20.30.40.50.6ly|max/D图6基于多体系统传递矩阵法和VanderPol的海洋热塑性增强管流固耦合仿真.（a）立管模型示意图，(b）剪切流情况下立管轴向涡激振动振幅分布（芮雪等2 0 2 0)abpressure/PaC0.6X10-2.第一发定向器导流环29203.4377.84651030.4第九发t=0.0125s第十八发0.200.50.50.50.2-0.4000-0.6发射管-0.8火箭弹-0.10-0.50-0.50-0.5000.20.40.60.81.0时间/spressure/Pa1720.504386129215103第一发(t-s.e)/回0.4第九发t=0.0410s0.3第十八发0.20.50.50.50.100.1000-0.20.3-0.4-0.500.500.5000.20.40.60.81.0时间/s图7基于多体系统传递矩阵法的多管火箭燃气射流-多体动力学单向耦合仿真.（a）多管火箭发射系统动力学模型，(b）不同时刻迎气面压强，(c)定向管口y方向位移、角速度（李书田2 0 2 0)壤耦合模型的降阶程序，并用于车辆越野机动仿真.模型自适应是指在多场耦合计算中利用移动网格、局部加密或合并网格等手段不断调节网格，在变化剧烈区域（或仿真阶段）采用精细模型(网格）提高仿真精度，在变化平缓区域（或仿真阶段）采用粗模型（网格）确保仿真效率（Rongetal.2019，Ih l e n b u r g 2 0 0 0）.Ih l e n b u r g（2 0 0 0）联合有限元-边界元法求解流体-结构耦合问题，并研究了模型自适应和网格细化技术对数值稳定性的影响.Ryan（2 0 0 6）研究了多物理和多尺度问题的网格自适应技术.Faigle等（2 0 14）基于网格自适应技术和序贯求解方法模拟了多孔介质两相流问题.在数值积分与送代算法改进方面，传统分区求解显式算法效率高，但单场求解异步性导致场477王国平，陶玲学仿真方法研究进展成保第2 期间耦合参量不平衡，从而带来计算误差，稳定性较差；隐式算法计算稳定性较好，但迭代求解计算量大.Zienkiewicz等（198 8）研究了孔隙-流体-土壤耦合问题分区求解算法的稳定性.Deg-roote等（2 0 0 8）利用傅立叶误差分析得到了管内非定常流固耦合问题分区求解数值稳定性条件.Joosten等（2 0 0 9）采用块Gauss-Seidel方法交错求解流固耦合问题，研究了不同时间积分器对收敛性的影响.Zhang等（2 0 0 3）研究了差分格式对稠密气固两相流动模拟精度和收敛性的影响.Fernandez等（2 0 0 7）提出半隐式分区迭代算法用于流固耦合问题.Astorino等（2 0 0 9）提出基于Robin条件的半隐式耦合算法，并研究了其收敛性.Busch和Schweizer（2 0 12）提出仅利用耦合变量Jacobian信息的半隐式耦合求解技术，并用于高速涡轮机非线性转子/轴承混合模型多场耦合分析.何涛（2 0 13）提出了基于特征分裂算法的全新半隐式分区耦合方法，提高了数值稳定性和精度.半隐式耦合方法计算稳定性与隐式算法相当，但效率更高（2 0 13）为解决迭代求解收敛缓慢问题，Aitken动态松弛技术（Irons&Tuck1969）、基于Robin边界条件的加速耦合迭代(Badia etal.2009）等改进迭代方法不断出现.协同仿真、自适应时间步长、多速率积分等研究越发受到重视协同仿真允许不同物理场（子系统）选择各自适用的积分方法、时间步长和模型离散方法，以提高计算精度、稳定性和效率.自适应时间积分可根据计算条件自动调整时间步长，当系统在相对较小的时间尺度上经历较大变化时自动降低步长，在系统变化相对较慢时增大步长，从而在保证计算精度/稳定性的同时尽可能提高计算效率.实际多场耦合系统通常包含分离时间尺度，多速率积分为此引入了不同时间步长，以避免系统最小积分步长限制.Niesner等(2006）研究了流固耦合问题自适应时间积分和平衡选代方法.赵振（2 0 11）研究了自适应步长的电磁-机械顺序耦合分析技术.Zohdi（2 0 17）研究了光学-热-机械多场耦合问题递归、交错、时间自适应求解方案.Rong等（2 0 13）将多速率方法应用于弹炮耦合武器系统动力学分析（如图8所示）.Lunk和Simeon（2 0 0 6)提出电弓一接触网动力学分析的网格-时间步长自适应算法.Mas-sat等（2 0 14)研究了列车-电弓-接触网耦合动力学的有限元-多体协同计算问题.随着多核计算、分布式计算、云计算的发展，利用多处理器或分布式系统完成多场耦合问题并行求解是计算力学发展趋势之一（吴帆2 0 16）.由于不同物理场的物理模型、离散网格、区域分解方式等不尽相同，加之场间耦合作用复杂，多场耦合问题高性能并行计算环境和并行算法的构建远比单物理场并行计算复杂，是目前急需解决的瓶颈性难题.Kittler和Wall（2 0 0 6）提出不可压缩牛顿流体-结构耦合问题的并行有限元方法.Vardon等（2 0 0 9）提出非饱和土热-水-力耦合MPI/OpenMP混合并行有限元模拟.Wang和Kolditz（2 0 0 9）提出便于线性方程组划分和并行求解的稀疏矩阵和迭代线性求解器，用于多场耦合问题大规模并行顺序有限元分析.张兵(2011）采用共享内存技术开发多场耦合计算平台，实现了流固耦合分区传热计算.Zhang等(2013）研究了车辆-板式轨道耦合系统的并行非平稳随机振动分析.Negrut等（2 0 12）研究了流固耦合问题和摩擦/接触多体动力学问题的并行模拟.Dong等（2 0 16）在OPENFOAM和RE-MARK代码基础上实现了核电站中子分析与热工水力耦合的高效并行计算.游检卫（2 0 16）提出多物理场MPI-OpenMP-GPU（M O G）混合并行求解技术.Bungartz等（2 0 16）研究了多物理场双向表面耦合问题的并行库preCICE.Hu等（2 0 18）实现了可压缩流体-多体系统耦合动力学GPU并行模拟.力478展20233年第5 3 卷进学a俯仰部分41弹丸子系统-IIbC炮尾4342.640固接铰44身链身管星周接铰46 炮口制退器47C0481.0仿真车体子系统、0.5一试验0.3仿真38一试验0回转部分车体37弹丸子系统0.60.530车体子系统-I0.94-1.0400.40.81.2 1.60123Ati-1车轮13 2 4t/st/sa后座位移b炮口位移00图8弹炮耦合武器系统多速率协同仿真.(a）弹丸-火炮多体系统耦合模型，(b）多速率积分方法，(c）后座位移和炮口位移时间历程(Rong et al.2013)5不确定性多场耦合系统分析方法由于制造水平、测量误差、建模精度等原因，实际多场耦合工程问题动力学模型求解所需的初始/边界条件、几何/材料参数、外部载荷甚至计算模型本身难以避免存在不确定性（如：随机不确定性、认知不确定性、模型误差）（Rong etal.2019，T h u n n i s s e n 2 0 0 3），随机不确定性将不确定性因素视为随机变量（场或过程），利用概率方法建立随机模型处理.认知不确定性由主观认识不足或信息贫乏所致，难以精确概率分布描述，可采用区间/凸集理论、可能性理论、模糊集理论、证据理论等表述（Ben-Haim&Elishakoff 1990，T z v i e l i 1990，Z a d e h 196 5，Sh a f e r1976）.模型误差包括建模假设、近似模型及数值计算等引起的误差，其主要包含认知不确定性.确定性模型无法反映上述不确定性因素对系统性能的影响，特别是对多场耦合系统，由于学科间耦合作用，不确定性因素的小幅变异亦可能使系统响应产生较大偏差，导致系统实际性能偏离设计指标较远，甚至造成关键部件失效.多场耦合系统不确定性建模与分析是复杂的，系统响应约束）通常是高维、高次非线性函数，系统不确定性变量数量、非线性程度等因素均会对计算精度和效率带来不可忽视的影响.如何建立可行的多场耦合系统不确定性数学模型，并在此基础上展开高效精确不确定分析，科学评估不确定性对系统动态特性（响应、性能、可靠性等）的影响，是多场耦合系统分析与优化不得不关注的难题之一多场耦合系统不确定性分析，按不确定性描述方法可分为概率、模糊、凸集、区间以及上述方法的混合方法等；按实现途径可分为非侵入式和侵入式；按实际目的可分为系统响应的不确定量化、可靠性分析等.本节将对多场耦合系统不确定性分析方法给予详细总结，特别对可靠性分析现状给予进一步概述.5.1单一不确定性分析对单一不确定性问题，当系统不确定性可用精确概率分布描述时，基于抽样方法、近似积分法、展开/随机元模型方法以及分解协调法等方法的概率不确定性分析是不确定性量化的首选.采样方法通过反复采样量化参数或模型的不确定性，对小概率或多场耦合问题计算量大，采用改进抽样策略，如：拉丁超立方采样（Rajabietal.2015)、重要性采样（Angetal.1992)、分层采样(Giles2013)、自适应采样（Shayanfaretal.2017）等，减少抽样次数和基于代理模型抽样以有效降低单次抽样计算量是克服其计算量大缺陷的两种主要手段.近似积分法应用高斯积分等数学479王国平，陶玲学仿真方法研究进展第2 期工具简化多维积分（Leeetal.2009），低维问题一般精度和效率较高，但高维、强耦合、高非线性问题计算精度（特别是高阶矩的精度）和效率面临严峻挑战.为减轻维数灾难问题，诸如单变量降维（Rahman&Xu2004)、稀疏网格（Gerstner&Griebel1998)、广义降维（Xu&Rahman2004）等非全阶数值积分受到青睐.展开/随机元模型方法主要包括Taylor展开、摄动法、随机配置法、多项式响应面法、多项式混沌展开、径向基函数、Kriging模型、支持向量回归、神经网络等，且绝大多数有侵入式与非侵入式两种形式（Martowicz&Uhl2012，Ji a n g e t a l.2 0 18,Man et al.2018,Kieu et al.2019,Chaudhuri et al.2018,Nannapaneni&Mahadevan 2020,唐新姿等2 0 2 0，Wangetal.2019).非侵入式方法将系统视为黑箱，利用试验设计选择模型输入计算相应输出，并基于不同方法建立原模型的近似模型进行不确定性分析，方便直接；侵入式方法对原问题中不确定性参数近似展开，进而推导系统输出，需设计新求解器，但与非侵入式相比一般计算成本较低.Chaudhuri等（2 0 18）研究了基于自适应代理模型的多学科耦合系统多保真不确定性分析.Nannapa