含孔層合板非協調廣義混合元模型的層間應力分析

中圖分類號：0343.8 DOI：10.16579/j.issn.1001.9669.2025.07.016

0 引言

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

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

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

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

目前國內外利用廣義混合元對復合材料層合板層間應力研究得較少。為了準確分析含孔復合材料層合板孔邊處的層間應力，本文建立了含孔層合板廣義混合元模型，利用廣義混合元能夠同時引入應力和位移邊界條件的特點，分別從厚度方向和環向對孔邊的層間應力進行了分析和討論。

1非協調廣義部分混合變分原理

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

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

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

利用8節點六面體非協調單元[22]建立含孔復合材料層合板廣義混合元模型。非協調項的引入使位移場函數增加了二次項，有利于消除剪切自鎖現象。對位移場和應力場分別進行離散，位移 u 和應力 σ 可分別表示為

u=N_dq_e+N_rr_e

σ=N_sp_e

式中， N_d 和 N_r 分別為位移形函數的協調項和非協調項； N_s 為應力形函數； q_e 和 p_e 分別為節點位移和節點應力； r_e 為單元內部節點的非協調位移。

對式（4）中的節點應力 p_e ，節點位移 q_e ，非協調位移 r_e 分別進行變分，可得

-（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

聯立式（7）和式（8）可得

式中，

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

K_qq=K_qq-K_qrK_rr^-1K_qr^T

式（9）為全混合元列式，由于不同鋪層材料參數的區別，因此，平面內應力在不同鋪層之間并不連續。為了克服這一缺點，對式（9）進行處理

將式（9）中的節點應力 p_e 分為層間應力 p_o 和平面內應力 p_i ，分別為 則式（9）中第1個方程可表達為

已知式（13）中 A₂₂ 可逆，則平面內應力 p_i 可表示為

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

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

式中

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

2 算例分析

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

考慮石墨/環氧樹脂復合材料層合板，采用與文獻 [25]^763-777 一致的含孔層合板模型，如圖1所示。板的邊長為 Ψ_a ，寬為 b ，板厚為 h ，圓孔直徑為 ? 。其中， a= 20h=254mm b=16h=203.2mm，h=?=12.7mm 4層對稱層合板，各層厚度相同。受單位單軸拉伸載荷為 σ₀= 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 平面的網格模型，網格總數為576。另外，對于只有 0^° 和 90^° 鋪層的層合板而言，可以采用1/4模型分析該問題，降低計算成本。

圖3所示為受單軸拉伸載荷的含孔復合材料層合板模型。由于厚度方向的材料屬性不連續，層間應力會在厚度方向出現劇烈變化，所以，分別對不同鋪層狀態、厚度方向網格數量不同的層間應力分布結果進行了分析。

2.2 圓孔厚度方向的層間應力分布

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

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

圖9＼～圖13所示分別為 [45^°/-45^°]_s 層合板在 θ=0^° 、30^°，45^°，60^°，90^° 這5個位置的層間應力沿厚度方向分布情況。層間剪應力 σ_xz 和 σ_yz 的大小、方向關于材料對稱面反對稱；層間正應力 σ_zz 大小、方向關于材料對稱面對稱。層間剪應力 σ_xz 和 σ_yz 的值在厚度方向發生了劇烈變化，隨著轉角 θ 變化， σ_xz 的值逐漸增大，大約在 θ=90^° 位置達到峰值。在角度變化過程中，非協調廣義部分混合元的層間剪應力 σ_xz 和 σ_yz 的結果始終與非協調位移元有較好的一致性，在 n=16 時就保證了與 n=32 的非協調位移元結果有相當的精度。另外，在計算層間正應力 σ_z 時， σ_zz 的值不斷減小直至反向，非協調廣義部分混合元始終保證了邊界應力值與實際值的一致性，有效預測了厚度方向層間正應力隨轉角θ 的變化過程。然而，非協調位移元難以保證邊界應力值的正確性。

2.3 圓孔環向的層間應力分布

選取如圖14所示的層合板界面1、界面2、界面3這3個層間界面位置，對 [0^°/90^°] 層合板在圓孔環向的應力分布情況進行分析。列出不同厚度方向網格數量 n 的非協調廣義部分混合元的計算結果，以及厚度方向劃分32個網格（ n=32 時，非協調位移元的計算結果。

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

3結論

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

1）本文方法符合層間應力在層間連續、平面內應力在層間不連續的客觀事實，能夠有效地捕捉孔邊的高應力梯度特性，并且精度優于位移元。

2）非協調廣義混合有限元能夠引入應力邊界條件，在稀網格和密網格下，層合板上、下表面的層間應力都始終與實際情況保持一致，更有利于分析層間應力的分布規律。

3）結合圓孔厚度方向和環向的應力分布，與非協調位移有限元相比，非協調廣義混合元分析得到的層間應力分布更加集中于層合板材料屬性的變化位置，能更加顯著、準確地反映應力奇異性情況，為層合板的優化設計提供了新的思路。

參考文獻（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］秦建兵，趙璽，梁榮娜．一種失效準則在面外載荷下含孔層合板 的應用[J].復合材料科學與工程，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]李姝，李宴，韓飛．開孔復合材料層合板強度不確定性分析[J]. 力學季刊，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］楊潔．含孔復合材料層合板結構應力分析[D]．鄭州：鄭州大學， 2009：48-52. YANG Jie. Stress analysis of composite laminates with holes[D]. Zhengzhou：Zhengzhou University，2009：48-52.（In Chinese）

[8］曾紅燕．低溫環境下復合材料層合板的應力分析[D].大連：大 連理工大學，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］王聿航，卿光輝．基于廣義混合有限元的壓電復合材料層合板 的數值分析[J].復合材料學報，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］王燮，卿光輝．基于非協調廣義部分混合元的層合板邊界效應 分析[J].復合材料學報，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］楊立洲，卿光輝．基于混合元的變剛度層合板靜力學分析[J]. 工程力學，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]田宗漱，卞學．多變量變分原理與多變量有限元方法[M]．北 京：科學出版社，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