本發(fā)明是關于復合材料多層厚板損傷預測領域，特別涉及預測低速沖擊下復合材料多層厚板漸進失效的有限元方法。

背景技術：

當前，復合材料多層厚板正廣泛應用于航空航天、風力發(fā)電、壓力容器、汽車等高新技術領域。但是低速沖擊損傷對復合材料多層厚板的強度、剛度及使用壽命都有較大的影響，因此必須清楚的了解在低速沖擊下復合材料多層厚板的漸進損傷演化過程。

現(xiàn)有的國內(nèi)外復合材料多層厚板低速沖擊研究都是基于二維的平板理論和三維逐層模擬方法。復合材料多層厚板面內(nèi)跨度與厚度尺寸相對較小，面外影響不可忽略，經(jīng)典的二維平板理論已經(jīng)無法適用。然而復合材料多層厚板往往因為層數(shù)較多且每層的單元尺寸很小，即使采用三維逐層模擬方法模擬低速沖擊載荷下的漸進失效無論是前處理時間還是有限元計算時間都較為龐大。目前能高效且較為準確的預測復合材料多層厚板的低速沖擊的方法幾乎空白，因此提出一種基于多尺度模型預測低速沖擊下復合材料多層厚板的漸進失效的方法尤為重要。

復合材料層合板強度計算分析一般可采用通用的有限元軟件ABAQUS進行，但是通過ABAQUS軟件自帶的模塊僅可以基于二維HASHIN失效準則采用逐層模擬方法對復合材料多層厚板進行低速沖擊下的漸進失效研究，但無法直接運用較有優(yōu)勢的三維PUCK失效準則，更無法直接采用多尺度模型預測低速沖擊下復合材料多層厚板的漸進失效。

技術實現(xiàn)要素：

本發(fā)明要解決的技術問題是，克服現(xiàn)有技術中的不足，提供一種預測低速沖擊下復合材料多層厚板漸進失效的有限元方法。

為解決上述技術問題，本發(fā)明的解決方案是：提供預測低速沖擊下復合材料多層厚板漸進失效的有限元方法，包括下述過程：

一、建立含沖錘、復合材料多層厚板以及支撐板的低速沖擊有限元模型；

二、建立復合材料多層厚板多尺度模型；

三、利用FORTRAN語言編寫的ABAQUS-VUMAT即ABAQUS用戶動態(tài)材料子程序模塊，基于多尺度模型求解低速沖擊載荷下應力、應變和損傷；

四、對低速沖擊進行計算，進一步獲得沖擊力、位移、速度和加速度；

所述過程一中，建立含沖錘、復合材料多層厚板以及支撐板的低速沖擊有限元模型：基于ABAQUS建立沖錘、復合材料多層厚板、支撐板部件。分別設置材料屬性和劃分網(wǎng)格，其中復合材料多層厚板在厚度方向網(wǎng)格可根據(jù)自身確定的子層壓板的個數(shù)劃分相應的網(wǎng)格，并不需要建立所有單層網(wǎng)格，不同子層壓板之間設置ABAQUS自帶的雙線性內(nèi)聚力單元。再用ABAQUS-ASSEMBLY模塊對其組裝之后設置分析步和通用接觸屬性；

所述過程二具體包括下述步驟：

步驟(1)：基于三維等效理論，建立子層壓板的三維等效剛度和應力-應變本構關系：

(a)子層壓板的三維等效剛度及其各個分量為：

其中，為子層壓板在整體坐標系下的等效剛度矩陣，

為子層壓板在整體坐標系下等效剛度矩陣分量，k為子層壓板的第k個單層，l為子層壓板的第l個單層，n為子層壓板的單層總數(shù)，V^k為第k個單層厚度占整個子層壓板的厚度百分比，V^l為第l個單層厚度占整個子層壓板 的厚度百分比，為第k個單層在整體坐標系下的剛度矩陣分量，為第l個單層在整體坐標系下的剛度矩陣分量。

(b)子層壓板的應力應變本構關系：

其中，為子層壓板在整體坐標系下的等效應力，為子層壓板在整體坐標系的等效應變。

步驟(2)：子層壓板應力應變分解為子層壓板內(nèi)各個單層的應力應變：

其中，

為子層壓板中第k個單層在整體坐標系下的應變分量， 為子層壓板在整體坐標系的應變分量，為子層壓板中第k個單層在整體坐標系下的應力分量，為子層壓板在整體坐標系的應力分量，為第k個單層在整體坐標系下的剛度矩陣分量。

步驟(3)：建立基于應變描述的PUCK失效初始判據(jù)和損傷演化準則，具體建立方式為：

(c)對于纖維拉伸和壓縮，損傷初始判據(jù)為：

其中，和分別為纖維拉伸和壓縮的初始失效應變；所述T,C分別指拉伸和壓縮；所述ε₁₁是指纖維方向應變；所述分別指纖維拉伸和壓縮失效判斷因子；

纖維拉伸和壓縮的損傷演化準則為：

其中，所述是指纖維拉伸和壓縮損傷變量；所述ε₁₁是指纖維方向應變；所述 是指纖維損傷變量達到1的纖維臨界拉伸和壓縮失效應變；所述是指纖維損傷變量為零的纖維初始拉伸和壓縮失效應變；

(d)對于基體拉伸損傷失效初始判據(jù)為：

其中，所述指基體拉伸失效判斷因子；所述ε₂₂是指基體方向應變；所述是指基體損傷變量為零的基體初始拉伸失效應變；

基體拉伸損傷演化準則為:

其中，為基體損傷變量達到1時基體臨界拉伸失效應變；所述是指基體拉伸損傷變量；

(e)對于基體壓縮損傷初始判據(jù)為：

其中，N是關于失效斷裂面的法向方向，T和L是關于失效斷裂面的切向方向；Y^C是橫向壓縮強度，斷裂平面上的應力σ_ij(i,j＝L,T,N)由笛卡氏坐標系下的Piola-Kirchhoff應力σ_ij(i,j＝1,2,3)通過旋轉矩陣T(α)旋轉獲得，T(α)為笛卡爾坐標系到斷裂面坐標系的旋轉矩陣，是斷裂面上的橫向剪切強度；所述σ_NN是指斷裂面的法向應力；σ_NT,σ_NL是指斷裂面的切向應力，μ_NL,μ_TN為斷裂面面內(nèi)兩個切向方向摩擦系數(shù)，θ_f為斷裂面的斷裂角；所述指基體壓縮失效判斷因子；所述σ₁₂是指笛卡爾坐標系下面內(nèi)剪切強度；所述是指在斷裂平面內(nèi)的橫向剪切強度；所述σ₁₂₃是指在笛卡爾坐標系下的六個Piola-Kirchhoff應力σ_ij(i,j＝1,2,3)；所述σ_LTN是指在斷裂面坐標系下的六個Piola-Kirchhoff應力σ_ij(i,j＝L,T,N)；所述T(α)^T是指T(α)的轉置矩陣；所述90°是指采用角度制計量的90度；

基體壓縮損傷演化準則為：

其中，所述γ_γ是指斷裂面聯(lián)合剪切應變；和是聯(lián)合剪切應變的初始和最大應變，γ_NT和γ_NL是斷裂面的剪切應變；是指基體壓縮損傷變量；

所述過程三具體包括下述步驟：

步驟(4)：通過ABAQUS用戶子程序VUMAT實現(xiàn)上述過程，ABAQUS-VUMAT是由應變驅動，隨著載荷增加，首先讀取各單層的材料信息和鋪層信息的參數(shù)代碼文件，計算各子層壓板內(nèi)各個單層在局部材料坐標系下的剛度矩陣，再轉換到整體坐標系下的單層剛度矩陣；

步驟(5)：再通過步驟(1)求取各子層壓板的三維等效剛度和子層壓板的應力應變；將子層壓板的應力應變通過步驟(2)分解為子層壓板內(nèi)各個單層應力應變，并轉換為局部材料坐標系下的單層應力應變。再根據(jù)步驟(3)判斷每個單層是否出現(xiàn)損傷以及求取損傷變量。如果沒有出現(xiàn)損傷，則更新應力應變后進入下一個增量步；如果出現(xiàn)損傷，則返回到步驟(4)通過損傷變量更新每個單層在整體坐標系下的剛度矩陣，再通過步驟(1)求取折減后的子層壓板的三維等效剛度，計算子層壓板三維等效剛度折減系數(shù)，并存儲為SDV狀態(tài)變量，更新應力應變后再進入下一個增量步；

所述過程四具體為：將過程一建立的模型主文件和過程三建立的ABAQUS-VUMAT用戶子程序聯(lián)合，使用ABAQUS/EXPLICT顯示方法對低速沖擊進行計算，進一步獲得沖擊力、位移、速度和加速度；即完成低速沖擊載荷下復合材料多層厚板漸進失效特性的預測。

與現(xiàn)有技術相比，本發(fā)明的有益效果是：

本發(fā)明利用ABAQUS用戶子程序界面的數(shù)值實現(xiàn)所建立的基于多尺度模型預測復合材料多層厚板低速沖擊下漸進失效的方法，該方法將復合材料多層厚板劃分為少量若干子層壓板，建立子層壓板和該子層壓板內(nèi)每個單層的轉化關系，相比于逐層模擬方法，本發(fā)明將大大減少復合材料多層厚板在低速沖擊下的計算時間，能高效且較為準確的預測復合材料多層厚板在低速沖擊載荷下的漸進失效行為。

附圖說明

圖1為本發(fā)明實施例復合材料多層厚板在低速沖擊下有限元模型圖。

圖2為本發(fā)明對所提出的基于多尺度模型VUMAT數(shù)值實現(xiàn)流程圖。

圖3為實施例中復合材料多層厚板在25J能量低速沖擊下沖擊力-時間數(shù)值模擬結果與實驗結果對比示意圖。

圖4為實施例中復合材料多層厚板在25J能量低速沖擊載荷下沖擊力-位移數(shù)值模擬結果與實驗結果對比示意圖。

圖5為實施例中復合材料多層厚板在25J能量低速沖擊載荷下沖擊力-時間數(shù)值模擬結果與實驗結果對比示意圖。

圖6為實施例中復合材料多層厚板在25J能量低速沖擊載荷下沖擊力-位移數(shù)值模擬結果與實驗結果對比示意圖。

具體實施方式

首先需要說明的是，本發(fā)明是計算機技術在復合材料多層厚板損傷預測領域的一種應用。在本發(fā)明的實現(xiàn)過程中，會涉及到多個軟件功能模塊的應用。申請人認為，如在仔細閱讀申請文件、準確理解本發(fā)明的實現(xiàn)原理和發(fā)明目的以后，在結合現(xiàn)有公知技術的情況下，本領域技術人員完全可以運用其掌握的軟件編程技能實現(xiàn)本發(fā)明。凡本發(fā)明申請文件提及的均屬此范疇，申請人不再一一列舉。

本發(fā)明中，預測低速沖擊下復合材料多層厚板漸進失效的有限元方法，包括下述過程：

一、建立含沖錘、復合材料多層厚板以及支撐板的低速沖擊有限元模型；

二、建立復合材料多層厚板多尺度模型；

三、利用FORTRAN語言編寫的ABAQUS-VUMAT即ABAQUS用戶動態(tài)材料子程序模塊，基于多尺度模型求解低速沖擊載荷下應力、應變和損傷；

四、對低速沖擊進行計算，進一步獲得沖擊力、位移、速度和加速度；

所述過程一中，包括基于ABAQUS建立沖錘、復合材料多層厚板、支撐板部件；分別設置材料屬性和劃分網(wǎng)格，其中復合材料多層厚板在厚度方向網(wǎng)格根據(jù)自身確定的子層壓板的個數(shù)劃分相應的網(wǎng)格，不需要建立所有單層網(wǎng)格，不同子層壓板之間設置ABAQUS自帶的雙線性內(nèi)聚力單元；再用ABAQUS-ASSEMBLY模塊對其組裝之后設置分析步和通用接觸屬性；

所述過程二具體包括下述步驟：

步驟(1)：基于三維等效理論，建立子層壓板的三維等效剛度和應力-應變本構關系：

(a)子層壓板的三維等效剛度及其各個分量為：

其中，為子層壓板在整體坐標系下的等效剛度矩陣，

為子層壓板在整體坐標系下等效剛度矩陣分量，k為子層壓板的第k個單層，l為子層壓板的第l個單層，n為子層壓板的單層總數(shù)，V^k為第k個單層厚度占整個子層壓板的厚度百分比，V^l為第l個單層厚度占整個子層壓板的厚度百分比，為第k個單層在整體坐標系下的剛度矩陣分量，為第l個單層在整體坐標系下的剛度矩陣分量；

(b)子層壓板的應力應變本構關系：

其中，為子層壓板在整體坐標系下的等效應力，為子層壓板在整體坐標系的等效應變；

步驟(2)：子層壓板應力應變分解為子層壓板內(nèi)各個單層的應力應變：

其中，

為子層壓板中第k個單層在整體坐標系下的應變分量， 為子層壓板在整體坐標系的應變分量，為子層壓板中第k個單層在整體坐標系下的應力分量，為子層壓板在整體坐標系的應力分量，為第k個單層在整體坐標系下的剛度矩陣分量；

步驟(3)：建立基于應變描述的PUCK失效初始判據(jù)和損傷演化準則，具體建立方式為：

(c)對于纖維拉伸和壓縮，損傷初始判據(jù)為：

其中，和分別為纖維拉伸和壓縮的初始失效應變；所述T,C分別指拉伸和壓縮；所述ε₁₁是指纖維方向應變；所述分別指纖維拉伸和壓縮失效判斷因子；

纖維拉伸和壓縮的損傷演化準則為：

其中，所述是指纖維拉伸和壓縮損傷變量；所述ε₁₁是指纖維方向應變；所述 是指纖維損傷變量達到1的纖維臨界拉伸和壓縮失效應變；所述是指纖維損傷變量為零的纖維初始拉伸和壓縮失效應變；

(d)對于基體拉伸損傷失效初始判據(jù)為：

其中，所述指基體拉伸失效判斷因子；所述ε₂₂是指基體方向應變；所述是指基體損傷變量為零的基體初始拉伸失效應變；

基體拉伸損傷演化準則為：

其中，為基體損傷變量達到1時基體臨界拉伸失效應變；所述是指基體拉伸損傷變量；

(e)對于基體壓縮損傷初始判據(jù)為：

其中，N是關于失效斷裂面的法向方向，T和L是關于失效斷裂面的切向方向；Y^C是橫向壓縮強度，斷裂平面上的應力σ_ij(i,j＝L,T,N)由笛卡氏坐標系下的Piola-Kirchhoff應力σ_ij(i,j＝1,2,3)通過旋轉矩陣T(α)旋轉獲得，T(α)為笛卡爾坐標系到斷裂面坐標系的旋轉矩陣，是斷裂面上的橫向剪切強度；所述σ_NN是指斷裂面的法向應力；σ_NT,σ_NL是指斷裂面的切向應力，μ_NL,μ_TN為斷裂面面內(nèi)兩個切向方向摩擦系數(shù)，θ_f為斷裂面的斷裂角；所述指基體壓縮失效判斷因子；所述σ₁₂是指笛卡爾坐標系下面內(nèi)剪切強度；所述是指在斷裂平面內(nèi)的橫向剪切強度；所述σ₁₂₃是指在笛卡爾坐標系下的六個Piola-Kirchhoff應力σ_ij(i,j＝1,2,3)；所述σ_LTN是指在斷裂面坐標系下的六個Piola-Kirchhoff應力σ_ij(i,j＝L,T,N)；所述T(α)^T是指T(α)的轉置矩陣；所述90°是指采用角度制計量的90度；

基體壓縮損傷演化準則為：

其中，所述γ_γ是指斷裂面聯(lián)合剪切應變；和是聯(lián)合剪切應變的初始和最大應變，γ_NT和γ_NL是斷裂面的剪切應變；是指基體壓縮損傷變量；

所述過程三具體包括下述步驟：

步驟(4)：通過ABAQUS用戶子程序VUMAT實現(xiàn)上述過程，ABAQUS-VUMAT是由應變驅動，隨著載荷增加，首先讀取各單層的材料信息和鋪層信息的參數(shù)代碼文件， 計算各子層壓板內(nèi)各個單層在局部材料坐標系下的剛度矩陣，再轉換到整體坐標系下的單層剛度矩陣；

步驟(5)：再通過步驟(1)求取各子層壓板的三維等效剛度和子層壓板的應力應變；將子層壓板的應力應變通過步驟(2)分解為子層壓板內(nèi)各個單層應力應變，并轉換為局部材料坐標系下的單層應力應變；再根據(jù)步驟(3)判斷每個單層是否出現(xiàn)損傷以及求取損傷變量；如果沒有出現(xiàn)損傷，則更新應力應變后進入下一個增量步；如果出現(xiàn)損傷，則返回到步驟(4)通過損傷變量更新每個單層在整體坐標系下的剛度矩陣，再通過步驟(1)求取折減后的子層壓板的三維等效剛度，計算子層壓板三維等效剛度折減系數(shù)，并存儲為SDV狀態(tài)變量，更新應力應變后再進入下一個增量步；

所述過程四具體為：將過程一建立的模型主文件和過程三建立的ABAQUS-VUMAT用戶子程序聯(lián)合，使用ABAQUS/EXPLICT方法對低速沖擊進行計算，進一步獲得沖擊力、位移、速度和加速度；即完成低速沖擊載荷下復合材料多層厚板漸進失效特性的預測。

以下用基于上述過程應用的具體實施例子來說明本發(fā)明的實現(xiàn)效果：

在ABAQUS/CAE中建立包括碳纖維/環(huán)氧樹脂復合材料多層厚板，沖錘以及支撐板的低速沖擊有限元模型，如圖1所示。其中復合材料多層厚板有兩種鋪層形式：

和平板大小為150×100×4mm，密度為

ρ＝1600kg/m³，總共分為16個單層，依次將每四個單層劃分為一個子層壓板，所以在厚度方向共有四個子層壓板，只需要在厚度方向劃分四個單元(每個單元1mm)。每個子層壓板均用減縮積分三維實體單元C3D8R來仿真，內(nèi)聚力單元層鋪設在纖維方向不同的子層之間研究分層損傷，沖錘前端建成半球形，沖錘的質(zhì)量為2kg,直徑為1.6cm。沖錘的沖擊能量分別為25J。

層合板下方有一個長方形撐框架模型，框架外尺寸與平板相同，并帶有125×75mm的長方形開孔,沖錘和支撐框架的剛度和泊松比分別為

207.6GPa和0.26，密度為ρ＝7830kg/m³，框架的下表面固定。

利用ABAQUS/EXPLICT顯示計算模擬沖錘沖擊多層厚板的過程，如圖2利用用戶子程序VUMAT首先求取子層壓板的等效剛度，得到平板在整體坐標系下的應力應變，并分解到每個單層判別損傷并求取相應的剛度比折減系數(shù)，并存儲到狀態(tài)變量SDV中，再更新應力應變。基于上述過程求取多層厚板漸進失效過程。

圖3和圖4分別為在25J沖擊能量下鋪層的沖擊力-時間和沖擊力-位移的曲線圖，圖5和圖6分別為在25J沖擊能量下的沖擊力-時間和沖擊力-位移的曲線圖，均與試驗值準確的吻合。所以本發(fā)明提出的多尺度模型可以較為準確的預測復合材料多層厚板在低速沖擊下的漸進損傷失效特征。

本發(fā)明在ABAQUS軟件的基礎上進行用戶子程序的開發(fā)，提出的多尺度模型建立子層壓板和該子層壓板內(nèi)每個單層的轉化關系，大大減少了計算量，能高效且較為準確預測復合材料多層厚板在低速沖擊下的漸進損傷失效，為深入闡明復合材料結構的損傷失效特性、提升輕量化強度設計水平提供了技術支撐。

最后，需要注意的是，以上列舉的僅是本發(fā)明的具體實施例。顯然，本發(fā)明不限于以上實施例，還可以有很多變形。本領域的普通技術人員能從本發(fā)明公開的內(nèi)容中直接導出或聯(lián)想到的所有變形，均應認為是本發(fā)明的保護范圍。