1 引言
    眾所周知，纖維增強復合材料板殼結構對橫向沖擊作用十分敏感，較小的橫向沖擊作用可能導致面內壓縮強度較大的損失^[1~3]，因此，這類結構的抗沖擊能力一直是人們關注的焦點。求解層合板橫向沖擊下的瞬態響應和沖擊破壞，一般都采用有限元法，并選用合適的本構理論^[4]。
    迄今為止，研究者們采用了不同的復合材料有限元模型，從不同角度去研究層合板的低速沖擊性能。Hou^[5]利用LS-DYNA3D有限元程序，選用實體單元，單元彈性矩陣由單向板彈性常數確定，選用修正的Chang-Chang^[6]破壞準則判斷纖維和基體的初始破壞以及Brewer-Lagace^[7]建議的破壞準則辨別分層,滿足破壞判據的單元采用相應的應力更新方案,以此研究沖擊過程中層合板的損傷問題。程小全等^[8]選用三維20節點等參元，單元彈性矩陣由單向板彈性常數確定，破壞判據選用Chang和Springer分層判據^[9]、Hashin基體開裂破壞判據^[10]，對層合板沖擊后壓縮全過程進行了分析。文獻_^[11]利用3DIMPACT有限元程序，選用實體單元，單元彈性矩陣由單向板彈性常數確定、Choi基體開裂破壞判據^[12] 、Choi和Chang分層破壞判據^[13]，分析低速沖擊荷載作用層合板的響應。文獻[14]利用Abaqus有限元程序，選用殼單元，因殼單元不計平面外應力，在層合板模型中添加純基體層，以該層破壞等效分層破壞，單元彈性矩陣由單向板彈性常數確定，破壞判據采用能量破壞判據、Hashin破壞判據。如同上述文獻，其它大多數有限元模型的單元彈性剛度矩陣由單向板彈性常數確定，雖然選用的破壞判據和單元應力更新方案存在一些區別，但這都需要事先實驗確定單向板動態彈性參數和強度參數。隨著層合板鋪設構造、纖維體積含量、沖擊速度等因素的變化，實驗將會變得相當煩瑣。細觀力學有限元模型直接通過兩種組分的物理性能參數確定層合板沖擊承載能力，僅需實驗確定纖維和基體的動態彈性參數和強度參數，能夠減少實驗的工作量。然而，目前的有限元模型較少采用細觀力學有限元模型，主要原因可能是細觀力學有限元模型^[15~17]用于解決層合板沖擊問題時存在計算量過大的難題。
    本文將復合材料細觀力學橋聯模型^[18]與有限元軟件ABAQUS結合，分析層合板受低速沖擊作用的極限沖擊響應。即將細觀力學本構模型、針對組分材料的破壞判據以及剛度衰減模式編制成用戶子程序VUMAT,為ABAQUS求解提供一種自定義材料模型。該模型僅需輸入纖維和基體的材料參數、纖維體積含量等有限數據并且無須單層板的實驗數據，通過確定組分材料是否破壞來判斷單層板是否破壞，對破壞后的單層實施一種常系數剛度衰減，就能使得復合材料層合板結構的沖擊承載能力分析得以順利實施。
2 有限元模擬
   纖維增強復合材料層合板受低速橫向沖擊作用，該沖擊系統包含沖擊頭、層合板和兩者之間的相互作用。層合板的單元選擇與低速沖擊的研究目的密切相關。如研究目的為沖擊后剩余強度,分層則是主要研究對象之一，選擇三維實體單元以及三維破壞判據是相當必要的。本文的研究目的為層合板低速沖擊作用下的極限承載能力,根據基體和纖維的平均應力判斷基體開裂、纖維斷裂，對破壞后的單  元組分材料的物理參數進行折減，而不在于識別分層破壞。一方面分層往往受面內基體開裂的誘發^[19],如果已經對發生基體開裂的層進行剛度衰減，就不需要因分層而再次折減該層剛度。另一方面分層往往出現在沖擊背面的受拉區^[11]，分層對復合材料面內受拉強度影響并不明顯,這種做法的實質是把分層吸收的能量包含在基體開裂中。基于上述理由并減少計算量，本文選擇殼單元。
    ABAQUS提供了4節點有限應變殼單元S4R^[20]，每個節點3個平移自由度、3個轉角自由度。為了獲得單層的應力和應變信息，每個單層采用Simpson 3點積分法，如圖1所示，在層合板中面布置節點并劃分網格，殼截面廣義力沿厚度積分得到。沖擊錘假定為剛性體，選擇R3D4剛性單元模擬。

              
    圖1中，節點1、2、3、4在殼參考平面上，層合板(5層)與S4R復合材料截面定義平面內1點積分；“x”表示厚度上Simpson積分點位置。[-page-] 
    纖維復合材料層合板雖屬于薄殼結構，仍需要適當考慮橫向剪切，而且對于一階剪切變形的殼單元，橫向剪切作用對動量守恒或虛功方程以及離散結點的動力方程有一定影響。不正確設置剪切剛度將導致收斂困難和計算時間漫長。對稱布置的層合板可視為正交各向異性材料，其橫向剪切剛度通過下面的簡式計算：

         
    其中，下標1,2,3為殼單元坐標方向；D₁₁，D₂₂，D₁₂為復合材料殼截面的抗彎剛度。
    復合材料層合板與沖擊頭之間的相互作用通過接觸對算法實現。接觸對中，沖擊頭為主面，受沖擊的層合板為從面，并指明與沖擊頭表面發生接觸的方向。另外從面網格細化須能滿足解的穩定性。由于不考慮兩者間的切向摩擦，選擇主從面動態強化約束接觸算法^[20]。
3 本構模型
3.1 應力更新方案
    橋聯模型細觀力學方法給出了準靜態荷載下基體、纖維與單向板的應力-應變之間的增量關系式，但在瞬態動力分析中需要材料率形式的本構方程。本文假定纖維、基體及復合材料為率無關材料。將準靜態性能簡單推廣到動態性能，纖維和基體以及單向復合材料在小變形下率形式的本構方程分別為：

       
    其中，“?”表示對時間求導；V_f和V_m分別為纖維和基體的體積含量；[S]、[S^f]、[S^m]依次為單向板、纖維和基體的柔度矩陣。類似文獻[18]的推導，單向復合材料應力-應變率形式的本構方程為：

         
    式中，[A_ij]為橋聯矩陣；[I_ij]為單位矩陣。由于在VUMAT中{△εi}已知，材料點在當前積分步的應力增量和總應力為：

              
3.2 破壞判據與剛度衰減
    采用大正應力破壞準則，當纖維或基體滿足下列不等式：

             
    則認為與之對應的單層板發生了破壞。其中,上標f,m分別為纖維和基體;下標1、2分別為分量坐標方向;u、c分別為極限拉伸和壓縮強度。
    單層板破壞后須進行剛度衰減，后續調用該單元時使用衰減的剛度矩陣。剛度衰減格式及應力更新策略采用經驗方法^[5,21]，針對組分材料進行。一旦纖維發生破壞，則認為該層完全喪失承擔后續荷載的能力，有：

           
    若基體發生破壞，認為基體開裂，但由于纖維和上下層的約束，該層還能夠繼續承受荷載作用，泊松比保持不變，只是纖維和基體的當前模量變為：

            
4 算例
    考慮一個周邊夾支的[0_n/90_m/0_n]層合板,2n＋m＝10層，密度為1.678kg?m^-3，有效直徑為160mm，單層厚為0.18mm,纖維體積含量為V_f＝0.43，中心承受橫向低速沖擊作用。文獻[22]給出的實驗數據：各向同性纖維，E_f＝74GPa、G_f＝30GPa、μ_f＝0.25；各向同性基體，E_m＝3GPa、G_m＝1.1GPa、μ_m＝0.4；沖擊頭質量2.3Kg，頭部直徑25mm的半圓球型。
    有限元模型見圖2，中心區域網格加密見圖3，圖中數字為單元編號。模型相對幾何中心對稱，材料性能則相對x、y軸對稱，取半圓進行分析，徑向設置為對稱邊界條件。

          
    給定沖擊能量1J，即沖擊速度為0.933m/s。從圖4可以看出，實驗值和預報值吻合較好，說明本文采用的細觀力學有限元模型計算有效。圖5給出了三種不同鋪層的層合板在沖擊載荷下的中點撓度，彼此差別不大，反映了目前工況下層合板承受沖擊作用與層合板的鋪設角影響不大，該結論與實驗數據一致。[-page-] 
    為了清晰反映基體應力沿徑向的變化趨勢，選取有代表性的單元3,1,80,144,192,具體位置見圖3。上述單元的纖維和基體平均應力變化曲線如圖6、7所示。

    圖6(a)、(b)表明,在沖擊背面0°鋪層內基體受雙向拉應力，沖擊作用點附近的應力峰值遠高于其它部分；隨著遠離沖擊中心，基體的應力峰值迅速變小；根據單元3,1,80,144,192的基體應力逐步變小的變化趨勢，該層距離中心比單元192單元更遠的區域，基體應力峰值較小。圖6(e)、(f)表明，沖擊正面內的基體受雙向壓應力，且應力峰值小于沖擊背面0°鋪層中的應力的峰值；隨著遠離沖擊中心，基體的應力峰值迅速變??；該層距離中心比192單元更遠的區域，基體應力峰值小。另外，根據基體拉壓性能的差異，容易判斷沖擊背面將先出現基體拉伸破壞。圖6(c)、(d)表明90°鋪層內的應力值峰值小且隨徑向變化平緩、應力梯度?。黄渌鼌^域基體應力峰值小。
    纖維平均應力隨時間變化曲線呈現的規律與基體類似。圖7(a)、(b)表明，在沖擊背面0°鋪層內纖維受雙向拉應力，沖擊作用點附近的應力峰值遠高于其它部分；隨著遠離沖擊中心，纖維的應力峰值迅速變小；根據單元3,1,80,144,192的基體應力逐步變小的變化趨勢，該層距離中心比單元192單元更遠的區域，纖維應力峰值較小。圖7(e)、(f)表明，單元3,1,80沖擊正面內的纖維受雙向壓應力，且應力峰值小于沖擊背面0°鋪層中的應力的峰值；144、192單元應力隨時間變化的過程體現不同特點，先是壓應力而后是拉應力。如將沖擊作用等效為橫力彎曲，一般來說反彎點靠近支座，該點應力應該為壓應力，這從側面反映了沖擊作用與橫力彎曲的區別。圖7(c)、(d)表明90°鋪層內纖維的應力值峰值小且隨徑向變化平緩、應力梯度??；該層距離中心比單元192單元更遠的區域，纖維應力峰值很小。
    計算結果(圖4、5)表明，三種鋪層[0₄/90₂/0₄],[0₃/90₄/0₃],[0₂/90₆/0₂]橫向反作用合力曲線差別很小，與實驗結論一致。將沖擊能量提高到27J，組分材料的強度取值，基體拉伸強度為80MPa、壓縮強度120MPa，纖維拉伸強度為2150MPa、壓縮強度1450MPa，其余參數同前。達到破壞后，纖維和基體性能按(12)式進行折減，結果見圖8，可知計算值與實驗值比較吻合。

  
    隨著沖擊能量不斷提高，沖擊過程持續的時間越來越短，大沖擊合力越來越大，如圖9所示。當沖擊速度為10m/s時，沖擊作用合力-時間曲線發生劇變(圖10)，合力的峰值很小。究其原因，是因為大量單元的基體發生破壞，按(9)、(10)式衰減剛度，導致整個結構剛度急劇變小，計算結果并不合理，至少沖擊作用合力-時間曲線組沒有體現出連續性，這需要對剛度衰減和應力更新方案做進一步研究。盡管如此，從中還是可以判斷層合板承受9~10m/s沖擊力作用將發生較大程度的破壞。
5 結語
    本文利用復合材料細觀力學有限元模型對層合板低速沖擊進行了數值模擬，并與實驗數據進行對比。分析結果中，除了得到沖擊荷載時間曲線、中點位移時間曲線、鋪層對抗沖擊能力的影響和變形特點等常規結論之外，還直接得到纖維、基體平均應力隨時間變化及其在空間上的分布。綜上所述，本文所述細觀力學方法，能夠反應層合板受橫向低速沖擊響應主要特點，計算過程簡潔，具有一定的可行性。
         參考文獻
[1] Abrate S. Impact on laminated composites：recent Advance[J]．Ap-pl Mech Rev，1991,44（4）：155-190.
[2] Abrate S. Impact on laminated composites：rencent Advance[J]．Appl Mech Rev，1994,47(11);517-544.
[3] 陳浩然,白瑞祥.低速沖擊后含損復合材料夾層板剩余強度研究進展[J].力學進展,2002,32(32):402-413．
[4] 黃爭鳴，張華山．纖維增強復合材料強度理論的研究現狀與發展趨勢-“破壞分析奧運會”評估綜述[J]．力學進展,2007,37(1):80-98.
[5] Hou JP，Petrinic N，Ruiz R. Prediction of impact damage in com-posite plates[J]．Composite Science and Technology，2006,60:273-281．
[6] Chang F，Chang K. A progressive damage model for laminated com-posites containing stress concentrations[J]．Journal of Composite Material，1987，21：834-55．
[7] Brewer JC,Lagace PA. Quadratic stress criterion for initiation of delami-nation[J].Journal of Composite Materials，1988,22(12):1141-55.
[8] 程小全,酈正能.復合材料層合板低速沖擊后壓縮的損傷積累模型[J].應用數學和力學,2005,26(5):569-576.
[9] Chang FK，Springer GS. The strength of fiber reinforced composite bends[J]．Journal of Composite Materials，1986,20：30-45.
[10] Hashin Z. failure ctiteria for unidirectional fiber composites[J].Jourual of Applied Mechanics，1980,47:329-334.
[11] Aslan A，Karakuzu R，Okutan B．The response of laminated composite plates under low-velocity impact loading[J]．Compos-ite Structures, 2003,59：119-127.
[12] Choi HY，Wu HYT，Chang FK. A new approach toward understand-ing damage mechanisms and mechanics of laminated composites due to low-velocity impact：part II -analysis[J]．J Compos Mater，1991,25:1012-38．
[13] Choi HY,Chang FK. A mode for predicting damage in graphite/ep-oxy laminated composites resulting from low velocity point impact[J].Journal of Composite Materials,1992,26(14):2134-69．
[14] Johnson HE，Louca LA，Mouring SE. Current research into model-ling of shock damage to large scale composite panels[J]．Journal of Materials Science，2006，41：1573-1590.
[15] Mayes SJ，Hansen AC. Composite laminate failure analysis using multicontinuum theory[J]．Compos Sci Tech，2004，64;379-394.
[16] Gotsis PK，Chamis CC，Minnetyan L. Prediction of composite lami-nate fracture：micromechanics and progressive frature[J]．Compos Sci Technol，1998，58：1137-1150.
[17] 田峰,胡更開,楊改英.塑性細觀力學模型的比較分析[J].北京理工大學學報,1996,16(5):477-482.
[18] 黃爭鳴．復合材料細觀力學引論[M]．北京：科學出版社，2004；49-60.
[19] 徐穎，溫衛東，崔海坡．復合材料層合板低速沖擊逐漸積累損傷預測方法[J]．材料科學與工程學報,2006,24(1):77-81．
[20] Abaqus Analysis User Manual，version 6.4.
[21] Hou JP，Petrinic N，RuiZ C. A delamination criterion for laminated composites under low-velocity impact[J]．Composites Science and Technology，2001，61：2069-2074.
[22] Collombet F，Lalbin X,Latailade JL. Impact of laminated compos-ites：physical basis for finite element analysis[J]．Composi Sci and Tech，1998,58:463-478．