二次有限元格式高階數(shù)值模擬矢量型彈性方程stress邊值問題

劉 雪, 丁 曉, 江 山

(南通大學(xué)理學(xué)院, 江蘇 南通 226019)

彈性力學(xué)主要研究彈性體受外界作用、環(huán)境變化或邊界約束等因素影響所導(dǎo)致的應(yīng)力、應(yīng)變和位移變化情況,被廣泛應(yīng)用于道路橋梁、水利工程和航空航天等領(lǐng)域. 近年來, 處理彈性力學(xué)方程的數(shù)值方法主要有有限元法[1-3]、混合有限元法[4-6]、弱有限元法[7-9]和多尺度方法[10-12]等, 其中有限元法因具有諸多優(yōu)點(diǎn)而備受關(guān)注: 可針對(duì)不同方程建立不同的變分形式; 等參單元與實(shí)際單元之間易于實(shí)現(xiàn)仿射變換;能通過獨(dú)立構(gòu)造分片基函數(shù)有效耦合局部與全局的數(shù)據(jù)信息; 數(shù)學(xué)理論完備等. 然而, 已有的工作主要采用線性基函數(shù)構(gòu)造多項(xiàng)式, 且多為處理彈性力學(xué)方程的本征邊界類型,抑或是格式較為復(fù)雜,計(jì)算成本偏高. 本文擬采用有限元法的高次基函數(shù)格式, 精確模擬帶有復(fù)雜stress邊界條件的矢量型彈性方程,通過法向量和切向量有效刻畫對(duì)應(yīng)的分量, 在更精細(xì)的數(shù)據(jù)結(jié)構(gòu)下組裝成總剛度矩陣并求解代數(shù)方程組, 以期利用有限的計(jì)算資源得到高階精度且穩(wěn)定收斂的數(shù)值模擬.

1 彈性方程及其變分

考慮二維彈性力學(xué)方程

(1)

(2)

(3)

2 有限元基函數(shù)的構(gòu)建

B(uh,vh)=(f,vh), ?vh∈Uh0×Uh0,

(4)

其中左端的雙線性形式

(5)

右端項(xiàng)(f,vh)需要根據(jù)具體的邊界條件進(jìn)行相應(yīng)的調(diào)整, 這里Uh0表示邊界為零的緊支集函數(shù)空間.

分別選用Lagrange型線性基函數(shù)和二次基函數(shù)構(gòu)造有限維泛函空間.將區(qū)域離散化剖分形成三角形網(wǎng)格，然后以每個(gè)三角形單元的3個(gè)頂點(diǎn)作為節(jié)點(diǎn)構(gòu)造局部線性基函數(shù).在此基礎(chǔ)上, 再增設(shè)3條邊的中點(diǎn)作為單元新節(jié)點(diǎn),進(jìn)而構(gòu)造局部二次基函數(shù).

1) 線性基函數(shù)的構(gòu)造.選取等參三角形單元的3個(gè)頂點(diǎn)Α1=(0,0),A2=(1,0),A3=(0,1),A1,A2,A3須滿足局部基函數(shù)ψj的構(gòu)造條件:

2) 二次基函數(shù)的構(gòu)造.除選取等參三角形單元的3個(gè)頂點(diǎn)Α1=(0,0),A2=(1,0),A3=(0,1)之外, 再增加3條邊的中點(diǎn)A4=(1/2,0),A5=(1/2,1/2),A6=(0,1/2).各節(jié)點(diǎn)均須滿足局部基函數(shù)ψj的構(gòu)造條件:

3 有限元法處理stress邊界條件

(6)

進(jìn)一步地,變分形式(3)為

(7)

(8)

根據(jù)式(7)可知, 在處理stress邊值時(shí)須先對(duì)其右端標(biāo)量積項(xiàng)進(jìn)行細(xì)節(jié)處理.對(duì)于局部邊界Γs??Ω, 設(shè)定二維問題邊界Γs的4條邊的法向量依次為nt=(0,1)T,nb=(0,-1)T,nl=(-1,0)T,nr=(1,0)T, 切向量依次為τt=(-1,0)T,τb=(1,0)T,τl=(0,-1)T,τr=(0,1)T. 將其代入具體邊界條件, 可得出pn和pτ.由式(8)可知, 右端的載荷向量還存在切線方向和法線方向的stress邊值.故在法線方向和切線方向進(jìn)行曲線積分, 分別對(duì)應(yīng)添加到f的兩個(gè)分量f1,f2上, 從而組裝為新的總載荷向量, 即全局區(qū)域的stress邊界總載荷向量.

4 誤差分析

引理1[13]若問題(1)的真解u和右端f具備必要的光滑性, 則線性有限元解uh的誤差估計(jì)為

‖u-uh‖∞≤Ch2‖f‖0,

(9)

‖u-uh‖0≤Ch2‖f‖0,

(10)

其中C為與單元步長h無關(guān)的常數(shù).

引理2在引理1條件下, 線性有限元解uh的H1半范數(shù)誤差估計(jì)為

|u-uh|1≤Ch‖f‖0.

(11)

定理1在引理1條件下, 二次有限元解uh的誤差估計(jì)如下:

‖u-uh‖∞≤Ch3‖f‖0,

(12)

‖u-uh‖0≤Ch3‖f‖0,

(13)

|u-uh|1≤Ch2‖f‖0.

(14)

證明 限于篇幅, 這里僅針對(duì)L2誤差范數(shù)證明式(13).二次有限元格式取k=2的分片二次多項(xiàng)式, 解向量u的分量ui(i=1,2)均有‖ui-uih‖0≤Ch|ui-uih|1, 而|ui-uih|1≤Chk|ui|1, 則‖ui-uih‖0≤Chk+1|ui|1.進(jìn)一步根據(jù)引理2, 式(13)得證.

5 數(shù)值實(shí)驗(yàn)

針對(duì)帶stress邊界的矢量型彈性方程, 通過MATLAB編程運(yùn)行算例的實(shí)際結(jié)果, 利用誤差圖示與數(shù)值分析,分別采用線性有限元格式和二次有限元格式驗(yàn)證本文有限元方法的模擬效果.

當(dāng)二維問題各方向取剖分?jǐn)?shù)為26時(shí), 其x方向等距步長h=2-6, 分別應(yīng)用線性有限元格式和二次有限元格式計(jì)算得到解向量的兩個(gè)分量u1,u2的絕對(duì)誤差如圖1～2所示.由圖1～2可見: 線性有限元格式下u1,u2的最大誤差上限為4×10-5, 而二次有限元格式下的最大誤差上限改進(jìn)為4×10-7,表明本文方法的精度更高;相較線性有限元格式,二次有限元格式的誤差的最大值、寬度和高度都明顯減小,故二次有限元格式求解帶stress邊界的矢量型彈性方程時(shí)精確性更高且收斂效果更優(yōu)越.

圖1 h=2-6時(shí)線性有限元格式下分量u1(a), u2(b)的絕對(duì)誤差Fig.1 Absolute errors of components u1(a), u2(b) in linear finite element scheme, with mesh size h=2-6

圖2 h=2-6時(shí)二次有限元格式下分量u1(a), u2(b)的絕對(duì)誤差Fig.2 Absolute errors of u1(a), u2(b) in quadratic finite element scheme, with mesh size h=2-6

為了突出對(duì)比線性有限元格式與二次有限元格式的具體精度和收斂效果, 對(duì)初始粗網(wǎng)格步長h=2-2進(jìn)行折半遞減的剖分加密.在各自粗細(xì)網(wǎng)格分別計(jì)算矢量型真解與有限元解向量之間誤差的L∞范數(shù)、L2范數(shù)和H1半范數(shù), 進(jìn)而計(jì)算對(duì)應(yīng)的收斂階,結(jié)果如表1～2所示.由表1～2可見: 在相同的步長剖分下,二次有限元格式較線性有限元格式的精度高達(dá)103倍.隨著區(qū)域離散三角形網(wǎng)格的加密,線性有限元解與真解之間誤差的L∞范數(shù)和L2范數(shù)為二階收斂,H1半范數(shù)為一階收斂;而二次有限元解與真解之間誤差的L∞范數(shù)和L2范數(shù)均實(shí)現(xiàn)了三階收斂,H1半范數(shù)為二階收斂,與定理1的誤差分析結(jié)果完全一致.這充分驗(yàn)證了二次有限元格式的誤差范數(shù)值顯著降低,收斂階和收斂速度顯著提升.

6 結(jié)語

本文基于有限維空間逼近無限維空間的數(shù)值思想,利用有限元法處理具有stress邊值的二維矢量型彈性力學(xué)方程. 有限元格式在每個(gè)三角形單元分別構(gòu)造線性基函數(shù)和二次基函數(shù), 對(duì)應(yīng)形成稀疏程度和數(shù)據(jù)結(jié)構(gòu)不同的代數(shù)方程組, 利用法向量和切向量刻畫stress邊值對(duì)應(yīng)的函數(shù)細(xì)節(jié). 理論證明與數(shù)值實(shí)驗(yàn)驗(yàn)證了本文方法的高階穩(wěn)定收斂,充分展現(xiàn)了高次有限元計(jì)算格式的高階模擬優(yōu)勢(shì),為有效求解彈性邊值問題提供了高效數(shù)值逼近的參考方案.

表1 各網(wǎng)格步長下線性有限元格式的誤差范數(shù)值及其收斂階

表2 各網(wǎng)格步長下二次有限元格式的誤差范數(shù)值及其收斂階