方思明 王賢明 郭書豪

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

0 引言

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

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

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

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

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

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

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

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

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

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

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

指數(shù)時間差分法之所以能夠有效處理剛性問題,關(guān)鍵在于通過解析地求得了包含快速變化指數(shù)項的積分,從而避免了對方程的剛性部分進(jìn)行直接的數(shù)值積分。而一般顯式格式數(shù)值求解剛性問題之所以不穩(wěn)定,就是由于對于快變剛性部分的直接數(shù)值積分帶來了數(shù)值誤差的快速放大。對于“擬非齊次項”采用不同的函數(shù)逼近方式,研究者發(fā)展出了不同的ETD 方案。一方面,產(chǎn)生了眾多的指數(shù)時間差分線性多步法[17-18];另一方面,進(jìn)一步發(fā)展出了指數(shù)時間差分Runge-Kutta(Runge-Kutta ETD,ETDRK)方法等[19]。其中,二階和四階指數(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)驗性的,一般直接選取方程線性項的系數(shù)矩陣作為線性算子[7,14]。當(dāng)線性項不是引起方程剛性的主要部分,甚至方程中并不存在顯式的線性項時,這樣選取可能會出現(xiàn)問題,影響指數(shù)積分法的精度和穩(wěn)定性。

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

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

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

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

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

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

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

P點在慣性基下的矢徑為

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

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

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

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

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

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

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

3 線性算子的選取

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

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

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

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

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

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

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

4 數(shù)值仿真

本節(jié)通過數(shù)值實驗來驗證上述的線性算子選取策略。以系統(tǒng)大范圍轉(zhuǎn)動已知的旋轉(zhuǎn)柔性梁這一典型問題為數(shù)值算例[24],從仿真精度、穩(wěn)定性以及CPU計算時間3 個方面進(jìn)行對比。梁的物理參數(shù)以及大范圍轉(zhuǎn)動的運動規(guī)律均根據(jù)文獻(xiàn)[24]給定。其中,長度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)動的運動規(guī)律為

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

4.1 線性算子對比

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

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

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

表1 采用ETD1 求解的CPU 計算時長(仿真時長50 s)

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

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

圖3 參考解,以四階Runge-Kutta 法求解,時間步長h=10-5

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

圖4 末端位移誤差圖

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

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

4.2 不同階ETD 格式的比較

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

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

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

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

圖6 末端位移絕對誤差圖,步長h=10 -2

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

圖7 末端位移絕對誤差圖,步長h=10 -3

5 結(jié)論

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

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