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

1、大规模动力系统改进的快速精细积分方法高强，吴锋，张洪武林家浩，钟万勰(大连理工大学，工业装备结构分析国家重点实验室，工程力学系，辽宁大连116024)摘要：木文提出一种针对大规模动力系统的改进的快速精细积分方法(fpim)。以精细积分方 法为基础，利用大规模动力系统矩阵的稀疏性和动力问题的物理特性，分析了矩阵指数的特 殊结构，并基于此给出一种计算大规模动力系统矩阵指数及其动力响应的高效率方法。关键词：动力系统；稀疏矩阵；精细积分；矩阵指数；快速算法1引言对于线弹性结构的动力学问题，比较成熟和常用的时域积分方法是newmark 方法111和runge-kutta方法|2_4,由于计算稳定性和精度

2、的要求，这两种方法的积 分步长一般都取得比较小，计算量比较大。钟万勰提出了精细积分方法ls7j，这 种方法计算精度非常高，稳定性好，允许采用很大的积分步长，特别是在处理刚 性问题时具有明显优势。精细积分方法提出后，得到了很多发展但是这种 方法的一个缺点是在处理规模很大的系统时，由于计算矩阵指数的计算量比较 大，效率是一个主要问题。本文针对大规模动力系统提出一种计算其动力响应的高效率算法。本文以精 细积分方法为基础15_71,利用动力问题的物理特性，利用一个关键思想，也就是 结构动力问题能量传递的有限性，来降低矩阵指数的计算量。利用上述思想，本 文分析了大规模动力结构对应矩阵指数的稀疏性质，并棊

3、于此给出一种计算矩阵 指数的高效率方法。在高效和精确计算矩阵指数的基础上，给出了人规模动力系 统响应的高效率和高精度算法。2动力系统的精细积分方法假设系统的刚度矩阵、阻尼矩阵和质量矩阵分别为k, c和m,那么结构 动力学方程为a 然科学基金(10902020, 10632030, 10721062, 2009aa04450i),辽宁省博士启动基金(20081091),辽宁省重 点实验室项h(2009s018)，铁道部科技研宂开发课题(2010t001-c)，大连理工大学青年教师培养拈金。通讯作者：张洪成，人连理工人学工程力学系，电i5: 86 411 84706249; e-mail addr

4、ess: zhanghw12a22a '(ie o其中f(r)为外力。方程(1)可以写为状态空间形式，即(2)hv +其中(3)其中p=为动量。数值积分时，将积分区间等分成步长为a的积分区间，即r0 = 0，t= h，l，tk= kh，tk+=(+ l)/z，l(4)若记v=v(，j(5)则方程(2)的解可以表示为v=tv d:exp鍵g- x)f(tk + x)d.x(6)其中t = exp(h/z)(7)通过方程(6)计算结构的动力响应，需要解决两个主要问题，一是要准确高 效的计算方程(7)定义的矩阵指数，二是要计算方程(6)中的积分。对于常见的外 载荷，方程(6)中积分大多可解析

5、积分，例如(1)如果外力在一个积分时间步长内为多项式，即f(z)= f()+f,z+l+f,/w，则tv, + tofo 4-t,f, 4-l 4-mm其屮，(=0, 1，l，m)cti2- t22),t(;2 = mti2(9)平卜膽匕，(k= 1, 2, l , m)上式中丁|7和1；7分别为矩阵指数t对应的分块矩阵，即t12(2)如果外力在一个积分时间步长内为简谐变化，即f (z)= rexp(/vvr)，则v,+ l=tv,+ exp(zw,x-i -扪 1 篆二)! - i: r 因此，上述计算过程的关键是计算矩阵指数。目前，计算矩阵指数最好的方 法是精细积分方法|5"71

6、。但是，精细积分法计算矩阵指数的计算量比较大，即使 对于h为稀疏矩阵的情况，也很难利用其稀疏性，得到的矩阵指数可能是满阵。3改进的快速精细积分方法精细积分方法基于两个要点，一个是加法定理，另一个是增量计算和存储。 对于给定矩阵h,它对应的矩阵指数exp(h)有如下性质，即exp(h) = expj2arexp(h92'，hzdi)如果tv足够大，则矩阵it比较小，那么矩阵it的矩阵指数可用如下的阶taylor 级数近似，即q'.即(12)<d = exp(h)-i + h +精细积分方法s7将h'的矩阵指数分为两部分，中-i + r， r=h +，一(13)然后对

7、增量部分应用加法定理，即r = r2+2r(14)通过执行tv次方程(14),然后将r加上单位矩阵，则得到h对应的矩阵指数。上面简要地介绍了精细积分方法，此方法编写程序，计算精度非常高。但是, 对于大规模系统，由于系统的自由度数目很大，计算矩阵指数将非常耗费时间和 内存。虽然，对于大规模动力系统，其刚度、阻尼和质量矩阵是稀疏矩阵，从而 h也是稀疏矩阵，但是直接通过以上的精细积分方法计算其矩阵指数，在计算过程，特别是方程(14)的加法定理的合并过程中，矩阵将逐渐变为满矩阵或非常稠 密的矩阵，很难利用矩阵的稀疏性质，从而计算量很大。本文利用结构动力响应的物理特点，从物理上说明，大规模动力系统对应的

8、 矩阵指数理论上应该是稀疏矩阵，这样在计算过程中可大大节约计算量，从而给 出一种计算矩阵指数的高效算法，且可节省存储要求。而原始精细积分方法之所 以得到非常稠密的矩阵是因为计算误差造成的。众所周知的事实是能量在一维杆屮的传播速度是有限值同理，在离散的结构中，虽然其能量的传播速度很难确切得到，但其动力学响应的能量传递 速度肯定是有限的，这将对矩阵指数的结构产生很大影响。根据方程(6)，如果 外力为零，贝ij方程(6)可写为如下分块形式，即y+l = t.y<- + ti2p(15p+l= t2iya- + t22pa-由上述方程可以得到矩阵指数元素的物理意义，g卩：tu的第f行第列元素表示

9、 第y个自由度上给定单位位移，而其他&由度位移为零且所有自由度动量为零 时，经过一个积分步长77后，第/个自由度的位移响应；t12的第/行第j列元素 表示第/个自由度上给定单位动量，而艽他自由度动量为零且所有自由度位移为 零时，经过一个积分步长7后，第/个自由度的位移响应；t2i的第第y列元 素表示第./个自由度上给定单位位移，而其他自由度位移为零且所有自由度动量 为零时，经过一个积分步长7?后，第/个自由度的动量响应；t22的第/行第y列元素表示第j个自由度上给定单位动量，而其他自由度动量为零且所有自由度位移为零吋，经过一个积分步长z/后，第/个自由度的动量响应。由于能量在结构屮的传

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

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

12、可将ri2，r2i，r22稀疏化。以上过程定义为一个矩阵的稀疏化。同理，方程(14)给出的加法定理同样有上述 类似的现象，因此，若直接应用方程(14),儿次合并后，矩阵将变为满阵或稠密 矩阵，同样可进行上述的矩阵稀疏化。对于大规模动力系统，对于一个合理的时间步长，某一处的扰动经过一个时 间步长后，其影响只是局部化的，不会传播到整个结构，因此系统对应的矩阵指 数一定是稀疏矩阵，则可将精细积分方法作如下改进，即对于给定矩阵h和积分 步长/7,我们可给出如下的快速精细积分算法(fpim)，来计算矩阵指数。1. 由于h是稀疏矩阵，将h按照稀疏矩阵存储；2. 确定/v，y8,91;3. 计算 h、h|;

13、4. 计算 r = h' +;5. 将矩阵r稀疏化；6. 执行如下语句，即for iter=l:nr=r*( r + 2*1);将r稀疏化；end7. 得到矩阵r后，将其与单位矩阵相加，即得到h和a对应的矩阵指数。上述计算过程与原始精细积分方法相比，只是增加了矩阵r的稀疏化过程， 但是这样的处理将极大地提高计算效率，具体比较请见数值算例。得到矩阵指数 后，还须考虑外力作用的影响。根据上面分析矩阵指数结构的相同思想，可以对 外力部分采用完全相同的处理方法，由于基本思想完全相同，本文不作详细介绍。4数值算例算例1:考虑如图1所示的弹簧质量系统，若取 =2000,则系统包贪2000个 质点，

14、本文随机选择每个质点的质量和弹簧的刚度，结果分别如图2和3。此结 构很容易写出其刚度矩阵k和质量矩阵m，而阻尼矩阵取为c=0.05k。假设 每个质点上都作用相同的外力cos(20,则系统的运动方程为mjjfe- c和 ky = f 爵i(2/)+ cos(2z)(16)其中f = 1, 1, l , if(17)分别采用木文方法(fpim)、原始精细积分方法(pim)、4阶r-k方法和 newmark方法积分此问题，稅分区间为0- 1000(s)o本文方法和原始精细稅分方法的积分步长为l(s), 4阶r-k方法采用matlab的ode45函数计算，并将绝对 误差和相对误差均设置为10_ 13,

15、 newmark方法的积分步长为0.001(s)。以r-k方法为参考解，分别计算木文方法、原始精细积分方法和newmark方法的相对 误差。位移和动量的相对误差分别定义为edl|y- yrklkiirkl|p- prk|2tpt(18)其屮yrk和分别表示r-k方法给出的位移和动量，而y和p可取本文方法、原始精细积分方法或newmark方法积分得到的位移和动量。本文方法积分到1000(s)时，各个质点的位移和动暈如图4所示。newmark方法、本文方法和原始精细积分方法的相对误差如图5中。图5表明，本文方法 和原始精细积分方法的精度都非常高，而本文方法和原始积分方法给出的结果非 常接近，并没有

16、损失计算精度，这说明了本文方法的正确性。图5同时表明 newmark方法采用0.001(s)步长积分，精度也没右达到本文方法的精度。四种方法的计算时间如表1所示，可以看到，本文方法的计算效率要大大优 于原始的精细积分方法、matlab的ode45和newmark方法。对于matlab的 ode45和newmark方法要达到较高的计算精度，必须采用小的多的积分步长， 因此效率降低，newmark方法要达到更高的精度，还必须减小步长，从而效率 还耍降低。按照前文的分析，由于矩阵指数中存在大量的零元素，原始精细积分 方法效率较低的原因在于浪费了大量零元素的乘法运算。图1弹簧质量系统图2各质点的质量

17、luauomidst;卜 -1ur. fi*，k-l»-?l:7 6 5 4 3 2 ? z?c.s»601801 # 10011201140116011801*2001node图3各弹簧的刚度b)5 0 5e3bcoe60180110011201140116011801node60180110011201140116011801node图4本文方法给出的(=1000(s)时各质点的位移和动量iou10叫neumarkio4110叫neumark10义图5算例1相对误差比较表1算例1计算效率比较fpimode45pimnewmarkcpu time(s) 18.31162.

18、7412.5440.4算例2:考虑如图6所示的平而应力的动力学问题，结构由两种材料组成，两种 材料的密度和泊松比相同，而杨氏模量不同，它们分别为r = 8? 103(kg/m3), m 0.2, £, = 3? 1011 (n/m), £22? 10n(n/m) (19)有限元网格和节点编号如图7所示，采用三节点三角形单元，质量矩阵采用集中 质量矩阵，共有3200个单元，1681个节点，施加边界条件后，共有3280个自 由度。初始位移为零，各自由度上具有初始动量1,无外力作用。分别采用木文方法(fpim)、原始精细积分方法(pim)、4阶r-k方法和 newmark方法积、

19、分此问题，积分区间为0- 4? 10' 3(s)o本文方法和原始精细釈 分方法的积分步长为10_6(s)，4阶r-k方法采用matlab的ode45函数计算，并将绝对误差和相对误差均设置为10_ 13, newmark方法的积分步长为10_ 8(s)。以r-k方法为参考解，分别计算木文方法、原始精细积分方法和newmark方法 的相对误差。位移和动量的相对误差的定义与算例1相同。本文方法积分到4' i(y30o时，各个质点x方向的位移和动量如图8所示，而方向的位移和动量如图9所示。newmark方法、本文方法和原始精细积分方法的相对误差如图10。图10表明，本文方法和原始精细积

20、分方法的精度都很 高，并且本文方法和原始积分方法给出的结果非常接近，这再次说明了本文方法 的正确性。图10同时表明newmark方法采用10_ 8(s)步长积分，相对误差的精 度仅达到约10"3 (s)量级。四种方法的计算时间如表2所示，可以看到，本文方法的计算效率要优于 matlab的ode45和newmark方法，比原始精细积分方法更是快了约220倍，与 原始精细积分方法相比效率得到了巨人的提高。lme'，r，me2, r, m0.5m 0.5mml <141'.*；*mvh2li4 - a 4« -m< 略“ ktjbtib4b461*2h

21、mmobwkmh5?ku嶙门 i »u蟮<> 图6两种材料组成的平囬应力m题hb w>«ui4 *>（«««>吻 ho參 ufl uhm«4»1»m4k，叫im1i >赛 >肩*|»«7*«鲁 1km 叫< * *«瓤4*mh*mi 鴒 m“詹<«）和鴇 m k嗱 m磷 m 5b *m» * *钃 * » *«*|ha> »wmi m« umma似籲

22、1;縐吻 m».4 鷀林霣 mb m* mmr >» w» m静籲如和礞>« w»| mv ml h k林*mfl m v4»k«« ?*<鉍摩 “h*m<<«hb5b*<?* mu him3i耗篇”鲁？'mrmm 勿h 科鴒卿<» *鑛 ?<>，> k碥肇mt hu u>笮*蝙”鲰 uh吻 *«詹，wi4-m，mmx w414-w m«4e« * mp，廖m« >, ;ai

23、w如” , *#< 一 螬 ihi3«i|i we*-iu纗鴦*j»*h»-钃 ?蠊？wk9ir 蟾r k m霸mo*?瓣> 镶 m > m« hh<om na m*m籌图7有限元网格和节点编号图84? l（rs）时刻各节点x方向的位移和动量 0123隹 s!»oj;pxjo ilemuolu图94? l（rs）时刻各节点y方向的位移和动量图10算例2相对误差比较表2算例2计算效率比较fpimode45pimnewmarkcpu time(s) 105.51070.923252.15319.2结论本文提出了一种针对大规模

24、动力系统的改进的快速精细积分方法(fpim),利 用大规模动力系统矩阵的稀疏性和动力问题能量传递有限的物理特性，表明对于 合理的积分步长，矩阵指数具有稀疏结构，并基于此给出一种计算大规模动力系 统动力响应的高效率方法。数值算例表明，本文方法仅需在原始精细积分方法基 础上进行简单修改，即可将原始精细积分方法的效率大幅度提高，并且效率也比 常用的4阶r-k方法和newmark方法高。参考文献1. n. m. newmark. a method of computation for structural dynamics. ascejournal of engineering mechanics，8

25、5(3): 67-94, 1959.2. hairer e， lubich c， wanner g. geometric numerical integration:structure-preserving algorithm for ordinaiy differential equations (secondedition). new york: springer, 2006.3. hairer e, n0rsett s p，wanner g. solving ordinary differential equationsi-no ns tiff problems (second edit

26、ion). berlin: springer, 1993.4. hairer e，wanner g. solving ordinary differential equations ll-stiff anddifferential- algebraic problems (second edition). berlin: springer，1996.5. 钟万勰。结构动力方程的精细时程积分法。大连理工大学学报，34(2):131-136， 1994。6. zhong wx, williams fw. a precise time step integration method. proceed

27、ings 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 adaptive algorithm of precise integration for transient analysis. acta mechannica solida sinica 2001; 14(3): 215-224.9. 谭述君，吴志刚，钟万勰。矩阵指数精细积分方法中参数的自适应选择。力学学报，41(6)：961-966, 2009。10. 孙雁。奇异系统矩阵的精细积