指數(shù)時(shí)間差分法在旋轉(zhuǎn)柔性梁動(dòng)力學(xué)仿真中的應(yīng)用策略①

方思明 王賢明 郭書豪

(浙江工業(yè)大學(xué)機(jī)械工程學(xué)院 杭州 310014)

0 引言

在諸如機(jī)器人機(jī)械臂、直升機(jī)旋翼、渦輪機(jī)葉片等的實(shí)際工程應(yīng)用中,其力學(xué)模型往往涉及旋轉(zhuǎn)柔性梁[1]。當(dāng)使用有限單元法等數(shù)值方法離散其動(dòng)力學(xué)方程后,其離散動(dòng)力學(xué)方程通常會(huì)具有數(shù)值剛性,且剛性病態(tài)程度與空間離散化的網(wǎng)格大小相關(guān)[2]。對(duì)剛性動(dòng)力學(xué)方程進(jìn)行高效、精確、穩(wěn)定的數(shù)值求解,是柔性體系統(tǒng)動(dòng)力學(xué)研究的重要內(nèi)容,也是動(dòng)力學(xué)仿真中的難點(diǎn)之一[3]。

剛性動(dòng)力學(xué)方程數(shù)值求解的關(guān)鍵,在于選擇一個(gè)有效的微分方程時(shí)間積分算法。常用的動(dòng)力學(xué)時(shí)間積分算法有四階經(jīng)典Runge-Kutta 法、Gear 法、Adams 法、Newmark 法,等等[4]。時(shí)間積分算法一般分為顯式算法和隱式算法[5]。一般地,這兩類算法求解剛性問題時(shí)都可能存在計(jì)算效率的問題,其中顯式算法單步計(jì)算效率高。但是,當(dāng)求解剛性問題時(shí),為了避免數(shù)值發(fā)散,一般的顯式算法必須使用非常小的時(shí)間步長(zhǎng)才能穩(wěn)定地求解,而問題本身的解總體變化緩慢,并不需要很高時(shí)間分辨率的小步長(zhǎng)。這是數(shù)值剛性問題的典型特征之一。隱式算法比顯式算法更穩(wěn)定,時(shí)間步長(zhǎng)一般只需要滿足時(shí)間分辨率的要求,沒有穩(wěn)定性的額外要求;但隱式算法一般每一步都需要求解代數(shù)方程組,對(duì)于大型非線性動(dòng)力學(xué)系統(tǒng),計(jì)算代價(jià)較高甚至無法計(jì)算。因此,對(duì)于大型剛性系統(tǒng),目前主要采用顯式算法進(jìn)行求解[6]。于是,人們期望發(fā)展出一種顯式算法,在滿足時(shí)間分辨率的基本要求下,能盡可能地放大時(shí)間步長(zhǎng)以提高整體求解效率。其中這樣一類顯式算法稱為指數(shù)積分法[7],指數(shù)積分法是利用了微分方程Jacobian 矩陣或其近似矩陣的矩陣指數(shù)和相關(guān)函數(shù)的一類數(shù)值格式,能夠有效地處理一些數(shù)值剛性問題。通過矩陣指數(shù)運(yùn)算,指數(shù)積分法避免了對(duì)引起方程剛性的部分直接進(jìn)行數(shù)值積分,從而可以使用較大的步長(zhǎng)進(jìn)行穩(wěn)定的數(shù)值求解。

指數(shù)時(shí)間差分法(exponential time differencing,ETD)[8]是一種典型的指數(shù)積分法。ETD 方法根據(jù)常數(shù)變易法,得到微分方程的一個(gè)積分形式解析解公式,通過多項(xiàng)式逼近積分中包含的“非齊次項(xiàng)”函數(shù),再半解析地解出積分而得到數(shù)值解。Cox 和Matthews[8]給出了任意階ETD 格式的推導(dǎo)以及Runge-Katta 類型的ETD 格式。指數(shù)積分法中,指數(shù)矩陣的計(jì)算是限制該方法的因素之一。隨著矩陣指數(shù)計(jì)算方法[9]的發(fā)展,ETD 格式在各個(gè)領(lǐng)域得到廣泛應(yīng)用[10-13]。ETD 格式在運(yùn)用過程中,其線性算子的選取會(huì)很大程度地影響其格式的形態(tài)[7,14],然而針對(duì)線性算子選取的相關(guān)研究目前并不是十分充分。

本文針對(duì)非慣性系下的旋轉(zhuǎn)柔性梁這一典型動(dòng)力學(xué)系統(tǒng)的數(shù)值仿真,研究采用ETD 求解時(shí)其中線性算子的選取策略。旋轉(zhuǎn)柔性梁存在動(dòng)力剛化現(xiàn)象,梁的剛度會(huì)隨角速度變化。一次近似理論[15-16]建立的動(dòng)力學(xué)模型可以反映動(dòng)力剛化現(xiàn)象,此時(shí)動(dòng)力學(xué)方程的剛度矩陣隨系統(tǒng)大范圍擺動(dòng)的情況而變化?？紤]到采用ETD 求解剛性問題時(shí),線性算子選取對(duì)格式的穩(wěn)定性會(huì)有重要的影響,因此,本文在應(yīng)用ETD 格式求解旋轉(zhuǎn)柔性梁?jiǎn)栴}中,考慮將系統(tǒng)大范圍擺動(dòng)的信息納入線性算子的構(gòu)造,并通過數(shù)值實(shí)驗(yàn)驗(yàn)證了方法的有效性。

具體地,本文的后續(xù)內(nèi)容安排如下。第1 節(jié)介紹指數(shù)時(shí)間差分法(ETD)。第2 節(jié)給出了旋轉(zhuǎn)柔性梁的動(dòng)力學(xué)方程。第3 節(jié)給出了相應(yīng)的線性算子的選取策略。第4 節(jié)是數(shù)值算例,用來驗(yàn)證上述的選取策略。最后是總結(jié)與討論。

1 指數(shù)時(shí)間差分法

本節(jié)介紹指數(shù)時(shí)間差分法的基本思想以及具體的幾個(gè)格式。本文將用該格式求解非慣性系下旋轉(zhuǎn)柔性梁的動(dòng)力學(xué)方程?？紤]如下一般形式的微分方程:

其中,y=y(t)是維度為n×1 的未知變量。當(dāng)采用指數(shù)積分法求解時(shí),將式(1)右端項(xiàng)改寫f(t,y)成如下形式:

其中,n×n的常系數(shù)矩陣A稱為指數(shù)積分法的線性算子,文中簡(jiǎn)稱為線性算子。A取定之后,可以確定擬非齊次項(xiàng)(t,y)=f(t,y)-Ay。在給定y(tk+h)=yk下,求得下一時(shí)間步tk+h時(shí)刻的解,就是求解微分方程的逐步法。根據(jù)常數(shù)變易法,可以得到式(2)如下積分形式的解析解公式:

式(3)是精確成立的,但是,關(guān)于下一時(shí)刻解y(tk+1)本身是隱式的。對(duì)于一般問題,這并不方便直接求解,需要將式中的積分項(xiàng)進(jìn)行顯式化計(jì)算。ETD 方法一般采用只含時(shí)間t的函數(shù),在tk≤ζ≤tk+1內(nèi)對(duì)(ζ,y(ζ)) 進(jìn)行逼近,通過分部積分就可以解析地求得積分。例如:若在tk≤ζ≤tk+1上采用零階多項(xiàng)式的常數(shù)逼近,即(ζ,y(ζ))=(tk,yk),再解析地積出積分項(xiàng),就得到了指數(shù)歐拉法,用ETD1 表示[8]:

其中,I 為n×n單位矩陣。

指數(shù)時(shí)間差分法之所以能夠有效處理剛性問題,關(guān)鍵在于通過解析地求得了包含快速變化指數(shù)項(xiàng)的積分,從而避免了對(duì)方程的剛性部分進(jìn)行直接的數(shù)值積分。而一般顯式格式數(shù)值求解剛性問題之所以不穩(wěn)定,就是由于對(duì)于快變剛性部分的直接數(shù)值積分帶來了數(shù)值誤差的快速放大。對(duì)于“擬非齊次項(xiàng)”采用不同的函數(shù)逼近方式,研究者發(fā)展出了不同的ETD 方案。一方面,產(chǎn)生了眾多的指數(shù)時(shí)間差分線性多步法[17-18];另一方面,進(jìn)一步發(fā)展出了指數(shù)時(shí)間差分Runge-Kutta(Runge-Kutta ETD,ETDRK)方法等[19]。其中,二階和四階指數(shù)時(shí)間差分Runge-Kutta 格式分別如下。

(1) 指數(shù)Runge-Kutta 的二階格式,用ETDRK2表示:

其中,r1、ak、φ1均為中間變量,公式推導(dǎo)見文獻(xiàn)[8],具體形式如下。

(2) 指數(shù)Runge-Kutta 的四階格式,用ETDRK4表示:

其中,r1、r2、r3、ak、bk、ck、φ1均為中間變量,公式推導(dǎo)見文獻(xiàn)[8],具體形式如下所示。

另外,為保證ETD 格式求解剛性微分方程,需要選取合理的線性算子。通常線性算子的選取更多是經(jīng)驗(yàn)性的,一般直接選取方程線性項(xiàng)的系數(shù)矩陣作為線性算子[7,14]。當(dāng)線性項(xiàng)不是引起方程剛性的主要部分,甚至方程中并不存在顯式的線性項(xiàng)時(shí),這樣選取可能會(huì)出現(xiàn)問題,影響指數(shù)積分法的精度和穩(wěn)定性。

2 旋轉(zhuǎn)柔性梁動(dòng)力學(xué)方程

本節(jié)中,針對(duì)旋轉(zhuǎn)柔性梁這一典型系統(tǒng),采用一次近似理論[16]建立其動(dòng)力學(xué)方程,并使用有限單元法進(jìn)行空間離散。

考慮圖1 所示的旋轉(zhuǎn)柔性梁系統(tǒng)。平面柔性梁繞O點(diǎn)做定軸轉(zhuǎn)動(dòng)。假設(shè)柔性梁的材料分布均勻且各向同性,變形時(shí)梁的橫截面仍然保持為平面且與中軸線垂直。于是,可以采用平面細(xì)長(zhǎng)Euler-Bernoulli 梁模型。梁的靜長(zhǎng)度為L(zhǎng),橫截面積為S,截面慣性矩為I,材料的質(zhì)量體積密度為ρ,彈性模量為E。OX0Y0為慣性坐標(biāo)系,固定不動(dòng),OXY為浮動(dòng)坐標(biāo)系,隨柔性梁轉(zhuǎn)動(dòng)。位于旋轉(zhuǎn)運(yùn)動(dòng)軸點(diǎn)的柔性梁端部的軸線與慣性坐標(biāo)系橫軸OX0的夾角為θ。對(duì)于柔性梁的變形位移描述如圖2 所示。P點(diǎn)為梁上任意一點(diǎn),P0點(diǎn)為P點(diǎn)未變形前的位置。

圖1 旋轉(zhuǎn)柔性梁

圖2 柔性梁在浮動(dòng)坐標(biāo)系下的變形位移描述

柔性梁上的P點(diǎn)的變形可以用2 個(gè)變量來描述:縱向變形量w1(x,t) 和橫向彎曲撓度w2(x,t)。那么,柔性梁上P點(diǎn)的變形位移可表示為

其中,wc1(x,t) 是橫向彎曲引起的梁縱向變形量:

P點(diǎn)在慣性基下的矢徑為

Θ浮動(dòng)坐標(biāo)系相對(duì)于固定坐標(biāo)系的方向余弦陣:

P0點(diǎn)關(guān)于浮動(dòng)坐標(biāo)系OXY的坐標(biāo)列向量r0=[x0]T。于是,系統(tǒng)的動(dòng)能和勢(shì)能表達(dá)式為

根據(jù)Hamilton 原理得到連續(xù)動(dòng)力學(xué)方程,并使用有限單元法以n個(gè)單元進(jìn)行空間離散,得到離散動(dòng)力學(xué)方程[20]:

其中,v=[θ,q]T,q是維度3(n+1)×1 的有限元節(jié)點(diǎn)廣義坐標(biāo)向量?？傮w質(zhì)量矩陣、總體陀螺效應(yīng)矩陣、總體剛度矩陣分別為

右端QB向量和FB向量分別為慣性力陣和廣義外力陣。

為了采用ETD 格式求解,首先將式(18)寫成標(biāo)準(zhǔn)的一階常微分方程組形式:

式中,系數(shù)矩陣中M-1K和2M-1G都與系統(tǒng)大范圍擺動(dòng)的信息相關(guān)。當(dāng)角速度隨時(shí)間變化時(shí),動(dòng)力學(xué)方程為時(shí)變常微分方程組。指數(shù)積分法求解時(shí),需要仔細(xì)考慮線性算子的選取。

3 線性算子的選取

當(dāng)使用ETD 求解時(shí)變或非線性問題時(shí),為了保證求解穩(wěn)定,需要合理地選取線性算子。針對(duì)系統(tǒng)大范圍運(yùn)動(dòng)已知的旋轉(zhuǎn)柔性梁?jiǎn)栴},考慮到其動(dòng)力學(xué)特征,根據(jù)系統(tǒng)大范圍轉(zhuǎn)動(dòng)的角速度來構(gòu)造ETD的線性算子。

在指數(shù)積分法的實(shí)際應(yīng)用中,通常經(jīng)驗(yàn)性地取方程線性項(xiàng)的系數(shù)矩陣作為格式中的線性算子。當(dāng)方程的剛性主要來源于其線性項(xiàng)時(shí),這種選取線性算子的方式是有效的。但是,當(dāng)方程的剛性并不完全來自于顯式線性項(xiàng)時(shí),若直接取顯式線性項(xiàng)作為線性算子,相應(yīng)的擬非齊次項(xiàng)仍會(huì)有較大的剛性。而這將弱化ETD 格式求解剛性問題的優(yōu)勢(shì),甚至?xí)钇渫嘶癁槌Ｒ?guī)的顯式格式。因此,對(duì)于剛性并不完全簡(jiǎn)單地來自于顯式線性項(xiàng)的問題,指數(shù)積分法中線性算子的選取,需要有進(jìn)一步的考慮。

文獻(xiàn)[21]研究了求解動(dòng)力學(xué)問題時(shí),ETD 格式中線性算子的選取對(duì)格式穩(wěn)定性的影響,并從穩(wěn)定性的角度給出了線性算子的選取建議。對(duì)于變系數(shù)或非線性動(dòng)力學(xué)問題,在無法完全動(dòng)態(tài)選取Jacobian 矩陣,所選線性算子偏差始終存在的一般情形下,選取對(duì)應(yīng)剛度略偏大而非偏小的線性算子是相對(duì)安全的。在這樣選取下,指數(shù)積分法底層方程解的變化速率大于所求解原始方程解的變化速率;相較于對(duì)應(yīng)偏差相當(dāng)而剛度偏小的線性算子選取方式,此時(shí)的指數(shù)積分法展現(xiàn)出了更好的穩(wěn)定性。這里,所謂的底層方程是指以指數(shù)積分法的線性算子作為線性項(xiàng)系數(shù)矩陣的齊次微分方程[22]。

對(duì)于旋轉(zhuǎn)柔性梁系統(tǒng),大范圍的旋轉(zhuǎn)運(yùn)動(dòng)會(huì)影響柔性梁的變形位移響應(yīng)[22]。由于離心力的存在,當(dāng)轉(zhuǎn)動(dòng)的角速度變大,柔性梁的剛度會(huì)隨之變大,出現(xiàn)動(dòng)力剛化現(xiàn)象[23]。這在高速旋轉(zhuǎn)中尤為明顯。對(duì)于大范圍運(yùn)動(dòng)已知的旋轉(zhuǎn)柔性梁系統(tǒng),根據(jù)前述以一次近似理論建立的動(dòng)力學(xué)方程式(18),其總體剛度矩陣與大范圍擺動(dòng)的信息相關(guān):

其中,Ki、Mi、D1、GT都是常數(shù)矩陣。GT矩陣元素的數(shù)量級(jí)很小,對(duì)于角加速度不太大的一般情況,其對(duì)系統(tǒng)總體動(dòng)力特性的影響較小?？傮w剛度矩陣中的第2 項(xiàng),與角速度的平方相關(guān)。當(dāng)角速度達(dá)到最大值時(shí),系統(tǒng)的各階固有頻率一般將會(huì)達(dá)到最大值。此時(shí),方程各個(gè)解分量的變化率是最快的。于是,在應(yīng)用指數(shù)積分法求解時(shí),需要考慮到旋轉(zhuǎn)梁系統(tǒng)的這一動(dòng)力學(xué)特征。如果僅以方程的線性常系數(shù)項(xiàng)系數(shù)矩陣作為線性算子,不能完整地涵蓋系統(tǒng)的快變“剛性”以及“振蕩”信息,擬非齊次項(xiàng)仍會(huì)含有較大的剛性,這將會(huì)導(dǎo)致ETD 格式求解時(shí)出現(xiàn)數(shù)值發(fā)散。

考慮到旋轉(zhuǎn)柔性梁的動(dòng)力學(xué)特征,本文將系統(tǒng)大范圍擺動(dòng)的信息納入指數(shù)積分法線性算子的構(gòu)造。總體剛度矩陣式(25)相應(yīng)的Jacobian 矩陣作為線性算子:

由式(25)可知,總體剛度矩陣與角速度和角加速度相關(guān)?？紤]到高速旋轉(zhuǎn)但角加速度較小的多數(shù)情況,如前所述,項(xiàng)對(duì)系統(tǒng)的固有特性影響很小,而角速度對(duì)柔性梁剛度和系統(tǒng)動(dòng)力學(xué)特性的影響較大。因此,本文在線性算子的構(gòu)造中只考慮角速度的變化。如果需要?jiǎng)討B(tài)選取線性算子,這樣依據(jù)單一的角速度指標(biāo)進(jìn)行調(diào)整是比較方便的。而對(duì)于實(shí)際上角速度變化范圍并不是特別大的大多數(shù)算例來說,可以一個(gè)固定的線性算子貫穿始終,這在應(yīng)用簡(jiǎn)便性和計(jì)算經(jīng)濟(jì)性上是比較理想的。結(jié)合旋轉(zhuǎn)柔性梁的動(dòng)力學(xué)特征,以最大角速度時(shí)對(duì)應(yīng)的Jacobian 矩陣作為線性算子。此時(shí),指數(shù)積分法底層方程的剛度略偏大而非偏小,從穩(wěn)定性的角度來看,一般是較為合適的。

4 數(shù)值仿真

本節(jié)通過數(shù)值實(shí)驗(yàn)來驗(yàn)證上述的線性算子選取策略。以系統(tǒng)大范圍轉(zhuǎn)動(dòng)已知的旋轉(zhuǎn)柔性梁這一典型問題為數(shù)值算例[24],從仿真精度、穩(wěn)定性以及CPU計(jì)算時(shí)間3 個(gè)方面進(jìn)行對(duì)比。梁的物理參數(shù)以及大范圍轉(zhuǎn)動(dòng)的運(yùn)動(dòng)規(guī)律均根據(jù)文獻(xiàn)[24]給定。其中,長(zhǎng)度L為8 m,梁的橫截面積S為7.2968 ×10-5m2,慣性矩I為8.2189 ×10-9m4,密度ρ為2.7667 ×103kg/m3,彈性模量E為6.8952 ×1010N/m2。系統(tǒng)大范圍轉(zhuǎn)動(dòng)的運(yùn)動(dòng)規(guī)律為

該運(yùn)動(dòng)規(guī)律表示旋轉(zhuǎn)柔性梁從靜止開始旋轉(zhuǎn)加速至穩(wěn)態(tài)轉(zhuǎn)速的情況。在下面的算例中,取加速過程時(shí)間T=15 s,最大轉(zhuǎn)速參數(shù)Ω分別取4 rad/s 和15 rad/s。對(duì)此均勻梁,以5 個(gè)單元離散,此時(shí)已有很好的空間精度。該算例是典型的數(shù)值剛性問題。

4.1 線性算子對(duì)比

為了驗(yàn)證本文線性算子選取方式的合理性,用ETD 格式對(duì)上述問題進(jìn)行數(shù)值求解,并比較采用不同的線性算子所獲得的數(shù)值結(jié)果。

對(duì)于Ω=4 的數(shù)值算例,角速度∈[0,4]。分別取不同旋轉(zhuǎn)角速度對(duì)應(yīng)的線性算子A(0)、A(1)、A(2)、A(3)、A(4)、A(5)、A(6)。其中,A(0)表示以=0 來構(gòu)造的線性算子,即取方程的線性項(xiàng)矩陣作為線性算子。A(4)表示以最大角速度=4 來構(gòu)造的線性算子,以此類推。A(1)、A(2)、A(3)是以中間角速度,而A(5)和A(6)是以更大角速度分別構(gòu)造的線性算子。本算例中,角加速度數(shù)值只在0～0.533 之間,構(gòu)造以上線性算子時(shí)角加速度均按初始時(shí)取為0,而只考慮角速度對(duì)系統(tǒng)動(dòng)力特性的影響來構(gòu)造線性算子是合理的。

首先,采用ETD1 分別根據(jù)上述的線性算子進(jìn)行數(shù)值計(jì)算,考察各參數(shù)下的ETD1 的穩(wěn)定性。表1中對(duì)比表中列出的4 個(gè)時(shí)間步長(zhǎng)下的結(jié)果?？梢钥吹?在求解時(shí)段內(nèi),用較小角速度的線性算子A(0)～A(3),只有以步長(zhǎng)h=10-4計(jì)算時(shí)的數(shù)值結(jié)果沒有出現(xiàn)發(fā)散。作為對(duì)比,較大角速度的A(4)、A(5)和A(6),在所有考察的步長(zhǎng)下都沒有出現(xiàn)數(shù)值發(fā)散。這是由于當(dāng)線性算子A(0)～A(3)沒有包含系統(tǒng)主要的剛性、振蕩信息,ETD 格式求解就容易發(fā)散。而當(dāng)線性算子包含了全部或者略為偏多的剛性、振蕩信息時(shí),格式相對(duì)不容易發(fā)散。

表1 采用ETD1 求解的CPU 計(jì)算時(shí)長(zhǎng)(仿真時(shí)長(zhǎng)50 s)

由上可知,通過選取合適的線性算子,指數(shù)積分法就可以放大步長(zhǎng)進(jìn)行數(shù)值計(jì)算,從而提高求解效率。當(dāng)然,求解的步長(zhǎng)能夠放大到多少,不僅取決于穩(wěn)定性條件,也取決于解本身的變化速度。

圖3 是梁末端位移的參考解?？梢?相對(duì)于求解參考解所用的時(shí)間步長(zhǎng)h=10-5來說,解本身的變化是極為緩慢的。從時(shí)間分辨率來說,并不需要如此小的步長(zhǎng),h=10-2大概已經(jīng)足夠。下面,考察以步長(zhǎng)h=10-2求解時(shí)的數(shù)值誤差情況,來比較ETD1 用以上各線性算子得到的數(shù)值結(jié)果。

圖3 參考解,以四階Runge-Kutta 法求解,時(shí)間步長(zhǎng)h=10-5

圖4(a)是以A(0)～A(3)求解時(shí)的數(shù)值誤差。從圖中看數(shù)值解都是發(fā)散的,發(fā)散的趨勢(shì)隨著線性算子呈現(xiàn)規(guī)律性,即線性算子中選用的角速度越接近最大角速度,數(shù)值解的發(fā)散時(shí)刻逐漸推遲?？梢?線性算子的選取對(duì)ETD 格式的發(fā)散速度會(huì)產(chǎn)生影響,較接近最大角速度的線性算子表現(xiàn)相對(duì)較好。圖4(b)是以A(4)～A(6)求解時(shí)的數(shù)值誤差?？梢钥吹?A(4)的誤差最小,而A(6)的誤差最大。單從誤差上看,A(5)和A(6)的數(shù)值誤差隨時(shí)間變大,似乎有發(fā)散的趨勢(shì)。但結(jié)合圖5(a)的解曲線可知,誤差增大是由于A(5)、A(6)的數(shù)值結(jié)果呈現(xiàn)衰減趨勢(shì)。實(shí)際求解到2000 s,兩者仍未出現(xiàn)數(shù)值發(fā)散。在使用更大步長(zhǎng)h=10-1后,數(shù)值衰減現(xiàn)象會(huì)更明顯,如圖5(b)所示。這是相應(yīng)ETD 格式的數(shù)值阻尼效應(yīng)?？梢?當(dāng)采用大于最大角速度的A(5)和A(6)時(shí),在整個(gè)求解時(shí)間內(nèi)線性算子都包含著過強(qiáng)的剛性、振蕩信息,雖然在有限時(shí)間內(nèi)可能不會(huì)出現(xiàn)數(shù)值發(fā)散,但長(zhǎng)時(shí)間的誤差累積會(huì)影響數(shù)值解的精度。

圖4 末端位移誤差圖

圖5 旋轉(zhuǎn)柔性梁的末端位移

綜上可知,采用ETD 格式數(shù)值求解非慣性系下的旋轉(zhuǎn)柔性梁系統(tǒng),線性算子需要根據(jù)角速度來構(gòu)造。從穩(wěn)定性的角度考慮,以最大角速度時(shí)的Jacobian 矩陣作為線性算子是相對(duì)合適的。這樣選取的線性算子可以兼顧穩(wěn)定性以及精度。若無法精確取到相應(yīng)線性算子,與偏小角速度相比,偏大角速度的線性算子穩(wěn)定性相對(duì)來說會(huì)更好一些,雖然在精度上會(huì)由于格式的數(shù)值阻尼而損失。

4.2 不同階ETD 格式的比較

本節(jié)中,對(duì)比不同ETD1、ETDRK2 以及ETDRK4 的數(shù)值結(jié)果。數(shù)值算例取Ω=15,此時(shí),∈[0,15]。線性算子分別取A(0)、A(15) 2 種來進(jìn)行計(jì)算。仿真時(shí)長(zhǎng)為50 s。其中,A(0)表示以=0 來構(gòu)造的線性算子,即取方程的線性項(xiàng)矩陣作為線性算子;A(15)是以最大角度構(gòu)造的線性算子。

表2 是各階ETD 格式在時(shí)間步長(zhǎng)h=10-2下計(jì)算所花的CPU 時(shí)間,其中,出現(xiàn)發(fā)散的情況用*表示。由表2 可知,就CPU 計(jì)算時(shí)間而言,ETD1 的計(jì)算效率幾乎是ETDRK2 的2 倍,是ETDRK4 的4 倍。另外,在以A(0)作為線性算子時(shí),各階ETD 格式均出現(xiàn)了數(shù)值發(fā)散;以A(15)作為線性算子時(shí),各階ETD 格式在表中所列的計(jì)算時(shí)間內(nèi)都沒有出現(xiàn)數(shù)值發(fā)散。這一結(jié)果,與本文的線性算子選取建議和前一算例的結(jié)論相符。

表2 仿真時(shí)間50 s,各格式的數(shù)值結(jié)果,步長(zhǎng)h=10 -2

圖6 是旋轉(zhuǎn)柔性梁末端位移的數(shù)值誤差隨時(shí)間的變化情況。圖6(a)是以A(0)作為線性算子時(shí)的數(shù)值誤差,可以清楚地看到ETD1、ETDRK2 以及ETDRK4 先后出現(xiàn)數(shù)值發(fā)散。圖6(b)是A(15)作為線性算子時(shí)的數(shù)值誤差,可以看到各階ETD 格式均沒有出現(xiàn)數(shù)值發(fā)散。從誤差的數(shù)值上來看,ETD1格式的誤差最大,而ETDRK2 和ETDRK4 的誤差相差不大。

圖6 末端位移絕對(duì)誤差圖,步長(zhǎng)h=10 -2

圖6 中出現(xiàn)一個(gè)“反?！钡默F(xiàn)象。ETDRK4 的精度階次比ETDRK2 高,但在圖6(b)中,后者的誤差反而比前者的小。當(dāng)0≤t≤15 時(shí),旋轉(zhuǎn)柔性梁處于加速旋轉(zhuǎn)狀態(tài),此時(shí),固定的線性算子始終與包含實(shí)時(shí)剛性信息的Jacobian 矩陣會(huì)存在偏差。作為高階格式,ETDRK4 相對(duì)于ETDRK2,對(duì)線性算子上的偏差應(yīng)該更為敏感。因此,在加速階段0≤t≤15,有可能在一定的時(shí)間段內(nèi),ETDRK4 格式是弱不穩(wěn)定的。由于時(shí)間段不長(zhǎng),不會(huì)引起數(shù)值發(fā)散,但導(dǎo)致數(shù)值解在此階段累積了較多誤差。于是,出現(xiàn)了二階格式ETDRK2 比四階格式ETDRK4 的誤差更小的現(xiàn)象。隨著時(shí)間步長(zhǎng)的縮小,ETDRK4 在加速階段的弱不穩(wěn)定情況逐漸消失,逐漸呈現(xiàn)收斂性。例如,圖7 中的ETDRK4 的精度遠(yuǎn)高于ETDRK2。對(duì)于本問題,采用ETDRK2 格式的CPU 計(jì)算時(shí)間,僅是采用ETDRK4 時(shí)的一半。綜合考慮,對(duì)實(shí)用性仿真來說,ETDRK2 是更適合用來求解該問題的ETD 格式。

圖7 末端位移絕對(duì)誤差圖,步長(zhǎng)h=10 -3

5 結(jié)論

本文針對(duì)非慣性系下旋轉(zhuǎn)柔性梁的數(shù)值仿真,給出了指數(shù)時(shí)間差分法求解該問題時(shí)的應(yīng)用策略。結(jié)合旋轉(zhuǎn)柔性梁的動(dòng)力學(xué)特征,將其大范圍轉(zhuǎn)動(dòng)的角速度納入指數(shù)積分法線性算子的構(gòu)造。由于動(dòng)力剛化現(xiàn)象,系統(tǒng)的剛度會(huì)隨角速度變化,線性算子需要根據(jù)大范圍轉(zhuǎn)動(dòng)的角速度來選取。對(duì)于本文中的數(shù)值算例,以上述的線性算子選取方式,指數(shù)積分法的穩(wěn)定性得到了改善,使得各階ETD 格式能夠采用更大的時(shí)間步長(zhǎng)來求解問題。這能夠提高旋轉(zhuǎn)柔性梁數(shù)值仿真的求解效率。其次,本文對(duì)比了不同階的ETD 格式,從精度、計(jì)算效率等方面綜合考慮,指出對(duì)于一般問題采用ETDRK2 進(jìn)行數(shù)值求解可能更適合。

值得指出的是,本文中所考察的旋轉(zhuǎn)柔性梁系統(tǒng),其大范圍運(yùn)動(dòng)相對(duì)簡(jiǎn)單。如果采用本文的線性算子構(gòu)造策略,單一固定的線性算子就能貫穿整個(gè)計(jì)算過程而獲得較好的數(shù)值結(jié)果。但是,對(duì)于系統(tǒng)大范圍運(yùn)動(dòng)未知、角速度變化范圍較大等更復(fù)雜的情況,更好的做法是階段性地根據(jù)角速度來動(dòng)態(tài)選取線性算子。而本文的相關(guān)結(jié)論可以為進(jìn)一步的研究提供基礎(chǔ)。