姚玉卿 光輝 王紹波

摘?要：首先，根據(jù)廣義HR變分原理建立了用于分析二維平面問題的非協(xié)調(diào)廣義混合元。然后，應用該混合元對中心含孔平板的應力及應力集中系數(shù)等問題進行了詳細的分析，并與二維的非協(xié)調(diào)位移元的結(jié)果進行了比較。在有限元模型相同的情況下，非協(xié)調(diào)廣義混合元的結(jié)果精度明顯優(yōu)于常規(guī)的非協(xié)調(diào)位移元。另一方面，為了回避混合元法要求計算機內(nèi)存空間較大的問題，提出了結(jié)構(gòu)的局部區(qū)域非協(xié)調(diào)廣義混合元方法。數(shù)值實例表明，局部區(qū)域非協(xié)調(diào)廣義混合元方法是可靠的。

關(guān)鍵詞：有限元方法;廣義混合元;非協(xié)調(diào)廣義混合元;含孔平板;應力集中系數(shù)

中圖分類號：V257;O343.4文獻標識碼：A

Stress Analysis of Plate witha hole Based

on Incompatible Generalized Mixed Element

Yao Yuzhe?Qing Guanghui*?Wang Shaobo

College of Aeronautical Engineering，CAUC?Tianjin?300300

Abstract：Firstly，based on the generalized HR variational principle，a nonconforming generalized mixed element is established for the analysis of twodimensional problems.Then，the mixed element is applied to analyze the stress and stress concentration factor of the central plate with holes in detail，and the results are compared with those of the twodimensional incompatible displacement element.In the case of the same finite element model，the result accuracy of the incompatible generalized mixed element is obviously better than that of the conventional incompatible displacement element.On the other hand，in order to avoid the problem that the mixed element method requires large memory space，a local incompatible generalized mixed element method is proposed.The numerical examples show that the local incompatible generalized mixed element method is reliable.

Key words：finite Element method;generalized mixed element;nonconforming generalized mixed element;plate with a hole;stress concentration factor

為了滿足裝配或服役過程中的維護修理等方面的要求，在現(xiàn)代工程結(jié)構(gòu)，尤其是飛機器這類復雜的工程結(jié)構(gòu)，常常要求在結(jié)構(gòu)的零部件上打孔或開口。這些孔或開口的出現(xiàn)勢必導致對結(jié)構(gòu)壽命有著重要影響的應力集中問題。對于含有多個孔或復雜開口的零部件，人們通常采用試驗或有限元數(shù)值方法分析其應力集中問題。在網(wǎng)格較稀疏的情況下，最常見的位移有限元法[1]的結(jié)果往往不是很理想。實質(zhì)上，通常的位移有限元法求得的位移精度高，但是應力精度比位移精度差一個階次，這是由于節(jié)點的應力結(jié)果是依據(jù)本構(gòu)關(guān)系逐一單元計算的，因而造成連接多個單元的同一節(jié)點上出現(xiàn)了不同的應力結(jié)果，所以必須對節(jié)點應力進行恢復或磨平的操作。另一方面，這種依據(jù)本構(gòu)關(guān)系計算應力的方案，不便引入應力邊界條件。因此，研究探索更為有效和可靠的新方法具有重要的工程意義。

榮廷玉、黨發(fā)寧和Felippa等人[25]最先提出了用于分析二維問題的廣義混合元。用于分析三維問題的非協(xié)調(diào)廣義混合元[69]是在三維廣義混合元的基礎(chǔ)上對位移變量的插值函數(shù)增加非協(xié)調(diào)項而建立的。一方面，在計算量增加不多的情況下，單元的精度得到了進一步提高。另一方面，非協(xié)調(diào)廣義混合元可克服剪切鎖死問題。不論是二維的、三維的廣義混合元還是非協(xié)調(diào)廣義混合元，其共同優(yōu)點是，與位移元相比較，在網(wǎng)格較稀疏的情況下，能給出更為精確的應力結(jié)果，位移值和應力值是同時求得的，且穩(wěn)定性好，收斂快[29]。但是，由于基于混合元的線性系統(tǒng)方程包含了位移和應力兩類變量，所以這類方法需要的計算機內(nèi)存空間相對較大。

實際上，對于一般的工程結(jié)構(gòu)問題，設(shè)計人員并不一定要分析結(jié)構(gòu)所有位置的應力結(jié)果。所以，針對結(jié)構(gòu)的開口、孔和裂紋尖端等高應力區(qū)域，可以進一步探索節(jié)省計算機內(nèi)存的局部區(qū)域的混合元法。

1 基礎(chǔ)理論

1.1 二維問題的廣義HR變分原理

忽略體積力，并且假設(shè)位移邊界條件u-u-=0，則廣義HR變分原理為[2]：

ΠGHR=V[-12σTSσ+σT（SymbolQC@

u）]dV+vλ（12σTSσ+12（SymbolQC@

u）TC（SymbolQC@

u）-σT（SymbolQC@

u））dV-SσT-TudS（1）

式中，u=[u1?u2]T為位移向量，SymbolQC@

為微分操作矩陣，C為材料剛度矩陣，S=C1，V代表體積，S表示平板類結(jié)構(gòu)的面積，T-=[T-1?T-2]T是作用平板類結(jié)構(gòu)在側(cè)面的已知面力。

1.2 非協(xié)調(diào)廣義混合元

四邊形單元的非協(xié)調(diào)位移場可以表示為[10，11]

u=Nqe+Nrre（2）

其中，qe為單元的位移向量，re為在單元內(nèi)部增加的結(jié)點位移向量，N和Nr分別為單元節(jié)點和內(nèi)部節(jié)點的插值函數(shù)矩陣。

假設(shè)應力變量的插值函數(shù)與位移變量插值函數(shù)的協(xié)調(diào)部分相同，則單元應力向量的近似表達式為

σ=Mpe（3）

式中，M是用于表示應力向量（σ11，σ22，σ12）的插值函數(shù)矩陣，pe為單元的應力向量。

根據(jù)文獻[5]，當式（1）中的λ=0.25時，數(shù)值結(jié)果精度最好。將式（2）和式（3）代入式（1），并取λ=0.25，可得到廣義HR變分原理的離散形式：

∏GHR（pe，qe，re）=∑ni=1-38pTeKpppe+34pTeKpqqe

+34pTeKprre+34qTeKqrre

+18qTeKqqqe+18rTeKrrre

-fTeqe（4）

式中，

Kpp=ViMTC-1MdV，

Kpq=ViMT（SymbolQC@

N）dV，

Kpr=ViMT（SymbolQC@

Nr）dV，

Kqr=Vi（SymbolQC@

N）TC（SymbolQC@

Nr）dV，

Kqq=Vi（SymbolQC@

N）TC（SymbolQC@

N）dV，

Krr=Vi（SymbolQC@

Nr）TC（SymbolQC@

Nr）dV，

fe=SσiNTT-dS。

對上式pe，qe，re進行變分，化簡后可得二維問題的非協(xié)調(diào)廣義混合元列式：

-6Rpp?6Rpq

6RTpq?Rqqpe

qe=0

8fq（5）

式中，Rpp=Kpp，

Rpq=Kpq-KqrK1rrKTpr，

Rqq=Kqq-KqrK1rrKTqr。

2 數(shù)值實例

例題1 中心含圓孔的方形有限平板，板的邊長為2a，兩端受均布拉伸載荷σ0作用，如圖1所示。圓孔的半徑r=a/50，材料楊氏模量E=210.0，泊松比μ=0.3。

根據(jù)文獻[12]，在極坐標系下，受均勻拉伸作用的含小圓孔無限大平板的應力解析解為：

σ11=σ02（1-r2ρ2）+σ02（1+3r4ρ4-4r2ρ2）cos2θ

σ22=σ02（1+r2ρ2）-σ02（1+3r4ρ4）cos2θ

σ12=σ21=-σ02（1-3r4ρ4+2r2ρ2）sin2θ（6）

其中，r為小圓孔的半徑。

考慮結(jié)構(gòu)的對稱性，四分之一的網(wǎng)格模型如圖2所示?？拷鼒A孔邊界上的單元于x1方向和x2方向的長度是其他單元相應長度的一半。

這里考慮平面應變問題。表1列出了商用軟件ABAQUS中二維非協(xié)調(diào)位移元（NCE4）和本文的二維非協(xié)調(diào)廣義混合元（NCME4）兩種方法的α點的應力集中系數(shù)。β點的應力結(jié)果列于表2。表1和表2中的精確解由式（6）得到。必須說明，二維的非協(xié)調(diào)位移元法（NCE4）的節(jié)點應力值是對高斯點應力采用外推法得到的。

從表1和表2中的結(jié)果可以看出，非協(xié)調(diào)位移元在α點的應力集中系數(shù)和β點的應力結(jié)果均都偏小，這主要是因為位移元的剛度大。表1表明，非協(xié)調(diào)廣義混合元的結(jié)果大于理論應力集中系數(shù)，這是由于本例分析的是有限尺寸的平板，而理論上的應力集中系數(shù)是以假設(shè)無限大平板為前提的。如果考慮孔的影響，NCME4的應力集中系數(shù)結(jié)果（參見表1的第3列）和β點的應力結(jié)果（參見表2的第3列）均比NCE4更接近精確解。

例題2中心含圓孔的無限大平板受單軸均勻拉伸載荷作用，如圖3（a）所示。材料楊氏模量E=1000，泊松比μ=0.3。這里僅對如圖3（b）所示的靠近孔周圍的子區(qū)域進行數(shù)值計算。

圖4給出了三種網(wǎng)格模型以及相應的位移邊界條件和應力邊界條件。通過精確解公式（6）可以計算出邊界BC和CD各個節(jié)點的應力結(jié)果。這些結(jié)果作為模型邊界BC和CD上節(jié)點的載荷。

這里采用NCME4對圖4（a）、（b）和（c）等三種網(wǎng)格密度不同的模型進行分析。

根據(jù)網(wǎng)格模型4（c）的應力結(jié)果所繪制出的應力云圖見圖5至圖7。

圖5至圖7的三個應力云圖直觀地表明：NCME4的結(jié)果比較真實可靠，尤其是模型邊界處，應力云圖更真實。這是由于基于NCME4的有限元混合模型可方便地引入應力邊界條件。

為了進一步分析NCME4的優(yōu)越性，在表3至表5中詳細列出三種網(wǎng)格模型P點（x1=1.0606602，x2=1.0606602或ρ=1.5，θ=45°）處，NCE4和NCME4兩種方法的σ11、σ22和σ12的數(shù)值結(jié)果。

表3至表5中，σG、σE和σP分別表示非協(xié)調(diào)位移元法求結(jié)構(gòu)節(jié)點位移之后，采用總體磨平、單元磨平和分片磨平等方案所得P點的應力結(jié)果。

表3至表5中的數(shù)據(jù)表明，與非協(xié)調(diào)位移元的結(jié)果相比，本文的非協(xié)調(diào)廣義混合元的結(jié)果更精確。

例題3 中心含孔板的材料以及三種網(wǎng)格模型、位移和應力邊界條件與例題2相同。

下面考慮對局部區(qū)域采用非協(xié)調(diào)廣義混合元計算應力的方法。即先用非協(xié)調(diào)位移元NCE4計算圖4中三個模型的節(jié)點位移，然后利用NCME4分析圖8所示模型中的粗實線所包圍區(qū)域的節(jié)點應力。這里稱這種局部區(qū)域采用非協(xié)調(diào)廣義混合元的方案為局部非協(xié)調(diào)廣義混合元法（LNCME4）。必須指出：圖8中各模型的粗實線是局部模型的邊界線，這些邊界線上節(jié)點的位移是已知的，被作為局部模型的位移邊界條件處理。

表6和表7分別給出了非協(xié)調(diào)位移元、非協(xié)調(diào)廣義混合元和局部非協(xié)調(diào)廣義混合元（LNCME4）等方法在α點σ11和β點σ22的數(shù)值結(jié)果。

表6和表7表明，非協(xié)調(diào)廣義混合元方法和局部非協(xié)調(diào)廣義混合元方法在α和β兩點處的σ11和σ22的精度相差不大。事實上，就圖8（a）這個網(wǎng)格較稀疏的模型而言，非協(xié)調(diào)廣義混合元法和局部非協(xié)調(diào)廣義混合元法的結(jié)果與精確解的最大誤差分別為1.3%和2.3%。但是局部非協(xié)調(diào)廣義混合元法在節(jié)省計算機資源方面有一定的優(yōu)勢。

3 結(jié)論

本文首先根據(jù)廣義HR變分原理建立了用于分析二維問題的非協(xié)調(diào)廣義混合元，并應用該單元對中心含孔平板的應力進行了詳細的分析。具體結(jié)論包括：

（1）有限元網(wǎng)模型相同的情況下，二維問題的非協(xié)調(diào)廣義混合元的應力結(jié)果精度明顯優(yōu)于常規(guī)的非協(xié)調(diào)位移元;

（2）數(shù)值實例表明：本文提出的局部區(qū)域采用非協(xié)調(diào)廣義混合元的方法的應力結(jié)果精度也高于非協(xié)調(diào)位移元。這種局部區(qū)域方法可以很簡便地擴展用于分析具有其他形式的局部高應力區(qū)域結(jié)構(gòu)的應力問題。

參考文獻：

[1]Zienkiewicz O C，Taylor R L，Zhu J Z.The finite element method：Its basis and fundamentals，sixth edition[M].The finite element method：its basis and fundamentals（Seventh Edition）.2005.

[2]榮廷玉.彈性力學廣義混合變分原理及有限元廣義混合法[J].力學進展，1987，17（2）：186199.

[3]黨發(fā)寧，榮廷玉，孫訓方.Reissner板問題的有限元廣義混合法[J].計算力學學報，2000，17（2）：218222.

[4]黨發(fā)寧，李寧，榮廷玉.廣義混合變分原理及其有限元法[J].自然科學進展，2002，12（3）：250254.

[5]Felippa C A.Parametrized multifield variational principles in elasticity：I.Mixed functionals[J].Communications in Applied Numerical Methods，1989，5（2）：7988.

[6]Qing G，Tian J.Highly accurate symplectic element based on two variational principles[J].Acta Mechanica Sinica，2018，34（1）：151161.

[7]劉艷紅，李銳.含參數(shù)辛元與熱彈性復合材料層合板分析[J].復合材料學報，2019，36（5）：13061312.

[8]Guanghui Qing，Junhui Mao，Yanhong Lui.Generalized mixed finite element method for 3D elasticity problems[J].Acta Mechanica Sinica，2018：34（2）：371380.

[9]趙直欽，卿光輝.改進的非協(xié)調(diào)廣義混合單元及性能分析[J].應用數(shù)學與力學，2019，40（5）：518526.

[10]Theodore H.H.Pian，ChangChun Wu，hybrid and incompatible finite element methods[M].Chapman and Hall CRC，2005.

[11]陳萬吉.一個高精度八結(jié)點六面體單元[J].力學學報，1982，18（3）：5867.

[12]吳家龍.彈性力學[M].北京：高等教育出版社，2001.

作者簡介：姚玉喆（1993—），男，漢族，山東德州人，碩士研究生，研究方向：從事民機復合材設(shè)計與研究。

*通訊作者：卿光輝（1968—），男，漢族，湖南新化人，博士，教授，研究方向：復合材料結(jié)構(gòu)力學。