郭玉杰，吳晗浪，李 薇，李迎港，吳劍晨

(南京航空航天大學(xué)飛行器先進(jìn)設(shè)計(jì)技術(shù)國防重點(diǎn)學(xué)科實(shí)驗(yàn)室，南京 210016)

現(xiàn)代工業(yè)及軍事裝備結(jié)構(gòu)往往較為復(fù)雜，因而對(duì)設(shè)計(jì)與分析工具提出較高要求。采用傳統(tǒng)數(shù)值手段(如有限元法等)進(jìn)行結(jié)構(gòu)分析時(shí)，常需要細(xì)密的網(wǎng)格以獲得較高精度解，往往會(huì)導(dǎo)致維數(shù)災(zāi)難問題(curse of dimensionality)，特別對(duì)于需要反復(fù)迭代求解的復(fù)雜問題，比如非線性動(dòng)力學(xué)、結(jié)構(gòu)多學(xué)科優(yōu)化、多場(chǎng)耦合等仍過于耗時(shí)[1]。

模型降階(model order reduction)方法通過構(gòu)造全階模型(full order model)的低階近似模型有效降低了求解問題的維度，同時(shí)也保留了原階模型的主要信息從而保證了較高的計(jì)算精度。常見的模型降階方法包括模態(tài)疊加法[2]、模態(tài)導(dǎo)數(shù)法[3-4]以及正交分解法[5-6](proper orthogonal decomposition, POD)等，但對(duì)于非線性動(dòng)力學(xué)問題，上述常規(guī)方法還存在一定問題，比如模態(tài)疊加法需要不斷更新模態(tài)空間，POD 方法需要不斷更新全階模型的內(nèi)力、剛度矩陣等，這顯著增加了計(jì)算復(fù)雜度并降低了計(jì)算效率[7]。針對(duì)這一問題，Astrid等[8]提出了一種MPE (missing point estimation，缺失點(diǎn)估計(jì))方法，其通過計(jì)算在預(yù)先選定的空間點(diǎn)處的函數(shù)值，再利用gappy POD 方法將其投影到降維空間以實(shí)現(xiàn)計(jì)算效率的提升。Rewienski 等[9]提出了TPWLMOR (trajectory piecewise linear model order reduction，軌跡分片線性模型降階)方法，其在一系列采樣點(diǎn)上采用Krylov 投射對(duì)非線性函數(shù)進(jìn)行降階，再對(duì)降階模型加權(quán)求和從而可以實(shí)現(xiàn)非線性函數(shù)的逼近，但對(duì)于一些復(fù)雜非線性函數(shù)，采用低階分片多項(xiàng)式逼近效果不佳。Gaonkar等[10-11]在此基礎(chǔ)上提出了一種改進(jìn)的TPWLMOR方法，有效提升了降階模型在非線性動(dòng)力學(xué)計(jì)算方面的精度與效率。滿興博等[12]提出了一種非線性Galerkin 方法，有效提升了降階模型精度。Barrault[13]及Grepl 等[14]提出了一種經(jīng)驗(yàn)插值法(empirical interpolation method, EIM)，其采用一組經(jīng)驗(yàn)基函數(shù)的線性組合來近似非線性項(xiàng)，組合系數(shù)可通過插值求解得到，在此基礎(chǔ)上可通過greedy算法選取合適的snapshots 作為降階基函數(shù)構(gòu)造降階模型[15]，但EIM 不能擴(kuò)展到任意非線性函數(shù)。Chaturanabut 等[7]基于EIM 提出了一種離散形式的經(jīng)驗(yàn)插值方法(DEIM)，其采用POD 基底作為輸入并適用于任意類型的常微分方程，相比于EIM，其離散的變體應(yīng)用范圍更廣，且無需施加額外的基底變換，同時(shí)兩者的計(jì)算結(jié)果在數(shù)學(xué)上是等價(jià)的，因此受到了廣泛關(guān)注。

在參數(shù)化降階模型研究方面，傳統(tǒng)方法存在與非線性降階模型類似的困難，即難以實(shí)現(xiàn)模型的輸入?yún)?shù)與微分算子的解耦，從而對(duì)于每個(gè)未知輸入?yún)?shù)都需要重新計(jì)算系統(tǒng)的剛度矩陣等，計(jì)算效率不高[16]。針對(duì)這一問題，許多學(xué)者嘗試將有限元模型與降階模型進(jìn)行解耦，以實(shí)現(xiàn)線下/線上的計(jì)算模式(offline/online procedure)[17-19]。在線下模式中，其計(jì)算復(fù)雜度及計(jì)算量與所求問題的有限元模型相關(guān)，但可以預(yù)先計(jì)算并存儲(chǔ)結(jié)果，在線上模式中，其計(jì)算量以及計(jì)算復(fù)雜度只與降階模型有關(guān)，結(jié)合后驗(yàn)誤差估計(jì)方法，即可以實(shí)現(xiàn)實(shí)時(shí)、高精度的降階模型求解。當(dāng)所求解的模型與輸入?yún)?shù)具有仿射關(guān)系時(shí)，可以將輸入?yún)?shù)與微分算子分離，再結(jié)合降階模型方法即可實(shí)現(xiàn)線下/線上的計(jì)算模式。但通常情況下這種仿射關(guān)系并不存在或者難以建立。許多學(xué)者[7,13-14,20-21]通過POD 或者EIM/DEIM 方法實(shí)現(xiàn)了參數(shù)化模型的降階，從而可以進(jìn)行快速設(shè)計(jì)與優(yōu)化。此外，在模型的參數(shù)化描述方面(特別是幾何形狀)，還缺乏統(tǒng)一、高效的手段，并且對(duì)于不同的幾何形狀輸入?yún)?shù)，需要重復(fù)建立有限元模型，當(dāng)幾何模型較為復(fù)雜時(shí)，其耗費(fèi)的網(wǎng)格劃分等前處理時(shí)間則較為顯著[22-23]，造成這一問題的根本原因在于CAD 與CAE 中幾何模型的描述方式存在不一致[23]。在CAD 中，幾何模型一般是以非均勻有理B 樣條(NURBS)或者構(gòu)造實(shí)體幾何(CSG)的方式來描述，而在CAE 中(比如有限元)，其模型一般是以線性或二次插值的方式來構(gòu)造，是真實(shí)幾何的一種近似。

等幾何分析[23]的基本思想是將CAD 中描述幾何形狀的NURBS 引入到等參有限元中，通過統(tǒng)一的幾何描述可以消除產(chǎn)品設(shè)計(jì)過程中CAD 與CAE 之間反復(fù)的數(shù)據(jù)轉(zhuǎn)換過程，從而節(jié)省了大量的前處理時(shí)間。等幾何分析自2005 年提出以來迅速成為計(jì)算力學(xué)以及計(jì)算幾何領(lǐng)域的研究熱點(diǎn)[24-27]，相比于傳統(tǒng)有限元方法，等幾何分析其具有幾何精確、光順以及高階連續(xù)等優(yōu)點(diǎn)，因此特別適合于薄壁類結(jié)構(gòu)的分析、優(yōu)化等[28-38]。在等幾何降階模型研究方面，一些學(xué)者考慮了等幾何分析在模型幾何形狀描述、高效的CAD/CAE 模型轉(zhuǎn)換等方面的優(yōu)勢(shì)，將其嵌入到了降階模型中。宋敏[39]研究了帶參數(shù)橢圓形方程的等幾何分析及POD 模型降階問題，驗(yàn)證了模型降階的有效性。Manzoni 等[40]結(jié)合了等幾何邊界元以及降階模型實(shí)現(xiàn)了實(shí)時(shí)的計(jì)算流體動(dòng)力學(xué)分析，其對(duì)NACA 翼型采用了參數(shù)化描述并利用EIM 方法近似實(shí)現(xiàn)了參數(shù)的仿射關(guān)系。Salmoiraghi 等[22]結(jié)合了自由形狀變換(free form deformation，F(xiàn)FD)以及等幾何分析進(jìn)行參數(shù)化建模，實(shí)現(xiàn)了以較少參數(shù)描述復(fù)雜形狀，并且保證了設(shè)計(jì)模型與分析模型的無縫融合。Ding 等[41]采用等幾何分析、POD 以及DEIM 研究了結(jié)構(gòu)動(dòng)力學(xué)問題，其結(jié)合了兩種時(shí)間積分方法，有效抑制了高頻振蕩。Rinaldi[20]采用了EIM/DEIM 以及等幾何分析研究了薄殼結(jié)構(gòu)的參數(shù)化模型降階問題，結(jié)果顯示只需要少數(shù)DEIM 插值單元以及降階基函數(shù)即可實(shí)現(xiàn)高精度、快速的參數(shù)化線性靜力學(xué)分析。此外，一些學(xué)者還將等幾何降階模型應(yīng)用到了結(jié)構(gòu)優(yōu)化[42]、裂紋擴(kuò)展[43]、結(jié)構(gòu)屈曲[44]等方面。

綜上所述，將等幾何分析的思想與降階模型相結(jié)合具有廣闊的應(yīng)用前景，因此許多學(xué)者在這方面做了初步探索，但針對(duì)參數(shù)化的薄壁結(jié)構(gòu)的非線性動(dòng)力學(xué)降階模型方面的研究還未有涉及。本文擬選取平面曲梁結(jié)構(gòu)作為研究對(duì)象，采用等幾何分析、POD-DEIM 方法研究其參數(shù)化的非線性動(dòng)力學(xué)模型降階問題。

1 NURBS 及幾何描述

非均勻有理B 樣條(NURBS)是計(jì)算機(jī)輔助設(shè)計(jì)領(lǐng)域用來表示曲線曲面的標(biāo)準(zhǔn)函數(shù)之一，其能夠精確描述圓形等圓錐曲線，其定義為：

式中：p為NURBS 基函數(shù)的階數(shù)；n為基函數(shù)的個(gè)數(shù)；Ni,p(ξ) 為p階B 樣條函數(shù)；wi為對(duì)應(yīng)的權(quán)重。階數(shù)為p的B 樣條函數(shù)Ni,p(ξ)可由一給定的單調(diào)不減的節(jié)點(diǎn)序列 Ξ={ξ1,ξ2,···,ξn+p+1}并結(jié)合Cox-de-Boor 迭代公式求得：

由式(2)和式(3)可知，B 樣條函數(shù)Ni,p(ξ)只在區(qū)間 ξ ∈[ξi,ξi+p+1)上非零，即每個(gè)B 樣條函數(shù)均跨越p+1 個(gè)節(jié)點(diǎn)區(qū)間。同時(shí)，每個(gè)節(jié)點(diǎn)區(qū)間上均包含p+1 個(gè)非零的B 樣條函數(shù)。除此之外，所有B 樣條函數(shù)均為非負(fù)，滿足單位分解性質(zhì)。B 樣條函數(shù)在節(jié)點(diǎn)處(即單元交界處)具有Cp-k次連續(xù)的特性，其中k為節(jié)點(diǎn)的重復(fù)度，因此可以通過插入與刪除節(jié)點(diǎn)操作來調(diào)節(jié)單元間的連續(xù)性。當(dāng)節(jié)點(diǎn)矢量兩端的節(jié)點(diǎn)重復(fù)度均為p+1 次時(shí)，其被稱為開放型節(jié)點(diǎn)向量(open knot vector)。

一條p次NURBS 曲線的定義為：

圖1 NURBS 曲線及NURBS 基函數(shù)Fig. 1 NURBS curve and NURBS basis functions

二維NURBS 曲面以及三維實(shí)體可通過張量積的形式構(gòu)造，具體可參見文獻(xiàn)[45]。

2 等幾何平面曲梁非線性動(dòng)力學(xué)模型

2.1 平面曲梁運(yùn)動(dòng)學(xué)及本構(gòu)關(guān)系

本文所討論的曲梁模型基于Euler-Bernoulli 假設(shè)，即忽略剪切變形和橫截面轉(zhuǎn)動(dòng)的影響。本文中所有與初始構(gòu)型(reference configuration)相關(guān)的參數(shù)均用大寫字母表示，與變形狀態(tài)(deformed configuration)相關(guān)的參數(shù)均采用小寫字母表示。

圖2 給出了平面曲梁模型示意圖，其可由一中性軸和一系列橫截面所表示。圖中X以及x為曲梁中性軸的位置向量，且x=X+u，u為曲梁的位移向量。為了便于計(jì)算，需要在中性軸上定義一個(gè)活動(dòng)標(biāo)架 {a1,a2,a3} ， 其中向量a1表示沿中性軸的單位切向量，向量a2及a3是位于曲梁的橫截面上并且互相正交的單位向量。單位向量a1可由下式計(jì)算得到：

圖2 初始構(gòu)型與當(dāng)前構(gòu)型Fig. 2 Reference and current configurations

式中，ds為中性軸的弧長微分。不失一般性，選取a3為垂直于xy平面的單位向量，則a2可表示為a1×a3。

對(duì)于二維平面曲梁，其應(yīng)變關(guān)系可由薄膜應(yīng)變?chǔ)?以及彎曲應(yīng)變?chǔ)?表示，其表達(dá)式分別為：

其中，弧長微分以及曲率可分別表示為：

其中，上標(biāo) (·)′以及 (·)′′表示對(duì)參數(shù)坐標(biāo)的一次和二次導(dǎo)數(shù)。根據(jù)式(4)可得：

對(duì)于彈性問題，基于圣維南原理，薄膜力與薄膜應(yīng)變、彎矩與彎曲應(yīng)變有如下關(guān)系：

式中：E為彈性模量；A為橫截面積；I為截面慣性矩。

2.2 控制方程及等幾何離散

對(duì)于幾何非線性問題，平面曲梁的動(dòng)力學(xué)控制方程可表示為：

式中：S為曲梁中性軸的長度； δ為泛函變分；ρ0為曲梁的密度；N=N A1，M=MA3為向量化的薄膜力以及彎矩。此外，q0為分布力，p0為集中力。

基于等幾何分析的思想，采用描述中性軸的NURBS 基函數(shù)插值位移場(chǎng)函數(shù)及其變分：

式中：Ui為第i個(gè)控制點(diǎn)的位移未知量；R和U分別為NURBS 基函數(shù)和位移未知量的向量集合。

將式(17)代入式(15)，并利用 δUi的任意性，得到半離散的平衡方程：

式中，M0為平面曲梁的質(zhì)量矩陣，可表示為：

式(18)中，fI和fE分別表示曲梁的內(nèi)力和外力向量，其表達(dá)式為：

式中， ξ0為集中載荷施加位置所對(duì)應(yīng)的參數(shù)值。

2.3 非線性動(dòng)力學(xué)求解

平面曲梁的半離散平衡方程式(18)為時(shí)變、非線性方程，一般可通過時(shí)間積分并結(jié)合牛頓迭代進(jìn)行求解。為便于說明，式(18)可改寫為：

式中：t為當(dāng)前時(shí)間步；Δt為時(shí)間增量。

為了尋求t+Δt時(shí)刻的平衡狀態(tài)，可采用HHT-α?xí)r間積分方法，其需要定義一個(gè)動(dòng)態(tài)殘余力向量：

式中，α 為時(shí)間積分常數(shù)，可以調(diào)節(jié)系統(tǒng)的數(shù)值耗散能力，通常取 α=-0.05[46]。此外，內(nèi)力的下標(biāo)0 表示當(dāng)前時(shí)刻的初始狀態(tài)。與Newmark 方法類似，HHT-α 方法假設(shè)加速度在時(shí)間步內(nèi)線性變化，因此，可以推導(dǎo)得出以下關(guān)系：

根據(jù)式(25)可求得t+Δt時(shí)刻的加速度：

式中， ΔU=Ut+Δt-Ut為位移增量。

式中，KT為平面曲梁的剛度矩陣，可由內(nèi)力向量fI的線性化得到，具體表達(dá)式為：

當(dāng)把平面曲梁位移場(chǎng)U設(shè)置為當(dāng)前載荷步的初始位移時(shí)，即可得到KT,0以及fI,0的表達(dá)式。

基于式(27)可得到初始載荷增量 ΔU1，從而可以更新內(nèi)力向量、加速度及剛度矩陣等變量，并進(jìn)行迭代求解。圖3 所示為采用HHT-α 時(shí)間積分方法進(jìn)行幾何非線性動(dòng)力學(xué)分析的算法流程。

圖3 HHT-α 時(shí)間積分算法流程圖Fig. 3 HHT-α time integration scheme for nonlinear dynamic analysis

3 POD-DEIM/MDEIM 模型降階方法

3.1 模型降階簡介

假設(shè)不考慮平面曲梁的阻尼，在某一時(shí)刻系統(tǒng)平衡方程的離散形式可以表示為：

其中，KT,M0∈Rn×n,U,fE∈Rn×1。

其中，rROM為代入近似公式(32)產(chǎn)生的殘差向量，其位于正交基底所張空間之外且正交于，因此如下關(guān)系式成立：

由于降階后系統(tǒng)維度遠(yuǎn)小于原階模型的系統(tǒng)維度(k?n)，因此，可以提升系統(tǒng)求解速度，但降階基底的選擇對(duì)于降階模型的求解精度至關(guān)重要。

3.2 POD：降階模型基底生成

特征正交分解(POD)通過搜集不同時(shí)刻、不同參數(shù)下結(jié)構(gòu)的響應(yīng)快照(snapshots)，再通過奇異值分解(singular value decomposition, SVD)可獲得一系列相互正交的基底向量以及對(duì)應(yīng)的奇異值。

假設(shè)結(jié)構(gòu)響應(yīng)快照集合為：式中，NS為所搜集的快照向量個(gè)數(shù)。對(duì)Usnap采用SVD 分解，得到：

式中：V∈Rn×n為left-hand 矩陣，包含n個(gè)互相正交的單位向量 [v1,v2,···,vn]；W∈RNS×NS為right-hand矩陣，包含NS個(gè)正交向量；Σ=diag(σ1,σ2,···σn)∈Rn×NS，其包含n個(gè)奇異值σi=(i=1,2,···,n)，且σ1≥σ2≥···≥σn＞0 每個(gè)奇異值 σi的大小代表對(duì)應(yīng)的基底正交向量vi占整個(gè)系統(tǒng)能量的比重。

假設(shè)選取前k(k?n)個(gè)正交向量作為模型降階基底，則其對(duì)于快照集在最小二乘意義下的逼近誤差可表示為：即可以通過剩余的n-k個(gè)奇異值的平方和來估計(jì)模型降階基底的逼近誤差。隨著所選取的基函數(shù)個(gè)數(shù)k的增大，其逼近誤差迅速減小，表明基函數(shù)所捕獲的能量占快照集總能量的比例迅速提升，其所捕獲的原系統(tǒng)的信息就越豐富。

3.3 DEIM/MDEIM 差值算法

對(duì)于非線性、參數(shù)化問題，在迭代求解過程中，往往需要更新結(jié)構(gòu)的剛度矩陣以及內(nèi)力向量，以準(zhǔn)確捕捉結(jié)構(gòu)的響應(yīng)，這顯著影響了降階模型的求解效率。

離散經(jīng)驗(yàn)插值方法(DEIM)提供了一種有效降低非線性函數(shù)計(jì)算規(guī)模的方案，采用DEIM 方法，只需計(jì)算少數(shù)給定位置的元素，再通過DEIM插值矩陣即可實(shí)現(xiàn)非線性項(xiàng)的近似表示，大大降低了非線性項(xiàng)的計(jì)算時(shí)間。其基本原理如下。

假設(shè)曲梁的內(nèi)力向量fI所在的n維空間可由一 組 正 交 單 位 向 量Z=[ζ1,ζ2,···,ζn]∈Rn×n所 描述，且 ζi(i=1,2,···,n)按照重要性從大到小排序。與模態(tài)疊加法類似，內(nèi)力向量fI可由Z的前m(m?n)階基向量近似得到：

式 中：Zm=[ζ1,ζ2,···,ζm]∈Rn×m為 前m階 基 向量；c=[c1,c2,···,cm]∈Rm×1為系數(shù)向量。

式(40)為過約束方程組，因此需要選擇基底Z中的m行進(jìn)行求解，其可以通過布爾矩陣P實(shí)現(xiàn)：

其中，eφi∈Rn×1為單位向量，且在位置φi(1≤φi≤n)處為1，其余位置均為0，式(40)左乘PT可得：

繼而可以得到：

將式(43)代入式(40)可得：

其中，矩陣D為DEIM 插值矩陣。

由式(44)可知，內(nèi)力向量fI可以由布爾矩陣P、基函數(shù)Zm以及特定位置φ=[φ1,φ2,···,φm]∈Rm×1處的內(nèi)力向量本身近似得到，因此，在進(jìn)行有限元計(jì)算時(shí)，只需要計(jì)算特定位置φ=[φ1,φ2,···,φm]處的內(nèi)力值，再通過插值矩陣D即可重構(gòu)得到完整的內(nèi)力向量，大大節(jié)省了計(jì)算時(shí)間。

在DEIM算法中，基底向量Zm以及位置坐標(biāo)φ=[φ1,φ2,···,φm]的選取對(duì)于算法的精度及穩(wěn)定性至關(guān)重要，通常Zm可以選為由式(38)得到的POD 基底向量V，而位置坐標(biāo)φ=[φ1,φ2,···,φm]則可以通過算法1 實(shí)現(xiàn)。

算法1 中函數(shù) max{·}用于選取向量的最大元素并返回其位置坐標(biāo)，其通過選取最大誤差元素的位置并在下一次迭代中精確插值該位置元素，從而能夠?qū)崿F(xiàn)誤差的可控，因此，殘差向量rD在所選取的位置φ=[φ1,φ2,···,φm]處為0。

算法1. DEIM 算法

算法1 可以用于非線性向量的插值近似，而對(duì)于剛度矩陣，其排列方式為矩陣形式，因此可以將其按照列向量的形式進(jìn)行排列，再利用DEIM方法進(jìn)行插值，最后再重新按照矩陣形式進(jìn)行排列即可得到插值近似的剛度矩陣，該方法可以稱為 MDEIM(matrix discrete empirical interpolation method)。

值得注意的是，對(duì)于求解帶有阻尼的結(jié)構(gòu)動(dòng)力學(xué)問題，當(dāng)其阻尼矩陣為時(shí)變或者非線性變化時(shí)，也可通過上述POD-DEIM/MDEIM 方法進(jìn)行模型降階，首先搜集阻尼矩陣的快照集，進(jìn)行奇異值分解得到一組基底向量Zm，再采用DEIM/MDEIM方法生成插值位置坐標(biāo)向量 φ，即可構(gòu)造阻尼矩陣的插值矩陣D，進(jìn)而實(shí)現(xiàn)阻尼矩陣的快速計(jì)算。

當(dāng)獲得插值位置向量 φ之后，需要計(jì)算這些位置處的元素值，而每個(gè)位置元素均有與之相對(duì)應(yīng)的控制點(diǎn)，根據(jù)NURBS 幾何以及NURBS 基函數(shù)的性質(zhì)，每個(gè)控點(diǎn)均對(duì)應(yīng)一個(gè)跨(p+1)個(gè)單元的基函數(shù)，因此需要在不同單元上計(jì)算該基函數(shù)對(duì)于非線性項(xiàng)的貢獻(xiàn)，基于這一規(guī)律，即可以快速地計(jì)算出位置 φ處的值，從而得到整個(gè)非線性項(xiàng)的近似表達(dá)。

4 數(shù)值算例

4.1 兩端固支半圓弧曲梁

如圖4 所示為兩端固支半圓弧曲梁模型，其橫截面寬度與高度分別為b=h= 25.4 mm，曲梁半徑R= 1705 mm，彈性模量E= 68 975 MPa，密度為 ρ0=2.6086×10-9t/mm3。

圖4 兩端固支半圓弧曲梁模型Fig. 4 Half-circle arc with clamped boundary conditions at both ends

首先，研究等幾何曲梁單元的收斂特性。該曲梁中部受集中力載荷F= 100 N，分別采用2、3、4 次NURBS 單元進(jìn)行離散，并進(jìn)行幾何非線性靜力分析，其中點(diǎn)位移的收斂曲線如圖5 所示。作為對(duì)比，圖5 中列出了ABAQUS 參考值，其采用了包含剪切自由度的B21 鐵木辛柯梁單元，并且采用了縮減積分來控制剪切和薄膜閉鎖。對(duì)比結(jié)果顯示，當(dāng)采用3 階和4 階NURBS 單元時(shí)，等幾何模型收斂速度快于ABAQUS 梁單元；當(dāng)階數(shù)較低為2 階時(shí)，收斂速度較慢。主要原因?yàn)椋罕疚牟捎玫牡葞缀吻簡卧獮闊o剪切自由度的歐拉-伯努利梁單元，當(dāng)階數(shù)較低時(shí)，存在薄膜閉鎖問題，當(dāng)提升單元階數(shù)時(shí)，薄膜閉鎖得到顯著緩解，因而提升了收斂速度。

圖5 兩端固支半圓弧曲梁位移收斂曲線(幾何非線性)Fig. 5 Convergence studies of the half-circle arc beam(geometrically nonlinear)

對(duì)于非線性動(dòng)力學(xué)分析，分別采用2、3、4 次NURBS 單元，單元數(shù)量為34，其中部受集中力載荷F= 3115 N，加載時(shí)間歷程為0.07 s，時(shí)間步長為 7×10-4s，共100 時(shí)間步。在線下(offline)計(jì)算過程中，首先對(duì)曲梁結(jié)構(gòu)進(jìn)行動(dòng)力學(xué)分析，搜集每個(gè)時(shí)間步所對(duì)應(yīng)的位移向量與剛度矩陣向量，再對(duì)位移以及剛度矩陣向量快照集合Usnap和KTsnap采用奇異值分解，得到基底向量ZU、ZK以及對(duì)應(yīng)的奇異值 ΣU、 ΣK。圖6 所示為不同階次NURBS 單元的位移以及剛度矩陣向量所對(duì)應(yīng)的奇異值分布圖，可以看出階次變化對(duì)奇異值影響不大。

圖6 兩端固支半圓弧曲梁奇異值Fig. 6 Singular values of displacement snapshots and stiffness matrix snapshots for polynomial degrees 2, 3, 4 for clamped half-circle arc

選取k= 20 階位移基底向量作為模型降階基底，選取m= 15 階剛度矩陣基底向量作為MDEIM基底，根據(jù)誤差估計(jì)式(39)計(jì)算可得逼近誤差為10-6%量級(jí)，具有足夠的精度。線上(online)計(jì)算采用與線下計(jì)算相同的載荷步長、加載時(shí)間歷程及載荷大小，計(jì)算結(jié)果如圖7所示。由圖可見，當(dāng)NURBS 單元為2 階時(shí)，全階模型與降階模型的位移載荷曲線與ABAQUS 參考值有一定偏差，這主要是由于2 階NURBS 單元存在一定的薄膜閉鎖現(xiàn)象，使其剛度偏大，從而導(dǎo)致了偏差，當(dāng)NURBS 單元階數(shù)升高為3 階、4 階時(shí)，其位移載荷曲線與ABAQUS 參考值吻合良好，體現(xiàn)了升階對(duì)于消除等幾何曲梁單元閉鎖的效果，避免了構(gòu)造無閉鎖單元的麻煩。因此本算例采用3 階NURBS 基函數(shù)進(jìn)行后續(xù)計(jì)算。

圖7 兩端固支半圓弧曲梁載荷-位移曲線Fig. 7 Load-displacement response of half-circle arc beam

圖8、圖9 所示為采用全階以及降階模型所獲得的曲梁中點(diǎn)的運(yùn)動(dòng)軌跡對(duì)比圖以及兩者差值的分布圖。由于所選模型是左、右對(duì)稱的，因此其x方向位移較小，可忽略不計(jì)，觀察y方向位移可以發(fā)現(xiàn)，降階模型與原階模型位移范圍基本一致(如圖8 所示)，誤差范圍在-1.25 mm～0.75 mm(如圖9 所示)，無量綱化后約為10-4量級(jí)，表明降階模型具有較高的精度。

圖8 兩端固支半圓弧曲梁中點(diǎn)軌跡對(duì)比圖Fig. 8 Comparison of the paths traced by the mid-point of the half-circle arc beam obtained using the FOM and ROM

圖9 兩端固支半圓弧曲梁降階與全階模型位移差值圖Fig. 9 Difference of mid-point displacement between FOM and ROM for half-circle arc beam

降階模型的計(jì)算時(shí)間對(duì)比見表1。由表可見，當(dāng)只采用POD 模型降階方法時(shí)，其計(jì)算時(shí)間相比全階模型提升有限；當(dāng)只選取前2 階基底作為模型降階基底時(shí)，其計(jì)算時(shí)間加速率(時(shí)間比)僅為1.4，效率提升較為有限。這主要是由于POD 降階模型需要重復(fù)計(jì)算曲梁的剛度矩陣所導(dǎo)致。當(dāng)采用POD-MDEIM 降階模型時(shí)，其計(jì)算時(shí)間加速率為8.2，顯然，MDEIM 方法顯著降低了剛度矩陣的計(jì)算時(shí)間，從而大大提升了降階模型的求解效率。

表1 兩端固支半圓弧曲梁降階模型計(jì)算時(shí)間對(duì)比Table 1 Comparisons of computational time for reduced order models of half-circle arc beam

表2 所示為采用POD-MDEIM 方法所得到的剛度矩陣列向量的插值位置以及包含該位置的單元編號(hào)。由表可見，對(duì)于兩端固支半圓弧曲梁模型，只需計(jì)算15 個(gè)位置處的剛度矩陣值即可通過插值矩陣D對(duì)剛度矩陣列向量進(jìn)行重構(gòu)，但由于每個(gè)NURBS 基函數(shù)的支撐區(qū)域包含p+1 個(gè)節(jié)點(diǎn)區(qū)間，因此所涉及的單元編號(hào)數(shù)量相比插值位置的個(gè)數(shù)更多。

表2 兩端簡支曲梁剛度矩陣MDEIM 插值位置Table 2 MDEIM Positions for stiffness matrix of simply supported arc beam

圖10 所示為不同單元數(shù)量對(duì)降階模型計(jì)算效率的影響?？梢?，降階模型與全階模型的計(jì)算時(shí)間隨著單元數(shù)量的增大呈線性增長趨勢(shì)，體現(xiàn)了該降階模型對(duì)于不同模型的適應(yīng)性。

圖10 兩端固支半圓弧曲梁降階、全階模型計(jì)算時(shí)間對(duì)比Fig. 10 Comparison of computational time for different reduced order models of half-circle arc beam

4.2 一端固支阿基米德螺旋曲梁

如圖11 所示為一端固支的阿基米德螺旋曲梁模型，幾何描述為r=5+2θ/m,θ ∈[0,2π]，另一端點(diǎn) ξ=1 處受集中力F=80 N 載荷，橫截面寬度和高度分別為b=h=0.5 m，彈性模量E=1.2×106Pa，密度 ρ0=1 kg/m-3，曲梁模型 采 用3 次NURBS單元描述，單元數(shù)為40。

圖11 一端固支阿基米德螺旋梁Fig. 11 Archimedes spiral shaped beam with one end clamped

線下計(jì)算時(shí)，載荷時(shí)間歷程為10 s，共300 個(gè)時(shí)間步。對(duì)位移及剛度矩陣快照集合Usnap和KTsnap進(jìn)行奇異值分解，得到奇異值分布(圖12)。

圖12 一端固支阿基米德螺旋曲梁奇異值Fig. 12 Singular values of displacement snapshots and stiffness matrix snapshots for clamped Archimedes spiral beam

選取k= 23 階位移基底向量作為模型降階基底，選取前m= 30 階剛度矩陣基底向量作為MDEIM 基底。根據(jù)誤差估計(jì)公式(39)可知，所選取的位移以及剛度矩陣基底的逼近誤差分別為1.85×10-6% 和6.81×10-5%，因此具有足夠的精度。

線上(online)計(jì)算過程中采用與線下計(jì)算相同的載荷步長、加載時(shí)間歷程及載荷大小，計(jì)算結(jié)果如圖13 所示。降階模型與全階模型吻合良好，表明降階模型具有較高的精度。在計(jì)算時(shí)間上，全階模型計(jì)算時(shí)間為460.25 s，而POD-MDEIM降階模型計(jì)算時(shí)間為64 s，速度提升率為13.8。

圖13 一端固支阿基米德螺旋曲梁載荷位移曲線Fig. 13 Load-displacement curve of clamped Archimedes spiral arc beam

圖14 所示為阿基米德螺旋曲梁端點(diǎn)軌跡對(duì)比圖，圖15 所示為全階模型與降階模型曲梁在端點(diǎn)處的位移差值。由圖可見，降階模型與全階模型的位移誤差較小，x方向位移誤差約為10-5m，y方向位移誤差約為10-6m 的水平，體現(xiàn)了降階模型的計(jì)算精度。

圖14 一端固支阿基米德螺旋曲梁端點(diǎn)的軌跡對(duì)比圖Fig. 14 Comparison of the paths traced by the end-point of the Archimedes spiral beam obtained using the FOM and ROM

圖15 阿基米德螺旋曲梁端點(diǎn)降階與全階模型位移差值圖Fig. 15 Difference of end-point displacement between FOM and ROM for the Archimedes spiral beam

4.3 兩端簡支曲梁

如圖16 所示為一兩端簡支曲梁模型，中點(diǎn)處ξ=0.5 受集中力F=800 N載荷，其跨度L= 10 m，橫截面寬度和高度分別為b=h= 0.5 m，曲梁采用二次NURBS 曲線描述，單元數(shù)為1，因此具有3 個(gè)控制點(diǎn)，其中間控制點(diǎn)坐標(biāo)為P2=(5,4)m。曲梁 彈 性 模 量E=1.2×106Pa ， 密 度 ρ0=1 kg/m3，為了獲得較高計(jì)算精度，將曲梁細(xì)化為37 個(gè)單元，階數(shù)提升為3 階，基于等幾何分析的特點(diǎn)，對(duì)曲梁進(jìn)行細(xì)化不會(huì)改變其幾何形狀。

圖16 兩端簡支曲梁模型Fig. 16 Arc beam with simply supported boundary conditions at both ends

線下計(jì)算時(shí)，載荷時(shí)間歷程為1 s，時(shí)間步長設(shè)為0.01 s，共100 個(gè)時(shí)間步。對(duì)位移以及剛度矩陣快照集合Usnap和KTsnap進(jìn)行奇異值分解，得到奇異值分布(如圖17 所示)。

圖17 兩端簡支曲梁奇異值Fig. 17 Singular values of displacement snapshots and stiffness matrix snapshots for simply supported arc beam

選取k= 20 階位移基底向量作為模型降階基底，選取前m= 15 階剛度矩陣基底向量作為MDEIM 基底。根據(jù)誤差估計(jì)公式(39)可知，選取的位移及剛度矩陣基底的逼近誤差分別為1.91×10-8% 和1.34×10-7%，因此具有足夠的精度。通過MDEIM 插值算法可以得到插值位置向量 φ，繼而得到位置向量中每個(gè)元素在剛度矩陣中對(duì)應(yīng)的行和列，從而可以得到包含這些特定位置元素的單元編號(hào)(見表3)。

觀察表3 可知，采用MDEIM 算法，選取前15 階插值基底，需要計(jì)算15 個(gè)位置處的剛度矩陣值，涉及的單元編號(hào)為：1, 8～25, 30～37，共27 個(gè)單元，但在每個(gè)單元中，無需計(jì)算所有的元素值，只需計(jì)算給定位置處的元素，因此，可以顯著提升計(jì)算效率。

表3 兩端簡支曲梁剛度矩陣MDEIM 插值位置Table 3 MDEIM Positions for stiffness matrix of simply supported arc beam

線上(online)計(jì)算過程中采用與線下計(jì)算相同的載荷步長、加載時(shí)間歷程及載荷大小，計(jì)算結(jié)果如圖18 所示，降階模型與全階模型吻合良好，表明降階模型具有較高的精度。在計(jì)算時(shí)間上，全階模型計(jì)算時(shí)間為281.6 s，而POD-MDEIM 降階模型計(jì)算時(shí)間為33.2 s，速度提升率為8.5。

圖18 兩端簡支曲梁載荷-位移曲線Fig. 18 Load-displacement curve of simply supported arc beam

圖19 所示為兩端簡支曲梁中點(diǎn)軌跡對(duì)比圖，圖20 所示為全階模型與降階模型曲梁在中點(diǎn)處的位移差值。由圖可見，降階模型與全階模型的位移誤差較小，y方向位移誤差約為10-6m 的水平，體現(xiàn)了降階模型的計(jì)算精度。

圖19 兩端簡支曲梁中點(diǎn)的軌跡對(duì)比圖Fig. 19 Comparison of the paths traced by the mid-point of the arc beam obtained using the FOM and ROM

圖20 兩端簡支曲梁降階與全階模型位移差值圖Fig. 20 Difference of mid-point displacement between FOM and ROM for the arc beam

圖21 所示為兩端簡支曲梁全階模型與降階模型的頻率響應(yīng)對(duì)比圖，頻率響應(yīng)可以通過曲梁模型中點(diǎn)位移的快速傅里葉變換得到。由圖可見，在頻率范圍內(nèi)，降階模型與原階模型亦吻合良好。

圖21 兩端簡支曲梁頻率響應(yīng)對(duì)比圖Fig. 21 Comparison of frequency response obtained using FOM and ROM

保持載荷時(shí)間歷程及時(shí)間步長不變，將外載荷(單位：N)變?yōu)镕=-800sin(36t)，通過在線模式進(jìn)行非線性動(dòng)力學(xué)計(jì)算，其載荷位移曲線如圖22所示。由圖可見，變換載荷后，降階模型所得結(jié)果與原階模型吻合良好，拓展了降階模型的適用范圍。在計(jì)算時(shí)間上，全階模型所需計(jì)算時(shí)間為252.9 s，POD-MDEIM 降階模型計(jì)算時(shí)間為33.5 s，速度提升率為7.5。

圖22 兩端簡支曲梁載荷位移曲線( F=-800sin(36t))Fig. 22 Load-displacement curve of simply supported arc beam ( F=-800sin(36t))

圖23、圖24 分別為曲梁受外載荷F=-800sin(36t)時(shí)，其中點(diǎn)軌跡對(duì)比圖以及位移差值圖，可以觀察到，降階模型與全階模型在y方向的位移誤差約為10-5m 的水平。如圖25 所示為改變外載荷后，曲梁全階模型與降階模型頻率響應(yīng)對(duì)比圖，兩者吻合良好。從圖23～圖25 中可以觀察到，本文提出的降階模型對(duì)于變載荷情況也能得到高精度的解。

圖23 兩端簡支曲梁中點(diǎn)的軌跡對(duì)比圖( F=-800sin(36t))Fig. 23 Comparison of the paths traced by the mid-point of the arc beam obtained using the FOM and ROM(F=-800sin(36t))

圖24 兩端簡支曲梁降階與全階模型位移差值圖(F=-800sin(36t))Fig. 24 Difference of mid-point displacement between FOM and ROM for the arc beam ( F=-800sin(36t))

圖25 兩端簡支曲梁頻率響應(yīng)對(duì)比圖( F=-800sin(36t))Fig. 25 Comparison of frequency response obtained using FOM and ROM ( F=-800sin(36t))

4.4 參數(shù)化兩端簡支曲梁

圖26 所示為兩端簡支曲梁的參數(shù)化模型，曲梁模型參數(shù)與算例4.3 相同。其由二階NURBS曲線描述，單元數(shù)為1，通過改變中間控制點(diǎn)P2的位置坐標(biāo)可以實(shí)現(xiàn)曲梁幾何模型的參數(shù)化。假設(shè)中間控制點(diǎn)P2可在x∈[2,8] m，y∈[0,10] m范圍內(nèi)取任意值，由于對(duì)稱性，故只需考慮范圍x∈[5,8] m，y∈[0,10] m (如圖26 中虛線方框所示)。

圖26 兩端簡支曲梁的參數(shù)化模型Fig. 26 Parametric model of simply supported curved beams

曲梁模型中點(diǎn)處 (ξ=0.5) 受集中力載荷F=-800 N，為了獲得較高計(jì)算精度，將曲梁細(xì)化為31 個(gè)單元，階數(shù)提升為3 階。由于等幾何分析模型與幾何模型保持一致，故無需進(jìn)行額外的模型轉(zhuǎn)換等工作，并且，模型的細(xì)化并不會(huì)改變其幾何形狀，這對(duì)模型的參數(shù)化提供了便利。

為了準(zhǔn)確捕捉不同參數(shù)下結(jié)構(gòu)的位移載荷響應(yīng)特征，本算例采用了兩種不同抽樣方法，即相對(duì)均勻抽樣(抽樣方法1)以及拉丁超立方抽樣(抽樣方法2)，在范圍x∈[5,8] m，y∈[0,10] m 中生成控制點(diǎn)P2的坐標(biāo)樣本，如圖27 中綠色控制點(diǎn)所示，其中，相對(duì)均勻抽樣點(diǎn)數(shù)為101 個(gè)，拉丁超立方抽樣點(diǎn)數(shù)為100。

圖27 控制點(diǎn)P2 抽樣分布圖Fig. 27 Sampling positions of the control point P2

線下計(jì)算時(shí)，分別對(duì)每個(gè)抽樣控制點(diǎn)所對(duì)應(yīng)的曲梁模型進(jìn)行非線性動(dòng)力學(xué)分析，并搜集位移以及剛度矩陣快照集合Usnap和KTsnap。對(duì)兩種抽樣方法所獲得的快照集合分別進(jìn)行奇異值分解，得到奇異值如圖28 所示。

圖28 參數(shù)化曲梁奇異值Fig. 28 Singular values of displacement snapshots and stiffness matrix snapshots for parameterized arc beam

選取k= 20 階位移基底向量作為模型降階基底，選取前m= 50 階剛度矩陣基底向量作為MDEIM基底。根據(jù)誤差估計(jì)公式(39)可知，所選取的位移以及剛度矩陣基底的逼近誤差分別為：1.23×10-5% 和9.78×10-6% (抽樣法1)， 9.06×10-6%和7.37×10-6%(抽樣法2)，因此具有足夠的精度。

在線上計(jì)算模式中，為了驗(yàn)證POD-MDEIM模型降階方法的精度及可靠性，在范圍x∈[5,8] m，y∈[0,10]m 中任意選擇與抽樣點(diǎn)不相重合的控制點(diǎn)，建立曲梁模型，進(jìn)行非線性動(dòng)力學(xué)降階模型計(jì)算。選取控制點(diǎn)P2坐標(biāo)為(7.5,9.5) m，(5.7,7.2) m?；谏鲜鰞煞N抽樣方法所構(gòu)建的降階模型，進(jìn)行非線性動(dòng)力學(xué)分析，其位移載荷響應(yīng)如圖29 所示，可以觀察到，對(duì)于不同抽樣方法，POD-MDEIM降階模型所得結(jié)果與原階模型吻合良好，驗(yàn)證了所提出的降階模型的可靠性以及準(zhǔn)確性。

圖29 不同抽樣方法所對(duì)應(yīng)的位移載荷響應(yīng)Fig. 29 Load-displacement response of half-circle arc beam

圖30、圖31 所 示分別為P2=(7.5,9.5) m 及P2= (5.7,7.2) m 時(shí)，降階模型與全階模型在ξ=0.5處的軌跡對(duì)比圖以及位移差值圖。由于參數(shù)化曲梁不再是對(duì)稱模型，因此需要考慮其x方向位移，可以觀察到，對(duì)于不同抽樣方法，降階模型與全階模型在 ξ=0.5處的運(yùn)動(dòng)軌跡均吻合良好，且位移誤差約為10-3m 的水平。

圖30 參數(shù)化曲梁 ξ=0.5處的軌跡對(duì)比圖Fig. 30 Comparison of the paths traced at the position of ξ=0.5for the parameterized arc beam obtained using the FOM and ROM

圖31 參數(shù)化曲梁降階與全階模型位移差值圖(ξ=0.5)Fig. 31 Difference of the displacement between FOM and ROM for the parameterized arc beam (ξ=0.5)

圖32 所示為參數(shù)化曲梁降階模型與全階模型的頻率響應(yīng)對(duì)比圖?？梢杂^察到，對(duì)于不同取樣方法以及不同控制點(diǎn)位置，降階模型與全階模型在頻域內(nèi)吻合良好。

圖32 參數(shù)化曲梁頻率響應(yīng)對(duì)比圖Fig. 32 Comparison of frequency response obtained using FOM and ROM for the parameterized arc beam

為了研究本文所提出的降階模型對(duì)于變載荷以及變時(shí)間步長的適應(yīng)性，以下分2 種情況進(jìn)行討論：首先，保持載荷時(shí)間歷程(1 s)和時(shí)間步長(0.01 s)不變，外力載荷(N)變?yōu)镕=-800sin(48t)；其次，保持上述載荷時(shí)間歷程和外載荷大小(N)不變(F=-800sin(48t))，減小時(shí)間步長為0.005 s，選擇控制點(diǎn)P2坐標(biāo)為(6, 5.5) m。針對(duì)上述兩種情況分別對(duì)兩種抽樣方法所建立的降階模型進(jìn)行非線性動(dòng)力學(xué)計(jì)算，所得載荷時(shí)間歷程曲線如圖33和圖34 所示，從圖中可以觀察到，對(duì)于不同載荷以及時(shí)間步長，采用不同抽樣方法所建立的降階模型與原階模型吻合良好，體現(xiàn)了降階模型的穩(wěn)定性。

圖33 不同抽樣方法所對(duì)應(yīng)的位移載荷響應(yīng)(ξ=0.5)(F=-800sin(48t))Fig. 33 Comparisons of load-displacement histories of the parameterized arc beam at (ξ=0.5) ( F=-800sin(48t))

圖34 不同抽樣方法所對(duì)應(yīng)的位移載荷響應(yīng)(ξ=0.5)(F=-800sin(48t), 時(shí)間步長：0.005 s)Fig. 34 Comparisons of load-displacement histories of the parameterized arc beam at(ξ=0.5)(F=-800sin(48t), time step: 0.005 s)

圖35 及圖36 所示為變載荷、變時(shí)間步長情況下參數(shù)化曲梁在 (ξ=0.5)處軌跡對(duì)比圖以及位移差值圖。由圖可見，對(duì)于兩種不同抽樣方法，降階模型與全階模型均吻合良好，位移誤差約為10-3m的水平。圖37 所示為降階模型與全階模型頻率響應(yīng)對(duì)比圖，兩者吻合良好。圖35～圖37 體現(xiàn)了參數(shù)化降階模型對(duì)于不同載荷以及不同抽樣方法的適應(yīng)性。

圖35 參數(shù)化曲梁 (ξ=0.5)處的軌跡對(duì)比圖(F=-800sin(48t), 時(shí)間步長：0.005 s)Fig. 35 Comparison of the paths traced at the position of(ξ=0.5)for the parameterized arc beam obtained using the FOM and ROM ( F=-800sin(48t), time step: 0.005 s)

圖36 參數(shù)化曲梁降階與全階模型位移差值圖(ξ=0.5) (F=-800sin(48t), 時(shí)間步長：0.005 s)Fig. 36 Difference of the displacement between FOM and ROM for the parameterized arc beam(ξ=0.5)(F=-800sin(48t), time step: 0.005 s)

圖37 參數(shù)化曲梁頻率響應(yīng)對(duì)比圖(F=-800sin(48t), 時(shí)間步長：0.005 s)Fig. 37 Comparison of frequency response obtained using FOM and ROM for the parameterized arc beam(F=-800sin(48t), time step: 0.005 s)

降階模型的計(jì)算時(shí)間對(duì)比如表4 所示。由表可見，POD-MDEIM 降階模型對(duì)于計(jì)算效率的提升較為顯著，且對(duì)于不同抽樣方法、不同計(jì)算樣本，以及不同的外載荷，均體現(xiàn)出了較高的速度提升效果，體現(xiàn)了降階模型的穩(wěn)定性及適應(yīng)性。

表4 參數(shù)化曲梁降階模型計(jì)算時(shí)間對(duì)比Table 4 Comparisons of computational time for reduced order models of parameterized arc beam

5 結(jié)論

本文結(jié)合等幾何分析、POD 以及離散經(jīng)驗(yàn)插值方法研究了參數(shù)化、平面曲梁結(jié)構(gòu)的非線性動(dòng)力學(xué)模型降階問題。等幾何分析具有幾何精確、高階連續(xù)等特點(diǎn)，較適合于薄壁類結(jié)構(gòu)(比如曲梁)的參數(shù)化描述以及力學(xué)性能分析，且無需網(wǎng)格生成等前處理過程，簡化了參數(shù)化模型的非線性動(dòng)力學(xué)分析。

POD-DEIM 模型降階方法通過尋找非線性項(xiàng)的插值位置向量來構(gòu)造插值矩陣，從而在非線性動(dòng)力學(xué)計(jì)算時(shí)，只需計(jì)算少數(shù)選定位置處的元素，顯著降低了非線性項(xiàng)的計(jì)算量，大大提升了計(jì)算效率。

在參數(shù)化、非線性動(dòng)力學(xué)分析方面，本文構(gòu)造的降階模型方法對(duì)于不同抽樣方法顯示出了良好的適應(yīng)性，其對(duì)于參數(shù)取樣空間中的幾何模型，展示了良好的計(jì)算精度與效率。

此外，本文所提出的POD-DEIM 降階模型也適用于變載荷以及變載荷步長的情形，因此，具有良好的適應(yīng)性與穩(wěn)定性。

綜上所述，基于等幾何分析的POD-DEIM 模型降階方法有效提升了平面曲梁結(jié)構(gòu)的非線性動(dòng)力學(xué)計(jì)算效率，可以推廣到復(fù)雜曲面類薄殼結(jié)構(gòu)的非線性動(dòng)力學(xué)分析中。