碳纖維增強復合材料沉頭螺栓連接失效分析

王 強,賈普榮,張 龍,王 剛,曾磊磊

（西北工業大學力學與土木建筑學院，西安710072）

復合材料因其比剛度高、比強度高、耐疲勞等優異的力學性能，在飛行器承載結構中得到了廣泛的應用。由于設計、制造、成本以及維護等方面的制約，飛行器結構中缺少不了一些工藝分離面、維護口蓋和多種外掛接口等[1]，因而結構件間的連接在飛行器中必不可少。螺栓連接因其承載能力強、連接可靠性高、便于安裝和拆卸等優點，被廣泛用于飛行器結構當中；但由于復合材料是各向異性材料，其損傷和破壞模式特別復雜，即使在單軸拉伸作用下，也會有多種損傷狀態同時產生，而在復雜的連接區域，其應力集中狀況會比金屬結構更為嚴重，尤其是對于沉頭螺栓連接，傳力過程中引入的二次彎矩使得連接結構的設計更加保守，進一步限制了復合材料減重優勢的發揮。雖然復合材料的損傷很難直接觀察到，但數值模擬方法有助于更好地分析復合材料的損傷演化過程。因此，建立精確的復合材料連接結構損傷演化模型十分重要。

在過去幾十年里，眾多學者對復合材料的損傷規律進行了深入的研究。Hill[2]和Tsai等[3]最先以金屬等各向同性材料為基礎，提出了Tsai-Wu和Tsai-Hill失效判據，這些簡單的失效準則不考慮復合材料失效模式的多樣化，但考慮了復合材料的正交各向異性特性。雖然這兩種準則都能在一定程度上預測復合材料的失效，但并不能分析損傷對材料剛度退化的影響。后來，學者們提出了考慮不同破壞模式（纖維、基體）的失效模型，其中，Hashin[4-5]的失效模型運用最為廣泛。與Tsai-Wu和Tsai-Hill準則不同的是，Hashin準則不僅區分失效模式，而且在判斷損傷起始后，根據計算得到的損傷變量值來降低損傷點的剛度。目前，仍然有眾多學者基于Hashin準則進行復合材料的損傷演化分析。Fang等[6]采用Hashin三維漸進損傷模型對三維編織復合材料代表性體積單元（RVC）的損傷及擴展進行了分析，并根據主要破壞模式討論了三維四向編織復合材料的非線性行為。Su等[7]對復合材料開孔板的拉伸失效進行了模擬，同時采用界面單元預測層合板的分層失效，得到了與實驗吻合良好的結果。

在復合材料連接結構的研究中，McCarthy等[8-9]建立了復合材料凸頭螺栓連接的三維有限元模型，并研究了螺栓孔間隙對連接結構力學性能的影響。Qin等[10]對比研究了沉頭和凸頭連接對復合材料雙搭接結構強度、剛度和破壞模式的影響，分析了沉頭連接結構弱于凸頭連接結構的主要原因，并通過有限元方法分析了螺栓預緊力和間隙配合對兩種搭接結構承載能力的影響。Joseph等[11]采用內裂紋帶模型對復合材料沉頭螺栓連接結構進行了詳細的漸進損傷和失效分析，驗證了該模型的有效性。Mandal等[12]建立了纖維增強復合材料多釘連接結構的三維漸進損傷模型，并研究了螺栓直徑以及螺栓預緊力對結構承載能力的影響。黃文俊等[13]計算了復合材料層合板單釘沉頭螺栓連接結構的條件擠壓載荷，分析了接觸面摩擦、螺栓剛度和連接金屬板剛度等因素對連接結構拉伸性能的影響。劉向東等[14]研究了層合板釘孔擠壓損傷后剩余剛度的變化規律，提出在模擬擠壓失效過程中采用模量先突減、然后逐漸提升的剛度修正方法。李星等[15]預測了機翼翼根后梁下緣條對接件模型的極限強度，發現材料非線性對計算結果的影響最為顯著，而網格密度、溫差等對計算結果的影響較小。何柏靈等[16]建立了T800級復合材料螺栓單剪連接的強度估算策略，并和實驗結果對比證明了該估算策略的可行性。姜曉偉等[17]研究了間隙與干涉兩種配合方式對接頭剛度的影響，結果表明間隙配合會降低復合材料單釘單剪螺栓連接結構的接頭剛度，而0.5%釘直徑的干涉量則可以有效提高連接接頭剛度。周龍偉等[18]建立了一種新的單向纖維增強樹脂基復合材料的剛度突降退化模型，并用該模型對典型復合材料螺栓連接結構的拉伸失效行為進行分析，驗證了該模型的計算精度和有效性。

本工作建立復合材料沉頭螺栓連接結構的漸進損傷模型，利用改進的三維Hashin失效起始判據和剛度退化方式對結構進行損傷演化分析，并通過編寫UMAT子程序對復合材料本構關系進行實現，同時采用黏性正則化方法改進數值算法的收斂性。最后將數值結果與實驗數據進行比較，為后續相關研究提供參考。

1 漸進損傷理論

1.1 材料本構模型

復合材料層合板是典型的正交各向異性材料，其主坐標系下單層板的應力應變關系為：

其中矩陣C表示剛度矩陣，如式（2）所示：

式中：C11=E11（1-ν23ν32）/?，C12=E22（ν12+ν13ν32）/?，C13=E33（ν13+ν12ν23）/?，C22=E22（1?ν13ν31）/?，C23=E3（ν23+ν21ν13）/?，C33=E33（1?ν12ν21）/?，C44=G12， C55=G13， C66=G23， ?= 1?ν12ν21-ν23ν32?ν13ν31?2ν21ν32ν13。E11、E22、E33分別為 1（軸向）、2（橫向）、3（面外）方向的彈性模量；G12、G13、G23分別為 12、13、23面的剪切模量；ν12、ν13、ν23分別為 12、13、23面的主泊松比。

由于復合材料層合板的橫向各向同性特性，通常認為 E22=E33，G12=G13，ν12=ν13。因此定義材料本構關系僅需材料參數 E11、E22、G12、G23、ν12和 ν23。

1.2 損傷起始判據

采用三維Hashin失效起始判據對復合材料進行損傷判斷，并將復合材料的損傷模式分為以下7種，一旦下列方程左邊的值等于1，即認為復合材料開始出現相應的損傷形式。

纖維拉伸損傷（σ11>0）：

纖維壓縮損傷（σ11<0）：

基體拉伸損傷（σ22>0）：

基體壓縮損傷（σ22<0）：

拉伸分層損傷（σ33>0）：

壓縮分層損傷（σ33<0）：

纖維基體剪切損傷（σ11<0）：

式中：XT和XC分別表示纖維方向的拉伸和壓縮強度；YT和YC分別表示基體方向的拉伸和壓縮強度；ZT和ZC分別表示面外方向的拉伸和壓縮強度；S12、S13和S23分別表示 12、13、23面的剪切強度。

1.3 損傷演化模型

當材料達到損傷起始點后，需要對材料的剛度矩陣進行退化。目前，運用比較廣泛的幾種退化模型有直接剛度降模型、指數型退化模型以及線性退化模型。線性退化模型基于Ba?ant等[19]提出的斷裂帶理論，當復合材料的成分失效時，單元耗散能等于其彈性能，材料的損傷均勻地分散在單元上，表現為材料剛度在單元體積上的均勻退化，其示意圖見圖1。但由于能量耗散隨著網格的細化而減小，所以數值解與單元尺寸密切相關。為了減輕模型的網格依賴性，學者們[20-23]在損傷演化表達式中引入了等效位移的概念來計算材料退化過程中的能量耗散：

式中：i=ft，fc，mt，mc，zt，zc，fms（下同）；表示不同破壞模式下的斷裂韌度；和表示幾種損傷模式下的等效位移和等效應力。表1列出了幾種損傷模式下等效位移的計算形式，其中LC表示單元的特征長度，對于三維實體單元，，其中Ve表示單元的體積。εjj，eq表示三個主方向的等效應變，定義如下式：

圖1 線性退化模型示意圖Fig.1 Schematic diagram of linear degenerate model

根據本構關系，也可將其表示為：

引入α1、α2、β1、β2、β3、β4這幾個參數的好處是可以將失效起始判據與損傷演化模型緊密聯系起來，當剛好滿足某一失效起始條件時，此時的等效失效起始應變表示為：

表1 不同失效模式下的等效位移Table1 Equivalence displacements corresponding to different failure modes

在ABAQUS隱式求解過程中，對單元剛度矩陣進行折減常常會導致計算收斂困難，故參照Lapczyk等[22]采用黏性正則化方法來提高計算的收斂性。加入黏性系數可以使軟化過程中材料的切線剛度矩陣在足夠小的時間增量下是正定的，因而有利于計算的收斂，但過大的黏性系數值往往會導致計算結果失真。因此，通常采用較小（相對于特征時間增量）的黏性系數值，這樣既有助于提高模型在退化狀態下的收斂速度，也不影響計算結果的精確性，采用黏性正則化后的損傷變量由式（18）計算：

對式（18）離散化，可得到t+Δt時刻的黏性損傷變量為：

引入黏性系數會使計算過程中產生黏性耗散能：

得到材料的損傷變量后，需對材料的剛度進行折減。參考文獻[6，25-26]后采用以下的剛度折減方案，折減后剛度矩陣表達式如下：

表2 各損傷狀態變量的最大值Table2 Maximum value of each damage status variable

2 實驗

2.1 試件

參照ASTM D5961設計試件，通過高鎖螺栓將兩塊復合材料層合板搭接，層合板材料為國產T800級碳纖維增強環氧樹脂復合材料（CCF800/AC531）和碳纖維斜紋織物，鋪層順序為[45*/?45/90/45/0/?45/0/45/0/?45/0/45/?45/]s。其中，45*鋪層為碳纖維斜紋織物層，其單層名義厚度為0.22mm；帶橫線上標的90°鋪層表示對稱鋪層的中心層，除織物外共25層，層合板單層名義厚度為0.14mm。沉頭抗剪高鎖螺栓的材料為Ti6Al4V，沉頭角為100°。螺栓和孔之間為理想緊配，并施加7.2N?m的預緊力將其緊固，并于兩板內側粘貼同屬性同厚度的墊板以保證實驗過程中試件的對中性。試件具體尺寸值及貼應變片位置如圖2所示。

2.2 拉伸實驗

拉伸實驗在DNS300液壓伺服電子萬能試驗機上進行。試件通過液壓夾頭夾持在試驗機上，上下夾持長度均為粘貼墊板的長度，采用位移控制加載，加載速率為1mm/min。用試驗機傳感器記錄拉伸載荷，載荷誤差在示值的±1%以內；使用與試驗機同一廠家的引伸計測量孔邊擠壓應變，其標距為50mm；用DH3820N靜態應變測試分析系統監控試件表面關鍵性部位的應變，重復測試3個試件。實驗現場照片見圖3（a），螺栓斷口照片見圖3（b），搭接面的損傷狀態見圖 3（c），外表面（非搭接面）的損傷態見圖 3（d）。

3 數值分析

用商業有限元軟件ABAQUS建立三維連續損傷有限元模型，預測其失效強度并分析損傷形式。復合材料層合板和鈦合金螺栓均采用實體單元建模，單元類型均選用C3D8R。參考文獻[16，27]，復合材料單層板和碳纖維斜紋織物的力學性能見表3；鈦合金螺栓的彈性模量和泊松比分別為110GPa和0.33，參照文獻[28-29]定義鈦合金螺栓的Johnson-Cook塑性模型和Ductile失效模型。

圖2 試件尺寸及貼應變片位置示意圖Fig.2 Schematic diagram of specimen size and position of strain gauge

圖3 實驗照片及斷口 （a）實驗設備；（b）螺栓斷口；（c）搭接面狀態；（d）外表面狀態Fig.3 Exprimental photographs and fractures （a）test device；（b）bolt fracture；（c）lap surface status；（d）outer surface status

網格質量嚴重影響計算結果，需對特定區域進行合適的網格劃分。“危險”區域為孔邊應力集中處，遠離孔邊處幾乎未觀察到損傷現象，所以在近孔處加密網格，并采用網格過渡的方式在遠孔處（夾持段）劃分較為稀疏的網格。經過計算驗證，這樣不僅可保證結果的精確性，還可以有效加快計算速度。在厚度方向上每層劃分一個單元，以保證得到每層單元的應力狀態和損傷情況。另外，對螺栓螺紋采用簡化方式建模，以間隔式的三角槽代替螺桿上的螺紋，并在螺母中建立與螺桿匹配的內螺紋。在螺栓與層合板接觸的區域盡可能保證接觸對之間的網格吻合度，網格劃分后的模型見圖4。

表3 碳纖維增強復合材料單層板和斜紋織物的力學性能Table3 Mechanical properties of carbon fiber reinforced composite laminate and fiber twill fabric

圖4 有限元模型及網格Fig.4 Finite element model and mesh

墊板與主板之間采用“tie”綁定約束，兩板之間、板與螺栓、板與螺母、螺栓和螺母之間的接觸面采用有限滑移的面-面接觸算法，接觸屬性為硬接觸，摩擦系數采用0.115。邊界條件為：模型夾持部分的表面限制橫向位移（U2）和面外方向的位移（U3），模型的左端固定，右端施加軸向位移（U1）。

通過軟件內置的“Bolt load”功能施加螺栓預緊力，需要注意的是一次施加較大預緊力會導致計算出現收斂性問題，因此在第一個分析步中加入較小的預緊力，并在后續分析步加入全部預緊力。螺栓預緊力Fbolt與擰緊力矩P之間的關系為：

式中：D為螺栓直徑；k為擰緊力矩系數；一般取k=0.2[16]。

4 結果與分析

4.1 黏性系數的影響

為了研究清楚黏性系數的取對數值分析結果的影響程度，采用三個常用的黏性系數值0.001、0.0005和0.0001進行數值分析并比較其結果。首先，對幾次計算結果的極限載荷值進行比較。三種黏性系數下的載荷位移曲線如圖5（a）所示，圖中縱軸的載荷表示固定端的軸向支反力，橫軸的位移表示加載端的軸向拉伸位移。從圖5（a）可以看出，三次數值分析的結果十分接近，極限載荷相差在0.1%以內。其次，對計算耗時進行比較。本次計算均是在Intel（R）Xeon（R）CPU E5-2630v4@2.20GHz（64GB RAM）處理器上運行，黏性系數為0.0001 時計算耗時180h，約是黏性系數為0.0005時的2.3倍，是黏性系數為0.001時的6.6倍。最后，為了驗證引入黏性系數對數值分析結果的影響程度，繪制三個黏性系數下計算的總應變能（ABAQUS輸出變量ALLSE）和黏性耗散能（ALLCD）隨軸向載荷的變化趨勢曲線（圖 5（b）），從圖 5（b）可以看出，整個計算過程中黏性耗散能的值相對于應變能都非常小，黏性系數取0.001時黏性耗散能最大，為0.032J，占應變能的0.43%。因此，可以認為黏性系數的加入對計算結果的影響可以忽略。為保證計算精確度的同時又提高計算效率，推薦取黏性系數η=0.001。

圖5 不同黏性系數下的對比 （a）載荷位移曲線；（b）應變能和黏性耗散能Fig.5 Comparison of different viscosity coefficients （a）load-displacement curves；（b）strain energy and viscous dissipation energy

4.2 剛度折減方案的影響

為了提高數值的收斂性，1.3節提出了改進折減剛度矩陣的方法。Fang等[6]在Lapczyk等[22]的基礎上將損傷變量的計算形式進行推廣，對剛度矩陣的主項（Cij，i=j且 i，j=1，6）進行平方次的折減。而筆者認為主項進行平方次的折減會“加速”材料的失效，尤其是對剪切項的折減過快，會使單元在計算過程中發生嚴重的扭曲，從而使計算難以收斂。經過多次計算驗證，筆者發現對剛度矩陣的次項（Cij，i≠j且 i，j=1，3）進行一次或二次的折減對計算結果影響甚微。在圓孔周圍必然存在一些剪應力主導的單元，根據損傷起始判據，剪應力過大會同時使兩個方向滿足損傷起始條件，如果采用Fang的折減方案，會導致剪切剛度折減過快，因此筆者對剪切項采用了式（21）中開方形式的折減因子。并分別采用改進前后的方案進行計算，得到如圖6所示的載荷位移曲線對比圖。可以看到，改進前后計算結果十分接近，且改進后模型的計算收斂性更好。從收斂速率上來看，改進前計算到停止共計301個分析步，而改進后僅用110個分析步便到同等分析量。因此，改進剛度折減方案可以明顯提高分析計算的效率。

4.3 數值分析與實驗結果對比

螺栓連接結構釘孔擠壓實驗共重復三次，表4列出了三次實驗的極限載荷值，其標準差和離散系數都比較小，說明該實驗的可重復性非常好。為了和實驗數據進行對比，在有限元模型上選取與試件貼應變片位置相同的4處單元，輸出其軸向應變，并取單板橫截面上的平均應力作為因變量，得到數值分析結果和實驗結果的應力-應變曲線對比圖（圖7）。

圖6 載荷-位移曲線對比Fig.6 Comparison of load-displacement curves

表4 各試樣的極限載荷值Table4 Ultimate load value of each specimen

圖7 應力-應變曲線 （a）1 號應變片和 2 號應變片；（b）：3 號應變片和 4 號應變片Fig.7 Stress-strain curves （a）strain gauge1#and strain gauge2#；（b）strain gauge3#and strain gauge4#

從圖7（a）看出，有限元預測的1號、2號應變片的軸向應變與實驗值吻合得非常好。由于單搭接產生的次彎曲效應，使得下搭接板出現向下“彎”的趨勢。單純的彎曲效應會使1號應變片呈現負值，2號應變片呈現正值。因此在此彎曲效應和拉伸載荷的共同作用下，1號應變片附近的應變在整個實驗過程中都非常小，實驗和有限元結果都在–150～600με之間變化，而即使在應變非常小的情況下，有限元模擬的1號應變值與實測值十分接近。2號應變片的應變相較于1號非常大，應變保持為正值并逐漸增大。對于2號應變，可以看到預測值和三次實驗值始終十分接近，載荷為7～8kN時誤差最大且不超過3%，說明該模型預測結果的精確性。從圖7（b）可以看出3號、4號的應變與實驗值相對1號、2號誤差略大，原因是3、4號應變片距離孔邊較近，應力狀態十分復雜，預測難度偏大，但從總體趨勢上來看，有限元結果和實驗結果仍比較一致。另外需指出的是，當有限元計算到圖6中所示的最高載荷時，計算結果不能收斂而停止，此時螺栓上的部分單元已經開始進入失效狀態。螺栓的Ductile狀態變量（ABAQUS輸出變量DUCTCRT）最終狀態如圖8（a）所示。因為實驗過程中一旦螺栓上出現少量裂痕，其在高應力水平下會快速擴展，從而出現螺栓突然斷裂的現象；而在ABAQUS隱式分析中，很難模擬該復雜模型的螺栓快速斷裂問題。因此，近似地認為有限元分析在出現收斂問題前的最終載荷值即為該連接結構的極限載荷預測值，為13.36kN。該數值結果比三次實驗結果的平均值低0.61kN，相對誤差為4.4%。

綜上可得，模型對復合材料沉頭螺栓連接結構失效過程預測較好，可用于分析結構損傷失效。

4.4 損傷失效分析

實驗過程中，加載初期載荷平穩上升，兩板搭接部分外側逐漸產生小間隙，板也由平直變得彎曲，說明單搭接的次彎曲效應十分明顯。當載荷達到約10kN時能聽到“滋滋”的聲音，說明此時已發生較為嚴重的擠壓損傷。繼續加載，板的彎曲程度不再明顯增加，此時孔的擠壓損傷進入擴展階段。直至載荷達到極限值時，隨著一聲巨響，螺釘斷裂，整個結構徹底失去承載能力。從實驗結果來看，圖3（c）中的搭接面發生了十分嚴重的擠壓損傷，圓孔被拉成了“橢圓狀”，從圖3（d）中可看出，上搭接板和螺栓頭、螺母接觸的表面上損傷并不明顯，這進一步說明了單搭接中次彎曲導致釘孔處應力不均的問題。

圖8（b）為復合材料極限載荷時螺栓的應力云圖。可以看出，沉頭螺栓單搭接較凸頭螺栓單搭接更容易產生應力集中，沉頭連接的應力集中主要發生在沉頭端和螺紋附近，釘頭錐面部分基本沒有應力集中，這將導致螺桿與沉頭之間應力梯度過大，螺桿頂端更容易發生較大變形，進而導致釘孔間距增大。而且在沉頭段斜面的擠壓作用下，復合材料板會出現分層損傷形式，這進一步加重了單搭接結構的次彎曲問題，也是沉頭連接弱于凸頭連接的原因之一。

根據輸出變量可以監控幾種損傷模式的擴展情況。首先，可以發現最早出現的損傷形式是基體拉伸損傷（圖9（a）），其起始位置在上搭接板（帶沉頭孔）離搭接面最近的90°層上，此時載荷為3.5kN。隨著繼續加載，該種損傷逐漸向上搭接板的內層擴展，并且下搭接板靠近搭接面的幾層也逐漸開始出

圖8 沉頭螺栓最終狀態 （a）DUCTCRT 狀態變量；（b）Mises應力Fig.8 Final status of countersunk bolt （a）DUCTCRT status variable；（b）Mises stress

圖9 損傷狀態 （a）基體初始拉伸損傷；（b）纖維最終拉伸損傷；（c）纖維最終壓縮損傷；（d）基體最終拉伸損傷；（e）基體最終壓縮損傷；（f）Z向最終拉伸損傷；（g）Z向最終壓縮損傷；（h）纖維基體剪切最終損傷Fig.9 Damage status （a）initial tension damage of matrix；（b）final tension damage of fiber；（c）final compression damage of fiber；（d）final tension damage of matrix；（e）final compression damage of matrix；（f）final tension damage ofZ-direction；（g）final compression damage ofZ-direction；（h）final shear damage of fiber-matrix

現基體拉伸損傷和纖維基體剪切損傷形式。然后隨著載荷繼續增大到5kN左右時，在上搭接板右側孔邊處（靠近加載端）和下搭接板左側孔邊處開始出現纖維拉伸和纖維壓縮的損傷形式，且失效都集中在上搭接板孔右側和下搭接板孔左側。載荷為7kN到11kN時為損傷快速擴展時間段，此時復合材料層合板的各損傷變量迅速達到最大值并向周邊擴展，同時螺栓的Ductile狀態變量也出現正值并逐漸上升。當載荷達到12kN時，復合材料板損傷擴展進入緩慢階段，而此時螺栓開始進入損傷積累階段，直至載荷達到極限值時計算出現收斂性問題而停止。在到達極限載荷時，整體的幾種損傷狀態如圖9（b）～（h）所示，可以看出損傷部位主要為上搭接板右側和下搭接板左側的擠壓區。另外，由于加載過程中導致螺栓偏斜，因此上搭接板左側局部也出現了Z向壓縮損傷的情況，而整個過程中并沒有Z向拉伸損傷形式的出現。對比前面的試件斷口照片，可以看出仿真得到的損傷區域和試驗結果吻合良好。

4.5 擠壓強度分析

圖10 擠壓應力-應變曲線Fig.10 Bearing stress-strain curve

5 結論

（1）基于漸進損傷理論，建立了復合材料沉頭螺栓連接的三維連續損傷有限元模型，采用三維Hashin失效準則進行損傷起始判斷，提出一種新的材料損傷后剛度折減方案。就破壞載荷而言，三次實驗結果十分接近，且實驗平均值與數值分析結果誤差為4.4%。而實驗與數值分析的載荷-應變曲線、擠壓應力-應變曲線也基本吻合，說明了基于連續損傷力學的有限元模型是可靠的。

（2）在漸進損傷理論中，提出了改進后的等效位移表達式和剛度矩陣退化方式，可以有效提高數值計算的收斂性；并考慮了引入黏性系數對數值計算結果的影響，發現當黏性系數為0.001時系統的黏性耗散能占應變能比值低于0.5%，因此推薦取該黏性系數值。

（3）采用數值模擬方法對復合材料沉頭螺栓連接結構的失效行為進行分析，發現復合材料板的損傷以孔兩側的擠壓損傷為主，同時伴隨著局部的其他損傷形式，其中帶沉頭孔的層合板損傷更為嚴重，結構最終破壞模式是沉頭螺栓在螺桿螺紋處被剪斷，而在螺栓斷裂之前，搭接板已經超過了其2%偏移擠壓強度值，因此，該搭接結構設計是合理的。