李 鋼，呂志超，余丁浩

(大連理工大學(xué)海岸和近海工程國家重點(diǎn)實(shí)驗(yàn)室，遼寧，大連 116024)

板殼單元作為數(shù)值計(jì)算分析中常見的單元類型，被廣泛應(yīng)用于建筑和橋梁等結(jié)構(gòu)的數(shù)值模擬中[1―3]。多年來，國內(nèi)外研究學(xué)者對板殼單元的計(jì)算精度和非線性求解效率兩方面做了大量研究[4―5]。目前，板單元主要有兩種分析理論，一種是Kirchhoff板理論[6]，該理論假定撓度和轉(zhuǎn)角相互獨(dú)立且要求位移為C1類連續(xù)，主要適用于由薄板組成的工程結(jié)構(gòu)中[7-9]，但當(dāng)板厚增加時(shí)，該理論模擬結(jié)果的精度降低甚至?xí)?dǎo)致錯(cuò)誤的計(jì)算結(jié)果；另一種理論是Mindlin-Reissner中厚板理論[10-11]，該理論克服了Kirchhoff板理論的不足，在解決中厚板問題上具有足夠的精度，被廣泛應(yīng)用于結(jié)構(gòu)有限元非線性分析中[12-15]。針對殼單元，其理論通常包括三種：由平面膜單元和板彎曲單元組成的平板型殼元[16]、經(jīng)典的曲面型殼元[17]以及三維實(shí)體退化型殼元[18]，其中，曲面型殼元和退化型殼元的有限元列式比較復(fù)雜且收斂性較差，在實(shí)際應(yīng)用中難以實(shí)施。而平板型殼元具有構(gòu)造簡單、計(jì)算精度高等優(yōu)點(diǎn)，被廣泛應(yīng)用于實(shí)際工程中[19]。近年，隨著復(fù)合材料層合板殼結(jié)構(gòu)的出現(xiàn)，單一材料的板殼單元已經(jīng)無法滿足工程應(yīng)用需要[20-21]。因此，研究人員以復(fù)合材料力學(xué)原理和平板型殼元基本理論為基礎(chǔ)，開發(fā)了一種分層殼單元模型，該模型將不同的材料本構(gòu)行為聯(lián)系起來，在模擬實(shí)際構(gòu)件或者結(jié)構(gòu)復(fù)雜受力方面具有明顯優(yōu)勢。陸新征等[22]和張有佳等[23]將分層殼模型應(yīng)用于鋼筋混凝土剪力墻和核電預(yù)應(yīng)力混凝土安全殼的非線性分析中，研究結(jié)果均表明該模型具有良好的計(jì)算精度。

在板、殼單元模型的非線性計(jì)算方面，高效的數(shù)值求解方法同樣占據(jù)重要地位。Farhat和Wilson[24]、楊玉明[25]將基于預(yù)處理共軛梯度法的并行計(jì)算方法應(yīng)用于板殼結(jié)構(gòu)的非線性分析中，大大提高了板殼結(jié)構(gòu)的并行計(jì)算效率。Li、Carrera等[26]提出一種分層殼有限元模型自適應(yīng)方法，提高了分層殼結(jié)構(gòu)的有限元非線性數(shù)值分析效率。盡管學(xué)者們提出的數(shù)值求解方法在很大程度上改善了板、殼單元的計(jì)算效率，但在非線性分析時(shí)仍需要對切線剛度矩陣進(jìn)行合成和分解，該過程將消耗大量的計(jì)算資源。

隔離非線性有限元法[27]作為一種結(jié)構(gòu)非線性分析新方法，將結(jié)構(gòu)的非線性部分隔離開來，并將結(jié)構(gòu)的切線剛度矩陣的逆矩陣表示成初始彈性剛度矩陣的逆與相應(yīng)修正矩陣的和，從而避免了切線剛度矩陣的直接分解，實(shí)現(xiàn)了非線性問題的高效求解。目前，該方法已被應(yīng)用于平面單元、纖維梁單元的材料非線性分析中[28]。而分層殼單元作為建筑結(jié)構(gòu)數(shù)值模擬的重要組成部分，有必要建立基于隔離非線性有限元法的分層殼單元分析模型，對提高板殼單元非線性分析效率具有較大意義。本文根據(jù)隔離非線性有限元法和分層殼單元理論，將單元截面變形分解為線彈性變形和非線性變形，以高斯積分點(diǎn)作為非線性變形插值結(jié)點(diǎn)，建立了非線性變形場。依據(jù)虛功原理和積分點(diǎn)處的內(nèi)力平衡條件建立了基于隔離非線性有限元法的分層殼單元控制方程，采用Woodbury公式和組合近似法聯(lián)合求解控制方程?；跁r(shí)間復(fù)雜度理論對模型的計(jì)算效率進(jìn)行評價(jià)，結(jié)果表明：隔離非線性有限元法的分層殼單元模型相比于傳統(tǒng)有限元方法的模型計(jì)算效率顯著提高，數(shù)值算例驗(yàn)證了本文單元模型的準(zhǔn)確性與高效性。

1 基本理論

1.1 分層殼單元

分層殼單元基于復(fù)合材料力學(xué)原理，將一個(gè)殼單元沿厚度方向上劃分成若干層，各層可根據(jù)需要設(shè)置不同的厚度和材料，每一層由平板型殼單元組成，如圖1所示。分層殼單元考慮了面內(nèi)應(yīng)變—面外彎曲之間的耦合作用，能較全面地反映殼體結(jié)構(gòu)的空間力學(xué)性能[29]。

圖1 分層殼單元示意圖Fig.1 Sketch of multi-layer shell elements

以常用的四節(jié)點(diǎn)平板型分層殼單元為例，對其基本理論和公式進(jìn)行闡述，如圖2所示。單元中面設(shè)置4個(gè)位移插值結(jié)點(diǎn)，每個(gè)位移插值結(jié)點(diǎn)有6個(gè)自由度，其結(jié)點(diǎn)位移增量可用列陣表示為[30-31]：

式中：Δqi={ΔuiΔviΔwiΔθxiΔθyiΔθzi}T，其中，Δui、Δvi、Δwi分別為沿x、y及z軸方向的平動(dòng)位移增量，Δθxi、Δθyi、Δθzi分別為繞x、y及z軸的轉(zhuǎn)角位移增量。平板型分層殼單元面內(nèi)和面外的變形行為可分別通過平面膜單元和板彎曲單元來進(jìn)行模擬，其中平面膜單元提供面內(nèi)位移以及面內(nèi)轉(zhuǎn)動(dòng)的分量，即板彎曲單元?jiǎng)t提供面外撓度及繞x、y軸轉(zhuǎn)角的分量，其節(jié)點(diǎn)位移分量表示為

圖2 四結(jié)點(diǎn)平板型分層殼單元Fig.2 Sketch of flat multi-layer shell element with four nodes

依據(jù)單元沿厚度方向的平截面假定，通過單元中面的任意一點(diǎn)的位移和橫截面上的轉(zhuǎn)角，可得單元內(nèi)任意一點(diǎn)的位移場增量方程，具體可表示為：

式中：Δu0、Δv0、Δw0為單元中面任意一點(diǎn)的位移增量；z為分層殼單元任意層與中面的距離。由彈性力學(xué)的幾何方程可知，對式(2)求一階偏導(dǎo)數(shù)可得到單元中面的應(yīng)變和曲率增量，即：

式中：Δεx、Δεy分別為x和y方向應(yīng)變增量；Δεxy為面內(nèi)剪切應(yīng)變增量；Δχx、Δχy分別為分層殼單元中面的x和y方向曲率增量；Δχxy為中面扭轉(zhuǎn)率增量。為下文表述方便，令Δdm={ΔεxΔεyΔεxy}T為膜單元變形場，Δdb={ΔχxΔχyΔχxy}T為板元彎曲變形場。

得到單元中面的應(yīng)變和曲率后，根據(jù)各層材料之間滿足平截面假定和各層在橫截面上的坐標(biāo)，可計(jì)算出各層的應(yīng)變，分層殼單元第j層的應(yīng)變增量可表示為：

在面外荷載作用下，分層殼單元在厚度方向上通常發(fā)生面外剪切變形，而單元厚度變薄將會(huì)發(fā)生剪切閉鎖現(xiàn)象，為解決該問題，學(xué)者采用假設(shè)剪應(yīng)變法計(jì)算面外剪切應(yīng)變，其基本原理是：單元內(nèi)的剪切應(yīng)變由插值點(diǎn)處的剪切應(yīng)變按線性插值表示，則單元的剪切應(yīng)變增量Δγxz、Δγyz可具體表示為[32]：

設(shè)單元截面變形增量分布函數(shù)為：

式中：Δd={ΔdmΔdb}T為截面變形向量；B=[BmBb]T為單元的變形矩陣；Bm和Bb分別為膜單元和板彎曲單元變形矩陣。

2.2 截面變形分解

隔離非線性方法作為一種高效的非線性求解方法[27]，是將材料本構(gòu)關(guān)系中的應(yīng)變分解為線彈性和非線性兩部分，如圖3所示。

圖3 應(yīng)變分解Fig.3 Strain decomposition of material

當(dāng)采用分層殼單元計(jì)算材料非線性問題時(shí)，可基于隔離非線性的思想，將分層殼單元截面變形分解為線彈性及非線性兩部分。由于面外剪切變形對單元非線性行為影響較小，在此不考慮面外剪切變形的非線性，認(rèn)為其一直處于彈性狀態(tài)，單元的膜單元變形和板元彎曲變形可分解為：

式中：Δde={Δdm,eΔdb,e}T為截面線彈性變形增量；Δdp={Δdm,pΔdb,p}T為截面非線性變形增量。

在截面內(nèi)力計(jì)算方面，分層殼第j層的應(yīng)力增量可通過該層的剛度矩陣和應(yīng)變增量相乘得到，對截面中各層的應(yīng)力增量進(jìn)行積分，得到膜內(nèi)力增量和彎矩增量，具體表達(dá)式為：

由圖3可知，材料的應(yīng)力等于彈性應(yīng)變和彈性模量的乘積，而分層殼單元中截面的內(nèi)力可通過截面彈性變形與截面彈性剛度矩陣相乘得到，結(jié)合式(8)截面變形的分解，截面內(nèi)力的表達(dá)式為：

式中：Δm={ΔNΔM}T為截面內(nèi)力；De=為彈性剛度矩陣；為膜單元應(yīng)力狀態(tài)下的第j層彈性剛度矩陣。

由式(11)可知，基于截面變形分解的方法求解截面內(nèi)力時(shí)，引入了未知量Δdp，依據(jù)有限元基本理論：單元的變形場增量Δd可通過變形插值函數(shù)和結(jié)點(diǎn)位移增量表示，同樣Δdp可通過非線性變形插值結(jié)點(diǎn)處的非線性變形增量Δdp進(jìn)行插值得到[28]，即：

在傳統(tǒng)變剛度有限元方法的非線性求解迭代過程中，截面內(nèi)力增量線性化表達(dá)式還可表示為[31]：

為切線剛度矩陣，為膜單元應(yīng)力狀態(tài)下的第j層的切線剛度矩陣。對比式(11)和式(13)可知：在非線性分析過程中，基于變形分解的截面內(nèi)力表達(dá)式可以保證截面剛度矩陣的始終不變。

3 隔離非線性分層殼控制方程

根據(jù)虛功原理，分層殼單元變分表達(dá)式為：

式中：δd為截面變形的變分；δq為單元節(jié)點(diǎn)位移的變分；Δf為單元節(jié)點(diǎn)力增量。將式(7)、式(8)、式(11)和式(12)代入虛功方程式(14)，經(jīng)整理可得：

此外，考慮塑性插值點(diǎn)處的內(nèi)力平衡條件，將式(11)和式(12)代入式(13)可得補(bǔ)充方程：

將式(16)代入式(15)，經(jīng)整理可得基于隔離非線性方法的分層殼單元控制方程：

式中：ke為單元初始彈性剛度矩陣；kmb表示非線性變形與單元節(jié)點(diǎn)力之間的關(guān)系；krr為單元的非線性矩陣，表示了單元的非線性信息。其具體表達(dá)式分別如下：

對單元控制方程集成，得到結(jié)構(gòu)的整體控制方程：

式中：ΔQ為結(jié)構(gòu)的位移向量；ΔF為結(jié)構(gòu)的荷載向量；為結(jié)構(gòu)中所有塑性插值點(diǎn)處的非線性變形向量。矩陣Krr為分塊對角矩陣，每個(gè)塊矩陣表示了每個(gè)高斯積分點(diǎn)的材料非線性信息，若該高斯積分點(diǎn)沒有進(jìn)入非線性，則該塊矩陣為零矩陣。在非線性分析過程中，只更新和分解矩陣Krr，矩陣Ke和Kmb均為常數(shù)矩陣且與材料非線性狀態(tài)無關(guān)。

4 控制方程求解

在結(jié)構(gòu)非線性分析任意迭代步內(nèi)，可利用Woodbury公式對結(jié)構(gòu)的整體控制方程式(18)進(jìn)行準(zhǔn)確和高效的求解。其求解高效性主要體現(xiàn)為：在局部非線性分析過程中，結(jié)構(gòu)中的大部分單元一般處于線彈性狀態(tài)，僅有小部分單元進(jìn)入非線性，因而，進(jìn)入非線性的塑性自由度數(shù)p遠(yuǎn)小于結(jié)構(gòu)位移自由度數(shù)n(p＜＜n)，在迭代步內(nèi)僅需對p維的矩陣合成和分解，避免了對n維整體剛度矩陣實(shí)時(shí)更新和分解，降低了結(jié)構(gòu)非線性分析的計(jì)算時(shí)間。而傳統(tǒng)變剛度有限元分析中，大部分計(jì)算時(shí)間消耗在n維切線剛度矩陣的分解上。

基于此，采用Woodbury公式求解整體控制方程式(18)，其展開形式為：

在數(shù)學(xué)上，矩陣Kpp稱為整體剛度矩陣Ke的舒爾補(bǔ)，其矩陣維數(shù)為p×p階，與塑性自由度數(shù)目一致[27]。從式(19)可知：Woodbury公式求解整體控制方程時(shí)，其計(jì)算量主要來自舒爾補(bǔ)矩陣Kpp的實(shí)時(shí)更新和分解，因而，非線性分析過程中產(chǎn)生塑性自由度數(shù)的多少直接影響著Woodbury公式的計(jì)算效率。以往的研究表明[33-34]：當(dāng)結(jié)構(gòu)局部進(jìn)入材料非線性時(shí)，Woodbury公式求解效率比傳統(tǒng)有限元方法求解效率高，對于一個(gè)位移自由度n=10000的結(jié)構(gòu)，當(dāng)進(jìn)入塑性的塑性自由度數(shù)p≥1000時(shí)，Woodbury公式求解效率低于傳統(tǒng)有限元方法。

4.1 Woodbury公式和CA法聯(lián)合求解

在板殼單元的非線性分析中，由于單元和其周圍單元聯(lián)系緊密，單元與單元之間的高斯積分點(diǎn)相互之間影響較大，即使在局部材料非線性階段，進(jìn)入非線性的塑性自由度數(shù)也較多，所以使用Woodbury公式求解整體控制方程在板殼單元非線性分析中存在一定的局限性。組合近似法(combined approximations method， CA法)通過構(gòu)造s個(gè)線性無關(guān)的基向量來線性逼近位移向量[35]，由于基向量個(gè)數(shù)s比結(jié)構(gòu)位移自由度數(shù)n小很多，所以大幅度的減小了非線性分析過程中的計(jì)算量。把式(19)的求解過程具體分解為從右至左的6個(gè)計(jì)算步驟，其計(jì)算過程如圖4(a)所示，由Woodbury公式求解分析可知，較大維數(shù)的矩陣Kpp直接分解計(jì)算是影響基于隔離非線性有限元法的板殼單元求解效率的關(guān)鍵因素，當(dāng)采用CA法求解圖4(a)中的步驟③時(shí)，在非線性求解過程中主要是矩陣回代與向量之間的運(yùn)算，避免了矩陣Kpp的直接分解計(jì)算，提高了基于隔離非線性有限元法的板殼單元的計(jì)算效率。CA法求解步驟③的具體計(jì)算過程如圖4(b)所示，圖4表示了結(jié)構(gòu)非線性分析過程中一個(gè)迭代步內(nèi)的完整計(jì)算。

圖4 Woodbury公式和CA法聯(lián)合求解過程示意圖Fig.4 Sketch of Woodbury formula and CA method solving process

4.2 效率分析

時(shí)間復(fù)雜度理論[33]是一種算法計(jì)算效率定量評價(jià)的方法，僅與算法本身相關(guān)，與硬件、軟件、編程語言以及程序優(yōu)化等因素?zé)o關(guān)。傳統(tǒng)有限元方法在增量步內(nèi)的每個(gè)迭代步的計(jì)算過程中都需要對規(guī)模為n×n階、帶寬為m階的整體剛度矩陣進(jìn)行一次LDLT分解，一次分解和回代的時(shí)間復(fù)雜度函數(shù)為[34]：

Woodbury公式和CA法聯(lián)合求解整體控制方程時(shí)，只需要在開始計(jì)算前對規(guī)模為n×n階、帶寬為m階的剛度矩陣Ke進(jìn)行一次LDLT分解。根據(jù)圖4所示的求解過程，通過對Woodbury公式和CA法聯(lián)合求解整體控制方程的計(jì)算量進(jìn)行統(tǒng)計(jì)，其時(shí)間復(fù)雜度函數(shù)約為[36]：

式中：m為剛度矩陣的帶寬，一般可取為s為基向量個(gè)數(shù)；p為進(jìn)入塑性的塑性自由度數(shù)。

5 數(shù)值算例

5.1 空心板

板的幾何尺寸及構(gòu)造如圖5所示，邊長為2000 mm×2000 mm，厚度為30 mm，板的四周為固定支座。將板模型劃分為400個(gè)單元，在厚度方向上每個(gè)單元?jiǎng)澐?0層，其中在單元中面的上、下兩層為空心層。板上承受均布荷載q=6.25 N/mm2，總計(jì)1251個(gè)增量步，力加載曲線如圖6所示。彈性模量為E=2×105MPa，泊松比為ν=0.3，屈服強(qiáng)度為262 MPa，屈服后剛度系數(shù)為0.1。模型位移自由度總數(shù)為1323，每個(gè)單元有12個(gè)塑性自由度，塑性自由度總數(shù)為4800。

分別采用有限元軟件ANSYS和隔離非線性有限元法建立分析模型，兩種方法的單元數(shù)量、截面劃分的層數(shù)和材料本構(gòu)關(guān)系均相同，其中，ANSYS模型中的單元類型為shell181單元，與隔離非線性有限元法中使用的分層殼單元類型相同。圖7(a)和圖7(b)分別給出了板中心點(diǎn)A的撓度-荷載曲線和第80個(gè)單元的1節(jié)點(diǎn)x方向的應(yīng)力-應(yīng)變曲線對比圖，可知：隔離非線性有限元法的計(jì)算結(jié)果與ANSYS的計(jì)算結(jié)果吻合較好，驗(yàn)證了隔離非線性分層殼有限元法的準(zhǔn)確性。

圖5 模型尺寸及分層 /mmFig.5 Size and layers of the model structure

圖6 力加載曲線Fig.6 Force loading curve

圖8給出了加載過程中板進(jìn)入非線性狀態(tài)的塑性自由度數(shù)量演變曲線，其中最大的塑性自由度數(shù)是4764，平均塑性自由度數(shù)是2498(位移自由度總數(shù)的187.96%)，表明了結(jié)構(gòu)發(fā)生了很大范圍的材料非線性變形，同時(shí)，板中進(jìn)入非線性狀態(tài)的塑性自由度數(shù)隨著外荷載和加卸載條件的改變而動(dòng)態(tài)變化。傳統(tǒng)變剛度有限元方法和隔離非線性有限元法的時(shí)間復(fù)雜度曲線如圖9所示，可知：有限元方法的時(shí)間復(fù)雜度保持不變，且僅與結(jié)構(gòu)位移自由度有關(guān)，而隔離非線性有限元法的時(shí)間復(fù)雜度與結(jié)構(gòu)自由度及產(chǎn)生的塑性自由度數(shù)均相關(guān)。對于大范圍的材料非線性問題，即板中所有單元全部進(jìn)入非線性狀態(tài)，隔離非線性有限元法的時(shí)間復(fù)雜度也是低于傳統(tǒng)有限元方法的，表明本文方法提高了分層殼單元在非線性分析過程中的數(shù)值計(jì)算效率。

圖7 結(jié)果曲線對比Fig.7 The comparison curve of the results

圖8 塑性自由度Fig.8 Plastic degree of freedoms

5.2 剪力墻

圖10為某一單片鋼板剪力墻，墻底部與基礎(chǔ)固結(jié)，墻的寬度為1000 mm、高度為2000 mm，將墻劃分為40個(gè)單元，每個(gè)單元截面在厚度方向上為100 mm并劃分為10層。在節(jié)點(diǎn)16施加x、y和z三個(gè)方向的荷載，大小均為Fx=Fy=Fz=336000 N，力控制加載，總計(jì)500個(gè)增量步。材料本構(gòu)關(guān)系采用理想彈塑性本構(gòu)模型，屈服強(qiáng)度為235 MPa，彈性模量為E=2.1×105MPa，泊松比為ν= 0.3。塑性自由度個(gè)數(shù)為960，位移自由度個(gè)數(shù)為330。

圖9 時(shí)間復(fù)雜度理論Fig.9 The comparison of time complexity theory

圖10 剪力墻模型Fig.10 The model of shear wall

分別采用有限元軟件ANSYS和隔離非線性有限元法建立分析模型。圖11給出了使用隔離非線性有限元法和傳統(tǒng)有限元方法計(jì)算得到的加載點(diǎn)A的x、y和z三個(gè)方向的位移-荷載(F-u)曲線，從圖中可知三個(gè)方向的位移吻合良好，因?yàn)榧袅Φ拿鎯?nèi)剛度大于面外剛度，所以z方向的位移進(jìn)入非線性程度較強(qiáng)，x、y方向的位移進(jìn)入非線性程度較弱，表明：剪力墻三個(gè)方向相互之間的耦合受力減低了剪力墻的承載能力。圖12給出了使用兩種方法得到的單元11中3節(jié)點(diǎn)的應(yīng)變和應(yīng)力曲線對比圖，可知：單元的應(yīng)變和應(yīng)力變化趨勢及計(jì)算精度與傳統(tǒng)有限元方法的基本相同，驗(yàn)證了本文提出分層殼單元模型的準(zhǔn)確性。

圖11 位移曲線Fig.11 The comparison of displacement curve

圖13給出了塑性自由度數(shù)隨增量步的變化曲線，非線性分析過程中產(chǎn)生的最大塑性自由度數(shù)為252，占總自由度數(shù)的76.36%，而每個(gè)增量步的平均塑性自由度數(shù)為129，占總自由度數(shù)的39.1%，表明結(jié)構(gòu)較大范圍的出現(xiàn)材料非線性。為論證其高效性，將傳統(tǒng)有限元方法和本文方法的時(shí)間復(fù)雜度進(jìn)行統(tǒng)計(jì)對比，如圖14所示，傳統(tǒng)有限元法和本文方法的平均時(shí)間復(fù)雜度分別為5.588×105和1.004×105，此外，傳統(tǒng)有限元法的時(shí)間復(fù)雜度為本文方法的5.57倍，表明本文方法能夠高效地求解結(jié)構(gòu)復(fù)雜受力狀態(tài)下的大范圍材料非線性問題。

圖12 應(yīng)變-應(yīng)力曲線Fig.12 The comparison of stress-strain

圖13 塑性自由度Fig.13 Plastic freedoms degree

圖14 時(shí)間復(fù)雜度理論Fig.14 The curve of time complexity theory

6 結(jié)論

本文基于隔離非線性有限元法和分層殼單元的基本理論，將單元截面變形分解為線彈性及非線性兩部分，推導(dǎo)了分層殼單元的隔離非線性控制方程，聯(lián)合Woodbury公式和CA法對控制方程進(jìn)行高效求解。最后，將該方法應(yīng)用于板、剪力墻結(jié)構(gòu)與構(gòu)件的非線性分析中，并得到如下結(jié)論：

(1) 與傳統(tǒng)的分層殼單元有限元模型相比，本文方法可保證單元初始剛度矩陣保持不變，將單元材料非線性集中于控制方程的右下角矩陣塊中，實(shí)現(xiàn)了單元非線性性狀態(tài)的隔離。

(2) 依據(jù)時(shí)間復(fù)雜度理論分析可知，隔離非線性有限元法的分層殼單元模型與傳統(tǒng)有限元方法相比，即使結(jié)構(gòu)大范圍出現(xiàn)材料非線性行為時(shí)，本文方法也可顯著提高分層殼單元的非線性求解效率。

(3) 通過與有限元軟件ANSYS的結(jié)果對比分析表明，本文方法與傳統(tǒng)有限元方法在求解分層殼單元的材料非線性問題時(shí)具有一致的計(jì)算精度。