楊慶平，王安平，田宗漱

（1.南京理工大學(xué) 能源與動(dòng)力工程學(xué)院，南京 210094；2.中國科學(xué)院空間應(yīng)用工程與技術(shù)中心，北京 100094；3.中國科學(xué)院大學(xué)，北京 100080）

0 引 言

工程結(jié)構(gòu)中大量遇到具有槽孔的層合材料構(gòu)件，由于槽孔的存在，可能在孔邊橫截面上產(chǎn)生高的局部應(yīng)力，導(dǎo)致構(gòu)件分層破壞。因此，研究槽孔附近應(yīng)力的正確分布，對(duì)分析構(gòu)件的破壞機(jī)理及確保其安全，具有重要意義。

對(duì)于槽孔邊沿的影響問題，正如Wang及Choi[1-2]所指出的，由于它們存在以下內(nèi)在的復(fù)雜性：各層力學(xué)性質(zhì)強(qiáng)烈的各向異性；通過板厚度方向材料的突然改變；沿板邊界幾何形狀的不連續(xù)性；靠近層合材料邊沿區(qū)域平面和橫向變形與應(yīng)力的耦合；準(zhǔn)確滿足自由邊無外力條件的困難等原因，致使這類問題的研究一直進(jìn)展緩慢。

鑒于有限元方法的通用性及簡(jiǎn)易性，對(duì)于具有一個(gè)中心圓孔層板承受拉伸時(shí)的應(yīng)力，許多學(xué)者應(yīng)用不同類型的有限元進(jìn)行過分析：Chauduri等[3-4]應(yīng)用新型三維板殼理論及勢(shì)能原理建立的三角形元；Temiz等[5]依據(jù)等參位移元；Hu等[6-7]采用三維等參元；Pian及Li[8]利用雜交應(yīng)力層合元；Nishioka等[9]依據(jù)一種特殊元；李及謝[10]則采用了一種擬三維有限元迭代法；Tian、Spilker等[11-14]利用基于一種修正余能原理的特殊元。

數(shù)值結(jié)果表明，層合材料靠近邊沿很小局部區(qū)域內(nèi)，呈現(xiàn)一種應(yīng)力梯度快速變化的復(fù)雜三維應(yīng)力狀態(tài)。由于有限元的收斂速度由鄰近高應(yīng)力梯度區(qū)域解的性質(zhì)所支配，用一般假定位移元及一般假定應(yīng)力元求解不僅精度不理想，而且收斂很慢，也就是根據(jù)高階多項(xiàng)式插值函數(shù)所構(gòu)造的一般高精度元，對(duì)改善這種工況下的收斂速度無益，除非應(yīng)用極密的計(jì)算網(wǎng)格，否則難以得到精確的自由孔邊應(yīng)力分布[15-18]。

本文依據(jù)田所導(dǎo)出的一種非勻質(zhì)材料的擴(kuò)展Hellinger-Reissner原理[19]，建立了一種新型具有一個(gè)無外力圓柱表面的三維雜交應(yīng)力元。對(duì)于具有圓柱形槽孔的層合構(gòu)件，它提供較一般假定位移元及一般雜交應(yīng)力元準(zhǔn)確的孔邊三維應(yīng)力分布。此種元亦可有效地分析具有倒圓角矩形孔、U-型孔、擬橢圓孔的多類圓柱型槽孔的層合板自由孔邊的三維應(yīng)力分布。

1 非勻質(zhì)材料的擴(kuò)展Hellinger-Reissner原理[19]

當(dāng)物體被離散為n個(gè)有限元時(shí)，其擴(kuò)展的Hellinger-Reissner原理為

約束條件

式中：B（σ ）為材料余能密度；σ 為應(yīng)力；D為微分算子陣；u為位移；及分別為已知的體積力、位移及表面力；Vnm為第 n 個(gè)元內(nèi)第 m 個(gè)子域的體積；?Vn為第 n 個(gè)元的外表面，?Vn=Sσn∪Sun∪Sab；Sun為第 n個(gè)元的位移邊界；Sσn為第n個(gè)元外力邊界；v為邊界外向法線方向余弦；m為子域數(shù)；）及）分別為第n個(gè)元內(nèi)第m個(gè)子域的應(yīng)力、位移及已知體積力；u~為子域c與d間引入的子域間位移；u~～為單元a與b間引入的元間位移；Sab及Scd為元a、b間及子域c、d間的邊界面。

當(dāng)在上式的層間及元間再引入條件

層間

元間

約束條件

2 層合元列式[19]

推導(dǎo)單元?jiǎng)偠染仃嚂r(shí)可以只考慮一個(gè)元，不計(jì)體積力及表面力。對(duì)線彈性體，（2）式給出的能量表達(dá)式為

對(duì)各層選取

根據(jù)Pian等所提出的雜交應(yīng)力元理性列式方法[20]，利用下式

作為應(yīng)力約束條件。同時(shí)，為了避免應(yīng)力場(chǎng)中的高階項(xiàng)與常數(shù)項(xiàng)發(fā)生耦合時(shí)，單元不能通過分片試驗(yàn)，根據(jù)田等建議[23]，將應(yīng)力）分為高階項(xiàng)）及常數(shù)項(xiàng)）兩部分，將平衡約束條件（5）也改為以下兩部分：

這樣，能量泛函（3）成為

同時(shí)，由于（4）式可知

代入（8）式，得

式中

再代入（10）式，即得到單元?jiǎng)偠汝?/p>

3 具有一個(gè)無外力圓柱表面的三維雜交應(yīng)力層合元

圖1 具有一個(gè)無外力圓柱表面的三維雜交應(yīng)力層合元Fig.1 Geometry of the three-dimensional multilayer hybrid stress element with a traction-free cylindrical surface

具有一個(gè)無外力圓柱表面的新型三維應(yīng)力雜交層合元，單元形狀如圖1所示。無外力表面，指單元外表面上無任何外載荷作用的自由表面，這里為圖中陰影部分的圓柱面ABCD。本特殊元嚴(yán)格滿足此圓柱表面上無外力作用的邊界條件。面ABCD和EIFGKH是兩個(gè)圓柱面，面ADHE和BCGF沿徑向并相交于z軸上。單元的下表面CDHKG和上表面BAEIF平行，并與z軸垂直。

3.1 單元位移場(chǎng)

由于可能發(fā)生彎曲/拉伸耦合現(xiàn)象，而且厚的層板也會(huì)產(chǎn)生嚴(yán)重的橫截面翹曲。而前述變分原理的有限元列式可用于對(duì)每層橫向剪應(yīng)變均獨(dú)立處置的厚層板，因此單元位移場(chǎng)選擇每層橫截面可以獨(dú)立的轉(zhuǎn)動(dòng)。

（1）單元協(xié)調(diào)位移uq

選取協(xié)調(diào)位移uq為

j

單元的結(jié)點(diǎn)位移

式中：k為總層數(shù)。

（2）單元內(nèi)位移 uλ

展開協(xié)調(diào)位移場(chǎng) uq，知其為 ξ，η，ζ的非完整二次多項(xiàng)式，缺 ξ2和 ζ2兩項(xiàng)；三次項(xiàng)缺 ξ2η，ξ2ζ，ξζ2，ηζ2，ξ3，η3，ζ3七項(xiàng)。 通過選取元內(nèi)位移 uλ補(bǔ)足所缺的項(xiàng)，得到

3.2 單元假設(shè)應(yīng)力場(chǎng)

初始應(yīng)力場(chǎng)選取為自然坐標(biāo)表示的非耦合完整的二次式。

利用理性約束平衡方程式（6）及（7），可消去35個(gè)β。

圓柱面無外力邊界條件

及層間應(yīng)力互等和底面無外力條件

（19a）式，表示單元下一層上表面的層間應(yīng)力與上一層下表面的層間應(yīng)力相等，滿足了此條件，在單元的層間界面上，應(yīng)力 σz，τrz，τθz保持連續(xù)。 （19b）式，表示單元的底面上無外力作用。 （19b）式不是擴(kuò)展變分原理（2）的必要條件，加入此條件，也不影響此原理的應(yīng)用。如元的底面上有外力作用，則代入相應(yīng)條件即可，無妨以下推導(dǎo)。

另外，單元引入了元內(nèi)位移uλ，因此在層間界面上位移u，v，w是不連續(xù)的。

利用（18）式、（19a）式及（19b）式，再消去40個(gè)β。這樣，總共消去 75個(gè)β。最終得到第i層的應(yīng)力場(chǎng)

式中

這個(gè)應(yīng)力場(chǎng)（20）及位移場(chǎng)（14）和（16）組成單元 Case A。

4 算 例

4.1 中心擬橢圓孔層板承受拉伸

圖2 中心擬橢圓孔層板承受拉伸Fig.2 Laminate plate with a quasi-elliptic hole under tension

圖3 1/8構(gòu)件有限元網(wǎng)格Fig.3 The mesh of 1/8 of structure

表1 層合板單層材料特性Tab.1 Material properties of a single laminar

兩端承受拉伸，中心擬橢圓孔[±45 ]s層板如圖2所示，每層厚度h=0.2a，材料特性見表1。取1/8板進(jìn)行有限元分析，有限元網(wǎng)格見圖3。

沿孔邊圓柱表面（圖中影線部分）用本文建立的特殊元；孔直邊表面用田、趙等[17]建立的具有一個(gè)無外力直表面的雜交應(yīng)力層合元，其余部分用Pian等建立的一般三維雜交應(yīng)力層合元[24]和10結(jié)點(diǎn)三維等參位移層合元聯(lián)合求解。

此問題沒有解析解，用8結(jié)點(diǎn)一般三維位移層合元(ANSYS中“Solid46”單元)，及三級(jí)子結(jié)構(gòu)所得結(jié)果作為參考解，其各級(jí)結(jié)構(gòu)的自由度分別為393、800和2 002，總自由度 3 195（圖 4）。

圖4三級(jí)子結(jié)構(gòu)技術(shù)Fig.4 Sub-structuring technique

θσ0值，圖6-7分別給出應(yīng)力σr及τrθ沿徑向的分布。可見，現(xiàn)在特殊元得到的解相當(dāng)接近參考解。計(jì)算得最大環(huán)向應(yīng)力為 σθ/σ0=6.106（θ=32.27°，[±45]s，z=h）；相應(yīng)參考解為 6.622（θ=32.40°），誤差為-7.8%。

而現(xiàn)在特殊元Case A求解所用自由度為1 000，僅為一般三維等參位移層合元方法的1/3。

表2 計(jì)算的層板孔邊環(huán)向應(yīng)力Tab.2 Computed circumferential stress

表2 計(jì)算的層板孔邊環(huán)向應(yīng)力Tab.2 Computed circumferential stress

單元類型 自由度數(shù)[+4 5/-4]5 s，z=h [9 0/]0 s，z=2 h σ θ max/σ 0 誤差（%） θ σ θmax/σ 0 誤差（%） θ 1 2 C a s e A參考解1 0 0 0 3 1 9 5 6.1 0 6 6.6 2 2-7.8 3 2.2 7°3 2.4 0°8.1 9 0 7.8 8 1 3.9 6 3.4 4°6 3.4 4°

圖5 孔邊環(huán)向應(yīng)力σθ分布Fig.5 σθdistribution along the rim of the hole

圖 6 應(yīng)力σr沿徑向分布（θ=46.31°）Fig.6 σrdistribution along the radial direction(θ=46.31°)

圖 7 剪應(yīng)力 τrθ沿徑向分布（θ=46.31°）Fig.7 τrθ distribution along the radial direction(θ=46.31°)

4.2 單側(cè)U型槽層板承受拉伸

考慮圖8所示單側(cè)U-型槽層板，各層厚0.2a，材料特性見表1。圖9給出1/4板有限元網(wǎng)格，圖中沿孔圓柱表面用本文建立的特殊元，孔直邊表面用田、趙等建立的元[17]，其余部分用Pian等建立的一般雜交元[24]和10結(jié)點(diǎn)等參元。計(jì)算所得孔邊正則化環(huán)向應(yīng)力σθ分布由圖10給出，圖11則給出沿徑向剪應(yīng)力 τrθ的分布。 同樣可見，對(duì)于 σθ及 τrθ，Case A的解也很接近參考解。 對(duì)于[90/0]s板為19.22，相應(yīng)參考解20.35，誤差為-5.6%。參考解同樣用子結(jié)構(gòu)法求得，所用自由度分別為840、1 950和5 162，總的自由度為7 952，而現(xiàn)在特殊元求解所用自由度為758（圖9），僅為位移層合元的1/10。

圖8 單側(cè)U型槽層板承受拉伸Fig.8 Laminate plate with a U-shaped groove under tension

圖9 1/4板的有限元網(wǎng)格Fig.9 The mesh of 1/4 of structure

圖10 孔邊環(huán)向應(yīng)力σθ分布Fig.10 σθdistribution along the rim of the hole

圖 11 剪應(yīng)力 τrθ沿徑向分布（θ=67.5°）Fig.11 τrθ distribution along the radial direction(θ=67.5°)

5 結(jié) 論

本文利用一種非勻質(zhì)材料擴(kuò)展的Hellinger-Reissner原理，及依據(jù)此原理當(dāng)應(yīng)力及位移分布不連續(xù)的非勻質(zhì)有限元?jiǎng)偠攘惺剑⒘司哂幸粋€(gè)給定無外力圓柱表面的層合雜交應(yīng)力元。

單元應(yīng)力場(chǎng)將Pian等發(fā)展起來的勻質(zhì)各向同性材料的雜交應(yīng)力理性列式擴(kuò)展至層合材料。它保存了此理性列式兩大基本點(diǎn)：（1）引入非協(xié)調(diào)位移，使所選應(yīng)力與位移相匹配，以提高單元求解精度；（2）引入自然坐標(biāo)進(jìn)行應(yīng)力及位移插值，從而明顯地降低了單元形狀歪斜對(duì)其精度的敏感性。

數(shù)值算例表明，在十分粗的網(wǎng)格下，此新的層合元不僅可提供較一般假定位移元及一般假定應(yīng)力元準(zhǔn)確的孔邊應(yīng)力，而且也可用于具有穿透圓孔、或半穿透圓孔層合構(gòu)件柱形自由邊的應(yīng)力分析。

[1]Wang S S,Choi I.Boundary-layer effects in composite laminates:Part 1:free-edge stress solutions and basic characteristics[J].Journal of Applied Mechanics-Transactions of the ASME,1982,49:541-548.

[2]Wang S S,Choi I.Boundary-layer effects in composite laminates:Part 2:free-edge stress solutions and basic characteristics[J].Journal Applied Mechanics-Transactions of the ASME,1982,49:549-560.

[3]Chauduri R A.A new three-dimensional shell theory in general(non-lines-of-curvature)coordinates for analysis of curved panels weakened by through/part-through holes[J].Composite Structures,2009,89:321-322.

[4]Wu Z,Chen W J.Stress analysis of laminated composite plates with a circular hole according to a single-layer higherorder model[J].Composite Structures,2009,90:122-129.

[5]Temiz S,?zel A,Aydinfe M D.FE stress analysis of thick composite laminates with a hole in bending[J].Applied Composite Materials,2003,10:103-117.

[6]Hu E Z,Soutis C,Edge E C.Interlaminar stresses in composite laminates with a circular hole[J].Composite Structures,1997,37:223-232.

[7]Yang Z,Kim C B,Cho C D,Beom H G.The concentration of stress and strain in finite thickness elastic plate containing a circular hole[J].International Journal of Solids and Structures,2008,45:713-731.

[8]Pian T H H,Li M.Stress analysis of laminated composites by hybrid finite elements[C]//Kuhn G,Mang H.Discretization Methods in Structural Mechanics,IUTAM/IACM Sym.Vienna,1989:539-548.

[9]Nishioka T,Atluri S N.Stress analysis of holes in angle-ply laminates:an efficient assumed stress‘special-hole-element’approach and a simple estimation method[J].Computers and Structures,1982,15:135-147.

[10]Li R F,Xie Z C.New finite element method for interlaminar stress analysis of composite laminates at free edge[C]//Cheung Y K,et al.Computational Mechanics,Proceeding of the Asian Pacific Conference Computational Mechanics Part 2(of 2).Rotterdam:Balkema A A,1991:1025-1030.

[11]Tian Z S,Liu J S,Ye L,Pian T H H.Studies of stress concentration by using special hybrid stress elements[J].International Journal for Numerical Methods in Engineering,1997,40:1399-1412.

[12]Tian Z S,Tian Z.Improved hybrid solid elements with a traction-free cylindrical surface[J].International Journal for Numerical Methods in Engineering,1990,29:801-809.

[13]Tian Z S,Zhao F D and Tian Z.Special hybrid stress element for stress analyses around circular cutouts in laminated composites[J].Science in China(Series E),2001,44(5):531-541.

[14]Spilker R L,Chou S C.Evalution of hybrid-stress formulation for thick multiplayer laminates[C]//Lenoe E M,Oplinger D W,Burke J J.Proceeding of 4th Conference on Fibrous Composites in Structural Design.New York:Plenum Press,1980:399-410.

[15]Spilker R L,Chou S C.Edge effects in symmetric composite laminates:importance of satisfying the traction-free-edge condition[J].Journal of Composite Materials,1980,14:2-20.

[16]Lakshminarayana H V.Stress distribution around a semi-circular edge-notch in a finite size laminated composite plate under uniaxial tension[J].Journal of Composite Materials,1983,17:357-3610.

[17]Tian Z S,Zhao F D,Yang Q P.Straight free-edge effects in laminated composites[J].Finite Elements in Analysis and Design,2004,41:1-14.

[18]Tian Z S.A study of stress concentration in solids with circular holes by 3-dimensional special hybrid stress finite elements[J].Journal of Strain Analysis,1990,25:29-35.

[19]田宗漱,卞學(xué)鐄.多變量變分原理與多變量有限元方法[M].北京:科學(xué)出版社,2011.Tian Z S,Pian T H H.Variational principles in multiple variables and finite element mothods in multiple variables[M].Beijing:Sciences Press,2011.

[20]Pian T H H,Sumihara K.Rational approach for assumed stress finite element[J].International Journal for Numerical Meth-ods in Engineering,1984,20:1685-1695.

[21]Pian T H H.Finite elements based on consistently assumed stresses and displacements[J].Finite Elements in Analysis and Design,1985,1:131-140.

[22]Pian T H H,Tong P.Relations between incompatible displacement model and hybrid stress model[J].International Journal for Numerical Methods in Engineering,1986,22:173-181.

[23]Tian Z S,Tian Z S.Further study of construction of axisymmetric finite element by hybrid stress method[C]//Proceeding of Third International Conference on Education,Practice and Promotion of Computational Methods in Engineering Using Small Computers.Macao,1990:549-558.

[24]Pian T H H,Mau S T.Some recent studies in assumed stress hybrid models[C]//Oden J T,Clough R W,Yamanoto Y.Advances in Computational Methods in Structural Mechanics and Design.Alabama:University of Alabama in Huntsville Press,1972:87-106.