[發明專利]預測低速沖擊下復合材料層合板漸進失效的有限元方法

有效

• 申請號: 201610833225.4
• 申請日: 2016-09-20
• 公開（公告）號: CN106503292B
• 公開（公告）日: 2018-04-24
• 發明（設計）人: 廖斌斌;劉鵬飛;鄭津洋
• 申請（專利權）人: 浙江大學
• 主分類號: G06F17/50
• 分類號: G06F17/50
• 代理公司: 杭州中成專利事務所有限公司33212
• 代理人: 周世駿

1.預測低速沖擊下復合材料層合板漸進失效的有限元方法，其特征在于，包括下述過程：

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

二、建立復合材料彈塑性損傷本構模型；

三、基于ABAQUS-VUMAT用戶動態材料子程序模塊，運用后向歐拉算法實現提出的彈塑性損傷本構模型，求解應力、應變和損傷；

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

所述過程一中，建立含沖錘、復合材料層合板以及支撐板的低速沖擊有限元模型：基于ABAQUS建立沖錘、復合材料層合板、支撐板部件，分別設置材料屬性和劃分網格，再進行ASSEMBLY對其組裝之后設置分析步和通用接觸屬性；

其中，復合材料層合板每層根據鋪層角度進行鋪層，不同鋪層之間設置ABAQUS自帶的雙線性內聚力單元；

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

步驟(1)：建立含損傷復合材料層合板層內本構關系；

復合材料應力-應變本構方程：S＝C_d：E^e；

其中，為未損傷材料的有效應力；S為損傷材料的名義應力；E^p為格林-拉格朗日應變張量的塑性部分，E^e＝E-E^p為格林-拉格朗日應變張量的彈性部分；C_d(C，d₁，d₂)為含損傷材料的四階彈性張量，C為未損傷材料的四階彈性張量，d₁和d₂是分別對應纖維斷裂和基體開裂的損傷變量；

步驟(2)：建立塑性模型，具體建立方式為：

(a)考慮面外應力的復合材料塑性流動準則為：

F(S‾,E~p)=3f(S‾)-k(E~p)=0;]]>

f(S‾ij)=12[(S‾22-S‾33)2+4S‾232+2a66(S‾132+S‾122)];]]>

其中，F為屈服方程；a₆₆為描述各向異性材料與塑性相關的常數，由偏軸拉伸測試確定；f為塑性勢函數；k和分別為塑性硬化應力和等效塑性應變；所述是指未損傷材料的有效應力，其中i，j用來確定應力方向；所述指層合板面內垂直于纖維方向的有效應力；所述指面外垂直于纖維的有效應力；所述指面外剪切有效應力；所述指面內剪切有效應力；

(b)塑性硬化應力k為：

k=β(E~p)n;]]>

其中，β和n是為滿足實驗硬化曲線的常數；塑性變形假定發生在損傷材料的未損傷區域，塑性流動準則和硬化準則表達在有效應力空間中；

(c)等效塑性應變率和塑形應變率分別為：

E~·p=λ·;]]>

E·p=λ·∂F/∂S‾;]]>

其中，為一致塑性因子；

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

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

F11T*=(E11E0,1T*)2-1≥0,F11C=(E11E0,1C)2-1≥0;]]>

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

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

d11T*(C)=Ef,1T*(C)Ef,1T*(C)-E0,1T*(C)(1-E0,1T*(C)E11);]]>

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

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

F22T*=(E22E0,2T*)2-1≥0;]]>

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

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

d22T*=Ef,2T*Ef,2T*-E0,2T*(1-E0,2T*E22);]]>

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

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

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

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

d22C(γNT,γNL)=γγmaxγγmax-γγf(1-γfγγ),γγ=γNT2+γNL2.;]]>

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

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

步驟(4)：通過用戶子程序VUMAT的SDV定義第n+1增量步的開始時的初始狀態變量值，同時也是第n增量步結束時的狀態變量值在第n+1增量步開始時，VUMAT讀入；

其中，所述n是指第n增量步，所述E_n是指第n增量步結束時的格林-拉格朗日總應變張量，所述是指第n增量步結束時的格林-拉格朗日塑性應變張量，所述是指第n增量步結束時的等效塑性應變，所述是指第n增量步結束時的未損傷材料的有效應力，所述k_n是指第n增量步結束時的塑性硬化應力，所述d_ij，n是指第n增量步結束時的損傷變量；

步驟(5)：VUMAT由應變增量驅動，計算試應力，將試應力代入到步驟(2)屈服方程中；

該公式是步驟(2)中的公式在第n+1增量步的特定計算；

如果F_n+1≤0，則處于彈性階段，將所有試應力和應變更新為n+1增量步狀態變量

如果F_n+1＞0，塑性加載出現，根據后向歐拉隱式算法，實現試應力到屈服面的最近點返回，將試應力賦予步驟(4)中的迭代初始條件；變量是塑性一致因子的增量Δλ_n+1的函數，運用牛頓-拉夫森算法求解再更新應變和有效應力，直到F_n+1≤0，結束迭代，得到n+1增量步狀態變量

步驟(6)：應變和有效應力得到更新后，根據有效應力和應變代入到步驟(3)的PUCK失效準則判斷是否出現損傷，如有損傷再通過損傷演化公式獲得損傷變量，再根據步驟(1)通過有效應力和損傷變量計算名義應力S_n+1；

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