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

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

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

?

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

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

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

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

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

有限元法是工程計算和結構分析中應用最為廣泛的數(shù)值計算方法。從解的結構上看，有限元法總體上可分為三種：h型有限元法，p型有限元法和h-p型有限元法。其中h-p型有限元法結合了h型有限元法和p型有限元法，即通過同時減少單元網(wǎng)格尺寸h(細分網(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型有限元法進行了細致的研究，得到了h-p型有限元法的基本逼近結果。Babu?ka[6]對文獻[1]的結果做了實質(zhì)性的改進，并把二維p-型有限元法的結果推廣到二維h-p型有限元法，為二維h-p型有限元法奠定了堅實的數(shù)學理論基礎。

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

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

1 H-P型有限單元解及形狀函數(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ù)c1和c2，使得

(7)

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

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

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

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

(8)

(9)

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

(10)

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

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

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

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

(11)

(12)

(13)

(14)

2. 標準正方形單元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. 標準等邊三角形單元T上的線性形狀函數(shù):

(25)

(26)

(27)

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

4. 標準等邊三角形單元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ù)：

標準正方形單元S上階譜類型的形狀函數(shù)：

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

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

(34)

(35)

(36)

(37)

b. p≥2時邊模式

邊模式指的是與邊相關聯(lián)的形狀函數(shù)，這里可以定義如下與邊Γk(2≤k≤4)相關聯(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)

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

2 L數(shù)值算例

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

2.1 相關參數(shù)

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

圖1 L型鋼板斷面尺寸

2.2數(shù)值計算結果

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

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

圖2 X方向最大位移Ux

圖3 角點處的第一主應力S1

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

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

圖5 角點處的第一主應力S1

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

表1 X方向最大位移Ux和角點處的第一主應力S1計算結果

表2 X方向最大位移Ux和角點處的第一主應力S1計算結果

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

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

表3 X方向最大位移Ux和角點處的第一主應力S1計算結果

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

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

表4 X方向最大位移Ux和角點處的第一主應力S1計算結果

3 結束語

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

(責任編輯：熊文濤)

2016-03-02

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

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

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

O242.21

A

2095-4824(2016)06-0087-06

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