动力学方程数值解法：精细时程积分法

weixin

　　1、问题的提出
　　工程结构在例如突加荷载或冲击荷载的作用下，往往要求作瞬态历程的分析。由于结构几何形状的复杂性，在对空间坐标采用有限元法离散后，得到的往往是数十、百、千自由度的系统。采用本征向量展开法普通只选用低频的若干个本征解，这是不够的。此时就要采用在时域中逐步积分的方法。

　　逐步积分法已经有了多年的研究，有许多种方法可供选用。从分类上说有显式积分与隐式积分两大类。显式积分对于每一时间步效率较高，但时间步必须取得很小方能保证其稳定性。隐式积分法则可以通过恰当的参数选择，保证积分的数值稳定性，因此时间步长delt长一些也可适用。

　　熟知的Newmark法及iWlso n-θ法都是隐式积分格式。然而这些积分法也有其弱点;由于delt选得大一些，一些高频振动的分量不能正确地反映出来。又如果将其用于一个保守体系，则经过若干步的计算之后，系统的能量等不能保持守恒。这种情况被称为“ 人工阻尼”或“ 算法阻尼” ;这是积分法选用所带来的，是无可奈何的。

　　对于最常见的动力体系方程：

　　其中： M、G、K为n x n矩阵，x待求，(r)t为给定外力。应当尽力设法寻求一类积分方法，对于定常体系，可使系统积分的结果高度准确。当M为对称正定，G 为反对称，K为对称且r=0时，即系统为保守时，积分结果保持系统的守恒量不变。这类积分方法虽不是提供解析解公式，但其数值结果却是高度准确的，因此可称之为精细分析;这也就是本文的目标。

　　2、方程的变换与时程积分
　　在研究偏微分方程与计算结构力学对最优控制理论的模拟关系时，对于哈密顿体系的积分推出了一套有效的算法。虽然方程(力已经不是哈密顿体系，但其基本方法与思路仍可用于当前课题。

　　仿效哈密顿体系对偶变量的引人，令

　　此时，式(1)成为：

　　以上方程可以写成为线性体系的一般形式：

　　其中各个量对比于方程(2)及(3)为：

　　当M阵对称正定，G反对称，K为对称，且rp=rq=0时，体系成为保守，即成为自由陀螺系统。此时B、D为对称矩阵，A = -Ct。保守体系可运用变分原理来推演，并充分利用哈密顿体系的已有成果。但当有阻尼时，方程(4)并不是保守体系;本征函数的频域法用起来就感到不便，因此，采用时程积分是很自然的。

　　方程组(4)是非齐次的。从线性常微分方程的理论知，应当先找出齐次方程组的通解，然后再用叠加原理或变量代换法找出其非齐次方程的特解。因此，应首先着重研讨齐次线性方程

　　常微分方程组常用的数值计算方法往往是将微分方程化成差分方程的积分法。然而差分方程往往会破坏保守体系的守恒性质，这是难于接受的。

　　将式(6)合并写成为：

　　对于定常系统，H是常矩阵，其通解可形式地写为：

　　令时间步长delt=r，则：

　　其中：

　　间题就归结到T阵的计算，只要很精细地算出了T阵，则时程积分就成为如下的一系列的矩阵乘法了。

　　3、指数矩阵的精细计算
　　指数矩阵的精细计算方法要点是利用加法定理，有：

　　其中可选用m=2ⁿ;例如N=20，则m=1048576。由问题看出r=delt本来是不大的时间区段，从而△t=r/m将是非常小的一个时间区段了。对于△t的时间区段，有

　　其中：

　　在数值计算时至关重要的一点是，矩阵的存储只能是式(13) 的Ta，而不可以是(12)的1+Ta。因为Ta是非常小的值;当它与1相加时，就成为尾数;在计算机的舍入操作时其精度将丧失殆尽。注意倍精度数的表示范围是10-308，所以Ta的值不会消失。另外应注意倍精度数的有效位数为16 位十进数，式(12)、(13) 的展开式中的项取得太多已无太大作用，然有益无害。究竟用几项还应当依课题而定。本文的例题计算都是用式(13)的。

　　为了计算T阵先应当将式(11) 作：

　　的分解。其次应当注意：

　　因此式(14)相当于执行语句：

　　当循环结束时有：

　　式(13)、(16)、(17)便是指数矩阵的精细计算公式。

　　还应当对于此算法的精度作一分析。设H阵有

　　其中： Y 为本征向量所组成的矩阵，μ为相应的本征值。于是就可以导出

　　于是式(12)、(13)的近似相当于：

　　因此计算精度的需要应当使μ×△t很小。当然应当选用最大的abs(μ)×△t来验证。由于N=20的选择，即使有delt×abs(μ)=100×2π，即delt有100个周期，仍有△t≈10∧-4的周期，划分仍是很细的。

　　式(16)、(17)对于迭代后有Ta≈-I的情况就丧失了其相对精度;出现这种情况意味着该体系的本征解是高度阻尼的，在dlet 的时间区段内已将近衰减掉了。这类解本来也是不太关心的，因此不必在意。如果关心很短时间内的系统反应，则delt本来就应选得更小的。

　　4、非齐次方程
　　可以认为非齐次项是在时间步(tk，tk+1)内线性的，即方程为

　　式中：r0、r1是2n维给定向量。该方程可以用叠加原理求解;令ψ(t-tk) 是齐次方程的解，即：

　　于是可以写出：

　　代入式(20)，利用(21) ，可知(20) 的微分方程与初始条件都已满足。

　　数值计算中虽没有ψ(t-tk)的解析式，然而逐步积分要求提供的是tk+1=tk+r时刻的向量vk+1此时

　　而T是已算得的。因此

　　这就是所要推导的时程积分公式。

　　5、精细时程积分应用实例

　　以下代码来自于钟万勰院士《应用力学的辛数学方法》。

　　以下代码是声振论坛会员编写的基于两自由度的车辆平顺性计算程序。

　　本文理论部分摘录自钟万勰院士撰写的《结构动力方程的精细时程积分法》，该文刊登于《大连理工大学学报》1994年第2期。