卿光輝，徐靖峪

（中國民航大學(xué)航空工程學(xué)院，天津 300300）

通常情況下，位移有限元法的節(jié)點(diǎn)位移精度高[1]，但由于導(dǎo)數(shù)運(yùn)算后應(yīng)變插值函數(shù)的階次降低，導(dǎo)致應(yīng)力結(jié)果精度低于位移結(jié)果精度。為了提高位移有限元法的節(jié)點(diǎn)應(yīng)力精度，Oden 等[2]提出了將相鄰單元在公共節(jié)點(diǎn)處的應(yīng)力取平均值來獲取改進(jìn)的應(yīng)力結(jié)果，但這種方法所得到的應(yīng)力結(jié)果誤差較大。Hinton 等[3]基于最小二乘法思想提出了總體應(yīng)力磨平法，但該方法的計(jì)算量大，一般不被人們所采用。單元應(yīng)力磨平法[3-4]可以方便地利用精度較高的高斯點(diǎn)（最佳應(yīng)力點(diǎn)）的應(yīng)力值，提高節(jié)點(diǎn)上的應(yīng)力結(jié)果精度，將高斯點(diǎn)應(yīng)力用節(jié)點(diǎn)應(yīng)力改進(jìn)值構(gòu)造的插值函數(shù)表示，然后4 個(gè)節(jié)點(diǎn)的應(yīng)力可由高斯點(diǎn)的應(yīng)力外推得到，最后可將相鄰單元在公共節(jié)點(diǎn)的不同數(shù)值結(jié)果取平均值作為該節(jié)點(diǎn)的應(yīng)力結(jié)果[5]。 該方法也稱為外推法，是目前最常用的一種應(yīng)力恢復(fù)方法。外推法具有普遍適用性，但主要缺點(diǎn)是節(jié)點(diǎn)上的應(yīng)力值是通過外推法得到的，未能充分利用單元最佳應(yīng)力點(diǎn)精度高的特性，在應(yīng)力梯度較大處或非均勻網(wǎng)格模型中，這種方法具有一定的局限性[6]。 文獻(xiàn)[5]也指出外推法所獲得的節(jié)點(diǎn)應(yīng)力精度相對于高斯點(diǎn)上的應(yīng)力而言并不十分理想。 文獻(xiàn)[7-9]提出的分片應(yīng)力磨平法和改進(jìn)的平衡分片磨平法能較大幅度地提高節(jié)點(diǎn)的應(yīng)力結(jié)果精度，但由于該方法較為復(fù)雜，且對規(guī)則的正交有限元網(wǎng)格無效，因此，包含這類應(yīng)力恢復(fù)法的商用軟件并不多見。Zhang 等[10]提出了一種對節(jié)點(diǎn)位移進(jìn)行最小二乘擬合的應(yīng)力恢復(fù)方法，可提高應(yīng)力結(jié)果精度；徐小明等[11]結(jié)合辛對偶體系下的解析辛本征展開解，提出了一種關(guān)于二維問題的應(yīng)力磨平方法，使得特定域內(nèi)應(yīng)力場具有與位移場相媲美的精度，但對于不同的問題，該方法中的應(yīng)力解析函數(shù)不確定。

隨著近年來研究的不斷深入，一些新的應(yīng)力恢復(fù)技術(shù)被提出。 Sharma 等[12]提出了一種低階有限元應(yīng)力恢復(fù)技術(shù)，通過平衡條件得到改進(jìn)的恢復(fù)應(yīng)力場；Moslemi 等[13]提出了一種基于每個(gè)高斯點(diǎn)應(yīng)力值分布統(tǒng)計(jì)的誤差估計(jì)方法，獲得了有利于應(yīng)力計(jì)算的最優(yōu)有限元網(wǎng)格；趙亞飛等[6]應(yīng)用高斯過程回歸方法對有限元應(yīng)力解進(jìn)行了改善研究。

隨著人工神經(jīng)網(wǎng)絡(luò)被廣泛應(yīng)用于有限元方法進(jìn)行結(jié)構(gòu)分析，Khoei 等[14]提出了一種采用前饋-反向傳播多層感知器（MLP，multi layer perceptron）神經(jīng)網(wǎng)絡(luò)的應(yīng)力恢復(fù)技術(shù)來提高節(jié)點(diǎn)處應(yīng)力場的精度。該方法基于后驗(yàn)Zienkiewicz-Zhu 誤差估計(jì)，實(shí)現(xiàn)了自動(dòng)自適應(yīng)網(wǎng)格細(xì)化， 適用于高應(yīng)力梯度區(qū)域的應(yīng)力場恢復(fù)，在應(yīng)力奇異性問題上也能有效地改善應(yīng)力場。

綜上，針對二維非協(xié)調(diào)4 節(jié)點(diǎn)四邊形單元模型，提出一種基于相鄰單元高斯點(diǎn)應(yīng)力精度高的特性應(yīng)力恢復(fù)方法，并通過數(shù)值算例驗(yàn)證該方法的精度和有效性。

1 基礎(chǔ)理論

1.1 4 節(jié)點(diǎn)四邊形單元的非協(xié)調(diào)形函數(shù)

根據(jù)非協(xié)調(diào)位移元理論[15]，4 節(jié)點(diǎn)四邊形單元的位移場u=[u1u2]T和應(yīng)力場σ=[σxσyσxy]T可表示為

式中：N 是單元形函數(shù)矩陣；Nλ是內(nèi)部點(diǎn)位移形函數(shù)矩陣；qe是單元節(jié)點(diǎn)位移列陣；λe是單元內(nèi)部點(diǎn)位移列陣；pe是單元節(jié)點(diǎn)應(yīng)力列陣。

1.2 二維問題非協(xié)調(diào)位移元列式

最小勢能原理的表達(dá)式[16]為

式中：Δ為微分算子矩陣；C 為材料彈性矩陣；V 為結(jié)構(gòu)體積；b 為物體所受的體積力；Sσ為結(jié)構(gòu)表面積；為作用在結(jié)構(gòu)表面積Sσ上的已知應(yīng)力。

將式（1）代入式（3），得到離散形式最小勢能原理表達(dá)式為

忽略內(nèi)部節(jié)點(diǎn)力fλ[15]，式（4）分別對節(jié)點(diǎn)位移qe和內(nèi)部節(jié)點(diǎn)位移λe變分，得到

由式（6），可將λe用qe表示為

將式（7）代入式（5）消去λe，得到非協(xié)調(diào)位移元的基本方程

2 基于相鄰單元高斯點(diǎn)的應(yīng)力恢復(fù)

2.1 單元高斯點(diǎn)的應(yīng)力計(jì)算

位移有限元法的變量精度問題可歸結(jié)為求解如下泛函

的極小值[3，17]，即求解

值得說明的是，以上理論的前提是單元的雅可比行列式是常數(shù)，且每個(gè)應(yīng)力分量是獨(dú)立假設(shè)的，所以該理論只適應(yīng)于一維單元。 但對于二維四邊形線性單元或三維六面體單元，實(shí)用結(jié)論是：在等參元中n+1階高斯點(diǎn)上的應(yīng)變或應(yīng)力近似解比其他部位具有較高的精度，這些點(diǎn)習(xí)慣上被稱為最佳應(yīng)力點(diǎn)。 因此，通常情況下，一般首先利用本構(gòu)關(guān)系求出單元具有高精度的高斯點(diǎn)應(yīng)力，然后以高斯點(diǎn)應(yīng)力值為基礎(chǔ)再設(shè)法求節(jié)點(diǎn)的應(yīng)力。

對于一般的等參元，由本構(gòu)關(guān)系可得到單元精度較高的高斯點(diǎn)應(yīng)力

式中：Bi為代入每個(gè)高斯點(diǎn)自然坐標(biāo)后的應(yīng)變矩陣；qe為該單元節(jié)點(diǎn)位移列陣。 由式（11）可得到待求應(yīng)力節(jié)點(diǎn)所有相鄰單元的高斯點(diǎn)應(yīng)力，在此基礎(chǔ)上，可充分利用這些高斯點(diǎn)的應(yīng)力，通過插值方法得到待求節(jié)點(diǎn)應(yīng)力，該方法可稱為相鄰單元法。對于模型中不同位置的節(jié)點(diǎn)，可采取不同的插值方案。

2.2 內(nèi)部節(jié)點(diǎn)的應(yīng)力恢復(fù)

當(dāng)待求節(jié)點(diǎn)為模型的內(nèi)部節(jié)點(diǎn)時(shí)，圍繞待求應(yīng)力的內(nèi)部節(jié)點(diǎn)的高斯點(diǎn)構(gòu)造等參單元，然后應(yīng)用內(nèi)推法計(jì)算包含在此單元內(nèi)部的節(jié)點(diǎn)應(yīng)力。

圖1 中，k 為一個(gè)內(nèi)部節(jié)點(diǎn)， 是單元A、B、C、D 的公共節(jié)點(diǎn)，這4 個(gè)單元為節(jié)點(diǎn)k 的相鄰單元，其中每個(gè)單元內(nèi)都有Ⅰ、Ⅱ、Ⅲ、Ⅳ4 個(gè)高斯點(diǎn)。 從相鄰單元的高斯點(diǎn)中選取距節(jié)點(diǎn)k 最近的4 個(gè)高斯點(diǎn)， 且這4個(gè)高斯點(diǎn)需要構(gòu)成一個(gè)凸四邊形。 這樣，這4 個(gè)高斯點(diǎn)可構(gòu)造出一個(gè)包圍節(jié)點(diǎn)k 的四邊形單元，如圖1 中虛線所示，4 個(gè)高斯點(diǎn)是此單元的4 個(gè)節(jié)點(diǎn)。

圖1 內(nèi)推法新單元的構(gòu)造方式Fig.1 The construction of the new unit of interpolation method

視新單元為等參元，令4 個(gè)高斯點(diǎn)的自然坐標(biāo)：A（Ⅲ）為（-1，1），B（Ⅰ）為（1，1），C（Ⅳ）為（-1，-1），D（Ⅱ）為（1，-1），新單元中任意一點(diǎn)的總坐標(biāo)可由等參元坐標(biāo)公式給出

式中：（xi，yi）是高斯點(diǎn)i 的總坐標(biāo)；Ni為高斯點(diǎn)i 對應(yīng)的插值函數(shù)，（1+ξiξ）（1+ηiη），ξi和ηi表示高斯點(diǎn)i 的自然坐標(biāo)（ξi，ηi）。

設(shè)（ξ，η）是節(jié)點(diǎn)k 的自然坐標(biāo)，由于ξ 和η 均在（-1，1）區(qū)間內(nèi)，則由式（12）可求出節(jié)點(diǎn)k 在新單元中的自然坐標(biāo)。 圖2 表示了等參元自然坐標(biāo)系下節(jié)點(diǎn)k與4 個(gè)高斯點(diǎn)的位置關(guān)系。

圖2 等參變換后的新單元Fig.2 The new unit after isoparametric transfer

基于等參元理論，單元內(nèi)任意一點(diǎn)的應(yīng)力可由4個(gè)節(jié)點(diǎn)的應(yīng)力表示

將通過式（12）求得的節(jié)點(diǎn)k 的自然坐標(biāo)代入式（13）中，即可得到節(jié)點(diǎn)k 的應(yīng)力結(jié)果。 顯然，相鄰單元法避免了常規(guī)的外推法中相鄰單元在公共節(jié)點(diǎn)處應(yīng)力結(jié)果取平均值的過程，計(jì)算程序更加簡便。

2.3 邊界節(jié)點(diǎn)的應(yīng)力計(jì)算

圖3 所示，邊界節(jié)點(diǎn)與最近的內(nèi)部節(jié)點(diǎn)之間也分布著單元內(nèi)的高斯點(diǎn)。如果采用線性拉格朗日插值法[19]求邊界節(jié)點(diǎn)應(yīng)力，則需要取兩處已知的具有高精度應(yīng)力值的點(diǎn)來構(gòu)造線性插值函數(shù)向外插值出邊界節(jié)點(diǎn)的應(yīng)力。圖3 中：c 點(diǎn)是距離待求邊界節(jié)點(diǎn)a 最近的內(nèi)部節(jié)點(diǎn)，假設(shè)內(nèi)部節(jié)點(diǎn)c 已采用2.2 節(jié)的方法求得了應(yīng)力值；b 點(diǎn)是內(nèi)部節(jié)點(diǎn)和邊界節(jié)點(diǎn)邊線上的中點(diǎn)，該點(diǎn)應(yīng)力的求解方法與2.2 節(jié)中內(nèi)部節(jié)點(diǎn)應(yīng)力的求解方法相同。

圖3 外推邊界節(jié)點(diǎn)應(yīng)力過程Fig.3 Process of extrapolating the stress

圖4 中，以σb表示中點(diǎn)b 的應(yīng)力，橫坐標(biāo)取1，以σc表示內(nèi)部節(jié)點(diǎn)c 的應(yīng)力，橫坐標(biāo)取2，構(gòu)造插值多項(xiàng)式，則滿足

圖4 外推法的插值函數(shù)Fig.4 Extrapolation function at external nodal

y=L（x）的幾何意義即通過兩點(diǎn)（1，σb）和（2，σc）的直線，L（x）的表達(dá)式為

由邊界節(jié)點(diǎn)a、中點(diǎn)b 和內(nèi)部節(jié)點(diǎn)c 的關(guān)系可知，L（0）即為所求邊界節(jié)點(diǎn)a 的應(yīng)力值。

若待求的邊界節(jié)點(diǎn)是角節(jié)點(diǎn)（如圖3 中的點(diǎn)e），同理，選取內(nèi)部節(jié)點(diǎn)c、兩點(diǎn)連線的中點(diǎn)d，應(yīng)用2.2 節(jié)中的方法可以計(jì)算包圍點(diǎn)d 的應(yīng)力，再由式（15）可給出角節(jié)點(diǎn)e 的應(yīng)力值。

值得說明的是，在使用相鄰單元法進(jìn)行應(yīng)力計(jì)算時(shí)：一方面，雖然邊界節(jié)點(diǎn)采用了外插法，但從理論上講，對于線性插值方法，內(nèi)部節(jié)點(diǎn)和外部節(jié)點(diǎn)的插值結(jié)果精度相同；另一方面，在實(shí)際問題中，不同應(yīng)力分量沿不同方向的變化規(guī)律不同，且在非正交網(wǎng)格中構(gòu)造插值函數(shù)時(shí)，邊界節(jié)點(diǎn)、中點(diǎn)和內(nèi)部節(jié)點(diǎn)連線方向具有任意性，且應(yīng)力可能并不呈線性變化，但若邊界的網(wǎng)格劃分較密，線性插值法得到的應(yīng)力值，通常不會(huì)出現(xiàn)較大誤差[20]。

3 數(shù)值算例及分析

3.1 經(jīng)典梁問題

算例1考慮平面應(yīng)力狀態(tài)下簡支矩形梁（如圖5所示）。 幾何參數(shù)：梁的半寬c=1，梁的長度L=10，厚度t=1。材料參數(shù)：彈性模量E=1 000，泊松比v=0.25。載荷P=80，位移邊界條件為：u1（L，0）=u2（L，0）=0，u1（L，c）=u2（L，-c）=0，應(yīng)力邊界條件[21]為：在x=0，-c ≤y ≤c 處

圖5 簡支矩形梁算例示意圖Fig.5 Diagram of a simply supported rectangular beam

在x=L，-c≤y≤c 處

算例2考慮平面應(yīng)力狀態(tài)下一端固支的矩形梁（如圖6 所示）。幾何參數(shù)：梁的半寬c = 1，梁的長度L =10，厚度t=1。 材料參數(shù)：彈性模量E=1 000，泊松比v = 0.25。載荷P = 100，位移邊界條件為：x = L，-c≤y≤c 時(shí)，u1=u2=0，應(yīng)力邊界條件[22]為：在x=0，-c ≤y ≤c 處

圖6 一端固支的矩形梁算例示意圖Fig.6 Diagram of rectangular beam fixed at one end

兩個(gè)算例的有限元網(wǎng)格模型如圖7 所示。 表1 和表2 分別列出了兩個(gè)算例中的點(diǎn)A（x=L，y=-c），點(diǎn)B（x=0.9L，y=-c），點(diǎn)C（x=0.8L，y=-c）和點(diǎn)D（x=0.95L，y=-0.75c）的應(yīng)力σx的值。表3 和表4 列出了算例1 中的點(diǎn)D（x = 0.95L，y = -0.75c），點(diǎn)E（x = L，y =0），點(diǎn)F（x=L，y=-0.75c）和算例2 中的點(diǎn)E（x=0.05L，y=-0.75c），點(diǎn)F（x=0，y=0），點(diǎn)G（x=0，y=-0.75c）的應(yīng)力σxy的值。 算例1 精確解參見文獻(xiàn)[21]，算例2精確解參見文獻(xiàn)[22]。 兩個(gè)算例中σy的精確解全域?yàn)?，故在此σy的計(jì)算結(jié)果不作對比。誤差率計(jì)算公式為

表1 網(wǎng)格模型1 A、B、C 和D 點(diǎn)的應(yīng)力σx（算例1）Tab.1 Stress σx at points A，B，C and D（Case 1）of mesh 1

表2 網(wǎng)格模型1 A、B、C 和D 點(diǎn)的應(yīng)力σx（算例2）Tab.2 Stress σx at points A，B，C and D（Case 2）of mesh 1

表3 網(wǎng)格模型1 D、E 和F 點(diǎn)的應(yīng)力σxy（算例1）Tab.3 Stress σxy at points D，E and F（Case 1）of mesh 1

圖7 網(wǎng)格模型1：20×4 正交四邊形網(wǎng)格Fig.7 Mesh 1：20×4 orthogonal quadrilateral mesh

式中：σexact為應(yīng)力精確解；σdata為應(yīng)力有限元解。

從表1～表4 中可看出，相鄰單元法的應(yīng)力和結(jié)果精確度比外推法的結(jié)果精確度明顯提高。

表4 網(wǎng)格模型1 E、F 和G 點(diǎn)的應(yīng)力σxy（算例2）Tab.4 Stress σxy at points E，F(xiàn) and G（Case 2）of mesh 1

采用圖8 所示的非正交網(wǎng)格模型時(shí)計(jì)算結(jié)果如表5～表8 所示。

圖8 網(wǎng)格模型2：非正交四邊形網(wǎng)格（157 單元）Fig.8 Mesh 2：non-orthogonal quadrilateral mesh（157 elements）

表5 網(wǎng)格模型2 A、B、C 和D 點(diǎn)的應(yīng)力σx（算例1）Tab.5 Stress σx at points A，B，C and D（Case 2）of mesh 2

表6 網(wǎng)格模型2 A、B、C 和D 點(diǎn)的應(yīng)力σx（算例2）Tab.6 Stress σx at points A，B，C and D（Case 2）of mesh 2

表7 網(wǎng)格模型2 D、E 和F 點(diǎn)的應(yīng)力σxy（算例1）Tab.7 Stress σxy at points D，E and F（Case 1）of mesh 2

表8 網(wǎng)格模型2 E、F 和G 點(diǎn)的應(yīng)力σxy（算例2）Tab.8 Stress σxy at points E，F(xiàn) and G（Case 2）of mesh 2

從表5～表8 中可以看出，對于不均勻的非正交網(wǎng)格模型，相鄰單元法的結(jié)果精確度也優(yōu)于外推法。

3.2 收斂性分析

考慮算例1，5 種均勻的正交網(wǎng)格模型分別為10×2，20×4，30×6，40×8，50×10。 圖9 和圖10 分別列出了5 種網(wǎng)格模型邊界上的應(yīng)力σx（L，-c）與σxy（L，0）的收斂情況，橫坐標(biāo)為單元x 方向尺寸的對數(shù)，縱坐標(biāo)為誤差率error。 圖9 和圖10 表明，相鄰單元法在計(jì)算應(yīng)力σx和σxy時(shí)均滿足收斂要求，且在相同密度的網(wǎng)格下計(jì)算結(jié)果誤差率要小于常規(guī)的外推法。

圖9 σx（L，-c）收斂性Fig.9 Convergence of σx（L，-c）

圖10 σxy（L，0）收斂性Fig.10 Convergence of σxy（L，0）

3.3 板中含孔應(yīng)力集中問題分析

算例3考慮如圖11 所示中心含圓孔的方形板。板長、寬均為a=8；孔半徑r=1；厚度t=1。材料參數(shù)：E=1 000，v=0.25。如圖12 所示，采用其1/4 模型劃分網(wǎng)格并進(jìn)行分析，位移邊界條件為：x=0 處，u1=0；y=0 處，u2=0；x=4 處，受均布載荷q=1。

圖11 含圓孔的方形平板幾何參數(shù)Fig.11 Geometric parameters of square plate with circular hole

圖12 算例3 網(wǎng)格模型與受載情況Fig.12 Mesh and loads of case 3

為研究模型應(yīng)力集中位置的應(yīng)力結(jié)果精確度，表9 列出了點(diǎn)A（0，1）處的應(yīng)力σx和點(diǎn)B（1，0）處的應(yīng)力σy的有限元解和精確解[22]及誤差率。 從表9 中可以看出，相鄰單元法在計(jì)算應(yīng)力集中位置的應(yīng)力結(jié)果精確度優(yōu)于常規(guī)的外推法。

表9 應(yīng)力集中位置σx 和σy（算例3）Tab.9 Stress σx and σy at position of stress concentration（Case 3）

4 結(jié)語

基于位移結(jié)果得到的單元高斯點(diǎn)的高精度應(yīng)力，提出了一種新的應(yīng)力恢復(fù)方法。 具體結(jié)論有：計(jì)算內(nèi)部節(jié)點(diǎn)應(yīng)力的內(nèi)推法避免了外推法中相鄰單元在公共節(jié)點(diǎn)的應(yīng)力結(jié)果取平均值的計(jì)算過程，計(jì)算邊界節(jié)點(diǎn)應(yīng)力采用線性插值方法，計(jì)算較為簡便。 數(shù)值算例表明，在相同的網(wǎng)格模型情況下，所提出的相鄰單元法比常規(guī)的外推法應(yīng)力計(jì)算結(jié)果精確度高，該方法對正交網(wǎng)格與非正交網(wǎng)格均適用。 為了進(jìn)一步驗(yàn)證方法的正確性，對邊界上應(yīng)力結(jié)果的收斂性情況進(jìn)行了分析，應(yīng)力結(jié)果滿足收斂性要求且相較于外推法精度更高。在對含圓孔方形板應(yīng)力集中問題的分析中，運(yùn)用該方法計(jì)算所得的應(yīng)力集中系數(shù)比外推法更接近理論解，可更有效地預(yù)測由應(yīng)力集中現(xiàn)象引起的含孔結(jié)構(gòu)纖維斷裂規(guī)律。與常規(guī)的外推法一樣，相鄰單元法也是一種具有普遍適用性的方法。