基于相鄰單元高斯點方案的應力恢復方法

卿光輝，徐靖峪

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

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

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

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

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

1 基礎理論

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

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

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

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

最小勢能原理的表達式[16]為

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

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

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

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

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

2 基于相鄰單元高斯點的應力恢復

2.1 單元高斯點的應力計算

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

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

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

對于一般的等參元，由本構關系可得到單元精度較高的高斯點應力

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

2.2 內(nèi)部節(jié)點的應力恢復

當待求節(jié)點為模型的內(nèi)部節(jié)點時，圍繞待求應力的內(nèi)部節(jié)點的高斯點構造等參單元，然后應用內(nèi)推法計算包含在此單元內(nèi)部的節(jié)點應力。

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

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

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

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

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

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

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

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

2.3 邊界節(jié)點的應力計算

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

圖3 外推邊界節(jié)點應力過程Fig.3 Process of extrapolating the stress

圖4 中，以σb表示中點b 的應力，橫坐標取1，以σc表示內(nèi)部節(jié)點c 的應力，橫坐標取2，構造插值多項式，則滿足

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

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

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

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

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

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

3.1 經(jīng)典梁問題

算例1考慮平面應力狀態(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，應力邊界條件[21]為：在x=0，-c ≤y ≤c 處

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

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

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

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

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

表1 網(wǎng)格模型1 A、B、C 和D 點的應力σ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 點的應力σ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 點的應力σ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為應力精確解；σdata為應力有限元解。

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

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

采用圖8 所示的非正交網(wǎng)格模型時計算結果如表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 點的應力σ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 點的應力σ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 點的應力σxy（算例1）Tab.7 Stress σxy at points D，E and F（Case 1）of mesh 2

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

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

3.2 收斂性分析

考慮算例1，5 種均勻的正交網(wǎng)格模型分別為10×2，20×4，30×6，40×8，50×10。 圖9 和圖10 分別列出了5 種網(wǎng)格模型邊界上的應力σx（L，-c）與σxy（L，0）的收斂情況，橫坐標為單元x 方向尺寸的對數(shù)，縱坐標為誤差率error。 圖9 和圖10 表明，相鄰單元法在計算應力σx和σxy時均滿足收斂要求，且在相同密度的網(wǎng)格下計算結果誤差率要小于常規(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 板中含孔應力集中問題分析

算例3考慮如圖11 所示中心含圓孔的方形板。板長、寬均為a=8；孔半徑r=1；厚度t=1。材料參數(shù)：E=1 000，v=0.25。如圖12 所示，采用其1/4 模型劃分網(wǎng)格并進行分析，位移邊界條件為：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

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

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

4 結語

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