显式和隐式力学计算方法对比研究

潘科琪

（上海核工程研究设计院，上海　200233）

显式和隐式力学计算方法对比研究

潘科琪

（上海核工程研究设计院，上海200233）

0　引言

柔性多体系统动力响应的计算可以通过建立多体系统的封闭的拉格朗日第一类方程，对于刚柔耦合系统，初始时刻的条件是已知，因此，动力学方程的计算归结于求解数学上的微分代数方程初值问题。

动力学方程的求解主要有显式和隐式两类算法。常用的显式算法主要有龙格库塔 （Runge-kutta Method）和中心差分法［1］，为了保证求解的精度，显式算法要使积分步长调整到很小，导致CPU的计算时间显著增加，尤其是存在几何非线性变形时，显式算法的计算效率较低，但是其程序编写较为简单、方便。与显式算法不同，隐式算法在计算时并不跟踪那些对计算结果没有较大影响的高频振荡成分，正因如此，步长稍微大一些时也不会对结果准确性有较大的影响，能够显著提高计算效率。常用的隐式算法有纽马克（Newmark）［2］、威尔逊（Wilson）-θ法［3］、HHT（Hilber-Hughes-Taylor）法［4］、广义（Generalized-α）法［5］。对于边界条件明确、规则形状的梁或者板构件的柔性多体系统动力学方程，可以优先考虑模态缩减法离散动力学方程，以减少节点坐标的数量，减少动力学方程的维数。此时，传统的龙格库塔法还是适用的。若采用的MATLAB软件编程计算，可以调用函数ODE直接求解。

本文采用两种计算方法，分别求解单摆和曲柄滑块多体系统方程，通过MATLAB仿真编程计算，分别比较了两种计算方法求解线性和非线性动力学方程的计算效率。

1　动力学方程

柔性多体系统的封闭的拉格朗日第一类方程

如（1），其中，M、Qe及Qd分别是多体系统的广义质量阵和广义力阵，Φe约束方程的雅可比阵，λ是约束方程的拉格朗日乘子列阵。上式是典型的微分代数方程（DAEs）。

2　计算方法

2.1纽马克方法计算步骤

（1）将公式（1）变为如下形式

根据纽马克教授于1959年提出来的纽马克方法［6］，在计算中首先作了如下假设：

其中，h是时间步长，γ和β分别是积分常数。

（2）将方程（2）修改为如下形式

2.2龙哥库塔方法的计算步骤

对于边界条件明确、规则形状的梁或者板构件的柔性多体系统动力学方程，可以优先考虑模态缩减法离散动力学方程，以减少节点坐标的数量，减少动力学方程的维数，此时，传统的龙格库塔法还是适用的。也可采用的MATLAB软件中的函数ODE命令直接求解，其求解的主要步骤为：

在MATLAB软件中，积分过程中的函数dy=f（t，q）可以表示为：

上面的求解过程中没有涉及到迭代，流程简单，对计算方法的高效性要求会有所降低。但是由于在大变形的条件下，模态阶数选取的多少对计算精度的影响很大，仿真程序的通用性和精确性较差。

3　仿真算例

3.1单摆模型

单摆梁的材料和几何参数为：密度=27667kg/m3，弹性模量E=6.8952×1010N/m2，横截面面积 A=8×10-5m2，横截面惯性矩 I=1.06667×10-8m2。定曲率曲梁，其坐标系上的初始形状可以用函数描述为：

其中，αs为曲梁的中心角。

表1所示为用不同计算方法曲梁单摆的刚-柔耦合动力学方程所需计算时间。从表中可以看出，在仿真计算曲梁单摆线性模型的动力学方程时，龙格库塔法耗时3个小时，纽马克方法耗时76.93秒，增量法仅需要18.03秒，龙格库塔法计算效率最低，增量法的计算效率最高；在计算非线性模型的动力学方程时，龙格库塔法的仿真耗时显著增多。

图1　曲梁单摆

表1　曲梁单摆的计算时间比较

3.2曲柄滑块多体系统模型

为了进一步比较各种计算方法在处理有铰约束的多体系统问题的计算效率，本文对曲梁多体系统的仿真时间做了比较。表2比较了变曲率连杆的曲柄滑块机构的动力学问题所耗费的计算机时间。

图2所示为连杆为曲梁的曲柄滑块机构。曲柄的长度为1m，横截面面积A=0.02×0.02m2，变?曲率曲梁的形状与公式（14）中的相同，其中参数d=1m，l=4m，滑块的质量0.5kg。施加在曲柄的驱动约束为：

其中，ω0=π是曲柄稳态运动阶段的角速度，ts=0.5为加速时间。

从表2中可以看出，在仿真计算曲柄滑块的线性模型时，龙格库塔法耗时要10多个小时，纽马克方法耗时957.50秒，但是在计算非线性模型时，龙格库塔法耗时要100多个小时，纽马克方法耗时120多个小时，在处理多柔性体系统的非线性动力学方程时，由于铰约束方程的增加，纽马克方法相对于龙格库塔法仍然具有显著优势。

图2　曲柄滑块结构

表2　曲柄滑块多体系统的计算时间比较

4　结语

本文分别采用龙哥库塔法和纽马克方法，采用MATLAB编程求解单摆和曲柄滑块多体系统的线性和非线性的刚-柔耦合动力学方程，在相同的结果精度下，纽马克方法的计算效率要明显高于龙哥库塔法，因此，在商用或者自编程序仿真计算时，采用隐式算法能够显著提高计算效率，节省CPU的计算时间。

［1］凌复华，殷学刚，何冶奇.常微分方程数值方法及其在力学中的应用，重庆大学出版社，1990.

［2］Newmark N.M.A method of computation for structural dynamics.Journal of the Engineering Mechanics Division ASCE，1959，85（3）：67-94.

［3］Wilson E L.A computer program for the dynamics stress analysis of underground structures，SESM Report No.68-1.Division of structural Engineering and Structural Mchanics.University of California，Berkeley，CA.

［4］Hilber H M，Hughes T J R，Taylor R L.Improved numerical dissipation for time integration algorithms in structural dynamics.Earthquake Engineering and Structural Dynamics，1977，5：283-292.

［5］Chung J，Hulbert G M.A time integration algorithm for structural dynamics with improved numerical dissipation：The Generalized-Method.Journal of Applied Mechanics，1993，60，371-375.

［6］王文亮.结构动力学.复旦大学出版社，1993.

Kinetic Equations；Explicit Algorithm；Implicit Algorithms

Comparative Research on the Explicit and Implicit Mechanical Calculation Method

PAN Ke-qi
（Shanghai Nuclear Engineering Research&Design Institute，Shanghai 200233）

1007-1423（2015）20-0020-04

10.3969/j.issn.1007-1423.2015.20.005

潘科琪（1984-），女，辽宁开原人，博士研究生，工程师，研究方向为反应堆结构力学2015-05-05

2015-07-03

阐述求解刚-柔耦合动力学方程的两种计算方法及其研究进展。基于显式算法即龙哥库塔法和隐式算法方法即纽马克方法结合牛顿迭代法，分别求解单摆和曲柄滑块系统的线性和非线性的动力学方程。通过结果对比研究发现，在保证相同的计算精度下，隐式算法求解线性和非线性动力学方程的计算效率都要高于显式算法。指出在工程计算建议采用隐式算法进行求解。

动力学方程；显式算法；隐式算法

Shows two kinds of numerical computation method for rigid-flexible dynamics equation and its developments.Based on the explicit numerical method，respectively solves i.e.Runge Kutta method and implicit numerical method，i.e.Newmark method，linear and nonlinear dynamics equations for curved beam pendulum and slider-crank mechanism.In condition of the same computation accuracy，results comparison show that implicit numerical method computation efficiency is obviously higher than the explicit numerical method no matter for solving linear nor the nonlinear dynamics equation.So，the conclusion that implicit numerical method，i.e.Newmark method is recommended for calculating the mechanics problem in engineering.