陳嵩濤， 段慶林， 馬今偉

(大連理工大學(xué) 工業(yè)裝備結(jié)構(gòu)分析國家重點(diǎn)實(shí)驗(yàn)室，大連116023)

1 引 言

許多工程問題的力學(xué)分析需考慮結(jié)構(gòu)構(gòu)形的變化，如橡膠材料的大變形等。這一類問題的剛體轉(zhuǎn)動(dòng)和應(yīng)變很劇烈，小變形假設(shè)失效，需要采用幾何非線性模型進(jìn)行計(jì)算。

傳統(tǒng)有限元法的形函數(shù)依賴于網(wǎng)格單元，在處理幾何非線性問題時(shí)容易因網(wǎng)格畸變而導(dǎo)致計(jì)算失敗。與此不同，無網(wǎng)格法[1]的形函數(shù)僅依賴于離散點(diǎn)，無需網(wǎng)格單元，因而處理大變形問題更具優(yōu)勢。目前，無網(wǎng)格法已在幾何非線性的數(shù)值分析中得到廣泛應(yīng)用[2-4]。

與線性分析不同，求解幾何非線性問題時(shí)參考構(gòu)形的選取關(guān)系到分析效率與結(jié)果的精確性，有限元通常選擇相鄰構(gòu)形作為參考構(gòu)形，即UL法，相比于有限元，無網(wǎng)格法求解形函數(shù)導(dǎo)數(shù)過程較為復(fù)雜，使用UL法不斷更新構(gòu)形，必然會(huì)降低計(jì)算效率。針對該問題，需要建立既能及時(shí)跟蹤變形又不產(chǎn)生極大計(jì)算代價(jià)的構(gòu)形更新方式。Léger等[5]發(fā)展了一種介于TL法和UL法之間的構(gòu)形更新方法，即在每個(gè)載荷步迭代平衡后將此構(gòu)形作為新的參考構(gòu)形，這種方法比較適用于形函數(shù)復(fù)雜、更新構(gòu)形代價(jià)較大的無網(wǎng)格法。

由于無網(wǎng)格形函數(shù)的豐富性，Galerkin弱形式的準(zhǔn)確積分要求較多的區(qū)域積分點(diǎn)，嚴(yán)重降低了計(jì)算效率。無網(wǎng)格領(lǐng)域已發(fā)展了一些高效積分方法，如Wang等[6]的子域積分、Beissel等[7]的節(jié)點(diǎn)積分和Duan等[8]的穩(wěn)定應(yīng)力點(diǎn)積分等。其中,Chen等[9]發(fā)展的基于應(yīng)變光順的穩(wěn)定相容節(jié)點(diǎn)積分SCNI(Stabilized Conforming Nodal Integration)方法極具吸引力，其突出優(yōu)點(diǎn)是不引入任何人工參數(shù)，且能精確通過線性分片試驗(yàn)。該方法已成功應(yīng)用于幾何非線性分析[10,11]，在高度不規(guī)則計(jì)算網(wǎng)格下仍展現(xiàn)出良好的計(jì)算穩(wěn)定性和效率。Wang等[12]進(jìn)一步將該方法拓展到板單元的大變形問題，成功克服了剪切自鎖。

Puso等[13]發(fā)現(xiàn)SCNI仍然可能導(dǎo)致邊界附近鋸齒形振蕩模式的出現(xiàn)。Duan等[14]也指出，即使對于二次無網(wǎng)格近似，SCNI也只能通過線性分片試驗(yàn)，而不能通過二次分片試驗(yàn)。本質(zhì)上，SCNI是一點(diǎn)積分方案，只能反映出常數(shù)應(yīng)變場。對于應(yīng)變分布可能較為復(fù)雜的非線性問題，高階方法更能夠精確地再現(xiàn)變形過程，避免體積自鎖。Ortiz等[15]指出，在一些大應(yīng)變算例中低階方法魯棒性較差，對這一類問題發(fā)展高效的高階方法是必要的。Duan等[14]針對線彈性小變形問題發(fā)展了一個(gè)三點(diǎn)積分方案，積分點(diǎn)上的形函數(shù)導(dǎo)數(shù)由其與形函數(shù)之間散度定理的離散形式計(jì)算。該方法可精確通過二次分片試驗(yàn)，因而稱之為二階一致三點(diǎn)積分方法QC3(Quadratically Consistent 3-point integration method)。相應(yīng)的，SCNI可稱為線性一致一點(diǎn)積分方法LC1(Linearly Consistent 1-point integration method)。目前，QC3積分方案已成功應(yīng)用于不可壓縮固體[15]、斷裂相場[16]和斯托克斯流[17]等問題，但這些研究仍主要集中在線性分析方面。

本文把QC3積分方案拓展到幾何非線性分析，并進(jìn)一步引入Léger等[5]的更新構(gòu)形方法，發(fā)展能夠高效分析大變形問題的高階無網(wǎng)格法，并通過數(shù)值結(jié)果考察和驗(yàn)證方法的計(jì)算精度與效率。

2 幾何非線性分析列式

如圖1所示，考慮當(dāng)前構(gòu)形為Ωn + 1，坐標(biāo)為xn + 1，自然邊界為Γn + 1，本質(zhì)邊界為Γun + 1，Γun + 1∪Γtn + 1=Γn + 1，Γun + 1∩Γtn + 1= ?。與標(biāo)準(zhǔn)的TL和UL方法不同，本文選取第n個(gè)載荷步迭代平衡后的構(gòu)形Ωn作為參考構(gòu)形，該構(gòu)形上的各物理量均標(biāo)注下標(biāo)n，如圖1的Γtn和Γun。相應(yīng)地，初始構(gòu)形Ω0的各物理量均標(biāo)注下標(biāo)0。圖1的Fn和fn + 1分別表示Ωn相對于Ω0以及Ωn + 1相對于Ωn的變形梯度。

當(dāng)前構(gòu)形下，平衡方程及邊界條件為

(1)

(2)

σi jnj=tn + 1ionΓtn + 1

(3)

(4)

式中左端為內(nèi)力虛功Wint，右端為外力虛功Wext。采用本文選取的參考構(gòu)形Ωn，式(4)可表示為

(5)

式中bn和tn為定義在參考構(gòu)形下的體力和面力載荷，Sn和En分別為

(6)

(7)

圖1 計(jì)算構(gòu)形

式中Jn= |Fn|，I為二階單位張量，S和E分別為參考初始構(gòu)形的第二Piola-Kirchhoff應(yīng)力(PK2應(yīng)力)和格林應(yīng)變。顯然，Sn和En分別為前推到參考構(gòu)形Ωn的PK2應(yīng)力和格林應(yīng)變[5]。

采用節(jié)點(diǎn)的無網(wǎng)格形函數(shù)建立位移近似u=NIuI，并代入式(5)，則該式左端可表示為

(8)

(9)

(10)

(11)

(12)

(13)

(14)

(15)

(16)

式中m為節(jié)點(diǎn)數(shù)目，un=xn + 1-xn為參考構(gòu)形到當(dāng)前構(gòu)形的位移增量。須著重指出的是，本文方法采用Ωn作為計(jì)算參考構(gòu)形，因而式(12，14，16)的空間導(dǎo)數(shù)均為對xn的導(dǎo)數(shù)，這是與參考初始構(gòu)形的傳統(tǒng)TL法以及參考相鄰構(gòu)形的傳統(tǒng)UL法的核心區(qū)別。

非線性方程的求解采用牛頓-拉普森迭代算法，計(jì)算流程簡述如下。

(17)

(3) 建立切線剛度陣K，并通過罰函數(shù)法施加本質(zhì)邊界條件求解增量平衡方程

(18)

(19)

3 二階一致三點(diǎn)積分格式

由以上列式可以看到，剛度陣的建立需要計(jì)算形函數(shù)導(dǎo)數(shù)的區(qū)域積分。由于無網(wǎng)格形函數(shù)是非多項(xiàng)式的有理函數(shù)，精確計(jì)算該區(qū)域積分需要較多的積分點(diǎn)，故計(jì)算效率下降。針對該問題，本文引入Duan等[14]針對線彈性小變形問題發(fā)展的三點(diǎn)積分格式QC3來減少積分點(diǎn)數(shù)目。

如圖2所示，QC3方法采用背景三角形單元進(jìn)行積分計(jì)算，每個(gè)單元僅使用三個(gè)積分點(diǎn)(五角星)，邊界積分在每條邊上使用兩個(gè)一維高斯點(diǎn)(三角形)。QC3方法的核心思想是導(dǎo)數(shù)修正，即積分點(diǎn)上形函數(shù)的空間導(dǎo)數(shù)不采用標(biāo)準(zhǔn)形式，而是由散度定理(20)確定

(20)

式中ΩSn為參考構(gòu)形下的背景積分子域，ΓSn為其邊界，q(x) =p,x(x)∪p,y(x)，p(x)為MLS形函數(shù)的基底函數(shù)向量。本文采用二次近似，即p(x) =[1xyx2xyy2]T，則q(x) = [1xy]T，因而

圖2 QC3積分方法

式(20)實(shí)際上表示了三個(gè)方程，剛好可以用于確定三個(gè)積分點(diǎn)上形函數(shù)的空間導(dǎo)數(shù)。采用數(shù)值積分，式(20)的離散形式為

(21)

4 數(shù)值算例

采用三個(gè)算例考察和驗(yàn)證本文發(fā)展的無網(wǎng)格方法處理幾何非線性問題的有效性。作為對比，對有限元(線性三節(jié)點(diǎn)三角形單元，標(biāo)識為FEM)、LC1以及積分點(diǎn)上直接使用標(biāo)準(zhǔn)導(dǎo)數(shù)(而非修正導(dǎo)數(shù))的無網(wǎng)格法(標(biāo)識為ST16，即每個(gè)背景三角形單元使用16個(gè)積分點(diǎn))也進(jìn)行了程序?qū)崿F(xiàn)。

4.1 淺拱

如圖3所示，兩端固支弧形淺拱在頂部中心點(diǎn)A處受向下的集中載荷作用[10]。淺拱中性軸半徑為100 mm，拱厚為2 mm，弧對應(yīng)圓心角為 28.06°，楊氏模量E=4.8×103N/mm2，泊松比υ=0.0。無網(wǎng)格離散采用105個(gè)節(jié)點(diǎn)，分20個(gè)加載步，每步載荷為1 N，各種方法得到的A點(diǎn)處的載荷-位移響應(yīng)如圖4所示?？梢钥闯觯疚牡腝C3方法以及ST16方法的結(jié)果與解析解吻合良好，但是FEM和LC1方法的計(jì)算結(jié)果明顯偏離了解析解。應(yīng)說明的是，F(xiàn)EM的結(jié)果通過加密網(wǎng)格和改變構(gòu)形更新方式(如采用UL法)可以得到顯著改善。

圖3 受集中載荷的淺拱

圖4 淺拱A點(diǎn)處的載荷-位移響應(yīng)

QC3與ST16的計(jì)算精度相當(dāng)，但QC3更為高效，表1比較了三種無網(wǎng)格方法在四種計(jì)算節(jié)點(diǎn)配置下消耗的CPU時(shí)間。顯然，QC3的計(jì)算時(shí)間大約僅為ST16的1/4。

4.2 懸臂梁

如圖5所示，懸臂梁一端固支，另一端部中心A處受向上的集中載荷作用[10]，梁長為20 mm，厚度為1 mm，楊氏模量E=4.8×103N/mm2，泊松比υ=0.0。采用405個(gè)節(jié)點(diǎn)離散，分100個(gè)加載步加載，每步載荷為0.1 N，QC3方法得到的變形過程如圖6所示，載荷-位移響應(yīng)的比較如圖7所示。QC3和ST16的結(jié)果與解析解吻合，LC1的精度相對較差，F(xiàn)EM則明顯偏離解析解。

圖8比較了各方法得到的應(yīng)力場。FEM和LC1的應(yīng)力場均出現(xiàn)了明顯的虛假數(shù)值振蕩。本文QC3方法使用的積分點(diǎn)明顯少于ST16，但兩種方法得到了類似的光滑無振蕩的應(yīng)力場。

表1 淺拱算例的CPU時(shí)間比較

圖5 受集中載荷的懸臂梁

圖6 QC3方法得到的不同載荷下變形后的懸臂梁

4.3 軸承座

如圖9所示，軸承座尺寸為10 cm×10 cm，大孔位于正中心，四角小孔對稱分布，h=1 cm，r=0.5 cm，R=4.5 cm。楊氏模量E=7.9×104N/cm2，泊松比υ=0.33。軸承座的上下表面分別施加向下和向上的不隨構(gòu)形變化的均勻面載荷。無網(wǎng)格計(jì)算采用了3830個(gè)節(jié)點(diǎn),如圖9所示，分50個(gè)加載步，前25步每步加載2 N直至達(dá)到極值載荷50 N，后25步每步卸載2 N直至完全卸載。

圖7 懸臂梁A點(diǎn)處的載荷-位移響應(yīng)

圖8 懸臂梁的σx x場

圖9 軸承座的幾何模型和計(jì)算節(jié)點(diǎn)

圖10為P=20 N時(shí)軸承座的y向正應(yīng)力場。以Abaqus結(jié)果為參照，QC3得到的應(yīng)力極值較為準(zhǔn)確，而LC1明顯偏離參考解。圖11展示了QC3方法得到的軸承座加卸載的變形過程。在完全卸載后軸承座與初始構(gòu)形吻合良好，無虛假殘余變形，這進(jìn)一步驗(yàn)證了本方法計(jì)算復(fù)雜結(jié)構(gòu)大變形的能力。

圖10 軸承座的σy y場

圖11 QC3得到的軸承座的變形過程

5 結(jié) 論

本文將QC3方法拓展到了幾何非線性分析，并采用上一載荷步的收斂構(gòu)形作為計(jì)算參考構(gòu)形，針對超彈性材料建立了大變形分析的高階無網(wǎng)格法。與標(biāo)準(zhǔn)的高階無網(wǎng)格法(ST16)相比，本文方法大幅度減少了積分點(diǎn)數(shù)目，節(jié)省了大量的計(jì)算時(shí)間。與一點(diǎn)積分的高階無網(wǎng)格法(LC1)相比，本文方法具有高得多的計(jì)算精度。因此，本文方法有效提高了高階無網(wǎng)格法幾何非線性分析的計(jì)算效率，是超彈性材料大變形問題值得考慮的計(jì)算方案。