基于微分求积法及V-变换的大规模动力系统快速数值计算方法

汪芳宗， 廖小兵

(三峡大学 电气与新能源学院，湖北 宜昌 443002)



基于微分求积法及V-变换的大规模动力系统快速数值计算方法

汪芳宗， 廖小兵

(三峡大学 电气与新能源学院，湖北 宜昌 443002)

摘要：针对大规模动力系统动态响应的数值计算，传统的微分求积法通常在时间域上逐步离散、整体求解，存在“维数灾”问题。在多级高阶时域微分求积法的基础上，提出了基于V-变换的大规模动力系统动态响应的快速数值计算方法。利用微分求积法的加权系数矩阵满足V-变换这一重要特性，将离散后的雅可比矩阵方程进行解耦分块，推导形成了多级分块递推计算方法。数值算例表明，即使采用相当于Newmark方法2s倍的步长，微分求积法的计算精度仍比Newmark方法要高出2~3个数量级。进一步对3个不同规模的算例系统进行了测试，结果表明：相对于传统的数值计算方法，多级分块递推计算方法可以获得较大的加速比，能够显著提高大规模动力系统动态响应的计算效率。

关键词：大规模动力系统；快速计算；微分求积法；V-变换；多级分块递推方法

大规模动力系统动态响应的数值仿真可以准确地描述动力系统的动态特性和每时每刻的运动状态，一直以来是人们研究的热点问题。迄今为止，所采用的数值计算方法主要有显式4阶Runge-Kutta方法[1-2]、Newmark类方法[3-5]、Wilson类方法[6-7]、精细积分法[8-10]等。前三类方法都是单级、低阶的数值积分方法，在实际应用中只能采用较小的积分步长。然而，采用较小积分步长不仅会增加积分步数，进而增加计算量，而且随着积分步数的增加也会产生更多的“累积误差”，从而损失计算精度。精细积分方法[8]自提出以来， 由于数值稳定性好、精度高，已得到了广泛的应用，但对大规模动力系统，其矩阵指数eAt的计算量很大[9-10]。

微分求积法(Differential Quadrature Method，DQM)[11-12]是一类求解微分方程(包括常微分方程和偏微分方程)数值解的新方法，具有数学原理简单，易于实现、计算效率较高等优点[13-15]。其基本原理是将数值积分思想类似地扩展为数值微分规则，具体是将一个函数关于某个坐标方向的导数表示为沿该方向所有离散点处的函数值的加权线性组合，以此作为对函数导数的离散规则[13]。

微分求积法自20世纪70年代提出以来，在微分动力学的数值计算中得到广泛应用[16-19]。相比单级、低阶的数值计算方法，多级高阶时域微分求积法[20]是一类s级2s阶、A稳定即无条件稳定的数值计算方法，比较适合于大规模动力系统动态响应计算。传统的微分求积法通常采用显式加权系数计算方式，在时间域上逐步离散、整体求解[17, 19]；当系统规模和离散网格点数较大时，需要求解高维的线性代数方程组，其巨大的计算量使人们难以接受。

本文在形成结构动力学方程一阶模型[7]的基础上，根据s级2s阶的时域微分求积法的加权系数矩阵满足V-变换(V-transformation)[20]这一重要特性，结合分块矩阵分解技巧，将离散后的雅可比矩阵方程进行解耦分块，推导形成了结构动力学方程的多级分块递推计算方法。最后，以多自由度的弹簧-质量系统为例，对比分析了基于微分求积法的多级分块递推计算方法与Newmark方法(γ=0.5，β=0.25)的计算精度，并进一步对3个不同规模的系统算例进行了时间测试和分析。

1多级高阶时域微分求积法

设函数f(x)在整个区间上足够光滑，将函数f(x)在网格点ci,i∈(1,s)处的一阶导数近似用函数f(x)在全部网格点处的函数值的线性加权和[10]来表示，即

(1)

式中，gij称为加权系数，s是对整个区间进行划分的网格数(不包括网格起始端点c0)。微分求积法的表达式(1)可以写成矩阵形式：

f(1)(c)=G0f(c0)+Gf(c)

(2)

其中

(3)

且满足G0≡-Ge，e=[1, 1, …, 1]T∈Rs×1。

微分求积法加权系数的计算方法通常分为两类：一类是选取拉格朗日插值基函数作为试函数代入(1)式求解，由于这种方法能求解出每个元素的具体表达式，因此称为“显式表达式”求解法[19]；另一类是选取一般多项式函数代入(1)式求解，由于这种方法的解由范德蒙德矩阵及其逆矩阵的乘积形式来整体计算加权系数，因此称为“隐式表达式”求解法[20]。

文献[20]利用隐式表达式求解法，详细推导了微分求积法的加权系数矩阵满足V-变换这一重要特性，并证明了采用切比雪夫(Chebyshev)网格点、切比雪夫-高斯-洛巴托(Chebyshev-Gauss-Lobatto)网格点以及均匀网格点时，传统的微分求积法均是s级s阶的、A-稳定的数值计算方法。在此基础上，文献[20]推导出了一类新的s级2s阶的、A-稳定的时域微分求积法。

依据文献[20]，s级的时域微分求积法可以转化为等值的隐式Runge-Kutta方法[2]，即

(4)

式(4)中

(5)

(6)

(7)

(8)

(9)

(10)

式(10)即是：

(11)

因此，利用等式(8)及相关等式即可求解出s级2s阶的时域微分求积法加权系数矩阵。

2基于微分求积法及V-变换的结构动力学方程计算方法

结构动力学方程的通用表达式如下：

(12)

(13)

再令z=[xTyT]，那么方程(12)可以重新写成：

(14)

其中

(15)

然后，运用s级的时域微分求积法求解等效的结构动力学方程组(14)：

(G⊗I2n)z+(G0⊗zn)=

h(Is⊗D)z+h(Is⊗I2n)R

(16)

即为

[(G⊗I2n)-h(Is⊗D)]z=

h(Is⊗I2n)R-(G0⊗zn)

(17)

式(17)中，I2n、Is分别表示2n维和s维的单位矩阵；h是积分步长；⊗表示矩阵的直积；R，z，zn分别定义如下：

(18)

由于A=G-1，G-1G0≡-e，方程(17)可等价为

[(Is⊗I2n)-h(A⊗D)]z=

h(A⊗I2n)R+(e⊗zn)

(19)

方程(19)称为雅可比矩阵方程。显然，雅可比矩阵方程是一组s×2n维的线性代数方程组，若直接LU分解求解，不考虑稀疏矩阵的影响，其计算量约为O(8s3n3)，当系统规模和离散网格点数较大时，其计算量之大使人们难以接受。因此需要“解耦”雅克比矩阵方程来降低求解方程的维数。

定义

(20)

J=(Vs⊗I2n)[(Is⊗I2n)-

(21)

令

(22)

依据As的表达式，很容易地写出其表达式：

(23)

(24)

(25)

(26)

式中

(27)

定义

(28)

雅可比矩阵方程等效转换为

(29)

方程(29)的求解主要是基于分块矩阵的前推回代过程。前推回代过程叙述如下：

前推方程：

(30)

(31)

回代方程：

(32)

其递推公式为：

(33)

在计算出步长内点值z后，步长末端的值zn+1由下式求解：

zn+1=zn+h(bT⊗I2n)

[(G⊗I2n)z+(G0⊗I2n)zn]

(34)

为了更好地阐述基于微分求积法及V-变换的结构动力学方程计算方法，可以将该方法在每一步积分中的主要实施步骤描述如下：

1)形成质量矩阵M，阻尼矩阵C，刚度矩阵K,确定外加力向量f；

3)选择网格点类型及网格点数，利用式(8)及相关等式计算出s级2s阶的时域微分求积法的加权系数矩阵：G，G0；

4)选择积分步长h，根据式(19)形成雅克比矩阵方程，即计算J，F；

从上述推导过程可以看出：基于微分求积法及V-变换的结构动力学计算方法不仅适用于线性结构动力学模型，而且对多级高阶时域微分求积法所使用的切比雪夫(Chebyshev)网格点、切比雪夫-高斯-洛巴托(Chebyshev-Gauss-Lobatto)网格点以及均匀网格点等[13]多类网格点均有效。在单步积分内，本文的算法只涉及一个2n维矩阵I2n+βs的LU分解，降低了LU分解的维数；其余的计算均由2n维分块矩阵的递推公式完成。因此，该方法称为“基于V-变换的多级分块递推计算方法(Multi-stage Block Recursive Method，MBRM)”。

3数值算例

图1　弹簧-质量系统Fig.1 The mass-spring system

图2　500个质点的质量分布Fig.2 Distribution of the mass of the nodes

图3　各弹簧的刚度系数分布Fig.3 Distribution of the stiffness of the springs

本文分别采用基于微分求积法及V-变换的多级分块递推计算方法和Newmark方法(γ=0.5，β=0.25)求解该系统。积分区间为0~1 200 s。Newmark方法的步长选取为0.05 s；基于微分求积法及V-变换的多级分块递推计算方法的步长选取为2s倍的Newmark方法，即为2s×0.05 s，s分别取5，10，20，微分求积法采用切比雪夫网格点。由于Newmark方法也是一类无条件稳定的数值方法，已经在结构动力学计算中取得了广泛的应用，因此，以Newmark方法采用小步长0.001 s所得的计算结果作为基准值。当积分到1 200 s时，分别跟踪各个质点的位移误差和速度误差曲线。图4、图5和图6分别是s=5、s=10、s=20时，多级分块递推计算方法(MBRM)与Newmark方法的位移误差和速度误差对比曲线。

图4 MBRM(s=5,h=0.5s)与Newmark方法误差对比曲线Fig.4ErrortrajectoriescomparisonofMBRM(s=5,h=0.5s)andNewmarkmethod图5 MBRM(s=10,h=1.0s)与Newmark方法误差对比曲线Fig.5ErrortrajectoriescomparisonofMBRM(s=10,h=1.0s)andNewmarkmethod图6 MBRM(s=20,h=2.0s)与Newmark方法误差对比曲线Fig.6ErrortrajectoriescomparisonofMBRM(s=20,h=2.0s)andNewmarkmethod

进一步，本文以图1所示的弹簧-质量系统为例，扩大系统规模，测试本文算法的计算效率。继续从1~10中随机选取2 000个质点的质量和2 001个弹簧的刚度组成2 000维系统；选取10 000个质点的质量和10 001个弹簧的刚度组成10 000维系统。分别采用传统微分求积法(traditional differential quadrature method，TDQM)[19]、基于V-变换的多级分块递推计算方法、Newmark方法三类方法求解500维系统、2 000维系统、10 000维系统这3个不同规模动力系统。传统微分求积法的计算过程在文献[19]中有详细的叙述，本文不再赘述。

积分区间同样为0~1 200 s，3类方法的计算时间如表1所示。将TDQM与MBRM的计算时间比值定义为加速比β1，Newmark方法与MBRM的计算时间的比值定义为加速比β2。加速比β1和β2随系统规模的变化趋势如图7、8所示。

从图7和图8可以看出：对同一个系统，随着级数s的增加，积分的步数减少，加速比β1和β2随之缓慢增加；随着系统规模的扩大，本文提出的算法的递推计算的规模也随之增加，所获得的加速比β1和β2亦愈大。实际上，无论是级数s的增加还是系统规模的扩大，都表明具有大步长积分、多级分块递推的优越性，能提高动力系统动态响应的计算效率。

表1　计算时间比较

图7　加速比β1Fig.7 Speedup β1

图8　加速比β2Fig.8 Speedup β2

4结论

本文利用s级2s阶时域微分求积法的V-变换特性，推导出了一类基于V-变换的多级分块递推计算方法。该算法不仅对时域微分求积法的多类网格点有效，而且对线性结构动力学方程普遍适用。通过本文的推导和测试表明：在保持高精度的情况下，微分求积法的级数s越大，采用的积分步长也能越大，多级分块递推计算方法所获得的加速比也随之增加，而且随着系统规模的扩大，其所获得的加速比亦愈大。因此，基于微分求积法及V-变换的多级分块递推计算方法适合于大规模动力系统的动态响应快速数值计算，能显著提高其计算效率，可以推广应用到实际工程中。

参 考 文 献

[1] Hairer E, Nørsett S P, Wanner G. Solving ordinary differential equations I: nonstiff problems[M]. Berlin: Springer, 1987.

[2] Butcher J C. Numerical methods for ordinary differential equations[M]. New York: Wiley, 2008.

[3] Newmark N M. A method of computation for structural dynamics[J]. Journal of the Engineering Mechanics Division, ASCE, 1959, 85: 67-94.

[4] 李鸿晶, 王通, 廖旭. 关于Newmark-β法机理的一种解释[J]. 地震工程与工程振动, 2011, 31(2): 55-62.

LI Hong-jing, WANG Tong, LIAO Xu. An interpretation on Newmark-βmethods in mechanism of numerical analysis[J]. Journal of Earthquake Engineering and Engineering Vibration, 2011, 31(2): 55-62.

[5] 邢誉峰, 郭静. 与结构动特性协同的自适应Newmark方法[J]. 力学学报, 2012, 44(5): 904-911.

XING Yu-feng, GUO Jing. A self-adaptive Newmark method with parameters dependent upon structural dynamic characteristics[J]. Chinese Journal of Theoretical and Applied Mechanics, 2012, 44(5): 904-911.

[6] 方德平,王全凤. 加速度修正的Wilson-θ法的精度分析[J]. 振动与冲击, 2010, 29(6): 216-218.

FANG De-ping, WANG Quan-feng.Accuracy analysis of Wilson-θmethod with modified acceleration[J]. Journal of Vibration and Shock, 2010, 29(6): 216-218.

[7] 张雄, 王天舒. 计算动力学[M]. 北京: 清华大学出版社, 2007.

[8] 钟万勰. 结构动力方程的精细时程积分法[J]. 大连理工大学学报, 1994, 34(2): 131-136.

ZHONG Wan-xie. On precise time-integration method for structural dynamics[J]. Journal of Dalian University of Technology, 1994, 34(2): 131-136.

[9] 高强,吴锋,张洪武,等.大规模动力系统改进的快速精细积分方法[J]. 计算力学学报, 2011, 28(4): 493-498.

GAO Qiang, WU Feng, ZHANG Hong-wu, et al. A fast precise integration method for large-scale dynamic structures[J]. Chinese Journal of Computational Mechanics, 2011, 28(4): 493-498.

[10] 吴泽艳, 王立峰, 武哲. 大规模动力系统高精度増维精细积分方法快速算法[J]. 振动与冲击, 2014, 33(2): 188-192.

WU Ze-yan, WANG Li-feng, WU Zhe. Fast algorithm for precise integration with high accuracy dimension expanding method with for large-scale dynamic systems[J]. Journal of Vibration and Shock, 2014, 33(2): 188-192.

[11] Bellman R E, Casti J. Differential quadrature and long term integration[J]. Journal of Mathematical Analysis and Application, 1971, 34 (2): 235-238.

[12] Bellman R E, Kashef B G, Casti J. Differential quadrature: a technique for the rapid solution of nonlinear partial differential equations[J]. Journal of Computational Physics, 1972, 10:40-52.

[13] Shu C. Differential quadrature and its application in engineering[M]. Berlin: Springer, 2000.

[14] 王鑫伟. 微分求积法在结构力学中的应用[J]. 力学进展,1995,25(2):232-240.

WANG Xin-wei. Differential quadrature in the analysis of structural components[J]. Advances in Mechanics, 1995,25(2):232-240.

[15] 程昌钧, 朱正佑. 微分求积法及其在力学应用中的若干新进展[J]. 上海大学学报:自然科学版, 2009, 15(6): 551-559.

CHENG Chang-Jun, ZHU Zheng-you. Recent advances in differential quadrature method with application to mechanics[J]. Journal of Shanghai University:Natural Science Edition, 2009, 15(6): 551-559.

[16] Fung T C. Solving initial value problems by differential quadrature method-Part 1: first order equations[J]. International Journal for Numerical Methods in Engineering, 2001, 50(6): 1411-1427.

[17] Fung T C. Solving initial value problems by differential quadrature method-Part 2: second and higher order equations[J]. International Journal for Numerical Methods in Engineering, 2001, 50(6): 1429-1454.

[18] Fung T C.Stability and accuracy of differential quadrature method in solving dynamic problems[J]. Computer Methods in Applied Mechanics and Engineering, 2002, 191(13/14): 1311-1331.

[19] 刘剑, 王鑫伟. 基于微分求积的逐步积分法与常用时间积分法的比较[J]. 力学季刊, 2008, 29(2): 304-309.

LIU Jian, WANG Xin-wei. Comparisons of successive integration method based on differential quadrature with commonly used time integration schemes[J]. Chinese Quarterly of Mechanics, 2008, 29(2): 304-309.

[20] Wang F Z, Liao X B, Xie X. Characteristics of the differential quadrature method and its improvement[J]. Mathematical Problems in Engineering, 2015. doi:10.1155/2015/657494.

Fast numerical calculation method for large-scale dynamic systems based on differential quadrature method andV-transformation

WANGFang-zong,LIAOXiao-bing

(College of Electrical Engineering & New Energy, China Three Gorges University, Yichang 443002, China)

Abstract:For numerical simulation of large-scale dynamic systems’ response, the traditional differential quadrature method (DQM) usually adopts successively discrete and global solution in time domain, where there is the problem of “curse of dimensionality”. On the basis of the multi-stage high-order time domain differential quadrature method, a fast numerical calculation method for large-scale dynamic systems’ response based on V-transformation was proposed. Using the V-transformation processed by the weighting coefficient matrix of DQM, the whole Jacobian matrix equations involved in the traditional approach of DQM were decoupled into blocks, thus a multi-stage block recursive method was achieved. The numerical examples show that, even using a step size of 2s times that as high as in the Newmark method, the calculation precision of the differential quadrature method is about 2~3 orders higher than that of the Newmark method. Furthermore, three different scale systems were used for computational efficiency test and the results show that the multi-stage block recursive method can obtain high speedup compared with the traditional numerical methods, which can significantly improve the computational efficiency of large-scale dynamic systems’ response.

Key words:large-scale dynamic systems; fast calculation; differential quadrature method; V-transformation; multi-stage block recursive method

中图分类号:O241.8；O313.3

文献标志码：A

DOI:10.13465/j.cnki.jvs.2016.03.012

收稿日期：2015-05-27修改稿收到日期：2015-08-23

基金项目：国家自然科学基金资助项目(51377098)

第一作者 汪芳宗 男，教授，1966年1月生