基于微分求積法及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月生