含孔層合板非協(xié)調(diào)廣義混合元模型的層間應(yīng)力分析

中圖分類號(hào)：0343.8 DOI：10.16579/j.issn.1001.9669.2025.07.016

0 引言

復(fù)合材料在機(jī)械制造、航空航天等許多領(lǐng)域中應(yīng)用廣泛，為了滿足工程結(jié)構(gòu)中的各種需要，復(fù)合材料結(jié)構(gòu)上通常會(huì)加工一個(gè)或多個(gè)孔洞，如螺栓孔、裝配孔等。在制孔過程中，孔邊容易發(fā)生應(yīng)力集中而引起疲勞問題。這是由于復(fù)合材料層合板各個(gè)單層之間的拉剪耦合系數(shù)和泊松比不匹配，其沿厚度方向的力學(xué)性能不連續(xù)。為了保證層合板整體變形的協(xié)調(diào)性，孔邊會(huì)產(chǎn)生嚴(yán)重的層間應(yīng)力，導(dǎo)致復(fù)合材料提前失效。因此，通過建立有限元模型來(lái)分析含孔復(fù)合材料板的最終失效強(qiáng)度具有重要價(jià)值，吸引了國(guó)內(nèi)外學(xué)者的廣泛關(guān)注[1-5]

CAMANHO等建立了一種基于有限斷裂力學(xué)的含孔層合板拉伸強(qiáng)度預(yù)測(cè)模型，與傳統(tǒng)的強(qiáng)度預(yù)測(cè)方法相比，該模型的預(yù)測(cè)精度有所提高。楊潔使用Ansys軟件對(duì)對(duì)稱鋪設(shè)的含孔層合板的層間應(yīng)力進(jìn)行了分析，得到了其在孔邊的分布規(guī)律。曾紅燕8研究了在低溫環(huán)境下復(fù)合材料孔邊的層間應(yīng)力，并分析了體積分?jǐn)?shù)和溫度對(duì)層間應(yīng)力分布特點(diǎn)的影響。SANTOS等研究了碳纖維復(fù)合材料開孔層合板在承受面外沖擊載荷時(shí)的損傷過程。JOSEPH等]建立了I2CBM有限元模型，對(duì)碳纖維復(fù)合材料開孔層合板在拉伸/壓縮過程中的漸進(jìn)損傷過程進(jìn)行了研究，分析了層合板的應(yīng)力分布，并與試驗(yàn)結(jié)果取得了良好的一致性，揭示了復(fù)合材料開孔層合板在拉伸過程中的漸進(jìn)損傷機(jī)制。GLIESCHE等[對(duì)碳纖維復(fù)合材料含孔層合板的主應(yīng)力分布進(jìn)行了分析，進(jìn)而設(shè)計(jì)了變角度牽引鋪縫（VariableAngleTowPlacement，VAT）補(bǔ)強(qiáng)結(jié)構(gòu)的纖維軌跡。ZHU等[12]開發(fā)了用橢圓形及圓形VAT補(bǔ)強(qiáng)結(jié)構(gòu)對(duì)開孔碳纖維層合板進(jìn)行粘接補(bǔ)強(qiáng)的方法，并根據(jù)孔周主應(yīng)力方向?qū)ρa(bǔ)強(qiáng)片纖維軌跡進(jìn)行了設(shè)計(jì)。VIDAL等[13]提出了一種變量分離方法，將x-y平面內(nèi)的變量和厚度z方向變量分開求解，在二維平面內(nèi)采用8節(jié)點(diǎn)位移元進(jìn)行離散，厚度方向采用4階分層理論求解，對(duì)不同位置、大小的孔的層合板的層間應(yīng)力分布進(jìn)行了研究，由于變量分離使得問題的維度降低，有效降低了計(jì)算成本。

然而，以上提到的數(shù)值方法均不能同時(shí)求解得到連續(xù)的位移和應(yīng)力結(jié)果，并且難以引入應(yīng)力邊界條件。這是由于基于最小勢(shì)能原理的位移有限元法僅包含位移場(chǎng)變量，求解應(yīng)力時(shí)，需要首先對(duì)最佳應(yīng)力點(diǎn)處的位移求微分，經(jīng)應(yīng)力外推和應(yīng)力磨平處理后得到的節(jié)點(diǎn)應(yīng)力精度差。

近年來(lái)，廣義混合元得到了進(jìn)一步的發(fā)展。QING等[14-16]提出了非協(xié)調(diào)廣義混合元，有效提高了應(yīng)力結(jié)果的收斂速度。王聿航等將非協(xié)調(diào)廣義混合元應(yīng)用于壓電層合板的靜力學(xué)分析中，并使用部分混合元來(lái)避免層間應(yīng)力的“超連續(xù)”問題。王燮等18利用部分混合元建立了復(fù)合材料層合板的有限元模型，對(duì)自由邊界附近的高應(yīng)力梯度特性進(jìn)行了有效分析。楊立洲等9將非協(xié)調(diào)廣義部分混合元模型應(yīng)用于變剛度層合板，分析了纖維鋪設(shè)角度與變剛度層合板平面內(nèi)位移場(chǎng)之間的關(guān)系。LEZGY-NAZARGAH等[20提出了一種低自由度的4節(jié)點(diǎn)四邊形部分混合元，用于Winker-Pasternak彈性地基上功能梯度材料板的靜力學(xué)和自由振動(dòng)分析，所得數(shù)值解的收斂速度快，且對(duì)網(wǎng)格畸變不敏感。

目前國(guó)內(nèi)外利用廣義混合元對(duì)復(fù)合材料層合板層間應(yīng)力研究得較少。為了準(zhǔn)確分析含孔復(fù)合材料層合板孔邊處的層間應(yīng)力，本文建立了含孔層合板廣義混合元模型，利用廣義混合元能夠同時(shí)引入應(yīng)力和位移邊界條件的特點(diǎn)，分別從厚度方向和環(huán)向?qū)走叺膶娱g應(yīng)力進(jìn)行了分析和討論。

1非協(xié)調(diào)廣義部分混合變分原理

設(shè)有限元模型體積為V，表面積為S，廣義混合變分原理[2]為

式中， σ 為應(yīng)力， σ=[σ₁₃σ₂₃σ₃₃σ₁₁σ₂₂σ₁₂]^T; （20c 為彈性材料的剛度系數(shù)矩陣； 為位移， u= [u₁u₂u₃]^T;α 為分裂因子； abla 為微分算子；T為 [Sσ^? 上施加的已知載荷， 。

Π_GHR 屬于多變量變分原理，控制方程中同時(shí)包含了位移場(chǎng)變量 u 和應(yīng)力場(chǎng)變量 σ 。一方面， Π_GHR 可以方便地引入位移邊界條件和應(yīng)力邊界條件；另一方面，應(yīng)力可以從控制方程中直接求解，且自然連續(xù)。這樣避免了位移元求解應(yīng)力時(shí)需先對(duì)位移求偏微分，再進(jìn)行應(yīng)力磨平的問題，減少了計(jì)算量以及對(duì)位移求偏微分所引入的誤差。

利用8節(jié)點(diǎn)六面體非協(xié)調(diào)單元[22]建立含孔復(fù)合材料層合板廣義混合元模型。非協(xié)調(diào)項(xiàng)的引入使位移場(chǎng)函數(shù)增加了二次項(xiàng)，有利于消除剪切自鎖現(xiàn)象。對(duì)位移場(chǎng)和應(yīng)力場(chǎng)分別進(jìn)行離散，位移 u 和應(yīng)力 σ 可分別表示為

u=N_dq_e+N_rr_e

σ=N_sp_e

式中， N_d 和 N_r 分別為位移形函數(shù)的協(xié)調(diào)項(xiàng)和非協(xié)調(diào)項(xiàng)； N_s 為應(yīng)力形函數(shù)； q_e 和 p_e 分別為節(jié)點(diǎn)位移和節(jié)點(diǎn)應(yīng)力； r_e 為單元內(nèi)部節(jié)點(diǎn)的非協(xié)調(diào)位移。

對(duì)式（4）中的節(jié)點(diǎn)應(yīng)力 p_e ，節(jié)點(diǎn)位移 q_e ，非協(xié)調(diào)位移 r_e 分別進(jìn)行變分，可得

-（1-α）K_ppp_e+（1-α）K_pqq_e+（1-α）K_prr_e=0

（1-α）K_pq^Tp_e+αK_qqq_e+αK_qrr_e-f_e=0

αK_rrr_e+（1-α）K_pr^Tp_e+αK_qr^Tq_e=0

由式（5c）可得

將式（6）代入式（5a）和式（5b）可得

（1-α）[K_pq^T-K_qrK_rr^-1K_pr^T]p_e+α[K_qq-K_qrK_rr^-1K_qr^T]q_e=f_e

聯(lián)立式（7）和式（8）可得

式中，

κ_pq=K_pq-K_prK_rr^-1K_qr^T

K_qq=K_qq-K_qrK_rr^-1K_qr^T

式（9）為全混合元列式，由于不同鋪層材料參數(shù)的區(qū)別，因此，平面內(nèi)應(yīng)力在不同鋪層之間并不連續(xù)。為了克服這一缺點(diǎn)，對(duì)式（9）進(jìn)行處理

將式（9）中的節(jié)點(diǎn)應(yīng)力 p_e 分為層間應(yīng)力 p_o 和平面內(nèi)應(yīng)力 p_i ，分別為 則式（9）中第1個(gè)方程可表達(dá)為

已知式（13）中 A₂₂ 可逆，則平面內(nèi)應(yīng)力 p_i 可表示為

p_i=-A₂₂^-1A₂₁p_o+A₂₂^-1κ_pqiq_e

將式（14）代入式（9），則有限元控制方程中只包含位移場(chǎng)變量和層間應(yīng)力場(chǎng)變量，即

式中

FELIPPA[23-24]就廣義混合元中分裂因子 α 的取值問題，提出了基于誤差分析理論的中值法，即假設(shè)精確解位于位移有限元法和應(yīng)力有限元法所求數(shù)值結(jié)果的正中間，推導(dǎo)出 α 值取0.75時(shí)，具有廣泛的適用性。因此，在以下算例中， α 均取值0.75。

2 算例分析

2.1 含孔層合板的廣義混合元模型

考慮石墨/環(huán)氧樹脂復(fù)合材料層合板，采用與文獻(xiàn) [25]^763-777 一致的含孔層合板模型，如圖1所示。板的邊長(zhǎng)為 Ψ_a ，寬為 b ，板厚為 h ，圓孔直徑為 ? 。其中， a= 20h=254mm b=16h=203.2mm，h=?=12.7mm 4層對(duì)稱層合板，各層厚度相同。受單位單軸拉伸載荷為 σ₀= 1MPa 。 0^° 層的材料屬性分別為彈性模量 E₁₁=145GPa E₂₂=10.7GPa E₃₃=10.7GPa ;剪切模量 G₁₂=4.5GPa;G₁₃= 4.5GPa ： G₂₃=3.6GPa ；泊松比 u₁₂=0.31 ： u₁₃=0.31 ：u₂₃=0.49 。

圖2所示為含圓孔層合板 x-y 平面的網(wǎng)格模型，網(wǎng)格總數(shù)為576。另外，對(duì)于只有 0^° 和 90^° 鋪層的層合板而言，可以采用1/4模型分析該問題，降低計(jì)算成本。

圖3所示為受單軸拉伸載荷的含孔復(fù)合材料層合板模型。由于厚度方向的材料屬性不連續(xù)，層間應(yīng)力會(huì)在厚度方向出現(xiàn)劇烈變化，所以，分別對(duì)不同鋪層狀態(tài)、厚度方向網(wǎng)格數(shù)量不同的層間應(yīng)力分布結(jié)果進(jìn)行了分析。

2.2 圓孔厚度方向的層間應(yīng)力分布

為了準(zhǔn)確獲得含孔層合板厚度方向的層間應(yīng)力分布，對(duì)圖1中 θ=0°，30°，45°，60°，90° 這5個(gè)位置的層間應(yīng)力結(jié)果進(jìn)行討論。利用非協(xié)調(diào)廣義部分混合元分別對(duì) [0^°/90^°]_s 和 [45^°/-45^°]_s 含圓孔層合板孔邊的應(yīng)力分布狀態(tài)進(jìn)行研究，分析沿厚度方向不同網(wǎng)格數(shù)量（圖4中的 n ）下的層間應(yīng)力結(jié)果，并與Abaqus軟件的8節(jié)點(diǎn)三維實(shí)體非協(xié)調(diào)位移元（NCSE8）厚度方向劃分32個(gè)網(wǎng)格（ n=32 的結(jié)果進(jìn)行對(duì)比。

圖4＼～圖8分別展示了 [0^°/90^°] 1層合板圓孔周向不同位置的層間應(yīng)力沿厚度方向（即 z 方向，從 0～ 12.7mm 變化，用 h_z 表示）的分布。由圖4＼～圖8可知，由于材料屬性的不連續(xù)，層間應(yīng)力的方向和大小不斷變化，與NCSE8計(jì)算結(jié)果的變化趨勢(shì)有較好的一致性，其中層間正應(yīng)力 σ_z 的大小、方向關(guān)于材料對(duì)稱面對(duì)稱;層間剪應(yīng)力 σ_xz，σ_yz 的大小、方向關(guān)于材料對(duì)稱面反對(duì)稱。所有層間應(yīng)力均先隨著轉(zhuǎn)角 θ 增加而逐漸增加，在 θ=45^°～60^° 時(shí)達(dá)到峰值，然后逐漸減小，至 θ=90^° 時(shí)最小。隨著厚度方向網(wǎng)格數(shù) n 的增加，廣義混合元對(duì)應(yīng)力峰值的預(yù)測(cè)結(jié)果有效改善。 n=8 時(shí)應(yīng)力結(jié)果仍存在小幅度的振蕩情況，且應(yīng)力峰值的預(yù)測(cè)并不準(zhǔn)確，但在 n=16 時(shí)的結(jié)果已和Abaqus軟件結(jié)果有較好的一致性，應(yīng)力曲線的變化情況得到了良好的改善。另外，由非協(xié)調(diào)廣義部分混合元分析含孔層合板的層間應(yīng)力時(shí)，所得層合板上、下表面的層間應(yīng)力值始終與實(shí)際情況一致。然而，傳統(tǒng)位移元不具備這樣的優(yōu)點(diǎn)。

圖9＼～圖13所示分別為 [45^°/-45^°]_s 層合板在 θ=0^° 、30^°，45^°，60^°，90^° 這5個(gè)位置的層間應(yīng)力沿厚度方向分布情況。層間剪應(yīng)力 σ_xz 和 σ_yz 的大小、方向關(guān)于材料對(duì)稱面反對(duì)稱；層間正應(yīng)力 σ_zz 大小、方向關(guān)于材料對(duì)稱面對(duì)稱。層間剪應(yīng)力 σ_xz 和 σ_yz 的值在厚度方向發(fā)生了劇烈變化，隨著轉(zhuǎn)角 θ 變化， σ_xz 的值逐漸增大，大約在 θ=90^° 位置達(dá)到峰值。在角度變化過程中，非協(xié)調(diào)廣義部分混合元的層間剪應(yīng)力 σ_xz 和 σ_yz 的結(jié)果始終與非協(xié)調(diào)位移元有較好的一致性，在 n=16 時(shí)就保證了與 n=32 的非協(xié)調(diào)位移元結(jié)果有相當(dāng)?shù)木?。另外，在?jì)算層間正應(yīng)力 σ_z 時(shí)， σ_zz 的值不斷減小直至反向，非協(xié)調(diào)廣義部分混合元始終保證了邊界應(yīng)力值與實(shí)際值的一致性，有效預(yù)測(cè)了厚度方向?qū)娱g正應(yīng)力隨轉(zhuǎn)角θ 的變化過程。然而，非協(xié)調(diào)位移元難以保證邊界應(yīng)力值的正確性。

2.3 圓孔環(huán)向的層間應(yīng)力分布

選取如圖14所示的層合板界面1、界面2、界面3這3個(gè)層間界面位置，對(duì) [0^°/90^°] 層合板在圓孔環(huán)向的應(yīng)力分布情況進(jìn)行分析。列出不同厚度方向網(wǎng)格數(shù)量 n 的非協(xié)調(diào)廣義部分混合元的計(jì)算結(jié)果，以及厚度方向劃分32個(gè)網(wǎng)格（ n=32 時(shí)，非協(xié)調(diào)位移元的計(jì)算結(jié)果。

圖15＼～圖17所示為層間應(yīng)力 σ_xz，σ_yz，σ_zz 在 [0^°/ 90^° 1層合板界面1、界面2、界面3的應(yīng)力分布結(jié)果。由于層合板形狀和材料屬性的對(duì)稱性，在界面1和界面3層間應(yīng)力值在相同角度的位置等大反向。另外，σ_xz 和 σ_z 沿圓孔環(huán)向 θ=180^° 位置對(duì)稱分布， σ_yz 沿 θ= 180^° 位置反對(duì)稱分布。在層合板對(duì)稱面（界面2）處，σ_xz，σ_yz 的值始終為0。非協(xié)調(diào)廣義部分混合元的應(yīng)力分布結(jié)果和非協(xié)調(diào)位移元一致性較好，層間正應(yīng)力 σ_z 在界面1和界面3的 θ=90^° 和 270^° 的位置達(dá)到峰值。層間剪應(yīng)力 σ_xz，σ_yz 在界面1和界面3的 θ=60^°?120^° 、240^° 和 300^° 達(dá)到峰值。當(dāng) n=8 時(shí)，隨著厚度方向網(wǎng)格數(shù)量 n 的提升，層間應(yīng)力峰值的計(jì)算結(jié)果逐漸趨近于真實(shí)情況。在 n=16 時(shí)，應(yīng)力變化情況已經(jīng)與NCSE8相當(dāng)一致。

3結(jié)論

本文建立了含孔復(fù)合材料層合板的非協(xié)調(diào)廣義混合有限元模型，利用廣義混合元控制方程包含位移場(chǎng)變量和應(yīng)力場(chǎng)變量、無(wú)需應(yīng)力恢復(fù)及額外的數(shù)值穩(wěn)定技術(shù)、求解應(yīng)力精度高的特點(diǎn)，分別從厚度方向和環(huán)向研究了 [0^°/90^°]_s 和 [45^°/-45^°] 含孔層合板的孔邊應(yīng)力分布狀態(tài)，通過與有限元軟件Abaqus的8節(jié)點(diǎn)三維實(shí)體非協(xié)調(diào)位移元結(jié)果做對(duì)比，驗(yàn)證了本文方法的有效性。得到主要結(jié)論如下：

1）本文方法符合層間應(yīng)力在層間連續(xù)、平面內(nèi)應(yīng)力在層間不連續(xù)的客觀事實(shí)，能夠有效地捕捉孔邊的高應(yīng)力梯度特性，并且精度優(yōu)于位移元。

2）非協(xié)調(diào)廣義混合有限元能夠引入應(yīng)力邊界條件，在稀網(wǎng)格和密網(wǎng)格下，層合板上、下表面的層間應(yīng)力都始終與實(shí)際情況保持一致，更有利于分析層間應(yīng)力的分布規(guī)律。

3）結(jié)合圓孔厚度方向和環(huán)向的應(yīng)力分布，與非協(xié)調(diào)位移有限元相比，非協(xié)調(diào)廣義混合元分析得到的層間應(yīng)力分布更加集中于層合板材料屬性的變化位置，能更加顯著、準(zhǔn)確地反映應(yīng)力奇異性情況，為層合板的優(yōu)化設(shè)計(jì)提供了新的思路。

參考文獻(xiàn)（References）

[1]VAN DER SYPT P，CHERIF M，BOIS C. Analysis of the fatigue behaviour of laminated composite holes subjected to pin-bearing loads[J]. International Journal ofFatigue，2017，103：86-98.

[2]KOSTOPOULOS V，KOTROTSOS A，SOUSANIS A，et al. Fatigue behaviour of open-hole carbon fibre/epoxy composites containing bis-maleimide based polymer blend interleaves as selfhealing agent[J]. Composites Science and Technology，2019，171： 86-93.

[3]GUO Q W，ZHANG Y F，LI D S，et al. Tensile properties and failure mechanism of 3D woven composites containing holes of different geometries [J].Thin-Walld Structures，2021，166： 108115.

[4］秦建兵，趙璽，梁榮娜．一種失效準(zhǔn)則在面外載荷下含孔層合板 的應(yīng)用[J].復(fù)合材料科學(xué)與工程，2021（6）：52-57. QIN Jianbing， ZHAO Xi， LIANG Rongna. Application of failure criteria in the composite laminate open hole for out-of-plane loads [J].Composites Science and Engineering，2021（6）：52-57.（In Chinese）

[5]李姝，李宴，韓飛．開孔復(fù)合材料層合板強(qiáng)度不確定性分析[J]. 力學(xué)季刊，2023，44（2）：293-305. LIShu，LI Yan，HANFei. Strength uncertainty analysis for notched composite laminates[J]. Chinese Quarterly of Mechanics，2023，44 （2）：293-305.（In Chinese）

[6]CAMANHO PP，ERCIN G H，CATALANOTTI G，et al. A finite fracture mechanics model for the prediction of the open-hole strength of composite laminates[J].Composites，Part A：Applied Science and Manufacturing，2012，43（8）：1219-122.

[7］楊潔．含孔復(fù)合材料層合板結(jié)構(gòu)應(yīng)力分析[D]．鄭州：鄭州大學(xué)， 2009：48-52. YANG Jie. Stress analysis of composite laminates with holes[D]. Zhengzhou：Zhengzhou University，2009：48-52.（In Chinese）

[8］曾紅燕．低溫環(huán)境下復(fù)合材料層合板的應(yīng)力分析[D].大連：大 連理工大學(xué)，2014：36-40. ZENG Hongyan.Stress analysis of compositelaminatesat cryogenic conditions[D].Dalian：Dalian University of Technology，2014：36-40.（In Chinese）

[9]SANTOS R A M，REIS P N B，SILVA F G A，et al. Influence of inclined holes on the impact strength of CFRP composites [J]. Composite Structures，2017，172（1）：130-136.

[10] JOSEPHAPK，DAVIDSONP，WAASA M.Open hole and filled holeprogressive damage and failure analysis of composite laminates with a countersunk hole [J]. Composite Structures， 2018，203：523-538.

[11]GLIESCHE K，HUBNER T，ORAWETZ H. Application of the tailored fibre placement（TFP）process fora local reinforcement on an“open-hole”tension plate from carbon/epoxy laminates[J]. Composites Science and Technology，2003，63（1）：81-88.

[12]ZHU Y D，QIN Y L，QI SJ，et al. Variable angle tow reinforcement design for locally reinforcing an open-hole composite plate[J]. Composite Structures，2018，202：162-169.

[13]VIDAL P，GALLIMARD L，POLIT O. Modeling of composite plates with an arbitrary hole location using the variable separation method[J]. Computers amp; Structures，2017，192：157-170.

[14]QING G H，MAO JH，LIU Y H. Highly accurate noncompatible generalized mixed finite element method for 3D elasticity problems [J].Journal of Mechanics of Materials and Structures，2017，12 （4）：505-519.

[15］QING G H，MAO JH，LIU Y H. Generalized mixed finite element method for3D elasticity problems[J].Acta Mechanica Sinica，2018，34（2）：371-380.

[16]QING G，TIAN J. Highly accurate symplectic element based on two variational principles[J].Acta Mechanica Sinica，2018，34 （1）：151-161.

[17］王聿航，卿光輝．基于廣義混合有限元的壓電復(fù)合材料層合板 的數(shù)值分析[J].復(fù)合材料學(xué)報(bào)，2022，39（6）：2987-2996. WANG Yuhang，QING Guanghui. Analysis of piezoelectric composite laminates based on generalized mixed finite element [J].Acta Materiae Compositae Sinica，2022，39（6）：2987-2996. （In Chinese）

[18］王燮，卿光輝．基于非協(xié)調(diào)廣義部分混合元的層合板邊界效應(yīng) 分析[J].復(fù)合材料學(xué)報(bào)，2022，39（8）：4172-4178. WANG Xie，QING Guanghui. Analysis of laminates free-edge effect based on noncompatible generalized partial mixed elements [J].Acta Materiae Compositae Sinica，2022，39（8）：4172-4178. （In Chinese）

[19］楊立洲，卿光輝．基于混合元的變剛度層合板靜力學(xué)分析[J]. 工程力學(xué)，2024，41（3）：19-25. YANG Lizhou，QING Guanghui. Static Analysisof variable stifess laminates based on mixed finite element[J].Engineering Mechanics，2024，41（3）：19-25.（In Chinese）

[20]LEZGY-NAZARGAH M，MESHKANI Z. An eficient partial mixed finite element model for static and free vibration analyses of FGMplates rested on two-parameter elastic foundations [J]. Structural Engineering and Mechanics，2018，66（5）：665-676.

[21]RONG TY，LU A Q.Generalized mixed variational principles and solutionsofill-conditioned problemsin computational mechanics： Part I. Volumetric locking[J].Computer Methods in Applied Mechanics and Engineering，2001，191（3/4/5）：407-422.

[22]田宗漱，卞學(xué)．多變量變分原理與多變量有限元方法[M]．北 京：科學(xué)出版社，2011：108-119. TIAN Zongshu，BIAN Xuehuang. Multivariablevariational principle and multivariable finite element method[M].Beijing： Science Press，2011：108-119.（In Chinese）

[23]FELIPPA C A. Parametrized multifield variational principles in elasticity： I. Mixed functionals[J].Communications in Applied NumericalMethods，1989，5（2）：79-88.

[24]FELIPPA C A.Parametrized multifield variational principles in elasticity：II. Hybrid functionals and the free formulation[J]. Communications in Applied Numerical Methods，1989，5（2）： 89-98.

[25]RAGHURAMPV，MURTY AV K.A high precision coupled bending-extension triangular finite element for laminated plates [J].Computers amp; Structures，1999，72（6）：763-777.

Abstract： Inorder to investigate the stressconcentration phenomenon and analyzethe distribution characteristics of the interlaminarstressintheholeedgeregionofcompositelaminates.Basedonthegeneralizedmixed variationalprinciple，the generalizedmixed finiteelement model forlaminated plates with various stackingmodes were established.The stress field variablesweredividedintotheinterlaminarstressandthein-planestress，withtheintroductionofstressboundaryconditiosto ensurethephysicalcontinuityof interlaminarstresses betweenlayersand thediscontinuityofin-plane stresesbetwenlayers. Theinterlaminarstresses attheedgeofthelaminatedplate hole wererespectivelyanalyzed throughthethickness directionand thecircumferentialdirection.Numericalexamplesdemonstrated thattheincompatiblegeneralized mixed elementcouldobtain more accurate stressingularityresults than the8-node thre-dimensional solid incompatible displacement element esults solved bythe finite element softwareAbaqus.Stressesonboth upperand lower surfaces of the laminated plateconsistently reflectedactualsituations.Theresearchindicatesthatcomparedwiththedisplacementelement，theincompatiblegeneralized mixedelementcanmoreefectivelycapture thehighstressgradientofthe interlaminarstresssattheedgeofthelaminated plate hole，which provides a new idea for the optimal design of the laminate.

KeyWords：Laminatedcomposites；Interlaminarstress；Incompatiblegeneralizedmixedelement；Stress boundary condition; Stress singularity