大规模动力系统改进的快速精细积分方法

1、大规模动力系统改进的快速精细积分方法高强，吴锋，张洪武*自然科学基金(10902020, 10632030, 10721062, 2009AA044501)，辽宁省博士启动基金(20081091)，辽宁省重点实验室项目(2009S018)，铁道部科技研究开发课题(2010T001-C)，大连理工大学青年教师培养基金。通讯作者：张洪武，大连理工大学工程力学系，电话：86 411 84706249; E-mail address: zhanghw，林家浩，钟万勰(大连理工大学，工业装备结构分析国家重点实验室，工程力学系，辽宁大连 116024）摘要：本文提出一种针对大规模动力系统的改进的快速精细积

2、分方法(FPIM)。以精细积分方法为基础，利用大规模动力系统矩阵的稀疏性和动力问题的物理特性，分析了矩阵指数的特殊结构，并基于此给出一种计算大规模动力系统矩阵指数及其动力响应的高效率方法。关键词：动力系统；稀疏矩阵；精细积分；矩阵指数；快速算法1 引 言对于线弹性结构的动力学问题，比较成熟和常用的时域积分方法是Newmark方法1和Runge-Kutta方法2-4，由于计算稳定性和精度的要求，这两种方法的积分步长一般都取得比较小，计算量比较大。钟万勰提出了精细积分方法5-7，这种方法计算精度非常高，稳定性好，允许采用很大的积分步长，特别是在处理刚性问题时具有明显优势。精细积分方法提出后，得到了

3、很多发展8-10，但是这种方法的一个缺点是在处理规模很大的系统时，由于计算矩阵指数的计算量比较大，效率是一个主要问题。本文针对大规模动力系统提出一种计算其动力响应的高效率算法。本文以精细积分方法为基础5-7，利用动力问题的物理特性，利用一个关键思想，也就是结构动力问题能量传递的有限性，来降低矩阵指数的计算量。利用上述思想，本文分析了大规模动力结构对应矩阵指数的稀疏性质，并基于此给出一种计算矩阵指数的高效率方法。在高效和精确计算矩阵指数的基础上，给出了大规模动力系统响应的高效率和高精度算法。2动力系统的精细积分方法 假设系统的刚度矩阵、阻尼矩阵和质量矩阵分别为，和，那么结构动力学方程为其中为外力

4、。方程可以写为状态空间形式，即其中其中为动量。 数值积分时，将积分区间等分成步长为的积分区间，即若记则方程的解可以表示为其中 通过方程计算结构的动力响应，需要解决两个主要问题，一是要准确高效的计算方程定义的矩阵指数，二是要计算方程中的积分。对于常见的外载荷，方程中积分大多可解析积分，例如(1) 如果外力在一个积分时间步长内为多项式，即，则其中上式中和分别为矩阵指数对应的分块矩阵，即(2) 如果外力在一个积分时间步长内为简谐变化，即，则 因此，上述计算过程的关键是计算矩阵指数。目前，计算矩阵指数最好的方法是精细积分方法5-7。但是，精细积分法计算矩阵指数的计算量比较大，即使对于为稀疏矩阵的情况，

5、也很难利用其稀疏性，得到的矩阵指数可能是满阵。3 改进的快速精细积分方法 精细积分方法基于两个要点，一个是加法定理，另一个是增量计算和存储。对于给定矩阵，它对应的矩阵指数有如下性质，即如果足够大，则矩阵比较小，那么矩阵的矩阵指数可用如下的阶Taylor级数近似，即精细积分方法5-7将的矩阵指数分为两部分，即然后对增量部分应用加法定理，即通过执行次方程，然后将加上单位矩阵，则得到对应的矩阵指数。上面简要地介绍了精细积分方法，此方法编写程序，计算精度非常高。但是，对于大规模系统，由于系统的自由度数目很大，计算矩阵指数将非常耗费时间和内存。虽然，对于大规模动力系统，其刚度、阻尼和质量矩阵是稀疏矩阵，

6、从而也是稀疏矩阵，但是直接通过以上的精细积分方法计算其矩阵指数，在计算过程，特别是方程的加法定理的合并过程中，矩阵将逐渐变为满矩阵或非常稠密的矩阵，很难利用矩阵的稀疏性质，从而计算量很大。本文利用结构动力响应的物理特点，从物理上说明，大规模动力系统对应的矩阵指数理论上应该是稀疏矩阵，这样在计算过程中可大大节约计算量，从而给出一种计算矩阵指数的高效算法，且可节省存储要求。而原始精细积分方法之所以得到非常稠密的矩阵是因为计算误差造成的。众所周知的事实是能量在一维杆中的传播速度是有限值，同理，在离散的结构中，虽然其能量的传播速度很难确切得到，但其动力学响应的能量传递速度肯定是有限的，这将对矩阵指数的

7、结构产生很大影响。根据方程，如果外力为零，则方程可写为如下分块形式，即由上述方程可以得到矩阵指数元素的物理意义，即：的第行第列元素表示第个自由度上给定单位位移，而其他自由度位移为零且所有自由度动量为零时，经过一个积分步长后，第个自由度的位移响应；的第行第列元素表示第个自由度上给定单位动量，而其他自由度动量为零且所有自由度位移为零时，经过一个积分步长后，第个自由度的位移响应；的第行第列元素表示第个自由度上给定单位位移，而其他自由度位移为零且所有自由度动量为零时，经过一个积分步长后，第个自由度的动量响应；的第行第列元素表示第个自由度上给定单位动量，而其他自由度动量为零且所有自由度位移为零时，经过一

8、个积分步长后，第个自由度的动量响应。由于能量在结构中的传递速度是有限的，假设某个自由度上有扰动，在较小的时间内，必然只能传播到有限的自由度，而不可能传播到所有自由度。根据上面给出的的物理含义，则它们一定是稀疏矩阵。这样，既可以将矩阵指数作为稀疏矩阵计算和存储，从而节约计算量和存储空间。对于给定矩阵和积分步长，原始的精细积分方法按如下过程计算矩阵指数：首先确定，然后令，其次按照方程计算，然后执行N次方程，最后将与单位矩阵相加得到矩阵指数。若采用上述的计算过程，虽然矩阵为稀疏矩阵，但是如果观察根据方程计算得到的矩阵，可以发现它可能是一个满矩阵或很不稀疏的矩阵，但仔细观察会发现，此矩阵有很多很小的元

9、素。另一方面，根据上面的分析，此时的矩阵相当于很小的时间步长对应的矩阵指数减去一个单位矩阵，因此按照上面分析的能量传播的性质，它一定是非常稀疏的矩阵。因此，此时矩阵中非常小的非零元素是计算中的误差，它们理论上应该为零，应该将它们判断出来，并令它们为零，从而将转换为稀疏矩阵存储。具体过程如下：从上面分析给出的矩阵指数的物理意义可知，矩阵的四块子矩阵的物理意义不同，因此我们将分为四块，假设为矩阵中绝对值最大的元素，并给定一个误差标准，如，则中的元素如果满足其绝对值小于，则认为此元素应该为零，根据此原则可将稀疏存储。同样，可将稀疏化。以上过程定义为一个矩阵的稀疏化。同理，方程给出的加法定理同样有上述

10、类似的现象，因此，若直接应用方程，几次合并后，矩阵将变为满阵或稠密矩阵，同样可进行上述的矩阵稀疏化。对于大规模动力系统，对于一个合理的时间步长，某一处的扰动经过一个时间步长后，其影响只是局部化的，不会传播到整个结构，因此系统对应的矩阵指数一定是稀疏矩阵，则可将精细积分方法作如下改进，即对于给定矩阵和积分步长，我们可给出如下的快速精细积分算法(FPIM)，来计算矩阵指数。1. 由于是稀疏矩阵，将按照稀疏矩阵存储；2. 确定8,9；3. 计算；4. 计算；5. 将矩阵稀疏化； 6. 执行如下语句，即For iter=1:NR=R*( R + 2*I );将R稀疏化;End7. 得到矩阵后，将其与单

11、位矩阵相加，即得到和对应的矩阵指数。上述计算过程与原始精细积分方法相比，只是增加了矩阵的稀疏化过程，但是这样的处理将极大地提高计算效率，具体比较请见数值算例。得到矩阵指数后，还须考虑外力作用的影响。根据上面分析矩阵指数结构的相同思想，可以对外力部分采用完全相同的处理方法，由于基本思想完全相同，本文不作详细介绍。4 数值算例算例1：考虑如图1所示的弹簧质量系统，若取，则系统包含2000个质点，本文随机选择每个质点的质量和弹簧的刚度，结果分别如图2和3。此结构很容易写出其刚度矩阵和质量矩阵，而阻尼矩阵取为。假设每个质点上都作用相同的外力，则系统的运动方程为其中 分别采用本文方法(FPIM)、原始精

12、细积分方法(PIM)、4阶R-K方法和Newmark方法积分此问题，积分区间为。本文方法和原始精细积分方法的积分步长为，4阶R-K方法采用Matlab的ODE45函数计算，并将绝对误差和相对误差均设置为，Newmark方法的积分步长为。以R-K方法为参考解，分别计算本文方法、原始精细积分方法和Newmark方法的相对误差。位移和动量的相对误差分别定义为其中和分别表示R-K方法给出的位移和动量，而和可取本文方法、原始精细积分方法或Newmark方法积分得到的位移和动量。本文方法积分到时，各个质点的位移和动量如图4所示。Newmark方法、本文方法和原始精细积分方法的相对误差如图5中。图5表明，本

13、文方法和原始精细积分方法的精度都非常高，而本文方法和原始积分方法给出的结果非常接近，并没有损失计算精度，这说明了本文方法的正确性。图5同时表明Newmark方法采用步长积分，精度也没有达到本文方法的精度。四种方法的计算时间如表1所示，可以看到，本文方法的计算效率要大大优于原始的精细积分方法、Matlab的ODE45和Newmark方法。对于Matlab的ODE45和Newmark方法要达到较高的计算精度，必须采用小的多的积分步长，因此效率降低，Newmark方法要达到更高的精度，还必须减小步长，从而效率还要降低。按照前文的分析，由于矩阵指数中存在大量的零元素，原始精细积分方法效率较低的原因在于

14、浪费了大量零元素的乘法运算。图 1 弹簧质量系统图 2 各质点的质量图 3 各弹簧的刚度图 4 本文方法给出的时各质点的位移和动量图 5算例1相对误差比较表 1 算例1计算效率比较FPIMODE45PIMNewmarkCPU time(s)18.31162.7412.5440.4算例2：考虑如图6所示的平面应力的动力学问题，结构由两种材料组成，两种材料的密度和泊松比相同，而杨氏模量不同，它们分别为有限元网格和节点编号如图7所示，采用三节点三角形单元，质量矩阵采用集中质量矩阵，共有3200个单元，1681个节点，施加边界条件后，共有3280个自由度。初始位移为零，各自由度上具有初始动量1，无外力

15、作用。分别采用本文方法(FPIM)、原始精细积分方法(PIM)、4阶R-K方法和Newmark方法积分此问题，积分区间为。本文方法和原始精细积分方法的积分步长为，4阶R-K方法采用Matlab的ODE45函数计算，并将绝对误差和相对误差均设置为，Newmark方法的积分步长为。以R-K方法为参考解，分别计算本文方法、原始精细积分方法和Newmark方法的相对误差。位移和动量的相对误差的定义与算例1相同。本文方法积分到时，各个质点方向的位移和动量如图8所示，而方向的位移和动量如图9所示。Newmark方法、本文方法和原始精细积分方法的相对误差如图10。图10表明，本文方法和原始精细积分方法的精度

16、都很高，并且本文方法和原始积分方法给出的结果非常接近，这再次说明了本文方法的正确性。图10同时表明Newmark方法采用步长积分，相对误差的精度仅达到约量级。四种方法的计算时间如表2所示，可以看到，本文方法的计算效率要优于Matlab的ODE45和Newmark方法，比原始精细积分方法更是快了约220倍，与原始精细积分方法相比效率得到了巨大的提高。图6 两种材料组成的平面应力问题图 7 有限元网格和节点编号图8 时刻各节点方向的位移和动量图9 时刻各节点方向的位移和动量图10算例2相对误差比较表 2 算例2计算效率比较FPIMODE45PIMNewmarkCPU time(s)105.5107

17、0.923252.15319.2结 论本文提出了一种针对大规模动力系统的改进的快速精细积分方法(FPIM)，利用大规模动力系统矩阵的稀疏性和动力问题能量传递有限的物理特性，表明对于合理的积分步长，矩阵指数具有稀疏结构，并基于此给出一种计算大规模动力系统动力响应的高效率方法。数值算例表明，本文方法仅需在原始精细积分方法基础上进行简单修改，即可将原始精细积分方法的效率大幅度提高，并且效率也比常用的4阶R-K方法和Newmark方法高。参考文献1. N. M. Newmark. A method of computation for structural dynamics. ASCE Journal

18、 of engineering mechanics, 85(3): 67-94, 1959.2. Hairer E, Lubich C, Wanner G. Geometric Numerical Integration: Structure-Preserving Algorithm for Ordinary Differential Equations (Second Edition). New York: Springer, 2006.3. Hairer E, Nrsett S P, Wanner G. Solving Ordinary Differential Equations I-N

19、onstiff Problems (Second Edition). Berlin: Springer, 1993.4. Hairer E, Wanner G. Solving Ordinary Differential Equations II-Stiff and Differential- Algebraic Problems (Second Edition). Berlin: Springer, 1996.5. 钟万勰。结构动力方程的精细时程积分法。大连理工大学学报，34(2)：131-136，1994。6. Zhong WX, Williams FW. A precise time s

20、tep integration method. Proceedings of the Institution of Mechanical Engineers. Part C. Journal of Mechanical Engineering Science 1994; 208(6): 427-430.7. Zhong WX. On precise integration method. Journal of Computational and Applied Mathematics 2004; 163(1): 59-78.8. Zhang HW, Chen BS, Gu YX. An ada

21、ptive algorithm of precise integration for transient analysis. ACTA Mechannica Solida Sinica 2001; 14(3): 215-224.9. 谭述君，吴志刚，钟万勰。矩阵指数精细积分方法中参数的自适应选择。力学学报，41(6)： 961-966，2009。 10. 孙雁。奇异系统矩阵的精细积分。上海交通大学学报，42(8)：1217-1225，2008。A Fast Precise Integration Method for Large-scale Dynamic StructuresGao Qiang, Wu Feng, Zhang Hongwu, Lin Jiahao, Zhong Wanxie (State Key Laboratory of Structural Analysis of Indu