基于非協(xié)調(diào)廣義混合元的含孔平板應(yīng)力分析

姚玉卿 光輝 王紹波

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

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

中圖分類號(hào)：V257;O343.4文獻(xiàn)標(biāo)識(shí)碼：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

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

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

實(shí)際上，對(duì)于一般的工程結(jié)構(gòu)問題，設(shè)計(jì)人員并不一定要分析結(jié)構(gòu)所有位置的應(yīng)力結(jié)果。所以，針對(duì)結(jié)構(gòu)的開口、孔和裂紋尖端等高應(yīng)力區(qū)域，可以進(jìn)一步探索節(jié)省計(jì)算機(jī)內(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)位移場(chǎng)可以表示為[10，11]

u=Nqe+Nrre（2）

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

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

σ=Mpe（3）

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

根據(jù)文獻(xiàn)[5]，當(dāng)式（1）中的λ=0.25時(shí)，數(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。

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

-6Rpp?6Rpq

6RTpq?Rqqpe

qe=0

8fq（5）

式中，Rpp=Kpp，

Rpq=Kpq-KqrK1rrKTpr，

Rqq=Kqq-KqrK1rrKTqr。

2 數(shù)值實(shí)例

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

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

σ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)的對(duì)稱性，四分之一的網(wǎng)格模型如圖2所示。靠近圓孔邊界上的單元于x1方向和x2方向的長(zhǎng)度是其他單元相應(yīng)長(zhǎng)度的一半。

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

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

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

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

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

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

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

為了進(jìn)一步分析NCME4的優(yōu)越性，在表3至表5中詳細(xì)列出三種網(wǎng)格模型P點(diǎn)（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é)點(diǎn)位移之后，采用總體磨平、單元磨平和分片磨平等方案所得P點(diǎn)的應(yīng)力結(jié)果。

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

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

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

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

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

3 結(jié)論

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

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

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

參考文獻(xiàn)：

[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]榮廷玉.彈性力學(xué)廣義混合變分原理及有限元廣義混合法[J].力學(xué)進(jìn)展，1987，17（2）：186199.

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

[4]黨發(fā)寧，李寧，榮廷玉.廣義混合變分原理及其有限元法[J].自然科學(xué)進(jìn)展，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ù)辛元與熱彈性復(fù)合材料層合板分析[J].復(fù)合材料學(xué)報(bào)，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]趙直欽，卿光輝.改進(jìn)的非協(xié)調(diào)廣義混合單元及性能分析[J].應(yīng)用數(shù)學(xué)與力學(xué)，2019，40（5）：518526.

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

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

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

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

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