CN108256211B - 基于abaqus的木材本构关系数值模拟方法 - Google Patents

Abstract

本发明提供一种基于ABAQUS的木材本构关系数值模拟方法，属于木材本构研究技术领域。本发明的模拟方法通过嵌入自定义的木材本构子程序VUMAT，能准确描述木材的正交各向异性，并预测木材在复杂受力状态下的破坏模式。

Description

基于ABAQUS的木材本构关系数值模拟方法

技术领域

本发明属于木材本构研究技术领域，提供了木材在复杂应力状态下应力-应变关系的数值模拟方法。

背景技术

木材是各向异性材料，其本构的复杂主要表现为在受压作用下发生塑性变形，而在拉、剪作用下发生脆性破坏，同时各向拉、压强度均不相等。在木结构设计中，除极少数情况外(如单轴受拉、受压或受剪)，木构件受力是复杂的，一个能准确描述木材各向异性本构关系的模型是木结构有限元模型的重要组成部分。

在有限元分析中，木材通常可以被简化为正交各向异性材料。目前，ABAQUS内嵌的材料本构中，针对正交各向异性材料发生受压塑性变形的只有Hill屈服准则，但该本构模型不能同时描述木材顺纹理想塑性和横纹塑性硬化。

此外，ABAQUS内嵌的渐进损伤和失效模型尽管可以描述拉、剪作用下发生的脆性破坏，但会引起顺纹和横纹方向同时失效，这与木材的真实受力情况不符。综上所述，在目前的ABAQUS软件中，其内嵌的材料本构还不能准确描述木材在不同方向受压的塑性变形和拉、剪的脆性破坏。

发明内容

本发明提供了一种基于ABAQUS的木材本构关系数值模拟方法，通过嵌入自定义的木材本构子程序VUMAT，能准确描述木材的正交各向异性，并预测木材在复杂受力状态下的破坏模式。

本发明的技术方案：

基于ABAQUS的木材本构关系数值模拟方法，步骤如下：

步骤一：建立弹性阶段本构方程，表达弹性阶段正交各向异性材料的应力-应变关系，如式(1)所示：

σ＝Dε (1)

式中：σ为应力矩阵、ε为应变矩阵、D为刚度矩阵；

展开式(1)得到式(2)：

式中，D₁₁＝E₁(1-v₂₃v₃₂)γ，D₂₂＝E₂(1-v₁₃v₃₁)γ，D₃₃＝E₃(1-v₂₁v₁₂)γ，D₁₂＝E₁(v₂₁+ν₂₃ν₃₁)γ，D₁₃＝E₁(ν₂₁+ν₂₃ν₃₁)γ，D₂₃＝E₂(ν₃₂+ν₁₂ν₃₁)γ，D₄₄＝G₁₂，D₅₅＝G₂₃，D₆₆＝G₁₃，γ＝(1-ν₂₁ν₁₂-ν₂₃ν₃₂-v₁₃v₃₁-2v₂₁v₃₂v₁₃)^-1；

E₁，E₂，E₃分别为L、R、T三个方向的弹性模量；G₁₂，G₁₃，G₂₃分别为L-R、L-T和R-T平面内的剪切模量；v_ij为泊松比；

步骤二：建立强度准则

采用Yamada-Sun强度准则建立顺纹方向L、横纹径向R和切向T三向独立的强度准则，描述木材的正交各向异性，当满足式(3)时，材料进入非弹性阶段。

式中，X，Y，Z分别为L、R、T方向的抗拉或抗压强度，当σ₁₁，σ₂₂和σ₃₃为拉或压应力时选用相应的抗拉或抗压强度；S_ij(i,j＝x,y,z)为i-j平面内的抗剪强度；

步骤三：当木材达到强度准则后，进行非弹性阶段木材的本构关系的模拟，具体如下：

(一)建立损伤模型

木材在拉、剪作用下时发生脆性破坏，通过引入连续损伤因子d_i(i＝L,R,T)，对应力进行折减来模拟材料的应变软化；当材料出现应变软化时，则会出现应变局域化；此时随着网格的细化，能量的耗散将减少，有限元计算结果将严重依赖网格尺寸的大小，所以采用能量释放率来减小网格的依赖性，损伤因子计算如式(4)所示

式中，

为i方向抗拉强度；L^C为网格特征长度，G_i(i＝L,R,T)为i方向断裂能，L^C与G_i的比值即为断裂能释放率；

应力更新得到式(5)：

式中，

为无损应力；

(二)塑性发展

设定木材横纹受压时，发生应变硬化，设定木材顺纹受压时，发生少许硬化后进入理想弹塑性；

(1)建立硬化法则：

采用Ziegler随动硬化模型，控制屈服面由初始屈服面向最终破坏面的转移来实现硬化，屈服面表达式如式(6)所示：

式中，n_i(i＝L,R,T)为屈服强度与极限强度的比值，0.9≤n_L≤1；b_ij为背应力，增量db_ij的表达式如式(7)所示：

式中，C_α,i(i＝L,R,T)为硬化参数；

为屈服面转移约束方程，保证最终的屈服面不超过破坏面；

其中，

为最终屈服应力，由式(3)可以求得；计算式如式(8)所示：

为等效塑性应变，其表达式如式(9)所示：

(2)流动法则：

当应力达到屈服面后，通过弹塑性径向返回算法，将应力点约束在屈服面上，如式(10)所示：

式中，dε_ij、

和

分别为总应变增量、弹性应变增量和塑性应变增量；dλ为塑性比例因子；g为塑性势能函数，由于在本发明中采用关联塑性，令g为屈服函数f；

(A)硬化阶段塑性因子求解：

类推各向同性材料塑性因子推导方法，推出正交各向异性材料的塑性因子表达式为：

式中，

其中，h_L、h_R、h_T为塑性模量，由一致性条件求得

(B)理想弹塑性阶段塑性因子求解：

基于式(3)，屈服函数f的表达式为：

由塑性理论的一致性条件，可推得木材顺纹(L)、横纹径向(R)以及切向(T)的塑性因子，三者的表达式如式(15)所示：

式中

为试探应力，

(3)塑性应变的求解：

塑性应变求解的公式如下：

屈服面交线上点的塑性应变增量取在该点相交各面法线方向所确定增量的线性组合，即

其中A_k为线性组合系数。

本发明的有益效果：本发明通过嵌入自定义的木材本构子程序VUMAT，能准确描述木材的正交各向异性，并预测木材在复杂受力状态下的破坏模式。

附图说明

图1是木材L，R和T三个方向的定义图示。

图2是本发明的本构关系模型。图中：1顺纹抗拉；2横纹抗拉；3顺纹抗压；4横纹抗压；5木材抗剪；6弹性阶段；7应变软化；8理想弹塑性；9应变硬化；10材料断裂；X_C，Y_C和Z_C分别为L，R和T三个方向的抗压强度；X_T，Y_T和Z_T分别为L，R和T三个方向的抗拉强度；S为木材平均抗剪强度；n_i(i＝L,R,T)为三个方向的屈服强度与极限强度的比值。

具体实施方式

以下结合附图和技术方案，进一步说明本发明的具体实施方式。

1.弹性阶段本构方程

弹性阶段正交各向异性材料应力-应变关系式如式(1)所示：

σ＝Dε (1)

式中：σ为应力矩阵、ε为应变矩阵、D为刚度矩阵。

展开表达式如式(2)：

式中，D₁₁＝E₁(1-v₂₃v₃₂)γ，D₂₂＝E₂(1-ν₁₃v₃₁)γ，D₃₃＝E₃(1-v₂₁v₁₂)γ，D₁₂＝E₁(v₂₁+v₂₃v₃₁)γ，D₁₃＝E₁(v₂₁+v₂₃v₃₁)γ，D₂₃＝E₂(v₃₂+v₁₂v₃₁)γ，D₄₄＝G₁₂，D₅₅＝G₂₃，D₆₆＝G₁₃，γ＝(1-v₂₁v₁₂-v₂₃v₃₂-v₁₃v₃₁-2v₂₁v₃₂v₁₃)^-1。

E₁，E₂，E₃分别为L、R、T三个方向的弹性模量；G₁₂，G₁₃，G₂₃分别为L-R、L-T和R-T平面内的剪切模量；v_ij为泊松比；

2.强度准则

采用Yamada-Sun强度准则描述木材的正交各向异性，具体表达式如下式(3)所示：

式中X，Y，Z分别为L,R,T方向的抗拉或抗压强度，当σ₁₁，σ₂₂和σ₃₃为拉或压应力时选用相应的抗拉或抗压强度；S_ij(i,j＝x,y,z)为i-j平面内的抗剪强度。

3.非弹性阶段木材本构关系描述

木材达到强度准则后(如图2所示)，对于拉、剪作用下的脆性破坏通过引入连续损伤因子以及断裂能释放率来模拟。对于压力作用下的塑性变形，其流动法则由Yamada-Sun屈服函数推导得到，而硬化法则采用Ziegler随动硬化模型。顺纹受压作用时的应力-应变关系为理想弹塑性；横纹受压作用时发生应变硬化，当压应变达到极限压应变时，材料在该方向发生断裂，同时该方向失去承载力。

3.1损伤演化

木材在拉、剪作用下时发生脆性破坏，通过引入连续损伤因子d_i(i＝L,R,T)，对应力进行折减来模拟材料的应变软化。当材料出现应变软化时，会出现应变局域化。这时，随着网格的细化，能量的耗散将减少，有限元计算结果将严重依赖网格尺寸的大小。基于此，本文中采用了能量释放率来减小网格的依赖性，损伤因子计算如式(4)所示

式中，

为i方向抗拉强度；

L^C为网格特征长度，G_i(i＝L,R,T)为i方向断裂能，两者的比值即为断裂能释放率；

应力更新如式(5)所示：

式中

为无损应力。

3.2塑性发展

木材在受压时，发生塑性变形。当顺纹受压时为理想弹塑性；横纹受压时发生硬化。

3.2.1硬化法则

木材横纹受压时，发生应变硬化，顺纹受压为理想弹塑性，本发明为了便于求解背应力，假定顺纹受压也是有少许硬化后进入理想弹塑性的。本发明通过采用Ziegler随动硬化模型，控制屈服面由初始屈服面向最终破坏面的转移来实现硬化，屈服面表达式如式(6)所示：

式中，n_i(i＝L,R,T)为屈服强度与极限强度的比值，0.9≤n_L≤1；b_ij为背应力，增量db_ij表达式如式(7)所示：

式中，C_α,i(i＝L,R,T)为硬化参数；

为屈服面转移约束方程，保证最终的屈服面不超过破坏面；其中，

为最终屈服应力，由式(3)可以求得。计算式如式(8)所示：

为等效塑性应变，其表达式如式(9)所示：

3.2.2流动法则

当应力达到屈服面后，通过弹塑性径向返回算法，将应力点约束在屈服面上，如式(10)所示：

式中，dε_ij、

和

分别为总应变增量、弹性应变增量和塑性应变增量；dλ为塑性比例因子；g为塑性势能函数，由于在本发明中采用关联塑性，可令g为屈服函数f。

3.2.2.1硬化阶段塑性因子求解

类推各向同性材料塑性因子推导方法可推出正交各向异性材料的塑性因子表达式为：

式中，

其中，h_L、h_R、h_T为塑性模量，由一致性条件求得

3.2.2.2理想弹塑性阶段塑性因子求解

基于式(3)，屈服函数f的表达式为：

由塑性理论的一致性条件，可推得木材顺纹(L)、横纹径向(R)以及切向(T)的塑性因子，三者的表达式如式(15)所示：

式中

为试探应力，

3.2.3塑性应变的求解

塑性应变求解的公式如下：

屈服面交线上点的塑性应变增量取在改点相交各面法线方向所确定增量的线性组合，即

其中A_k为线性组合系数。

4.编写材料本构子程序

利用Fortran编写VUMAT材料子程序，通过ABAQUS提供的二次开发接口将VUMAT嵌入主程序中进行木材构件的数值模拟。

Claims (1)

1.一种基于ABAQUS的木材本构关系数值模拟方法，其特征在于，步骤如下：

步骤一：建立弹性阶段本构方程，表达弹性阶段正交各向异性材料的应力-应变关系，如式(1)所示：

σ＝Dε (1)

式中：σ为应力矩阵、ε为应变矩阵、D为刚度矩阵；

展开式(1)得到式(2)：

式中，D₁₁＝E₁(1-v₂₃v₃₂)γ，D₂₂＝E₂(1-v₁₃v₃₁)γ，D₃₃＝E₃(1-v₂₁v₁₂)γ，D₁₂＝E₁(v₂₁+v₂₃v₃₁)γ，D₁₃＝E₁(v₂₁+v₂₃v₃₁)γ，D₂₃＝E₂(v₃₂+v₁₂v₃₁)γ，D₄₄＝G₁₂，D₅₅＝G₂₃，D₆₆＝G₁₃，γ＝(1-v₂₁v₁₂-v₂₃v₃₂-v₁₃v₃₁-2v₂₁v₃₂v₁₃)^-1；

E₁，E₂，E₃分别为L、R、T三个方向的弹性模量；G₁₂，G₁₃，G₂₃分别为L-R、L-T和R-T平面内的剪切模量；v_ij为泊松比；

步骤二：建立强度准则

采用Yamada-Sun强度准则建立顺纹方向L、横纹径向R和切向T三向独立的强度准则，描述木材的正交各向异性，当满足式(3)时，材料进入非弹性阶段；

式中，X，Y，Z分别为L、R、T方向的抗拉或抗压强度，当σ₁₁，σ₂₂和σ₃₃为拉或压应力时选用相应的抗拉或抗压强度；S_ij(i，j＝x，y，z)为i-j平面内的抗剪强度；

步骤三：当木材达到强度准则后，进行非弹性阶段木材的本构关系的模拟，具体如下：

(一)建立损伤模型

木材在拉、剪作用下时发生脆性破坏，通过引入连续损伤因子d_i(i＝L，R，T)，对应力进行折减来模拟材料的应变软化；当材料出现应变软化时，则会出现应变局域化；此时随着网格的细化，能量的耗散将减少，有限元计算结果将严重依赖网格尺寸的大小，所以采用能量释放率来减小网格的依赖性，损伤因子计算如式(4)所示

式中，

为i方向抗拉强度；L^C为网格特征长度，G_i(i＝L，R，T)为i方向断裂能，L^C与G_i的比值即为断裂能释放率；

应力更新得到式(5)：

式中，

为无损应力；

(二)塑性发展

设定木材横纹受压时，发生应变硬化，设定木材顺纹受压时，发生少许硬化后进入理想弹塑性；

(1)建立硬化法则：

采用Ziegler随动硬化模型，控制屈服面由初始屈服面向最终破坏面的转移来实现硬化，屈服面表达式如式(6)所示：

式中，n_i(i＝L，R，T)为屈服强度与极限强度的比值，0.9≤n_L≤1；b_ij为背应力，增量db_ij的表达式如式(7)所示：

式中，C_α，i(i＝L，R，T)为硬化参数；

为屈服面转移约束方程，保证最终的屈服面不超过破坏面；

其中，

为最终屈服应力，由式(3)求得；计算式如式(8)所示：

为等效塑性应变，其表达式如式(9)所示：

(2)流动法则：

当应力达到屈服面后，通过弹塑性径向返回算法，将应力点约束在屈服面上，如式(10)所示：

式中，dε_ij、

和

分别为总应变增量、弹性应变增量和塑性应变增量；dλ为塑性比例因子；g为塑性势能函数，采用关联塑性，令g为屈服函数f；

(A)硬化阶段塑性因子求解：

类推各向同性材料塑性因子推导方法，推出正交各向异性材料的塑性因子表达式为：

式中，

其中，h_L、h_R、h_T为塑性模量，由一致性条件求得

(B)理想弹塑性阶段塑性因子求解：

基于式(3)，屈服函数f的表达式为：

由塑性理论的一致性条件，推得木材顺纹L、横纹径向R以及切向T的塑性因子，三者的表达式如式(15)所示：

式中：

为试探应力，

(3)塑性应变的求解：

塑性应变求解的公式如下：

屈服面交线上点的塑性应变增量取在该点相交各面法线方向所确定增量的线性组合，即

其中A_k为线性组合系数。

Images