一種求解非線性模態(tài)的改進(jìn)Galerkin法

李 誠， 李鴻光

(1.上海衛(wèi)星工程研究所，上海 201109；2.上海交通大學(xué) 機(jī)械系統(tǒng)與振動(dòng)國家重點(diǎn)實(shí)驗(yàn)室，上海 200240)

非線性模態(tài)理論自Ronsenberg提出開始，作為線性模態(tài)理論在包含某些剛度非線性(如幾何非線性)的振動(dòng)系統(tǒng)中的一種自然拓展，一直以來受到國內(nèi)外學(xué)者的廣泛關(guān)注。在該理論發(fā)展期間，關(guān)于非線性模態(tài)的定義、與其相關(guān)的非線性系統(tǒng)的動(dòng)力學(xué)特性、以及非線性模態(tài)的數(shù)值計(jì)算及半解析求解等方面，不少學(xué)者都做出了重大的貢獻(xiàn)。Ronsenberg針對非線性自治保守系統(tǒng)，將系統(tǒng)各自由度之間的一種同步的周期運(yùn)動(dòng)模式首先定義為非線性模態(tài)[1]。該定義也適用于線性模態(tài)，所以可以認(rèn)為非線性模態(tài)是線性模態(tài)的自然推廣。Shaw等[2]將系統(tǒng)相空間中二維不變流形上的運(yùn)動(dòng)與非線性模態(tài)運(yùn)動(dòng)相聯(lián)系，用不變流形的約束關(guān)系來描述非線性模態(tài)，且根據(jù)不變流形定理，不變流形與對應(yīng)線性系統(tǒng)的不變子空間相切于平衡點(diǎn)，不變子空間即對應(yīng)線性系統(tǒng)的線性模態(tài)。吳志強(qiáng)等[3]進(jìn)一步推廣了該思想，將非線性模態(tài)定義為系統(tǒng)模態(tài)空間中偶數(shù)維不變流形上的運(yùn)動(dòng)，提出求解系統(tǒng)動(dòng)力學(xué)方程的規(guī)范型來描述非線性模態(tài)上的動(dòng)力學(xué)特性，將非線性模態(tài)進(jìn)一步分為非耦合模態(tài)和內(nèi)共振模態(tài)。在上述定義的基礎(chǔ)上，李欣業(yè)等[4]研究了內(nèi)共振非線性系統(tǒng)的內(nèi)共振模態(tài)的分岔特性。徐鑒等[5]研究了兩自由度非對稱系統(tǒng)在非奇異時(shí)的非線性模態(tài)，著重研究了疊加解的有效性與模態(tài)動(dòng)力學(xué)方程靜態(tài)分岔間的關(guān)系。Cirillo等[6]揭示了非線性模態(tài)與Koopman算子譜特性之間的聯(lián)系。

關(guān)于非線性模態(tài)的計(jì)算，Nayfeh等[7-9]應(yīng)用了多尺度法，King等[10]根據(jù)非線性模態(tài)的頻率與能量的關(guān)系，提出了一種基于能量的級數(shù)展開式求解方法，文獻(xiàn)[11-13]采用了規(guī)范型方法，求得級數(shù)展開形式的近似解。Shaw等[14-15]根據(jù)其關(guān)于非線性模態(tài)的不變流形的定義，在提出了相應(yīng)級數(shù)展開式的漸近解法之后、發(fā)展出了基于Galerkin法的半解析方法[16]。一系列Galekin法的求解精度一般情況下均優(yōu)于基于攝動(dòng)法的漸近展開解，尤其在非線性系統(tǒng)響應(yīng)幅值較大時(shí)，Galerkin法在精度上更具有顯著優(yōu)勢。對于非線性模態(tài)所對應(yīng)的一系列周期運(yùn)動(dòng)常采用打靶法等數(shù)值算法[17-18]，結(jié)合延拓連續(xù)法獲得數(shù)值解。

提高動(dòng)力系統(tǒng)非線性模態(tài)的求解精度，其目的之一在于它可以描述系統(tǒng)接近非線性模態(tài)運(yùn)動(dòng)時(shí)各自由度之間的耦合動(dòng)力學(xué)特性。保守系統(tǒng)非線性模態(tài)的一個(gè)重要特性，是當(dāng)其對應(yīng)的非保守系統(tǒng)發(fā)生主共振時(shí)，在相空間中該受迫振動(dòng)的響應(yīng)軌線發(fā)生在保守系統(tǒng)非線性模態(tài)的周圍。所以，非線性模態(tài)的曲面幾何結(jié)構(gòu)，可用來描述非線性系統(tǒng)發(fā)生主共振時(shí)的動(dòng)力學(xué)特性，以及系統(tǒng)各自由度之間的耦合關(guān)系[19-20]。本文以一種非線性雙級隔振器系統(tǒng)為例，說明在求解非線性模態(tài)時(shí)，經(jīng)典的兩種Galerkin方法均存在一些限制和不足；提出了一種改進(jìn)的Galerkin法，可在一定程度上彌補(bǔ)經(jīng)典Galekin法的限制和不足，獲得整體更高精度的非線性模態(tài)曲面解。

1 非線性模態(tài)的控制方程

本文主要針對的是無阻尼的非線性自治系統(tǒng)，方程形如式(1)

(1)

式中：ηi(i=1,2,…,n)為該系統(tǒng)在平衡位置處的派生線性系統(tǒng)的廣義模態(tài)坐標(biāo)；ωi為系統(tǒng)第i階線性模態(tài)的固有頻率；fi為作用在該階模態(tài)振子上的非線性回復(fù)力，假設(shè)其為廣義模態(tài)坐標(biāo)位移向量的非線性函數(shù)。暫不考慮自由度振子之間出現(xiàn)內(nèi)共振的情況，根據(jù)采用不變流形形式的非線性模態(tài)的定義，系統(tǒng)某一階非線性模態(tài)的運(yùn)動(dòng)中各自由度之間滿足如下的約束關(guān)系

(2)

(3)

主坐標(biāo)的控制方程變?yōu)?/p>

(4)

待求的非線性模態(tài)約束關(guān)系也變?yōu)榱?/p>

(5)

因該非線性模態(tài)的約束曲面一般是關(guān)于主坐標(biāo)相角φ在[0,2π]上的周期連續(xù)函數(shù)，所以在后續(xù)的求解中其半解析解可以假設(shè)為關(guān)于φ的諧波函數(shù)疊加形式。將式(5)代入式(1)從自由度對應(yīng)的微分方程，并采用鏈?zhǔn)椒▌t，得

(6)

再將式(4)代入式(6)，兩邊同乘以a得到該階非線性模態(tài)的控制方程

(7)

2 Galerkin系列方法求解非線性模態(tài)

假設(shè)控制方程式(7)的解為展開式(8)的形式

(8)

式中：Pi(a,φ),Qi(a,φ)分別為關(guān)于a和φ的未知函數(shù)；C和D為展開式中各項(xiàng)的待求系數(shù)；Na和Nφ分別為展開式中關(guān)于a和φ的基函數(shù)的個(gè)數(shù)。在求解域{(a,φ)|a∈[0,a0],φ∈[0,2π]上應(yīng)用Galerkin法，依次選取T(a,φ),U(a,φ)為試函數(shù)，得到式(9)、式(10)

(9)

i=1,2,…,n;i≠M(fèi);p=1,2,…,Na;q=1,2,…,Nφ

(10)

該式一般為待求系數(shù)C和D的非線性代數(shù)方程組，可用Powell hybrid算法[21]等非線性代數(shù)方程求根數(shù)值算法求解。

非線性模態(tài)假設(shè)解的展開式形式將決定其對應(yīng)的展開式系數(shù)代數(shù)方程組的形式。因假設(shè)非線性模態(tài)的約束曲面是關(guān)于主坐標(biāo)相角φ在[0,2π]上的周期連續(xù)函數(shù)，所以T(a,φ),U(a,φ)中關(guān)于φ的展開函數(shù)可選取一系列三角函數(shù)

Tk,l(a,φ)=Lk(a)cos[(l-1)φ],Uk,l(a,φ)=Lk(a)sin(lφ)

(11)

對于保守系統(tǒng)，因?yàn)樵谄浞蔷€性模態(tài)的周期運(yùn)動(dòng)中，各自由度的響應(yīng)均保持同步，所以式(11)中關(guān)于位移響應(yīng)的展開式選取余弦函數(shù)作為基函數(shù)，關(guān)于速度響應(yīng)的展開式選取正弦函數(shù)作為基函數(shù)[16]。關(guān)于式(11)中的Lk(a)，一種方法(方法1)可選擇一組定義在[0,a0]上的多項(xiàng)式，滿足如下定義的正交性，可使積分式(9)、式(10)中相應(yīng)的乘積項(xiàng)為0，簡化最終的代數(shù)方程組

(12)

同時(shí)，對于線性剛度項(xiàng)已完全解耦的式(1)，因其非線性模態(tài)曲面與派生線性系統(tǒng)的不變子空間相切于系統(tǒng)的平衡點(diǎn)，這里的不變子空間即為(a,φ)平面，所以上述的這組多項(xiàng)式的一次導(dǎo)數(shù)需在a=0處為0。圖1展示了Pesheck等研究中給出的前7個(gè)多項(xiàng)式函數(shù)。

圖1 非線性模態(tài)解展開式采用的正交多項(xiàng)式Fig.1 Orthogonal polynomials assumed in nonlinear normal mode (NNM) expansion

在該方法中將式(11)代入式(9)、式(10)，得到的代數(shù)方程組為

(13)

(14)

為了減小待求系數(shù)非線性方程組的復(fù)雜度，便于快速求解，另一種策略(方法2)是將不變流形曲面關(guān)于主坐標(biāo)(a,φ)的求解區(qū)域沿a方向離散，劃分為幾個(gè)環(huán)狀區(qū)域，再逐個(gè)獨(dú)立求解。每個(gè)環(huán)狀求解域內(nèi)，函數(shù)L(a)假設(shè)為線性函數(shù)。兩個(gè)相鄰的環(huán)域在分別求解后，在相接區(qū)域曲面的函數(shù)值一般并不相等，形成的不變流形曲面會(huì)出現(xiàn)階躍不連續(xù)，在接下來的例子中將體現(xiàn)這一點(diǎn)。一般可取相接處的節(jié)點(diǎn)在各自求解域中的函數(shù)值的均值，作為該節(jié)點(diǎn)的函數(shù)值，以保證解曲面的連續(xù)。

本文在這兩種方法的基礎(chǔ)上，采用線性的有限元格式作為假設(shè)函數(shù)L(a)，式(15)～式(17)示意圖，如圖2所示。

圖2 采用有限元格式的分段線性函數(shù)L(a)Fig.2 Piecewise linear function as finite element approximation L(a)

(15)

(16)

(17)

將式(15)～式(17)代入式(11)后再代入式(10)，采用Galerkin法求解。不同于方法2中每個(gè)積分方程的積分域?yàn)閍方向上兩個(gè)節(jié)點(diǎn)之間的環(huán)域，此方法(方法3)中每個(gè)方程的積分域?yàn)榘瑢?yīng)節(jié)點(diǎn)的相鄰環(huán)域，得到的代數(shù)方程即同時(shí)考慮了在相關(guān)鄰域內(nèi)非線性模態(tài)的控制方程對待求曲面的約束要求。因?yàn)楣?jié)點(diǎn)處函數(shù)的連續(xù)性，各個(gè)子域上經(jīng)過Galerkin積分得到的關(guān)于節(jié)點(diǎn)函數(shù)值的代數(shù)方程之間互相耦合，需要聯(lián)立求解。相比于方法2中求解一系列獨(dú)立的代數(shù)方程組再取平均的策略，本方法得到的節(jié)點(diǎn)值整體考慮了每個(gè)節(jié)點(diǎn)相關(guān)鄰域的影響，有可能獲得更精確的解，這將在后續(xù)的例子中展示。

采用方法3獲得的代數(shù)方程組由式(18)～式(23)給出

當(dāng)p=1;q=1,2,…,Nφ;i,j=1,2,…,n;i,j≠M(fèi)時(shí)

(18)

(19)

當(dāng)p=2,3,…,Na-1;q=1,2,…,Nφ;i,j=1,2,…,n;i,j≠M(fèi)時(shí)

(20)

(21)

當(dāng)p=Na,q=1,2,…,Nφ;i,j=1,2,…,n;i,j≠M(fèi)時(shí)

(22)

(23)

相比于方法1中式(13)、式(14)里的每個(gè)方程均耦合了全部的未知數(shù)，本方法生成的方程組中每個(gè)代數(shù)方程只含有與相關(guān)節(jié)點(diǎn)相鄰的節(jié)點(diǎn)值作為未知數(shù)，每個(gè)方程組的耦合項(xiàng)大大減少。因?yàn)榇鷶?shù)方程組的這個(gè)特點(diǎn)，其對應(yīng)的Jacobian矩陣也具有一定的稀疏性。后續(xù)的非線性代數(shù)方程組求根采用基于梯度的迭代數(shù)值算法，求解最小二乘意義下方程組的余量函數(shù)的最小值，其Jacobian矩陣的函數(shù)偏導(dǎo)數(shù)采用差分法進(jìn)行近似。稀疏Jacobian矩陣示意圖，如圖3所示。

圖3(d)、圖3(e)中深色的網(wǎng)格代表非零元素。為了簡化近似Jacobian矩陣的計(jì)算量，減少調(diào)用待求方程組的計(jì)算次數(shù)，本方法根據(jù)生成的代數(shù)方程組、即式(18)～式(23)所具有的稀疏性，設(shè)計(jì)了特定的Jacobian矩陣計(jì)算流程(見附錄A)。設(shè)方程組式(18)～式(23)中從坐標(biāo)的方程個(gè)數(shù)為Nsub=n-1，每個(gè)從坐標(biāo)的位移約束函數(shù)P(a,φ)和速度函數(shù)Q(a,φ)中，關(guān)于a和φ的基函數(shù)的個(gè)數(shù)分別為Na和Nφ，總的未知數(shù)個(gè)數(shù)Nx=2(n-1)NaNφ。若不考慮Jacobian矩陣的稀疏性，采用向前差分計(jì)算待求函數(shù)在當(dāng)前迭代點(diǎn)的近似Jacobian矩陣時(shí)，總共需調(diào)用計(jì)算函數(shù)值的次數(shù)為Nf=1+Nx次，即當(dāng)前迭代點(diǎn)處的每個(gè)分量都需要單獨(dú)計(jì)算其獲得微小增量后對應(yīng)的函數(shù)值以進(jìn)行差分。參考方法3中Jacobian稀疏矩陣的示意圖3，并采用附錄A中特定的稀疏Jacobian矩陣計(jì)算流程后，總共需調(diào)用計(jì)算函數(shù)值Nf=1+2(n-1)·3·Nφ次。

圖3 稀疏Jacobian矩陣示意圖Fig.3 Schematic diagramof sparse Jacobian matrix

3 Galerkin系列方法的算例結(jié)果比較

下文以一種非線性雙級隔振器為例，研究其對應(yīng)的動(dòng)力系統(tǒng)所包含的非線性模態(tài)的動(dòng)力學(xué)特性。

3.1 雙級非線性隔振器動(dòng)力學(xué)模型

圖4給出了該隔振器的構(gòu)型示意圖，相較于文獻(xiàn)[22]中引入非線性彈性元件的雙級隔振器，此例中的隔振器考慮了另一種相關(guān)構(gòu)型，其附加質(zhì)量只與負(fù)載所在的隔振平臺(tái)通過彈性及阻尼元件連接。

圖4 非線性雙級隔振器示意圖Fig.4 Schematic diagram of a nonlinear two-stage isolator

式(24)給出了該隔振器負(fù)載所在的隔振平臺(tái)受到諧波激勵(lì)時(shí)的控制方程。

(24)

式中：kv1和kv2分別為兩個(gè)縱向彈性元件的剛度；lo1和lo2分別為側(cè)向彈性元件kh1和kh2的原長。當(dāng)兩個(gè)平臺(tái)的位移響應(yīng)滿足x1<0.4l1和x2<0.4l2[23]時(shí)，方程中的剛度非線性項(xiàng)可以用泰勒展開式近似，式(24)簡化為式(25)

(25)

將系統(tǒng)方程投影到線性模態(tài)坐標(biāo)空間，考慮其對應(yīng)的保守系統(tǒng)，得到形如式(26)的控制方程，

(26)

圖5(b)同時(shí)給出了這一組主共振響應(yīng)所對應(yīng)的無阻尼式(26)的自由響應(yīng)，該組自由響應(yīng)即是發(fā)生在系統(tǒng)非線性模態(tài)上的周期運(yùn)動(dòng)響應(yīng)，即圖5(c)中的軌線，均由打靶法數(shù)值求解出。圖5(c)中的各條強(qiáng)迫響應(yīng)曲線均在非線性模態(tài)周期軌線的周圍，且兩者變化趨勢一致，說明此處無阻尼系統(tǒng)的非線性模態(tài)可以描述原系統(tǒng)的強(qiáng)迫響應(yīng)在達(dá)到主共振峰時(shí)各自由度之間的耦合關(guān)系。圖5(c)、圖5(d)分別給出了式(26)在原物理坐標(biāo)空間和廣義模態(tài)坐標(biāo)空間內(nèi)的非線性模態(tài)曲面，發(fā)生在曲面上的周期運(yùn)動(dòng)軌線與圖5(b)中的自由響應(yīng)相對應(yīng)。

圖5 主共振受迫響應(yīng)與非線性模態(tài)自由周期響應(yīng)Fig.5 Forced responses in primary resonance and the corresponding free periodic motions on the NNM

3.2 非線性模態(tài)半解析解的結(jié)果比較

為求解與原式(25)的主共振相對應(yīng)的保守式(26)的非線性模態(tài)，采用以上的3種方法求解。表1給出了3種方法在不同求解域上選取的參數(shù)。展開式中關(guān)于a的基函數(shù)的個(gè)數(shù)Na，在方法1中指的是多項(xiàng)式基函數(shù)的最高階數(shù)；在方法2中指的是劃分的獨(dú)立子域的數(shù)量，也是獨(dú)立的代數(shù)方程組的個(gè)數(shù)；在方法3中指的是沿a方向劃分的網(wǎng)格數(shù)。對于所有的非線性代數(shù)方程組的數(shù)值求解，待求系數(shù)的迭代初值均取0，精度取1×10-6。注意到在方法1中，當(dāng)Na≥5時(shí)，某些參數(shù)條件下求根的數(shù)值迭代中止(表1中用“無”表示)，無法求得滿足方程的數(shù)值解。

表1 求解參數(shù)和計(jì)算時(shí)間 {(a,φ)|a∈[0,0.5],φ∈[0,2π]}

由圖6(a)中方法2計(jì)算的解曲面的側(cè)視圖可以看出，在兩個(gè)相鄰的求解域上分別單獨(dú)計(jì)算非線性模態(tài)曲面解，在相接處不連續(xù)，非線性模態(tài)曲面存在間斷。對節(jié)點(diǎn)處的函數(shù)值取均值后，解曲面如圖6(a)中所示。由圖6(e)可以看出，在a∈[0,0.5]內(nèi)，3種方法求得的解曲面在P1(a,0)處的投影差別不大，均接近于各個(gè)精確解的投影點(diǎn)。

圖6 非線性模態(tài)解曲面及其投影 {(a,φ)|a∈[0,0.5],φ∈[0,2π]}Fig.6 Approximation surfaces and profiles for NNM in {(a,φ)|a∈[0,0.5],φ∈[0,2π]}

圖7、圖8分別給出了該非線性模態(tài)關(guān)于位移約束曲面在a∈[0,1.0]，a∈[0,2.0]內(nèi)采用上述3種方法求得的解曲面。圖7(a)中的解分別采用了：方法1，選取Na=5，Nφ=10的展開式求解；方法2，劃分了25個(gè)獨(dú)立的環(huán)域分別求解再在節(jié)點(diǎn)處取平均；方法3，選取a方向的網(wǎng)格數(shù)為25，作為采用了方法3求得的一系列解的代表。圖8(a)中的解曲面分別采用了：方法1，選取Na=5，Nφ=10的展開式；方法2，劃分了50個(gè)環(huán)域分別求解再在相接節(jié)點(diǎn)處取平均處理；方法3獲得的解曲面a方向劃分的網(wǎng)格數(shù)為50。圖7(b)～圖7(e)、圖8(b)～圖8(e)分別給出了在a∈[0,1.0]，a∈[0,2.0]內(nèi)采用3種方法獲得的各個(gè)解曲面在P1(a,0)處的截面，同時(shí)給出對應(yīng)求解域上部分精確周期解在該截面的投影點(diǎn)。

由圖7(a)、圖7(b)及圖8(a)、圖8(b)可以看出，當(dāng)積分求解域增大時(shí)，方法1的部分假設(shè)解在值較大處與精確解存在偏差。如表2、表3所示，若展開式中的多項(xiàng)式選取更高的階數(shù)，有可能在未知系數(shù)的初始值為0時(shí)方程組非線性迭代求根不收斂；還有可能，解曲面雖然能在整個(gè)求解域內(nèi)滿足加權(quán)積分平均意義上的控制方程，但實(shí)際得到的解在a=1處或a=2處已經(jīng)遠(yuǎn)遠(yuǎn)偏離精確解。

表2 求解參數(shù)和計(jì)算時(shí)間 {(a,φ)|a∈[0,1.0],φ∈[0,2π]}Tab.2 Solution parameters and computing time {(a,φ)|a∈[0,1.0],φ∈[0,2π]}

表3 求解參數(shù)和計(jì)算時(shí)間 {(a,φ)|a∈[0,2.0],φ∈[0,2π]}

由圖7、圖8可以看出，在a∈[0.4,1.8]附近，方法2求得的解曲面在P(a,0)處的投影與各個(gè)精確解的投影點(diǎn)存在一定的偏差，網(wǎng)格密度加密后，該偏差依然存在。采用方法3，取不同網(wǎng)格密度均能獲得較好的解，解曲面在P(a,0)處的投影相比于方法1和方法2的解，更為接近各精確解的投影點(diǎn)。

圖7 非線性模態(tài)解曲面及其投影 {(a,φ)|a∈[0,1.0],φ∈[0,2π]}Fig.7 Approximation surfaces and profiles for NNM in {(a,φ)|a∈[0,1.0],φ∈[0,2π]}

圖8 非線性模態(tài)解曲面及其投影 {(a,φ)|a∈[0,2.0],φ∈[0,2π]}Fig.8 Approximation surfaces and profiles for NNM in {(a,φ)|a∈[0,2.0],φ∈[0,2π]}

3.3 一種加速計(jì)算策略

當(dāng)網(wǎng)格數(shù)較大時(shí)，因?yàn)榉椒?中每個(gè)獨(dú)立子域求解的代數(shù)方程組未知數(shù)個(gè)數(shù)都遠(yuǎn)小于對應(yīng)方法3中未知數(shù)的個(gè)數(shù)，采用方法2的所有求解域的整體計(jì)算時(shí)間仍遠(yuǎn)小于方法3的計(jì)算時(shí)間。而在方法3的上述應(yīng)用中，未知數(shù)的迭代初值均取為0。為了利用方法3獲得更精確的解曲面，同時(shí)減少計(jì)算成本，我們進(jìn)而提出一種綜合了方法2和方法3的策略，即先采用方法2求解，獲得解曲面的一個(gè)初步的全貌，再將方法2得到的解作為方法3中未知數(shù)的初值，方法2中的網(wǎng)格數(shù)、諧波基函數(shù)保持不變，采用方法3獲得相較于方法2更精確的解。本例中，方法3采用了方法2的解作為初值后的計(jì)算時(shí)間在表4中給出，即使加上方法2的求解時(shí)間，整體計(jì)算時(shí)間也大大減少。該解與初值為0時(shí)求得的解基本相同，如圖9所示。

表4 求解參數(shù)和計(jì)算時(shí)間 {(a,φ)|a∈[0,2.0],φ∈[0,2π]}Tab.4 Solution parameters and computing time {(a,φ)|a∈[0,2.0],φ∈[0,2π]}

圖9 非線性模態(tài)解曲面及其投影 {(a,φ)|a∈[0,2.0],φ∈[0,2π]}Fig.9 Approximation surfaces and profiles for NNM in {(a,φ)|a∈[0,2.0],φ∈[0,2π]}

4 結(jié) 論

本文提出的一種改進(jìn)的 Galerkin方法可以較為準(zhǔn)確地計(jì)算一類無阻尼非線性系統(tǒng)的非線性模態(tài)在不變流形定義下的解曲面。相比于已有的第一種非線性模態(tài)Galerkin求解法，該方法在假設(shè)函數(shù)中引入了一次多項(xiàng)式的有限元格式的展開式，使得待求系數(shù)的非線性代數(shù)方程組中各方程的耦合項(xiàng)減少。在非線性方程組迭代求根的算法中，根據(jù)Jacobian矩陣在該問題中特定的稀疏性，本文采用了一種特定的Jacobian矩陣計(jì)算流程，可以有效減少待求函數(shù)的計(jì)算次數(shù)，從而加速求解。相比于已有的兩種非線性模態(tài)Galerkin解法，當(dāng)該方法在求解域較大時(shí)，仍能獲得較為準(zhǔn)確的解。當(dāng)該方法與采用分段獨(dú)立求解策略的第二種Galerkin法相結(jié)合，將其解作為本方法迭代計(jì)算的初值，不僅可以保證解的精度，更能夠進(jìn)一步減少計(jì)算時(shí)間。

附錄A

非線性模態(tài)求解中稀疏Jacobian矩陣的計(jì)算流程JACdf_Sparse，F(xiàn)ortran Free format格式。

“!”后為注釋。

SUBROUTINE JACdf_Sparse(X, FCN, FJAC)

!‘X’輸入向量，計(jì)算‘X’處的Jacobian矩陣‘FJAC’作為輸出

!‘FCN’待求根的函數(shù)向量，返回‘X’處的函數(shù)值‘FF’

FJAC = 0 !Jacobian全部元素先賦0值

h = eps! 自變量增量h; eps為獲取浮點(diǎn)數(shù)精度數(shù)值的腳本

call FCN (X, FF, N) !調(diào)用FCN，計(jì)算X點(diǎn)處的函數(shù)值FF

Band: do i_Bnd = 1,3 !非零元素帶帶寬方向含有的分塊矩陣的塊數(shù)為3

Sub: do i_Sub = 1,j !從坐標(biāo)的個(gè)數(shù)j

Harmonic: doi_Hrm = 1,N_phi !諧波基函數(shù)的總個(gè)數(shù)N_phi

Col_1: do i_X = i_Bnd, N_a, 3 !節(jié)點(diǎn)的總個(gè)數(shù)N_a；

每相隔3個(gè)節(jié)點(diǎn)所對應(yīng)的自變量獲得微小增量

! 需獲得增量的自變量的編號(hào)kk和jj

kk=(i_Sub-1)*N_a*N_phi+(i_X-1)*N_phi+i_Hrm

jj = (i_Sub-1)*N_a*N_phi + (i_X-1)*N_phi + i_Hrm+ j*N_a*N_phi

! 將這兩個(gè)自變量添加微小的增量，并賦給中間變量：

Temp_X = X(kk) + h

Temp_Y = X(jj) + h

end do Col

!以上的‘Col循環(huán)’獲得了一組添加了增量的自變量，自變量之間的節(jié)點(diǎn)序號(hào)間隔為3

! 調(diào)用子函數(shù)FCN，輸入以上這組獲得增量的自變量向量，得到新的函數(shù)值向量Wa1, Wa2：

call FCN (Temp_X, Wa1, N)

call FCN (Temp_Y, Wa2, N)

Col_2: do i_X = i_Bnd, N_a，3

Row:do i_Row = 1 to N_a! 沿'行方向'遍歷，查找非零元素

if (i_Row >= i_X - 1 .and. iRow<= i_X + 1) then

! 判斷當(dāng)前的元素是否在非零元素帶的分塊矩陣內(nèi)

! 利用Wa1,Wa2,FF,h計(jì)算差商，近似當(dāng)前元素的偏導(dǎo)FJAC(略)

end if

end do Row

end do Col_2

end do Harmonic

end do Sub

end do Band

end subroutine JACdf_Sparse