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

《大规模动力系统改进的快速精细积分方法》由会员分享，可在线阅读，更多相关《大规模动力系统改进的快速精细积分方法（14页珍藏版）》请在文档大全上搜索。

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)、原始精