多体系统动力学时程积分算法研究进展.pdf

20239Building Construction1922多体系统动力学时程积分算法研究进展林汛达 李忠学浙江大学建筑工程学院 浙江 杭州 310058摘要：发展精确且稳定的多体系统动力学计算方法一直是研究者们努力的方向。时程积分法已经有数十年的发展历程，这种算法的核心思想是将结构模型的运动过程在时间域上离散，在每个时间步内进行积分，逐步计算出结构的完整运动轨迹。其中，能量守恒时程积分法有精度高，稳定性好等优点，具有较大的应用潜力。通过对时程积分法相关文献的分析与总结，介绍了该算法的基本原理和发展现状。最后对该算法的后续研究提出展望。关键词：多体系统动力学；非线性分析；时程积分法；数值模拟中图分类号：TU13 文献标志码：A 文章编号：1004-1001(2023)09-1922-04 DOI：10.14144/ki.jzsg.2023.09.053Progress in Time-integration Algorithms for Multi-body System DynamicsLIN Xunda LI ZhongxueCollege of Civil Engineering and Architecture,Zhejiang University,Hangzhou 310058,Zhejiang,ChinaAbstract:The development of accurate and stable multi-body system dynamics calculation methods has been the direction of researchers efforts.The time integration method has been developed for decades.The core idea of this algorithm is to discretize the motion process of the structural model in the time domain,integrate in each time step,and gradually calculate the complete motion trajectory of the structure.Among them,the energy conservation time integration method has the advantages of high accuracy and good stability,and has great application potential.Based on the analysis and summary of the relevant literature on time integration method,the basic principle and development status of the algorithm are introduced.Finally,the future research of the algorithm is prospected.Keywords:multi-body system dynamics;nonlinear analysis;time-integration method;numerical simulation组成的多体系统运动规律1。如今，诸如航空航天、机械制造、建筑设计等工程领域，都在开展多体系统动力学的研究。在建筑领域，多体系统动力学常被用于计算大跨桥梁、悬索桥梁、高层建筑的动力响应，在抗震设计方面也有较大应用潜力。许多基于多体系统动力学的仿真模拟软件不断出现，例如美国的ADAMS、韩国的RecurDyn和德国的Simpack等。通过仿真模拟可以减少大型实验需耗费的资源，为工程设计提供可行性验证，提高设计效率，并为工程施工提供参考。而发展出精确且稳定的多体系统动力学计算理论则可以为数值模拟提供理论基础。为了计算多体系统的动力响应，许多学者都提出了各种不同的求解思路，例如Bellman等2提出了微分求积法（differential quadrature method），钟宏志等3将变分原理和微分求积法结合，提出了弱形式求积元法（weak-form-quadrature element method），钟万勰等4提出了微分代数方程的保辛方法。而本文将主要介绍时程积分法（time integration method）的基本原理和发展现状。2 时程积分法概述时程积分法的基本思想是将动力方程在时间域0,T上多体系统是指多个物体（包括刚体、弹性体、质点）相互联结形成的复杂结构，多体系统动力学（multi-body system dynamics,MSD）研究的是多体系统受动力荷载作用下产生的力学响应。发展出一套原理清晰、计算精度高、稳定性较好的多体系统动力学计算方法一直是学者们研究的方向。本文首先讨论了动力学计算理论的研究意义，随后简述了时程积分法的主要特点，第3节和第4节分别介绍了线性时程积分法和非线性时程积分法的发展现状，最后对时程积分法的发展进行了展望。1 多体系统动力计算理论的研究意义多体系统动力学是近代一般力学与基础力学的一个分支，它是古典的刚体力学、分析力学和计算机等学科融合而成的一门新兴学科，重点研究相互连接的刚体或柔性体作者简介：林汛达（1998），男，在读硕士。通信地址：浙江省杭州市浙江大学紫金港校区安中大楼A712（310058）。电子邮箱：L收稿日期：2023-04-19科学研究SCIENTIFIC RESEARCH建筑施工第45卷第9期1923离散成N个微小时间步t，假设系统的本构关系在离散的时间步内保持线性，在单个时间步tn,tn1内根据初始条件进行积分，逐步求出一系列离散时刻上的动力响应。根据对tn1时刻结构构型的不同描述方式，时程积分法又被划分为显式积分法（explicit integration method）和隐式积分法（implicit integration method）两类。其中，显式积分法认为结构在tn1时刻的构型仅与tn时刻的状态有关，于是利用tn1和tn时刻的物理量直接推导出 tn1时刻的物理量值，其中的代表算法为中心差分法（central difference method）。需要说明的是，几乎所有的显式方法都存在稳定条件。以中心差分法为例，它的稳定条件是：tTnGrD（1）其中Tn为结构的自振周期，对于多自由度体系则为结构的最小自振周期。只有当离散的时间步长满足稳定条件时，数值计算得到的结果才能保持稳定，否则计算结果会逐步发散至无穷大，最终发生失稳。当结构的自由度较多，显式积分法容许的时间步长往往比结构最大自振周期低数个量级，求解效率很低。隐式积分法认为tn1时刻的构型会同时受到tn和tn1时刻结构状态的影响，于是将tn1时刻的某些物理量作为未知数代入微分方程中进行求解，其中的代表算法有Newmark-法5，它的稳定条件为：tT2121nGrcbD-（2）当1/2和1/4时，Newmark-法的稳定条件为t，即为无条件稳定（unconditional stability）。与显式积分法相比，隐式积分法所容许的时间步长更大，计算效率更高。3 线性时程积分法的发展假设结构模型在动力荷载的作用下发生运动，运动过程中任意时刻系统的构型参数将满足下列运动方程：MuCuKuf+=po（3）其中，M、C、K分别代表系统的质量矩阵、阻尼矩阵以及切线刚度矩阵，f 代表系统所受外荷载，u代表节点位移，上标圆点代表对时间的导数（有几个圆点，即表示为几阶）。通过运动方程求解出结构在运动过程中任意时刻的节点位移，即可得到结构的动力响应。1959年，Newmark5提出了著名的Newmark-法，此算法是第一个被广泛用于求解多体系统动力问题的计算方法，在整个学术界都有着巨大的影响力。Newmark-法假设tn1时刻节点速度和加速度的表达式具体为：tuuuu1nnnn11ccD=+-+oopph6（4）（5）ttuuuuu11211nnnnn121bbbDD=-+popahk将式（4）和（5）代入tn1时刻的运动方程，此时的方程只有唯一的未知量un1，求解运动方程即可得到tn1时刻的节点位移。1968年，Wilson提出了Wilson-法6，该算法假定节点变量的加速度在tn,tnt时间段内保持线性变化：tuuuunnnixD=+-i+pppph （6）其中0,t，1。在时间段内将u对求一阶和二阶积分即可得到tnt时刻的节点速度un和位移un+的表达式：tuuuuu2nnnnn2xixD=+-ii+ooppph（7）6tuuuuuu21nnnnnn23xxixD=+-ii+oppph（8）联立上式，将加速度un和速度un+均由un+表示，再将它们带入tnt时刻的运动方程，即可解得un+，令t，将取适当的值并代入式（6）式（8），即可得到tn1的构型参数。经验证，当1.37时，Wilson-法可以实现无条件稳定。1993年，为了在时程积分过程中实现数值耗散以排除结构高频能量震荡的不利影响，Chung和Hulbert7提出了广义法（generalized-method）。这种方法用时间步首尾时刻的构型参数求插值来计算广义中点时刻的各项参数：uuu1nfnfn11faa=-+a+-+h（9）uuu1nfnfn11faa=-+a+-+oooh（10）fff11nfnfn1faa=-+a+-+h（11）uuu1nmnmn11maa=-+a+-+ppph（12）其中tn1时刻的速度和加速度依然选择Newmark-法的基本假定，见式（4）、式（5）。将广义中点时刻的各参数代入广义中点时刻的运动方程，再使用增量迭代算法求解，即可得到tn1时刻的节点位移un1。相比于其他具有数值耗散的时程积分法，广义法在排除结构高频响应的同时，还能最大程度地降低对结构低频响应的干扰。同时，通过f 和m的不同取值，广义法可以转化为其他耗散方法，例如当m0时，公式变成HHT-法8；当f 0时，公式变成WBZ-法9；当m和f 都为0时，公式则回归到Newmark-法。为了方便调整耗散系数，Chung引入了谱半径用于计算f 和m：林汛达、李忠学:多体动力学时程积分算法研究进展20239Building Construction1924 121matt=+-33（13）1fatt=+33（14）其中0,1，当1时，高频耗散不会被激活；而当0时，耗散程度达到最大，高频响应会在一个时间步内被完全排除。1994年，钟万勰10提出了精细时程积分法，此方法首先将系统的运动方程转化成一阶常微分方程组：uMuCuM CCM CKMCMuMuCuf2242201111+=-+-opoRTSSSSVXWWWWHHF（15）先利用精细积分法计算出齐次方程的通解，再利用叠加原理找出非齐次方程组的特解，即可求出结构的动力 响应。精细时程积分法能够得到趋近于精确解的数值结果，但是在计算非齐次方程组的特解时需要进行矩阵求逆，不仅计算过程比较烦琐，还可能遇到逆矩阵不存在的问题，令该算法的应用范围受到了限制。2000年，顾元宪等11提出了增维精细积分法，在精细时程积分法的基础上，将式（15）中的非齐次项看作方程的状态变量，纳入求解对象中，将非齐次方程转换为齐次方程，在进行精细积分的过程中避免了矩阵求逆，拓宽了精细积分法的应用范围。2004年，汪梦甫等12将精细时程积分法与高斯积分法相结合，用高斯积分来计算微分方程中的特解，同样避免了矩阵求逆。此算法的计算精度取决于所选高斯积分点的数量，如果选取的积分点足够多，理论上可以达到任意高的精度。4 非线性能量守恒时程积分法的发展为求解几何非线性动力问题，最早的方法是将原本适用于计算线性问题的算法直接应用在非线性问题上，例如Newmark-法5、HHT-法8和WBZ-9法等。上述算法在求解线性问题时可以保证无条件稳定且拥有极高的计算精度，但当它们被应用于非线性问题时，会产生严重的能量波动现象。1975年，Belytschko等13提出利用能量准则来衡量时间积分公式的稳定性。该准则要求自由运动阶段系统的总能量不会发生改变。此后，这种能量守恒判定准则得到了学界的广泛认可，许多基于能量守恒的时程积分算法不断被提出。1978年，Hughes等14在梯形法的基础上引入拉格朗日乘子来实现能量守恒，提出了最早的能量守恒方法：约束能量法（constraint energy method,CEM）。1996年，Kuhl等15又在CEM的基础上加入动量守恒作为约束条件，发展了约束能量动量法（constraint energy momentum algorithm,CEMA）。然而由于该法在迭代计算过程中存在收敛性问题，这种通过拉格朗日乘子来实现能量守恒的方法最后停止了发展。Bathe及其团队提出了复合时程积分法（composite time integration scheme）16-17，在计算过程中将一个时间步tn,tn1分为两个子步，在tn,tn1/2时间步内使用梯形法（trapezoidal rule）求解出tn1/2时刻的构型参数；第二个时间步tn1/2,tn1内将tn1/2时刻的节点位移作为初始条件，使用三点后向欧拉法（three-point Euler backward method）求解出tn1时刻的构型，成功地实现了能量守恒和动量守恒。相比于传统的梯形法，复合时程积分法的稳定性得到了明显提高。19911995年，Simo及其团队在中点法（mid-point method）的基础上进行修改，提出了能量-动量守恒法1821（energy-momentum method,EMM）。和中点法的内力计算方式不同，Simo及其团队将时间步中点时刻的应力张量替换为时间步首尾两时刻应力张量的平均值式（16），成功实现了系统的能量守恒、线动量守恒和角动量守恒。S tS tS t2/nnn1 21=+hhh（16）由于能量-动量法在计算非线性问题时能够实现无条件稳定，许多结合了此算法的单元公式被不断提出。Crisfield及其团队先后将能量-动量法与协同转动框架下的二维梁单元22、三维梁单元23相结合，计算了桁架结构在动力荷载下的动力响应。此外，还有许多研究者尝试将能量-动量守恒法进一步优化。1998年，Armero等24在能量-动量守恒法的基础上引入耗散系数，通过将耗散系数取不同的数值可以控制能量耗散的程度，当0时，该法便回归到能量-动量法。其应力计算公式如下：S tS tS t2121/nnn1 21pp=+-+aahkhkh（17）此后一年，Kuhl及其团队将广义法7和能量-动量守恒法结合起来，在时间步的广义中点时刻tn1建立平衡方程，提出了广义能量-动量法（generalized energy-momentum method,GEMM）25-26。该法既包含能量-动量守恒法的基本思想，同时也可以排除结构高频能量震荡的影响，实现能量耗散。其应力计算公式如下：S tS tS t1nnn11apap=-+-a+-+hhhhh （18）近年来，国内也开始对非线性时程积分法展开研究。2016年，清华大学的Zhang等27在退化壳单元（degenerated shell element）的基础上引入复合时程积分法16，实现了壳结构的非线性动力分析。林汛达、李忠学:多体动力学时程积分算法研究进展建筑施工第45卷第9期19252019到2021年，北京航空航天大学的Zhang及其团队将能量守恒作为约束条件，通过在计算内力时选取三个乃至更多的插值点来提升公式的精确性，计算框架与广义能量-动量法相似2830。华南理工大学的Zhang等31以Simo的几何精确壳模型32（geometrically exact shell model）为基础，提出了带自旋自由度（drilling rotation）的几何精确壳单元公式，在此公式的基础上引入能量-动量守恒算法19，成功计算了壳体结构在动力荷载作用下的力学响应。5 时程积分法的研究展望本文对时程积分法的原理进行了介绍，总结了该算法的发展现状。如今，越来越多的研究者对时程积分法展开研究，诸如复合积分法、能量动量守恒法等性能优良的动力计算公式不断被提出，同时该算法还有很大的发展潜力，以下提出2个发展方向：1）结构模型在动力荷载的作用下往往会产生较大的位移和旋转，此时许多传统的有限元计算理论在计算时往往会产生较大误差。而一些新兴的有限元描述方法例如协同转动法（co-rotational method）在计算几何非线性问题时具有较好的准确度和收敛性。将协同转动法与时程分析法结合，可以得到准确度更高的动力计算公式。2）在进行仿真模拟时，结构模型往往由多种不同的单元构成。可以将目前已经发展完善且性能优良的单元公式与时程积分法结合，发展出适用于计算组合结构的动力计算理论。在进行数值模拟时，可以选择最符合实际情况的组合单元进行分析，提高计算的准确性。1 阿哈姆德A.谢巴纳.多体系统动力学M.4版.张青斌,丰志伟,葛健 全,译.北京:国防工业出版社,2019.2 BELLMAN R,CASTI J.Differential quadrature and long-term integrationJ.Journal of Mathematical Analysis and Applications,1971,34(2):235-238.3 ZHONG H Z,YU T.Flexural vibration analysis of an eccentric annular Mindlin plateJ.Archive of Applied Mechanics,2007,77(4):185-195.4 钟万勰,高强,彭海军.经典力学辛讲M.大连:大连理工大学出版 社,2013.5 NEWMARK N M.A method of computation for structural dynamics J.Transactions of the American Society of Civil Engineers,1962,127(1):1406-1433.6 WILSON E.A computer program for the dynamic stress analysis of underground structuresR.California:California UNIV Berkeley Structural Engineering LAB,1968.7 CHUNG J,HULBERT G M.A time integration algorithm for structural dynamics with improved numerical dissipation:The generalized-methodJ.Journal of Applied Mechanics,1993,60(2):371-375.8 HILBER H M,HUGHES T J R,TAYLOR R L.Improved numerical dissipation for time integration algorithms in structural dynamicsJ.Earthquake Engineering&Structural Dynamics,1977,5(3):283-292.9 WOOD W L,BOSSAK M,ZIENKIEWICZ O C.An alpha modification of Newmarks methodJ.International Journal for Numerical Methods in Engineering,1980,15(10):1562-1566.10 钟万勰.结构动力方程的精细时程积分法J.大连理工大学学 报,1994,34(2):131-136.11 顾元宪,陈飚松,张洪武.结构动力方程的增维精细积分法J.力学 学报,2000,32(4):447-456.12 汪梦甫,周锡元.结构动力方程的更新精细积分方法J.力学学 报,2004,36(2):191-195.13 BELYTSCHKO T,SCHOEBERLE D F.On the unconditional stability of an implicit algorithm for nonlinear structural dynamicsJ.Journal of Applied Mechanics,1975,42(4):865-869.14 HUGHES T J R,CAUGHEY T K,LIU W K.Finite-element methods for nonlinear elastodynamics which conserve energyJ.Journal of Applied Mechanics,1978,45(2):366-370.15 KUHL D,RAMM E.Constraint energy momentum algorithm and its application to non-linear dynamics of shellsJ.Computer Methods in Applied Mechanics and Engineering,1996,136(3/4):293-315.16 BATHE K J,BAIG M M I.On a composite implicit time integration procedure for nonlinear dynamicsJ.Computers&Structures,2005,83(31/32):2513-2524.17 BATHE K J.Conserving energy and momentum in nonlinear dynamics:A simple implicit time integration schemeJ.Computers&Structures,2007,85(7/8):437-445.18 SIMO J C,WONG K K.Unconditionally stable algorithms for rigid body dynamics that exactly preserve energy and momentumJ.International Journal for Numerical Methods in Engineering,1991,31(1):19-52.19 SIMO J C,TARNOW N.The discrete energy-momentum method.conserving algorithms for nonlinear elastodynamicsJ.Zeitschrift Fr Angewandte Mathematik Und Physik ZAMP,1992,43(5):757-792.20 SIMO J C,TARNOW N.A new energy and momentum conserving algorithm for the non-linear dynamics of shellsJ.International Journal for Numerical Methods in Engineering,1994,37(15):2527-2549.21 SIMO J C,TARNOW N,DOBLARE M.Non-linear dynamics of three-dimensional rods:Exact energy and momentum conserving algorithmsJ.International Journal for Numerical Methods in Engineering,1995,38(9):1431-1473.（下转第1929页）林汛达、李忠学:多体动力学时程积分算法研究进展建筑施工第45卷第9期1929降承压水方案的可行性。模拟中可以精确调整降水管井的流量，但管井降水实际施工时，管井多采用固定流量的水泵，不能细微调整管井出水量，记录流量为分段均值，所以通过降水试验最终确定深坑开挖的管井开启方案为15口40 m管井和1口50 m管井，此时坑内水位已经超过深坑所需安全水位1 m，与模拟结果有细微出入。工程止水帷幕深度设计主要考虑普挖大底板施工，局部深坑开挖依赖深井降承压水。虽然工程造价有所降低，但是深井与止水帷幕的深度接近，深坑开挖降水时对坑外影响较大，有可能引起土体沉降，影响周边土工环境。施工过程中应仅在需要时开启深井，且加快施工进度，缩短深井降水时间，保护周边土工环境。降水过程是所有降水井共同作用形成沿平面和沿深度方向的降水漏斗，降水方案设计的降水井流量是最低要求，需保证所有降水井的流量满足要求。如果局部降水井出水量不足，有可能导致基坑局部区域水位升高，不满足抗承压水突涌安全要求。在承压水层之间越流补给充沛，承压水层竖向渗透系数较高时，局部降水井效能降低导致的风险更高，尤其是距观测井较远的降水井。所以，在基坑土方开挖过程中，除坑内承压水位外，降水井流量在运行期间也需监测，确保基坑施工安全。5 结语1）基于Python软件的Flopy自编程序和当前常用Visual Modflow软件结果一致，验证了Flopy的准确性。利用Python模拟基坑在普挖大底板和深坑开挖2阶段的降水方案，降水井布置满足基坑开挖抗承压水突涌安全要求。2）生产性抽水试验结果表明，开启相应数量降压井后，坑内水位满足基坑开挖抗承压水突涌要求，降水方案具有合理性，证明了基于Python模拟地下水可以用于基坑降承压水方案分析。3）降水分析表明，降水井流量是保证降水效果的关键。在施工过程中，要确保降水管井施工质量，保证降水管井的最大出水效能和长期运行效能。在运行期间，要同时监测坑内水位和降水管井流量，维持承压水安全水位。1 上海市住房和城乡建设管理委员会.基坑工程技术标准:DG/T J08-612018S.上海:同济大学出版社,2018.2 兰韡,王卫东,常林越.超大规模深基坑工程现场抽水试验及土层变 形规律研究J.岩土力学,2022,43(10):2898-2910.3 骆祖江,成磊,张兴旺,等.悬挂式止水帷幕深基坑降水方案模拟优 化J.吉林大学学报(地球科学版),2022,52(6):1946-1956.4 BAKKER M,POST V,LANGEVIN C D,et al.Scripting MODFLOW model development using python and FloPyJ.Ground Water,2016,54(5):733-7395 康燕楠,降亚楠,苏振辉.基于NSGA-和FloPy的灌区水资源多目 标模拟优化模型J.水利与建筑工程学报,2021,19(3):17-23.杨喻声:Python软件在基坑降承压水模拟中的应用（上接第1925页）22 CRISFIELD M A,SHI J.An energy conserving co-rotational procedure for non-linear dynamics with finite elementsJ.Nonlinear Dynamics,1996,9(1):37-52.23 CRISFIELD M A,GALVANETTO U,JELENI G.Dynamics of 3-D co-rotational beamsJ.Computational Mechanics,1997,20(6):507-519.24 ARMERO F,PETCZ E.Formulation and analysis of conserving algorithms for frictionless dynamic contact/impact problemsJ.Computer Methods in Applied Mechanics and Engineering,1998,158(3/4):269-300.25 KUHL D,CRISFIELD M A.Energy-conserving and decaying Algorithms in non-linear structural dynamicsJ.International Journal for Numerical Methods in Engineering,1999,45(5):569-599.26 KUHL D,RAMM E.Generalized energy-momentum method for non-linear adaptive shell dynamicsJ.Computer Methods in Applied Mechanics and Engineering,1999,178(3/4):343-366.27 ZHANG J E,LIU D H,LIU Y H.Degenerated shell element with composite implicit time integration scheme for geometric nonlinear analysisJ.International Journal for Numerical Methods in Engineering,2016,105(7):483-513.28 ZHANG H M,XING Y F.A framework of time integration methods for nonsmooth systems with unilateral constraintsJ.Applied Mathematics and Computation,2019,363:124590.29 ZHANG H M,XING Y F,JI Y.An energy-conserving and decaying time integration method for general nonlinear dynamicsJ.International Journal for Numerical Methods in Engineering,2020,121(5):925-944.30 ZHANG H M,XING Y F.A three-parameter single-step time integration method for structural dynamic analysisJ.Acta Mechanica Sinica,2019,35(1):112-128.31 ZHANG R,STANCIULESCU I,YAO X H,et al.An energy-momentum conserving scheme for geometrically exact shells with drilling DOFsJ.Computational Mechanics,2021,67(1):341-364.32 SIMO J C,FOX D D.On a stress resultant geometrically exact shell model.Part I:Formulation and optimal parametrizationJ.Computer Methods in Applied Mechanics and Engineering,1989,72(3):267-304.