章 敏，張大衛(wèi)，張建銘

(昆明理工大學(xué) 建筑工程學(xué)院, 云南 昆明 650500)

?

H-P型有限單元法在L型鋼模型優(yōu)化設(shè)計中的應(yīng)用

章 敏，張大衛(wèi)，張建銘

(昆明理工大學(xué) 建筑工程學(xué)院, 云南 昆明 650500)

提出了一種新的有限單元方法——H-P型有限單元方法，該方法主要是通過減小網(wǎng)格尺寸和提高插值多項式的階數(shù)，在相對較少網(wǎng)格的前提下實(shí)現(xiàn)較高的計算精度，而且該方法在高階計算時利用了低階的計算結(jié)果，避免了重復(fù)計算，大大減少了計算工作量。重點(diǎn)介紹了該方法的理論基礎(chǔ)，形狀函數(shù)的構(gòu)造，目前的主要研究成果，以及使用Stress Check軟件對L型鋼板模型進(jìn)行數(shù)值模擬分析，所得的結(jié)果驗證了所提算法的有效性。

h-p型有限單元方法；型函數(shù)；Stress Check；L型鋼板；數(shù)值模擬

有限元法是工程計算和結(jié)構(gòu)分析中應(yīng)用最為廣泛的數(shù)值計算方法。從解的結(jié)構(gòu)上看，有限元法總體上可分為三種：h型有限元法，p型有限元法和h-p型有限元法。其中h-p型有限元法結(jié)合了h型有限元法和p型有限元法，即通過同時減少單元網(wǎng)格尺寸h(細(xì)分網(wǎng)格)和增加單元上插值多項式的階數(shù)p來提高有限元解的精度。

1981年，Babu?ka[1]最先從理論上研究了p型有限元法，并證明了在擬一致網(wǎng)格上p型有限元法至少具有和h型有限元法一樣的收斂速度，如果解具有rγ類型的奇性，則p型有限元法的收斂速度是h型有限元法收斂速度的兩倍。Gui[2-4]從理論上系統(tǒng)地研究了一維h-p型有限元解的收斂性，并且給出了一維h-p型有限元解的誤差估計。Babu?ka[5]也對h-p型有限元法進(jìn)行了細(xì)致的研究，得到了h-p型有限元法的基本逼近結(jié)果。Babu?ka[6]對文獻(xiàn)[1]的結(jié)果做了實(shí)質(zhì)性的改進(jìn)，并把二維p-型有限元法的結(jié)果推廣到二維h-p型有限元法，為二維h-p型有限元法奠定了堅實(shí)的數(shù)學(xué)理論基礎(chǔ)。

近三十年來，一維和二維的p型和h-p型有限元法取得了很多重要的進(jìn)展，國內(nèi)外的許多專家和學(xué)者取得了大量的研究成果。例如，除前面所提到的一些重要結(jié)果外，Demokowicz[7]較早地研究了一維的h-p型自適應(yīng)有限元法；Szabó[8]進(jìn)行了有限單元法分析；Guo[9]對R3空間中的h-p型有限元法的理論和算法進(jìn)行了研究；Guo[10]研究了二維h-p型有限元法中的一類Schwarz方法；Schwab[11]分析了p和h-p型有限元方法；Sun[12]研究了二維h-p型有限元法的最佳收斂性，等等。

對于三維p型和h-p型有限元法，相應(yīng)的研究成果尚不多見。文獻(xiàn)[13]詳細(xì)研究了三維空間中標(biāo)準(zhǔn)四面體單元上的多項式延拓算子，并將其成功地應(yīng)用到三維h-p型有限元法中。Guo[14]研究了三維空間中標(biāo)準(zhǔn)五面體和六面體單元上的多項式延拓算子，并結(jié)合文獻(xiàn)[13]得到的標(biāo)準(zhǔn)四面體單元上的多項式延拓算子結(jié)果，將其成功地應(yīng)用到三維p型有限元法和h-p型有限元法中。2015年，Zhang[15]研究了擬一致網(wǎng)格上三維h-p型有限元法的收斂性。目前國內(nèi)外對三維p型有限元法和h-p型有限元法的研究還在進(jìn)行之中，p型有限元法和h-p型有限元法在工程計算中的應(yīng)用還有待進(jìn)一步深入研究。

1 H-P型有限單元解及形狀函數(shù)的構(gòu)造

假設(shè)S?U是有限維子空間，有限元法是根據(jù)變分公式

B(u，v)=f(v)，?v∈S.

(1)

尋找有限元解uN∈S，并且滿足

B(uN，v)=f(v)， ?v∈S.

(2)

選擇基函數(shù)φi(x)， 1≤i≤N 且

(3)

得到

(4)

若B是雙線性算子，則我們可得到

(5)

(6)

在有限元單元剖分時，如果定義hi為第i個單元Ωi的單元尺寸，即hi=diag(Ωi)，則若h=hi對所有i都成立，即所有單元尺寸都相同，則這樣的網(wǎng)格剖分稱為一致的網(wǎng)格剖分。若存在實(shí)數(shù)c1和c2，使得

(7)

那么隨著網(wǎng)格的不斷細(xì)分，hi變得越來越小，而實(shí)數(shù)c1和c2不變，則具有這種性質(zhì)的網(wǎng)格稱為擬一致的網(wǎng)格。本文主要討論擬一致網(wǎng)格上的h-p型有限元法。

有限元空間S上的全局基函數(shù)是定義在毎一個單元Ωi上的分段連續(xù)函數(shù)拼接為定義在整個求解區(qū)域Ω上的連續(xù)函數(shù)。為了計算和編程時的快速、準(zhǔn)確和高效，現(xiàn)代有限元法使用如下策略來構(gòu)造基函數(shù)：

(1) 選擇一個標(biāo)準(zhǔn)單元，也叫主單元，在標(biāo)準(zhǔn)單元上定義一系列標(biāo)準(zhǔn)函數(shù)Ni(ξ，η)，Ni(ξ，η)是關(guān)于變量ξ和η的階數(shù)為pj的多項式，這種多項式被稱為形狀函數(shù)。

(2) 以二維情形為例，定義一系列映射Mk，Mk把標(biāo)準(zhǔn)單元S或T映射到第k個單元Ωk，1≤k≤K，Mk: x=Xk(ξ，η)， y=Yk(ξ，η)

(8)

(9)

=Nj(Φk(x，y)，Ψk(x，y))

(10)

根據(jù)形狀函數(shù)的構(gòu)造，有兩種不同類型的形狀函數(shù)：拉格朗日插值類型的形狀函數(shù)和階譜類型的形狀函數(shù)。

下面以二維情形為例，對標(biāo)準(zhǔn)正方形單元S和標(biāo)準(zhǔn)等邊三角形單元T上的這兩類形狀函數(shù)的定義進(jìn)行簡要說明。

(A) 標(biāo)準(zhǔn)正方形單元S上拉格朗日插值類型的形狀函數(shù)：

1. 標(biāo)準(zhǔn)正方形單元S上的雙線性形狀函數(shù)(每一變量的多項式階數(shù)p=1)：

(11)

(12)

(13)

(14)

2. 標(biāo)準(zhǔn)正方形單元S上的雙二次形狀函數(shù):

這里列舉有限元計算分析軟件中常用的兩組形狀函數(shù)。第一組形狀函數(shù)包括以下8個函數(shù)：

(15)

(16)

(17)

(18)

(19)

(20)

(21)

(22)

第二組形狀函數(shù)包括9個函數(shù)，包括兩種選擇：

選擇1：N1至N8與前面相同，只增加一個新的函數(shù)N9=(1-ξ2)(1-η2)使得

N9(P9)=1，N9(Pj)=0，1≤j≤8.

(23)

(24)

3. 標(biāo)準(zhǔn)等邊三角形單元T上的線性形狀函數(shù):

(25)

(26)

(27)

式中，L1+L2+L3≡1。(L1，L2，L3)稱為面積坐標(biāo)，Li代表點(diǎn)(ξ，η)∈T到邊Γi的距離。

4. 標(biāo)準(zhǔn)等邊三角形單元T上的二次形狀函數(shù):

以Li表示有6個二次形狀函數(shù)如下：

N1=L1(2L1-1)，

(28)

N2=L2(2L2-1)，

(29)

N3=L3(2L3-1)，

(30)

N4=4L1L2，

(31)

N5=4L2L3，

(32)

N6=4L3L1.

(33)

式中，每一個Ni都是階數(shù)為2的多項式。

(B) 階譜類型的形狀函數(shù)：

標(biāo)準(zhǔn)正方形單元S上階譜類型的形狀函數(shù)：

a. p=1時節(jié)點(diǎn)模式

節(jié)點(diǎn)模式指的是與節(jié)點(diǎn)相關(guān)聯(lián)的形狀函數(shù)，節(jié)點(diǎn)模式是與拉格朗日插值類型一樣的雙線性函數(shù)：

(34)

(35)

(36)

(37)

b. p≥2時邊模式

邊模式指的是與邊相關(guān)聯(lián)的形狀函數(shù)，這里可以定義如下與邊Γk(2≤k≤4)相關(guān)聯(lián)的形狀函數(shù)：

(38)

(39)

(40)

c. p≥4時內(nèi)部模式

(41)

p=5時，增加兩個內(nèi)部模式形狀函數(shù)

(42)

(43)

p=6時，增加三個內(nèi)部模式形狀函數(shù)

(44)

(45)

(46)

p=7時，增加四個內(nèi)部模式形狀函數(shù)

(47)

(48)

(49)

(50)

p=8時，增加五個內(nèi)部模式形狀函數(shù)

(51)

(52)

(53)

(54)

(55)

以此類推，當(dāng)p>8時，我們可以增加新的內(nèi)部模式形狀函數(shù)，而且這些形狀函數(shù)的順序數(shù)關(guān)于多項式的階數(shù)p是升階譜的。

2 L數(shù)值算例

以h-p型有限元法作為理論和數(shù)值計算基礎(chǔ)，采用基于h-p型有限元法的Stress Check有限元計算分析軟件，分析并計算了L型鋼板在受到側(cè)面均布壓力的荷載作用下，位移和應(yīng)力的變化情況在受到水壓力、淤沙壓力、自重和場壓力等基本組合荷載作用下的應(yīng)力和位移。

2.1 相關(guān)參數(shù)

首先選取受側(cè)面均布荷載壓力的L型鋼板進(jìn)行數(shù)值模擬計算，L型鋼板所用的材料類型為線彈性，具體參數(shù)如下：鋼板長約5 m，彈性模量為200 GPa，泊松比為0.3，均勻分布的壓力為1×106N/m2，除底部采用完全約束外，其他部位沒有約束，圖形實(shí)例如圖1所示。

圖1 L型鋼板斷面尺寸

2.2數(shù)值計算結(jié)果

基于h-p型使用有限元法計算分析軟件Stress Check在不同網(wǎng)格劃分情況下和在取不同插值多項式階的情況下，計算和分析L型鋼板在受到水平向右的均勻分布的壓力的荷載作用下的應(yīng)力大小和位移。

當(dāng)網(wǎng)格數(shù)為588時，受側(cè)面均布荷載1×106N /m2作用下L型鋼板有限元計算模型的具體網(wǎng)格剖分、X方向最大位移Ux和角點(diǎn)處的第一主應(yīng)力S1如圖2和圖3所示。

圖2 X方向最大位移Ux

圖3 角點(diǎn)處的第一主應(yīng)力S1

通過計算我們分別得到單側(cè)受壓L型鋼板在X方向最大位移Ux(單位：mm)、角點(diǎn)處的第一主應(yīng)力S1(單位：Pa)和角點(diǎn)處的第一主應(yīng)力S1的誤差估計(%)(當(dāng)P=1到8時)如表1所示。

當(dāng)網(wǎng)格數(shù)為1728時，受側(cè)面均布荷載1×106N /m2作用下L型鋼板有限元計算模型的具體網(wǎng)格剖分、X方向最大位移Ux和角點(diǎn)處的第一主應(yīng)力S1如圖4和圖5所示。

圖5 角點(diǎn)處的第一主應(yīng)力S1

通過計算，可得單側(cè)受壓L型鋼板在X方向最大位移Ux、角點(diǎn)處的第一主應(yīng)力S1和角點(diǎn)處的第一主應(yīng)力S1的誤差估計(%)(當(dāng)P=1到8時)，結(jié)果見表2。

表1 X方向最大位移Ux和角點(diǎn)處的第一主應(yīng)力S1計算結(jié)果

表2 X方向最大位移Ux和角點(diǎn)處的第一主應(yīng)力S1計算結(jié)果

當(dāng)網(wǎng)格數(shù)為2700時，受側(cè)面均布荷載1×106N /m2作用下L型鋼板有限元計算模型的具體網(wǎng)格剖分、X方向最大位移Ux和角點(diǎn)處的第一主應(yīng)力S1如圖6和圖7所示。

通過計算我們分別得到單側(cè)受壓L型鋼板在X方向最大位移Ux、角點(diǎn)處的第一主應(yīng)力S1和角點(diǎn)處的第一主應(yīng)力S1的誤差估計(%)(當(dāng)P=1到8時)如表3所示。

表3 X方向最大位移Ux和角點(diǎn)處的第一主應(yīng)力S1計算結(jié)果

當(dāng)網(wǎng)格數(shù)為3267時，受側(cè)面均布荷載1×106N /m2作用下L型鋼板有限元計算模型的具體網(wǎng)格剖分、X方向最大位移Ux和角點(diǎn)處的第一主應(yīng)力S1如圖8和圖9所示。

通過計算我們分別得到單側(cè)受壓L型鋼板在X方向最大位移Ux、角點(diǎn)處的第一主應(yīng)力S1和角點(diǎn)處的第一主應(yīng)力S1的誤差估計(%)(當(dāng)P=1到8時)如表4所示。

表4 X方向最大位移Ux和角點(diǎn)處的第一主應(yīng)力S1計算結(jié)果

3 結(jié)束語

從以上計算結(jié)果可以看出，對于L型鋼板在側(cè)向受壓的均布荷載時，使用Stress Check有限元計算分析軟件在劃分3267個網(wǎng)格數(shù)的時候，X方向最大位移Ux=2.791，角點(diǎn)處的第一主應(yīng)力S1=4.419×8Pa，并且誤差在可接受的范圍內(nèi)，說明h-p型有限元法通過同時減小網(wǎng)格尺寸和提高插值多項式的階數(shù)而實(shí)現(xiàn)了較高的計算精度。同時從h-p型有限單元法的表上可以看出來，在P=1、P=2、P=3、P=4、P=5、P=6的情況下，應(yīng)力增大的值是呈跳躍狀的，說明h-p型有限單元法比h型有限元法具有網(wǎng)格劃分少、收斂速度快、計算精度高和誤差小的優(yōu)點(diǎn)，并且保證結(jié)果的正確性和合理性，大大節(jié)省了我們的計算工作量。

[1] Babu?ka I, Szabó M, Katz N. The p-version of the finite element method [J].SIAM J. Numer. Anal., 1981，20：515-545.

[2] Gui W, Babu? ka I. The h, p and h-p versions of the finite element method in 1 dimension, Part 1. The error analysis of the p-version[J].Numer. Math., 1986，49：577-612.

[3] Gui W, Babu? ka I. The h, p and h-p versions of the finite element method in 1 dimension, Part Ⅱ. The error analysis of the h and h-p version[J].Numer. Math., 1986，49：613-657.

[4] Gui W , Babu? ka I . The h, p and h-p versions of the finite element method in 1 dimension, Part Ⅲ. The adaptive h-p version[J].Numer. Math., 1986，49：659-683.

[5] Babu? ka I, Guo B Q. The h-p version of the finite element method, Part 1 : The basic approximation results[J].Comp. Mech., 1986，1：22-41.

[6] Babu? ka I, Suri M. The h-p version of the finite element method with quasiuniform meshes[J].RAIRO Model Math. Anal. Numer., 1987，21：199-238.

[7] Demokowicz L, Oden T, Rachowicz W，et al. Towards a universal h-p adaptive finite element method, I. Constrained approximation and data structure[J].Comp. Meth. Appl. Mech. Engg., 1989，77：79-112.

[8] Szabó B, Babu? ka I . Finite element analysis[M].John Wiley and Sons, Inc, 1991：75-218.

[9] Guo B Q. The h-p version of the finite element method in R3: theory and algorithm[C]// Proceeding of ICOSAHOM 1995, 1995，56：487-500.

[10] Guo B Q, Cao W. An additive Schwarz method for the h-p version of the finite element method in two dimensions[J].SIAM J. Sci. Comput., 1997，18：1267-1288.

[11] Schwab C. p-and h-p finite element methods[M].Oxford Science Publication, 1998：124-256.

[12] Guo B Q, Sun W. The optimal convergence of the h-p version of the finite element method with quasi-uniform meshes[J].SIAM J. Numer. Anal., 2007，45：698-730.

[14] Guo B Q, Zhang J M. Stable and compatible polynomial extensions in three dimensions and applications to the p and h-p finite element method[J].SIAM J. Numer. Anal.，2009，47： 1195-1225.

[15] Zhang J M,He Y,Qi X.The convergence of the h-p version of the finite element method with quasi-uniform meshes in three dimensions[C]//AIP Conference Proceedings 2015, 15：222-338.

(責(zé)任編輯：熊文濤)

2016-03-02

國家自然科學(xué)基金項目(11261026)

章 敏(1987- )，男，安徽省安慶人，昆明理工大學(xué)建筑工程學(xué)院碩士研究生。

張大衛(wèi)(1992- )，男，安徽省宿州人，昆明理工大學(xué)建筑工程學(xué)院碩士研究生。

O242.21

A

2095-4824(2016)06-0087-06

張建銘(1964- )，男，云南省昆明人，昆明理工大學(xué)建筑工程學(xué)院副教授，博士。