柔性多體系統(tǒng)動(dòng)力學(xué)的并行迭代算法

孫 衛(wèi), 周小平

(1. 西京學(xué)院 應(yīng)用統(tǒng)計(jì)與理學(xué)系, 陜西 西安 710123; 2. 四川師范大學(xué) 數(shù)學(xué)與軟件科學(xué)學(xué)院, 四川 成都 61006)

柔性多體系統(tǒng)動(dòng)力學(xué)的并行迭代算法

孫 衛(wèi)1, 周小平2

(1. 西京學(xué)院 應(yīng)用統(tǒng)計(jì)與理學(xué)系, 陜西 西安 710123; 2. 四川師范大學(xué) 數(shù)學(xué)與軟件科學(xué)學(xué)院, 四川 成都 61006)

為了提高多體系統(tǒng)動(dòng)力學(xué)數(shù)值計(jì)算的穩(wěn)定性和精度，針對(duì)柔性多體系統(tǒng)動(dòng)力學(xué)建模方式和求解算法問(wèn)題，采用分類算法、迭加算法完成了柔性多體系統(tǒng)動(dòng)力學(xué)模型求解，利用波形松弛技術(shù)和RK離散方法設(shè)計(jì)并實(shí)現(xiàn)了柔性多體系統(tǒng)的并行迭代算法，為了降低計(jì)算的復(fù)雜性改進(jìn)這些迭代算法的收斂速度，提出了內(nèi)外步并行迭代進(jìn)程，并進(jìn)一步證明了相應(yīng)迭代算法的收斂性，完成了誤差估計(jì).

柔性; 多體系統(tǒng); 波形松弛; 并行進(jìn)程; 優(yōu)化

多體系統(tǒng)動(dòng)力學(xué)的研究核心任務(wù)是高精度地求解在外力作用下系統(tǒng)整體動(dòng)力學(xué)行為,將多體系統(tǒng)動(dòng)力學(xué)問(wèn)題模型化,對(duì)該系統(tǒng)進(jìn)行數(shù)值分析求解.由于多體系統(tǒng)的相互耦合及位形空間的位數(shù)增大,導(dǎo)致多體動(dòng)力學(xué)方程產(chǎn)生嚴(yán)重的剛性問(wèn)題且維數(shù)增多,如何提高數(shù)值計(jì)算的穩(wěn)定性和精度是多體系統(tǒng)動(dòng)力學(xué)數(shù)值解法的重要課題[1-3]，亦是本文研究工作的目的.

隨著機(jī)械工程、航天技術(shù)應(yīng)用復(fù)雜程度日益增加,多剛體系統(tǒng)模型無(wú)法考慮柔性體的變形動(dòng)力耦合效應(yīng),更能真實(shí)地反映實(shí)際系統(tǒng)動(dòng)力行為的柔性多體動(dòng)力學(xué)應(yīng)運(yùn)而生.時(shí)至今日柔性多體系統(tǒng)動(dòng)力學(xué)已成為計(jì)算力學(xué)領(lǐng)域的重要發(fā)展方向,各種新興的研究方法層出不窮,建模和數(shù)值計(jì)算難題探討越發(fā)顯現(xiàn)其重要性[4-6].由于柔性系統(tǒng)是具有無(wú)限多個(gè)自由度的連續(xù)模型,其動(dòng)力學(xué)方程通常為維數(shù)十分龐大的微分-代數(shù)混合方程組,只能用各種離散方法求出其近似解.為提高動(dòng)力學(xué)方程的數(shù)值計(jì)算精度和穩(wěn)定性,近年來(lái)出現(xiàn)了一系列改進(jìn)的計(jì)算策略,然而微分-代數(shù)方程的求解仍是多體系統(tǒng)動(dòng)力學(xué)與控制中的一個(gè)難點(diǎn)[7-9],仍無(wú)非常通用和程式化的方法.

針對(duì)上述難題,將波形松弛技術(shù)應(yīng)用于柔性多體系統(tǒng)演化方程的計(jì)算剖分中,對(duì)數(shù)值計(jì)算的數(shù)據(jù)進(jìn)行采集解耦處理,采用分步算法、迭加算法完成方程的并行迭代;同時(shí)考慮到隱式Runge-Kutta方法具有單步方法的固有優(yōu)點(diǎn),特別是具有處理剛性問(wèn)題所必備的穩(wěn)定性和高精度，所以選用隱式Runge-Kutta方法完成模型的時(shí)間域離散.

基于上述迭代技術(shù)和離散方法,設(shè)計(jì)了基于柔性多體系統(tǒng)演化方程RK離散的波形松弛迭代技術(shù)并行進(jìn)程,關(guān)于修正牛頓(modified Newton)法的內(nèi)外步并行迭代進(jìn)程,論證了通過(guò)降低計(jì)算的復(fù)雜性可以改進(jìn)這些迭代算法的收斂速度.對(duì)于內(nèi)外步迭代算法證明了算法的收斂性,并分別做了相應(yīng)的誤差估計(jì).

1 柔性多體系統(tǒng)動(dòng)力學(xué)程式化模型的實(shí)現(xiàn)

由于系統(tǒng)部件大范圍運(yùn)動(dòng)和構(gòu)件本身變形的相互耦合作用,給柔性多體系統(tǒng)動(dòng)力學(xué)建模與分析計(jì)算帶來(lái)了許多困難.基于不同建模方法得到的多體系統(tǒng)動(dòng)力學(xué)方程形式不盡相同,運(yùn)動(dòng)方程可以寫在許多不同的方式,其中之一是著名的Lagrange方程第一形式[1,7]

(1)

dC(t,q)=G(t,q)dq+b(t,q)dt=0

成立.在定常約束條件下,非完整約束方程形式為

亦為速度極約束.

基于完整約束的Lagrange模式為非線性微分-代數(shù)方程,由表達(dá)式(1)和運(yùn)動(dòng)約束C(t,q)=0組成

(2)

若廣義質(zhì)量陣M(t,q)∈Rnq×nq和約束陣G(t,q)∈Rnλ×nq滿秩,方程(2)為指標(biāo)-3微分-代數(shù)方程.完整約束C(t,q)=0對(duì)t求二階微分

ψ=-(Ctt(t,q)+2Ctq(t,q)z+Cqq(t,q)z2),

設(shè)初始條件為ψ(t0,z0,q0)=0,z0=z(t0),q0=q(t0),當(dāng)且僅當(dāng)矩陣G滿秩時(shí)方程(2)有唯一解.

(3)

其中,qn和zn分別是q(tn-1+cih)和z(tn-1+cih)在網(wǎng)格點(diǎn)tn處的近似值,?為Kronecker乘積,I為s×s單位陣,A為nq×nq矩陣,向量b、d、e是nq維向量,且滿足條件

(d1,…,dnq)=(b1,…,bnq)A-1.

2 隱式RUNGE-KUTTA迭代進(jìn)程

為了求解微分方程(2)考慮動(dòng)力學(xué)迭代或波形松弛技術(shù).波形松弛技術(shù)作為一類重要的解耦算法具有內(nèi)在并行性的優(yōu)勢(shì),使實(shí)際系統(tǒng)靈活選擇適合的迭代算法,波形松弛方法的實(shí)現(xiàn)可以根本地改進(jìn)數(shù)值方法的穩(wěn)定性能,迭代技術(shù)具體實(shí)現(xiàn)形式為

其中q(0)(·)和z(0)(·)是給出的初始猜測(cè)函數(shù),在時(shí)間域[t0,tT]上假設(shè)q(0)(·)≡q0和z(0)(·)≡z0.剖分后的函數(shù)φ:φ(t,u,u,v,v)=φ(t,u,v),矩陣

波形松弛技術(shù)最大的優(yōu)點(diǎn)是將復(fù)雜的系統(tǒng)解耦,子系統(tǒng)由不耦合或弱耦合的組成并保持著原系統(tǒng)的某些特性,進(jìn)行并行處理時(shí)可節(jié)省大量運(yùn)算時(shí)間,當(dāng)系數(shù)矩陣為髙度稀疏矩陣,適宜用計(jì)算機(jī)計(jì)算.

基于迭代技術(shù)和離散方法設(shè)計(jì)了針對(duì)柔性多體系統(tǒng)演化方程的波形松弛迭代技術(shù)的RK離散并行進(jìn)程,關(guān)于修正牛頓(modified Newton)法的內(nèi)外步并行迭代進(jìn)程.

fork=1 tomPar-do;

forn=κω+1 to (κ+1)ωPar-do;

set

end for;

end for.

(5)

2.2 修正Newton法的內(nèi)外步迭代進(jìn)程 考慮Waveform Newton迭代法,求解strongly非線性難題時(shí)該迭代法收斂條件弱且速度較快,只需迭代幾次就能達(dá)到原問(wèn)題的解所要求的精度.應(yīng)用IRK離散后的迭代張量形式為

(6)

forj=1tomPar-do

forv=1torPar-do;

set

endfor;

endfor.

(7)

分解矩陣

i,j=1,…,s

共有σ塊對(duì)角子陣(矩陣B為下三角實(shí)陣),則矩陣N0算法復(fù)雜度O(sσp-1(m*)3).

3 收斂分析

綜合迭代張量形式(6)和內(nèi)外步迭代過(guò)程(7)得

上式也可表示為誤差的遞推式

(8)

其中

v=1,…,r.

進(jìn)一步地分析內(nèi)外步迭代序列的收斂性,給出下列的廣義特征方程

(9)

其中λ是方程的特征值，設(shè)μ和b?ω分別是矩陣M的特征值和特征向量,則M(b?ω)=μ(b?ω)成立,即有下列的特征方程成立

(10)

定義z:=λh2并將廣義特征方程(9)式代入(10)式,得到

z(A-B)(b?ω)=μ(I-zB)(b?ω). (11)

若ω≠0成立,則矩陣

z(z):=z(I-zB)-1(A-B)

的特征值為μ=μ(z),其中I-zB非奇異;對(duì)于大多數(shù)IRK矩陣A來(lái)說(shuō),滿足下列條件的對(duì)角及下三角陣B存在:當(dāng)Re(z)≤0時(shí)矩陣z(z)的特征值μ(z)在單位圓里.

如果矩陣B擁有上面的特性,定義矩陣集B(A)，與矩陣A有關(guān).

定義 1 設(shè)

z(z):=z(I-zB)-1(A-B).

(12)

(13)

成立,其中j=1,…,m.更進(jìn)一步地,記

H(Δ):=-(A?I)[Φ(qn+Δ)-Φ(qn)]+

(14)

綜合(8)、(13)和(14)式,可以得到誤差遞推式

由IRK數(shù)值離散張量形式(6)可得迭代差推導(dǎo)式

當(dāng)穩(wěn)態(tài)函數(shù)R(∞)滿足‖R(∞)‖≤1時(shí),IRK張量式(6)是A-stable.類似地有下面結(jié)論.

4 數(shù)值算例分析

通過(guò)一非線性微分-代數(shù)方程算例驗(yàn)證上述的系列離散并行迭代算法.設(shè)非線性微分-代數(shù)方程的形式如下

(15)

方程的波形松弛迭代形式如下

其中常數(shù)滿足A=A1+A2,B=B1+B2,D=D1+D2.對(duì)方程的波形松弛迭代式(15)利用四階LobattoIIIARK法進(jìn)行離散,并分別用波形松弛迭代并行算法(5)和內(nèi)外步迭代并行算法(7)在網(wǎng)格實(shí)驗(yàn)室4個(gè)節(jié)點(diǎn)并行機(jī)上完成該數(shù)值算例的運(yùn)算.設(shè)步長(zhǎng)h,并行算法(5)的迭代次數(shù)k和內(nèi)外步迭代次數(shù)r、m是固定的.

定義記號(hào)εzG,εzI,εqG,εqI具體如下:

通過(guò)比較表1(εzG波動(dòng)松弛/εzI內(nèi)外部迭代)和表2(εqG波動(dòng)松弛/εqI內(nèi)外部迭代)中并行迭代算法運(yùn)算誤差結(jié)果可以看到:內(nèi)外步迭代并行算法(7)的并行迭代模式效能是優(yōu)于波形松弛迭代并行算法(5).

表 1 z值運(yùn)算誤差

表 2 q值運(yùn)算誤差

5 結(jié)語(yǔ)

對(duì)于完整約束下的柔性多體運(yùn)動(dòng)系統(tǒng)動(dòng)力學(xué)模型，提出了隱式Runge-Kutta數(shù)值離散張量形式,實(shí)現(xiàn)了基于波形松弛迭代技術(shù)、修正Newton法內(nèi)外步并行迭代進(jìn)程，這些迭代算法具有在并行機(jī)上運(yùn)行的適用性，可以實(shí)時(shí)模擬和求解復(fù)雜的機(jī)械系統(tǒng)問(wèn)題，并適用于現(xiàn)有的多體系統(tǒng)動(dòng)力學(xué)仿真軟件的實(shí)現(xiàn).

對(duì)系統(tǒng)使用波形松弛迭代法時(shí)也存在著效率較低、收斂慢的主要缺點(diǎn),為了改善迭代算法的收斂性能,采用Krylov子空間預(yù)處理技術(shù)來(lái)改進(jìn)并行迭代算法的收斂速度,并使用當(dāng)前并行算法設(shè)計(jì)主流的粗粒度構(gòu)造并行算法,以增強(qiáng)軟件的維護(hù)性和可移植性.最后針對(duì)一個(gè)非線性微分-代數(shù)算例具體地測(cè)試這些迭代算法的實(shí)用有效性,在多粒度的并行程序設(shè)計(jì)和實(shí)現(xiàn)時(shí)有多種平臺(tái)支持,本算例選擇性能較好的MPI+OpenMP應(yīng)用編程接口作為并行計(jì)算平臺(tái).

[1] ARCZEWSKI K, BLAJER W, FRACZEK J, et al. Multibody Dynamics: Computational Methods and Applications[M]. New York:Springer-Verlag,2012:135-148.

[2] 王琪,莊方方,郭易圓,等. 非光滑多體系統(tǒng)動(dòng)力學(xué)數(shù)值算法的研究進(jìn)展[J]. 力學(xué)進(jìn)展,2013,43(1):101-111.

[3] 李偉東,胡永明,韓小強(qiáng),等. 基于多體動(dòng)力學(xué)遞推算法的電動(dòng)汽車建模與仿真分析[J]. 計(jì)算力學(xué)學(xué)報(bào),2012,29(3):357-362.

[4] 戎保,芮筱亭,王國(guó)平,等. 多體系統(tǒng)動(dòng)力學(xué)研究進(jìn)展[J]. 振動(dòng)與沖擊,2011,30(7):178-187.

[5] SHABANA A A. Dynamics of Multibody Systems[M]. 3rd Rev. Cambridge:Cambridge University Press,2010:112-130.

[6] BAUCHAU O A. Flexible Multibody Dynamics[M]. New York:Springer-Verlag,2011:146-186.

[7] 丁潔玉,潘振寬. 多體系統(tǒng)動(dòng)力學(xué)微分-代數(shù)方程廣義:α投影法[J]. 工程力學(xué),2013,30(4):380-384.

[8] 徐自然,吳明兒. 彈性桿剛性體系統(tǒng)動(dòng)力學(xué)分析的廣義逆矩陣方法[J]. 計(jì)算力學(xué)學(xué)報(bào),2013,30(S):74-78.

[9] 吳鋒,高強(qiáng),鐘萬(wàn)勰. 剛-柔體動(dòng)力學(xué)方程的保辛攝動(dòng)迭代法[J]. 應(yīng)用數(shù)學(xué)和力學(xué),2014,35(4):341-352.

[10] NIELSEN M B, KRENK S. Conservative integration of rigid body motion by quaternion parameters with implicit constraints[J]. International Journal for Numerical Methods in Engineering,2012,92(8):734-752.

[11] POURSINA M, ANDERSON K S. An extended divide-and-conquer algorithm for a generalized class of multibody constraints[J]. Multibody System Dynamics,2013,29(3):235-254.

[12] CRITCHLEY J H, ANDERSON K S. A parallel logarithmic order algorithm for general multibody system dynamics[J]. Multibody System Dynamics,2004,12(1):75-93.

[13] 崔蔚,曾憲雯,趙國(guó)偉. 線性方程組正交化行處理法并行算法[J]. 四川師范大學(xué)學(xué)報(bào)(自然科學(xué)版),2004,27(5):492-496.

[14] 貢玉軍,何中全,牛森. 一類廣義非線性多值變分包含的迭代算法[J]. 廣西師范學(xué)院學(xué)報(bào)(自然科學(xué)版),2008,70(2):16-20.

[15] YOO W S, KIM K N, KIM H W, et al. Development of multibody system dynamics: computer simulations and experiments[J]. Multibody System Dynamics,2007,18:35-58.

[16] TIAN Q, ZHANG Y Q, CHEN L P, et al. Dynamics of spatial flexible multibody systems with clearance and lubricated spherical joints[J]. Computers & Structures,2009,87(13/14):913-929.

[17] ZHANG C J, DU S Y, LIU J, et al. Robust iterative closest point algorithm with bounded rotation angle for 2D registration[J]. Neurocomputing,2016,195:172-180.

[18] XU J H. An iterative algorithm for computing mean first passage times of Markov chains[J]. Applied Mathematics and Computation,2015,250(2):372-389.

[19] POPOVICIU L. Parallel iterative algotitjm for solving poisson’s equation[C]//Proceedings of the Scientific Conference. AFASES,2013:371-374.

[20] MARTIN G, LUKAS K, JONAS T, et al. On the parallel iterative solution of linear systems arising in the FEAST algorithm for computing inner eigenvalues[J]. Parallel Computing,2015,49:153-163.

[21] 闞子云,彭海軍,陳飆松,等. 開(kāi)放式多體系統(tǒng)動(dòng)力學(xué)仿真算法軟件研發(fā)(II)DAEs求解算法對(duì)比[J]. 計(jì)算力學(xué)學(xué)報(bào),2015,32(6):707-715,721.

[22] 李安志,任繼念,崔蔚. 三對(duì)角方程組通用性迭代解法[J]. 四川師范大學(xué)學(xué)報(bào)(自然科學(xué)版),2013,36(1):57-60.

[23] 馬秀騰,翟彥博,羅書(shū)強(qiáng). 多體系統(tǒng)動(dòng)力學(xué)仿真向后微分公式算法[J]. 計(jì)算機(jī)集成制造系統(tǒng),2013,19(1):119-126.

[24] 萬(wàn)波,江曉濤. 求解多值廣義混合隱似平衡問(wèn)題的迭代算法[J]. 四川師范大學(xué)學(xué)報(bào)(自然科學(xué)版),2011,34(2):197-200.

2010 MSC：65L05

(編輯 周 俊)

Parallel Iterative Algorithm Studies of Flexible Multibody System Dynamics

SUN Wei1, ZHOU Xiaoping2

( 1. Department of Applied Statistics and Science, Xijing College, Xi’an 710123, Shaanxi; 2. College of Mathematics and Software Science, Sichuan Normal University, Chengdu 610066, Sichuan

In order to improve the stability and accuracy of multibody system dynamics calculation, the paper focuses on the modeling and algorithm for flexible multibody system dynamics. By using classification algorithm and iterative algorithm, the flexible multibody system dynamics model is solved. The parallel iterative algorithm is accomplished by the waveform relaxation technique and the Runge-Kutta method. In order to reduce the complexity of the calculation and improve the convergence speed of the iterative algorithm, the inner-outer step parallel iterative process has been proposed, and the convergence of the iterative algorithm is further proved, and the error estimate has been accomplished as well.

flexible; multibody system; waveform relaxation; parallel process; optimization

2015-12-13

國(guó)家自然科學(xué)基金(11371289)和陜西省自然科學(xué)基礎(chǔ)研究計(jì)劃項(xiàng)目(2016JM1032)

孫 衛(wèi)(1964—),副教授,主要從事非線性系統(tǒng)算法的研究,E-mail:smsunwei@163.com

O193

A

1001-8395(2016)06-0869-06

10.3969/j.issn.1001-8395.2016.06.018