FRACTURE BEHAVIOR OF PERIODIC POROUS STRUCTURES BY PHASE FIELD METHOD
-
摘要: 周期性多孔结构具有质量轻、比密度低、比强度高、隔音等优良特点, 同时也能很好地满足结构-功能一体化的需求, 在许多领域具有广泛的应用前景. 目前, 对周期性多孔结构在复杂载荷下的力学响应和断裂行为的研究较少. 采用细观力学和相场方法相结合, 基于二维代表性体积单元RVE模型, 施加能实现比例加载的周期性边界条件, 研究周期性多孔结构在复杂多轴比例加载状态下的裂纹萌生位置、断裂模式、承载极限及其变化规律. 本文的数值模拟结果表明: 首先, 周期性多孔结构在竖直方向拉伸载荷作用下, 裂纹均从孔边萌生并沿水平方向同步扩展; 其次, 在双轴载荷作用下, 随着水平载荷的增加, 结构在竖直方向的极限拉伸载荷逐渐增大; 当双轴拉伸载荷等值时, 结构的抗拉强度达到最大, 此时断裂模式呈现为十字正交型开裂; 最后, 面内剪切应力的引入会导致结构的拉伸强度极限降低, 孔边裂纹的萌生位置和扩展路径发生偏移, 裂纹模式从单S型转变为双弧线型, 裂纹向水平位置上相邻的孔洞扩展. 随着水平载荷的增加, 裂纹模式最终转变为斜裂纹, 从孔边对角线位置萌生并沿着45°方向扩展.Abstract: Periodic porous structures have excellent characteristics such as low mass, low density, high specific strength, sound insulation, and they are also well satisfy the needs for structural-functional integration, which have a wide range of applications in many fields. At present, the mechanical response and fracture behavior of periodic porous structures under complex loads have been poorly investigated. In this paper, we use a combination of micro-mechanics method and phase field method to investigate the crack initiation location, crack propagation path, fracture mode and the ultimate strength of periodic porous structures under combined multiaxial loading based on a two-dimensional representative volume element (RVE) model with the periodic boundary condition (PBC) that can implement multiaxial proportional loading. Numerical simulation results in this paper show that all cracks in the periodic porous structure established in this paper initiate from the edge of the holes and propagate consequently along the horizontal direction under the uniaxial tensile loading in the vertical direction. Secondly, under the biaxial loadings in both vertical and horizontal directions, the ultimate strength of the periodic porous structure gradually increases with the increase of the horizontal tensile loading. When the horizontal load is equal to the vertical load, the fracture pattern exhibits as orthogonal cross-type cracking and the ultimate strength reaches the maximum value. Thirdly, the in-plane shear stress simultaneously acted on the RVE model of the periodic porous structure results in a significant decrease of the ultimate strength and the variations of initiation location and propagation trajectory of the hole-edge cracks. Hence, the fracture pattern of periodic porous structure subjected to combined multiaxial loadings changes from the single S-type cracking to the double arc-type cracking, and cracks extend toward adjacent holes in horizontal direction. Finally, with the increase of the horizontal tensile loading, cracks initiate diagonally at the edge of the holes and propagate along the 45-degree direction which lead to the oblique cracking of periodic porous structure.
-
Keywords:
- periodic porous structure /
- phase field /
- fracture /
- RVE model /
- multiaxial loading
-
引 言
若材料在空间上由重复的且周期性的代表性单元(单胞)组成, 则称之为周期性结构[1]. 多孔结构是一类由随机或周期性的微观结构组成的轻质结构, 其材料主要包括金属、高分子以及陶瓷3大类[2]. 随机多孔结构又可称为泡沫结构, 周期性多孔结构常包括蜂窝结构、波纹结构等. 因多孔结构独特多样的性能特点, 逐渐成为了诸多领域学者研究的重点和热点. 多孔结构具有隔热、轻质、吸音、相对密度低、比强度高、可变形等优异的性能, 在航空航天、柔性电子器件、能量吸收、隔音和组装框架等领域有着广泛的应用[3].
断裂破坏问题的预测一直以来都是学术界和工程界关注的难题. 由于内部存在相互贯通或封闭的孔洞和不同的外部加载模式, 周期性多孔结构的断裂行为更加具有不确定性. Hayes等[4]探讨了单胞孔洞为正方形的周期性多孔合金在复杂加载模式下的断裂失效行为, 杜映洪[5]给出了分层蜂窝周期性多孔结构在复杂应力状态下的断裂模式, Jelitto等[6]探究了泡沫多孔材料断裂韧性、裂纹扩展与孔隙率、加载方式的关系, 均说明了多孔材料断裂时的裂纹扩展路径会受到孔洞影响, 且随着外部加载模式的改变, 裂纹扩展行为与裂纹数量也发生改变. 因此可以发现, 周期性多孔结构在复杂力学环境下, 会呈现出复杂的力学响应与失效模式[3]. 经典断裂力学理论难以解决裂纹形核以及扩展方向等问题[7-8], 因此迫切需要寻找合适的断裂问题数值模拟方法.
传统有限元框架下模拟裂纹扩展的数值分析方法主要有单元删除法[9]、界面单元法[10]、扩展有限元(XFEM)[11-12]等. 单元删除法将满足条件的单元应力置0, 但难以模拟裂纹分岔问题. 界面单元法通过设置内聚力单元模拟裂纹扩展, 但具有较强的网格依赖性. 扩展有限元法通过扩充形函数使裂纹可以在网格内扩展, 但较难处理三维问题[13]. 本文所采用的相场(phase-field)断裂模型是一种弥散式裂纹模型. 该方法基于传统Griffith[14]理论, 通过能量平衡理论研究裂纹的扩展行为. Francfort等[15]在该理论基础上提出了断裂变分准则, Bourdin等[16]通过引入弥散的相边界来表征裂纹的尖锐边界, 引入序参量(即相场)来得到空间中描述损伤的标量场, 并经过相场方程控制序参量的演化来显式追踪裂纹扩展路径, 极大地简化了算法实现的复杂度. 目前, 相场断裂方法已经广泛应用于脆性断裂[17]、塑性断裂[18-19]、动态断裂[20]、疲劳裂纹扩展[21]、各向异性材料断裂[22]、多场耦合[23]等问题. 该方法在三维模拟方面也体现出了优势[24]. 对于复合材料、周期性点阵、周期性多孔结构等可能出现的裂纹轨迹难以预测、多裂纹交汇问题, 相场法具有其独特的优势[25].
基于多尺度分析方法和均匀化理论[26], 周期性结构宏观模型的力学行为可以通过分析细观尺度的代表性体积单元(representative volume element, RVE)模型的力学性能来表征, 细观力学分析方法如图1所示. 在多尺度分析过程中, 需要满足基于均匀化理论的细−宏观能量等价条件(hill-mandel条件)[27]. 此外, 需要给RVE模型施加周期性边界条件, 来保证RVE边界处的变形、应力连续, 由宏观模型向微细观模型传递信息[28].
尽管多孔结构在单轴、双轴或剪切载荷下的力学响应已经被广泛研究[29-30], 但在实际情况中, 周期性多孔结构往往承受多轴复杂载荷, 目前对于该方面断裂问题的研究较少. 此外, 以往的大多数周期性边界条件只是控制RVE模型中对称边界的位移, 但不能保证作用在RVE模型上的多轴宏观应力的比值恒定不变.
本文基于ABAQUS有限元软件平台建立了周期性多孔结构的RVE模型, 施加以多轴加载比例为控制参数的周期性边界条件, 实现了对于RVE模型的比例加载, 结合相场断裂方法从而研究周期性多孔结构在复杂多轴加载状态下的断裂失效行为.
1. 相场断裂方法
1.1 断裂变分准则
首先在本文中定义
$ \varphi \in [0,1] $ 的标量来定义相场变量表示裂纹拓扑. 当$ \varphi = 0 $ 时表示材料完好,$ \varphi = 1 $ 时表示材料完全破坏, 引入一维杆件的弥散裂纹模型与相场函数[31], 如图2(a)和图2(b)所示, 其裂纹表面密度函数为$$ \varphi (x) = \exp ( - {{\left| x \right|} \mathord{\left/ {\vphantom {{\left| x \right|} {{l_0}}}} \right. } {{l_0}}}) $$ (1) 由此推导出二维和三维的表面密度函数
$$ \gamma (\varphi ,\nabla \varphi ) = \frac{{{\varphi ^2}(x)}}{{2{l_0}}} + \frac{{{l_0}}}{2}\left| {\nabla \varphi (x)} \right| $$ (2) 式中,
$ {l_0} $ 代表裂纹的弥散宽度.断裂变分原理认为弹性体
$ \varOmega $ 的势能$ {\varPi _\varOmega } $ 由弹性应变能$ {E_{{\rm{str}}}} $ 和断裂表面能$ {W_{{\rm{frac}}}} $ 组成, 其表达式为$$ \begin{split} {\varPi _\varOmega }({\boldsymbol{u}},\varGamma ) =& {E_{{\rm{str}}}} + {W_{{\rm{frac}}}} = \\ & \int_\varOmega {{\psi _e}} {\boldsymbol{\varepsilon (u)}}{\rm{d}}\varOmega + \int_\varGamma {{G_c}} {\rm{d}}\varGamma \end{split} $$ (3) 式中,
$ {\psi _e} $ 为弹性应变能密度,$ {\boldsymbol{\varepsilon (u)}} $ 为应变张量,$ {G_c} $ 为临界能量释放率.由于在损伤过程中应变能会减少, 故引入二次退化函数
$$ g(\varphi ) = {(1 - \varphi )^2} + k $$ (4) 式中,
$ k $ 是一个很小的参数, 可以防止$ \varphi = 1 $ 时刚度矩阵奇异. 本文取$ k = 1.0\;\times\;{10^{ - 7}} $ , 因此有$$ \begin{split} \int_\varGamma {{G_c}} {\rm{d}}\varGamma =& \int_\varOmega {{G_c}} \gamma (\varphi ,\nabla \varphi ){\rm{d}}V =\\ &\int_\varOmega {{G_c}} \left(\frac{{{\varphi ^2}}}{{2{l_0}}} + \frac{{{l_0}}}{2}{\left| {\nabla \varphi } \right|^2}\right){\rm{d}}\varOmega \end{split} $$ (5) 1.2 相场控制方程
假设弹性体受到面力
$ {\boldsymbol{\bar t}} $ 和体力$ {\boldsymbol{\bar b}} $ , 如图3所示, 则外部功$ {W_{{\rm{ext}}}} $ 可以表示为![]()
图 3 弥散裂纹与弹性体受力情况, 左侧为尖裂纹, 右侧为弥散裂纹Figure 3. Diffused crack and elastomeric stresses, with sharp cracks on the left and diffused crack on the right$$ {W_{{\rm{ext}}}} = \int_\varOmega {({\boldsymbol{\bar b}} \cdot {\boldsymbol{u}}){\rm{d}}\varOmega } + \int_{\partial \varOmega } {({\boldsymbol{\bar t}} \cdot {\boldsymbol{u}}){\rm{d}}S} $$ (6) 系统的总势能
$ \varPi = {\varPi _\varOmega } - {W_{{\rm{ext}}}} $ , 其对$ \varphi $ 和$ {\boldsymbol{u}} $ 求变分后可以得到$$ \begin{split} \delta \varPi =& \int_\varOmega {{\boldsymbol{\sigma }}\delta {\boldsymbol{u}} \cdot {\boldsymbol{n}}} {\rm{d}}S - \int_\varOmega {\nabla {\boldsymbol{\sigma }}\delta {\boldsymbol{u}}} {\rm{d}}V + \int_\varOmega {{\psi _0}\frac{{\partial g}}{{\partial \varphi }}} \delta \varphi {\rm{d}}V +\\ & {G_c}\int_\varOmega \left[{\frac{\varphi }{{{l_0}}}} - {l_0}{(\nabla \varphi )^{\rm{T}}}\nabla \varphi \right]\delta \varphi {\rm{d}}V + \\ &{G_c}{l_0}\int_{\partial \varOmega } {(\nabla \varphi } \cdot {\boldsymbol{n}})\delta \varphi {\rm{d}}S -\\ & \int_\varOmega {({\boldsymbol{\bar b}} \cdot \delta {\boldsymbol{u}}){\rm{d}}V - \int_{\partial \varOmega } {({\boldsymbol{\bar t}} \cdot \delta {\boldsymbol{u}}){\rm{d}}S} } \end{split} $$ (7) 式中,
$ {\boldsymbol{n}} $ 为法线向量,$ {\boldsymbol{\sigma }} $ 和$ {\boldsymbol{u}} $ 分别为柯西应力张量和位移张量. 根据变分原理, 式(7)应满足$ \delta \varPi = 0 $ , 因此得到相场与位移场耦合的形式为$$ \left. {\begin{array}{*{20}{l}} {\nabla {\boldsymbol{\sigma }} + {\boldsymbol{\bar b}} = {\boldsymbol{0}}} \\ {g'(\varphi ){\psi _0} + {G_c}\left(\dfrac{\varphi }{{{l_0}}} - {l_0}\Delta \varphi \right) = 0} \end{array}} \right\} $$ (8) 为了保证在相场值达到1时损伤不可逆, 引入历史场变量
$ \tilde H $ , 取值为损伤过程中弹性应变能密度最大值, 可以防止弹性体受压或卸载时裂纹愈合. 因此式(8)可以改写为$$ \left.{\begin{array}{*{20}{l}} {\nabla {\boldsymbol{\sigma }} + {\boldsymbol{\bar b}} = {\boldsymbol{0}}} \\ {g'(\varphi )\tilde H + {G_c}\left(\dfrac{\varphi }{{{l_0}}} - {l_0}\Delta \varphi \right) = 0} \end{array}} \right\} $$ (9) 1.3 弹性应变能的拉压分解
考虑到材料中只有拉伸载荷会引起裂纹扩展, 若不分离应变能中的拉应力和压应力部分, 会产生裂纹伪分岔现象[32]. 首先对应变张量进行谱分解
$$ {{\boldsymbol{\varepsilon }}_ \pm } = \sum\limits_{i = 1}^m {{{\left\langle {{{\boldsymbol{\varepsilon }}_{{i}}}} \right\rangle }_ \pm }} {{\boldsymbol{n}}_i} \otimes {{\boldsymbol{n}}_i} $$ (10) 式中,
$ {{\boldsymbol{\varepsilon }}_ + } $ 为拉应变张量,$ {{\boldsymbol{\varepsilon }}_ - } $ 为压应变张量,${{\boldsymbol{\varepsilon }}_{{i}}}$ 为主应变值,${{\boldsymbol{n}}_{{i}}}$ 为主应变值的方向.$ m $ 代表空间维数, 尖括号算子定义为:$ {\langle x\rangle }_{ + } = (x + \left|x\right|)/2 $ ,$ {\langle x\rangle }_{-} = (x-\left|x\right|)/2 $ . 由此, 弹性应变能密度分解为$$ \psi _e^ \pm ({\boldsymbol{\varepsilon }}) = \frac{\lambda }{2}\left\langle {{\rm{tr}}({\boldsymbol{\varepsilon }})} \right\rangle _ \pm ^2 + \mu {\rm{tr}}({\boldsymbol{\varepsilon }}_ \pm ^2) $$ (11) 式中,
$ \lambda $ 和$ \mu $ 为拉梅常数. 考虑到刚度退化仅作用于拉伸应变能, 则$ {\psi _e} $ 可表示为$$ \begin{split} {\psi _e}({\boldsymbol{\varepsilon }}) =& g(\varphi )\psi _e^ + ({\boldsymbol{\varepsilon }}) + \psi _e^ - ({\boldsymbol{\varepsilon }}) =\\ & [{(1 - \varphi )^2} + k]\psi _e^ + ({\boldsymbol{\varepsilon }}) + \psi _e^ - ({\boldsymbol{\varepsilon }}) \end{split} $$ (12) 1.4 相场断裂模型的分步算法
相场控制方程组往往是非线性的, 若通过Newton迭代法直接求解, 在ABAQUS隐式求解中容易造成不收敛. 因此, 可以对控制方程进行解耦, 使位移场方程和相场演化方程交错迭代求解. 本文用一层UEL计算相场, 一层UMAT计算位移场并完成可视化[33].
首先推导式(9)中方程的弱形式
$$ \left. \begin{array}{*{20}{l}} {\displaystyle\int_\varOmega {({\boldsymbol{\sigma }}\delta {\boldsymbol{\varepsilon }} - {\boldsymbol{\bar b}} \cdot \delta {\boldsymbol{u}}){\rm{d}}V - \displaystyle\int_{\partial {\varOmega _t}} {{\boldsymbol{\bar t}} \cdot \delta {\boldsymbol{u}}{\rm{d}}S = 0} } } \\ \displaystyle\int_\varOmega \left[g''(\varphi )\tilde H\delta \varphi + {G_c}\left({l_0}\nabla \varphi \cdot \nabla \delta \varphi + \dfrac{1}{{{l_0}}} \varphi \delta \varphi \right)\right]{\rm{d}}V = 0 \end{array} \right\} $$ (13) 对位移场和相场变量进行离散可得
$$ {\boldsymbol{u}} = {{\boldsymbol{N}}}^{u}{{\boldsymbol{u}}}^{e} = {\displaystyle \sum _{i = 1}^{n}{{\boldsymbol{N}}}_{i}^{u}{{\boldsymbol{u}}}_{i}} $$ (14) $$ \varphi = {{\boldsymbol{N}}^\varphi }{\varphi ^e} = \sum\limits_{i = 1}^n {N_i^\varphi } {\varphi _i} $$ (15) 式中,
$ n $ 是每个单元的结点数,$ {{\boldsymbol{N}}^u} $ 和$ {{\boldsymbol{N}}^\varphi } $ 分别为单元的位移场矩阵和相场形状函数,$ {\boldsymbol{N}}_{\boldsymbol{i}}^u $ 和$ N_i^\varphi $ 分别为结点$ i $ 的位移场形状函数矩阵和相场形状函数, 其表达式如下所示$$ {\boldsymbol{N}}_i^u = \left[ {\begin{array}{*{20}{c}} {{N_i}}&0 \\ 0&{{N_i}} \end{array}} \right] \text{, } N_i^\varphi = {N_i} $$ (16) 式中,
${{{N}}_i}$ 为结点$ i $ 的形状函数. 相应地可以得到$ {\boldsymbol{u}} $ 和$ \varphi $ 的梯度表达式为$$ \nabla {\boldsymbol{u}} = {{\boldsymbol{B}}^u}{{\boldsymbol{u}}^e} = \sum\limits_{i = 1}^n {{\boldsymbol{B}}_i^u{{\boldsymbol{u}}_i}} $$ (17) $$ \nabla \varphi = {{\boldsymbol{B}}^\varphi }{\varphi ^e} = \sum\limits_{i = 1}^n {{\boldsymbol{B}}_i^\varphi } {\varphi _i} $$ (18) 式中,
$ {{\boldsymbol{B}}^u} $ 和$ {{\boldsymbol{B}}^\varphi } $ 为单元的位移场和相场形状函数梯度,$ {\boldsymbol{B}}_i^u $ 和$ {\boldsymbol{B}}_i^\varphi $ 为结点$ i $ 的位移场和相场形状函数梯度, 其表达式为$$ {\boldsymbol{B}}_i^u = \left[ {\begin{array}{*{20}{c}} {{N_{i,x}}}&0 \\ 0&{{N_{i,y}}} \\ {{N_{i,y}}}&{{N_{i,x}}} \end{array}} \right] \text{, } {\boldsymbol{B}}_i^\varphi = \left[ {\begin{array}{*{20}{c}} {{N_{i,x}}} \\ {{N_{i,y}}} \end{array}} \right] $$ (19) 式中,
$ {N_{i,x}} $ 和$ {N_{i,y}} $ 为结点$ i $ 处的形状函数分别对$ x $ 和$ y $ 的导数. 将上述离散后的值代入弱形式平衡方程(13)中, 可以分别得到位移场和相场的右端残余向量为$$ \begin{split} {{\boldsymbol{r}}^u} = &\int_\varOmega {g(\varphi ){{({{\boldsymbol{B}}^u})}^{\rm{T}}}} {\boldsymbol{\sigma }}{\rm{d}}V -\\ & \int_\varOmega {{{({{\boldsymbol{N}}^u})}^{\rm{T}}}} {\boldsymbol{\bar b}}{\rm{d}}V - \int_{\partial {\varOmega _t}} {{{({{\boldsymbol{N}}^u})}^{\rm{T}}}} {\boldsymbol{\bar t}}{\rm{d}}S \end{split} $$ (20) $$\begin{split} {{\boldsymbol{r}}^\varphi } = &\int_\varOmega \Biggr\{ g'(\varphi ){({{\boldsymbol{N}}^\varphi })^{\rm{T}}}\tilde H +\\ & {G_c}\left[\frac{1}{{{l_0}}}{({{\boldsymbol{N}}^\varphi })^{\rm{T}}}\varphi + {l_0}{({{\boldsymbol{B}}^\varphi })^{\rm{T}}}\nabla \varphi \right] \Biggr\}{\rm{d}}V \end{split}$$ (21) 为了能够获得稳定的隐式求解公式, 采用Newton-Raphson迭代法, 对位移场和相场进行解耦求解, 表达式如下
$$ \left[ {\begin{array}{*{20}{c}} {{\boldsymbol{K}}_n^u}&0 \\ 0&{{\boldsymbol{K}}_n^\varphi } \end{array}} \right] \left[ {\begin{array}{*{20}{c}} {{u_{n + 1}}} \\ {{\varphi _{n + 1}}} \end{array}} \right] = - \left[ {\begin{array}{*{20}{c}} {{\boldsymbol{r}}_n^u} \\ {{\boldsymbol{r}}_n^\varphi } \end{array}} \right] $$ (22) 式中,
$ {\boldsymbol{K}}_n^u $ 和$ {\boldsymbol{K}}_n^\varphi $ 分别为$ {t_n} $ 时刻的位移场和相场刚度矩阵, 具体表达式如下$$ {\boldsymbol{K}}_n^u = \int_\varOmega {g(\varphi ){{({{\boldsymbol{B}}^u})}^{\rm{T}}}} {{\boldsymbol{C}}_0}{{\boldsymbol{B}}^u}{\rm{d}}V $$ (23) $$\qquad \begin{split} {\boldsymbol{K}}_n^\varphi =& \int_\varOmega {\Big[{{({{\boldsymbol{B}}^\varphi })}^{\rm{T}}}{G_c}{l_0}} ({{\boldsymbol{B}}^\varphi }) +\\ &{({{\boldsymbol{N}}^\varphi })^{\rm{T}}}\left(2\tilde H + \frac{{{G_c}}}{{{l_0}}}\right){{\boldsymbol{N}}^\varphi }\Big]{\rm{d}}\varOmega \end{split}$$ (24) 式中,
$ {{\boldsymbol{C}}_0} $ 为弹性体的初始刚度矩阵.2. 多轴比例加载的周期性边界条件
在尺度转换的分析过程中, 常用体积平均应变和体积平均应力来描述RVE模型的各项宏观特性. 在本文中, 以作用在RVE模型上的名义应力作为宏观量. 首先, RVE模型的宏观和细观第一Piola-Kirchhoff应力张量表示为
${{\boldsymbol{P}}_{{M}}}$ 和${{\boldsymbol{P}}_{{m}}}$ , 其宏观和细观变形张量表示为${{\boldsymbol{F}}_{{M}}}$ 和${{\boldsymbol{F}}_{{m}}}$ $$ \left. {\begin{array}{*{20}{l}} {{{\boldsymbol{P}}_{{M}}} = \dfrac{1}{{{V_0}}}\displaystyle\int_\varOmega {{{\boldsymbol{P}}_{{m}}}{\rm{d}}V} } \\ {{{\boldsymbol{F}}_{{M}}} = \dfrac{1}{{{V_0}}}\displaystyle\int_\varOmega {{{\boldsymbol{F}}_{{m}}}{\rm{d}}V} } \end{array}} \right\}$$ (25) 式中,
$ {V_0} $ 和$ \varOmega $ 为RVE模型的体积和域. 通过式(25), 可以实现尺度转换下对应参量的转换, 使RVE模型的宏观名义应力可用于表征宏观模型的承载能力.2.1 周期性边界条件
如图4所示, 在二维笛卡尔坐标系下,
$ {x_1} $ 和$ {x_2} $ 分别表示水平方向和竖直方向, RVE模型的尺寸为$ 2{l_1} \times 2{l_2} $ . 考虑RVE模型处于多轴加载状态[34],$ {P_{11}} $ ,$ {P_{22}} $ 和$ {P_{12}} $ 分别为作用在RVE模型上的水平、竖直方向的宏观名义应力和面内剪切名义应力.在施加双轴以及面内剪切载荷的工况下, RVE模型的宏观变形张量可以写为
$$ \begin{split} {{\boldsymbol{F}}_M} = &{\lambda _{11}}{{\boldsymbol{e}}_{{1}}} \otimes {{\boldsymbol{E}}_{\boldsymbol{1}}} + {\lambda _{22}}{{\boldsymbol{e}}_{{2}}} \otimes {{\boldsymbol{E}}_{\boldsymbol{2}}} + {\lambda _{12}}{{\boldsymbol{e}}_{{2}}} \otimes {{\boldsymbol{E}}_{{1}}} =\\ & \left[\begin{array}{cc} \lambda_{11} & \lambda_{12} \\ 0 & \lambda_{22} \end{array}\right] \end{split} $$ (26) 式中,
$ {\lambda _{11}} $ 和$ {\lambda _{22}} $ 分别是RVE模型水平、竖直方向的宏观变形量,$ {\lambda _{12}} $ 是面内剪切的宏观变形量.${{\boldsymbol{e}}_{{i}}}(i = 1,2,3)$ 是当前变形下的笛卡尔坐标系基向量,${{\boldsymbol{E}}_{{i}}}(i = 1,2,3)$ 是未变形时的笛卡尔坐标系基向量.RVE模型的广义周期性边界条件表示为
$$ {u_B} - {u_A} = ({{\boldsymbol{F}}_{{M}}} - {\boldsymbol{I}})({X_B} - {X_A}) = {{\boldsymbol{H}}_{{M}}}({X_B} - {X_A}) $$ (27) 式中,
$ {u_A} $ 和$ {u_B} $ 是RVE模型中任意两对称边界上节点对的当前位移, 即图4中用红色部分标出边界A和边界B, 或用蓝色标出的边界A和边界B; XA和XB是RVE模型中对称边界上节点对的初始位置;${{\boldsymbol{H}}_{{M}}}$ 是宏观应变张量,$ {\boldsymbol{I}} $ 则是单位张量.根据Hill-Mandel的宏、细观能量守恒条件, 整个RVE模型的宏观应变能密度应等于该模型中所有单元的平均应变能密度, 其数学表达式为
$$ \frac{1}{{{V_0}}}\int_\varOmega {{{\boldsymbol{P}}_{{m}}}} :\delta {{\boldsymbol{F}}_{{m}}}{\rm{d}}V = {{\boldsymbol{P}}_{{M}}}:\delta {{\boldsymbol{F}}_{{M}}} $$ (28) 静力加载下, 忽略体力和惯性的作用, RVE模型在多轴加载下的宏观总功变化率可以表示为
$$\qquad \begin{split} \dot W = & A{\boldsymbol{P}}:{{{\dot{\boldsymbol F}}}_M} = \\ & A({P_{11}}{{\dot \lambda }_{11}} + {P_{12}}{{\dot \lambda }_{12}} + {P_{22}}{{\dot \lambda }_{22}}) \end{split} $$ (29) $$ {\boldsymbol{P}} = \left[ {\begin{array}{*{20}{c}} {{P_{11}}} \\ {{P_{12}}} \\ {{P_{22}}} \end{array}} \right] \text{, } {{\dot{\boldsymbol F}}_M} = \left[ {\begin{array}{*{20}{c}} {{{\dot \lambda }_{11}}} \\ {{{\dot \lambda }_{12}}} \\ {{{\dot \lambda }_{22}}} \end{array}} \right] $$ (30) 式中,
$ A $ 是二维RVE模型的面积,$ {\boldsymbol{P}} $ 为RVE模型的宏观应力张量,$ {{\dot{\boldsymbol F}}_M} $ 是宏观变形变化率梯度张量,$ \dot \lambda $ 为宏观变形率.2.2 比例加载周期性边界条件的数值实现
引入虚拟节点, 在虚拟节点上施加相应的位移载荷, 可以将周期性边界条件施加到RVE模型上. 基于均匀化思想, 通过能量守恒将RVE模型的宏观应力、应变与虚拟节点的自由度、反力关联起来.
如式(26)所示, 在多轴加载条件下, RVE模型的宏观变形梯度张量
$ {{\boldsymbol{F}}_M} $ 存在3个非零分量. 在二维平面内, 每个节点可以提供两个自由度, 因此本文需要引入两个虚拟节点ghostnode1和ghostnode2.首先, 考虑RVE模型在水平和竖直方向的变形能与作用在虚拟节点ghostnode1的外力做功应满足能量相等的条件, 用变化率的形式可以表示为
$$ A{P_{11}}{\dot \lambda _{11}} + A{P_{22}}{\dot \lambda _{22}} = {p_1}{\dot q_1} + {p_2}{\dot q_2} $$ (31) 式中,
$ {q_1} $ ,$ {q_2} $ 和$ {p_1} $ ,$ {p_2} $ 分别是虚拟节点ghostnode1在水平(x1)方向和竖直(x2)方向的位移与反力. 上面的公式意味着变换前后的内积不变, 因此两者应满足正交变换$$ \left[ {\begin{array}{*{20}{c}} {{p_1}} \\ {{p_2}} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {\cos {\alpha _1}}&{ - \sin {\alpha _1}} \\ {\sin {\alpha _1}}&{\cos {\alpha _1}} \end{array}} \right] \left[ {\begin{array}{*{20}{c}} {A{P_{11}}} \\ {A{P_{22}}} \end{array}} \right] $$ (32) $$ \left[ {\begin{array}{*{20}{c}} {{{\dot q}_1}} \\ {{{\dot q}_2}} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {\cos {\alpha _1}}&{ - \sin {\alpha _1}} \\ {\sin {\alpha _1}}&{\cos {\alpha _1}} \end{array}} \right] \left[ {\begin{array}{*{20}{c}} {{{\dot \lambda }_{11}}} \\ {{{\dot \lambda }_{22}}} \end{array}} \right] $$ (33) 在式(32)中, 令
$ {p_1} = 0 $ , 可以得到$$ \left. {\begin{array}{*{20}{l}} {\cos {\alpha _1} = \dfrac{{{P_{22}}}}{{\sqrt {P_{11}^2 + P_{22}^2} }}} \\ {\sin {\alpha _1} = \dfrac{{{P_{11}}}}{{\sqrt {P_{11}^2 + P_{22}^2} }}} \end{array}} \right\} $$ (34) 从而可以得到虚拟节点ghostnode1上的非零反力
$ {p_2} $ 的表达式为$$ {p_2} = A\sqrt {P_{11}^2 + P_{22}^2} $$ (35) 然后引入剪切应变能, 考虑RVE模型的总变形能与作用在虚拟节点ghostnode2的外力做功应满足守恒的条件, 同样用变化率的形式可以表示为
$$ \dot W = A{P_{12}}{\dot \lambda _{12}} + {p_2}{\dot q_2} = {p_3}{\dot q_3} + {p_4}{\dot q_4} $$ (36) 式中,
$ {q_3} $ ,$ {q_4} $ 和$ {p_3} $ ,$ {p_4} $ 分别是虚拟节点ghostnode2在水平方向和竖直方向的位移与反力. 同理, 通过正交变换可以得到下式$$ \left[ {\begin{array}{*{20}{c}} {{p_3}} \\ {{p_4}} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {\cos {\alpha _2}}&{ - \sin {\alpha _2}} \\ {\sin {\alpha _2}}&{\cos {\alpha _2}} \end{array}} \right] \left[ {\begin{array}{*{20}{c}} {A{P_{12}}} \\ {{p_2}} \end{array}} \right] $$ (37) $$ \left[ {\begin{array}{*{20}{c}} {{{\dot q}_3}} \\ {{{\dot q}_4}} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {\cos {\alpha _2}}&{ - \sin {\alpha _2}} \\ {\sin {\alpha _2}}&{\cos {\alpha _2}} \end{array}} \right] \left[ {\begin{array}{*{20}{c}} {{{\dot \lambda }_{12}}} \\ {{{\dot q}_2}} \end{array}} \right] $$ (38) 在式(37)中, 令
$ {p_3} = 0 $ , 则有$$ \left. {\begin{array}{*{20}{l}} {\cos {\alpha _2} = \dfrac{{\sqrt {P_{11}^2 + P_{22}^2} }}{{\sqrt {P_{11}^2 + P_{22}^2 + P_{12}^2} }}} \\ {\sin {\alpha _2} = \dfrac{{{P_{12}}}}{{\sqrt {P_{11}^2 + P_{22}^2 + P_{12}^2} }}} \end{array}} \right\} $$ (39) 此时, 作用在虚拟节点ghostnode2的非零反力
$ {p_4} $ 为$$ {p_4} = A\sqrt {P_{11}^2 + P_{22}^2 + P_{12}^2} $$ (40) 通过式(33)和式(38)的逆变换, 可以得到
$$ \left[ {\begin{array}{*{20}{c}} {{{\dot \lambda }_{11}}} \\ {{{\dot \lambda }_{22}}} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {\cos {\alpha _2}}&{\sin {\alpha _2}} \\ { - \sin {\alpha _2}}&{\cos {\alpha _2}} \end{array}} \right] \left[ {\begin{array}{*{20}{c}} {{{\dot q}_1}} \\ {{{\dot q}_2}} \end{array}} \right] $$ (41) $$ \left[ {\begin{array}{*{20}{c}} {{{\dot \lambda }_{12}}} \\ {{{\dot q}_2}} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {\cos {\alpha _2}}&{\sin {\alpha _2}} \\ { - \sin {\alpha _2}}&{\cos {\alpha _2}} \end{array}} \right] \left[ {\begin{array}{*{20}{c}} {{{\dot q}_3}} \\ {{{\dot q}_4}} \end{array}} \right] $$ (42) 综合式(41)和式(42), 对虚拟节点自由度重新编号, 可以得到关联RVE模型宏观变形率与虚拟节点位移变化率的转换矩阵
$ {\boldsymbol{T}} $ $$ \begin{split} \left[ {\begin{array}{*{20}{c}} {{{\dot \lambda }_{11}}} \\ {{{\dot \lambda }_{12}}} \\ {{{\dot \lambda }_{22}}} \end{array}} \right] =& \left[ {\begin{array}{*{20}{c}} {\cos {\alpha _1}}&{ - \sin {\alpha _1}\sin {\alpha _2}}&{\sin {\alpha _1}\cos {\alpha _2}} \\ 0&{\cos {\alpha _2}}&{\sin {\alpha _2}} \\ { - \sin {\alpha _1}}&{ - \cos {\alpha _1}\sin {\alpha _2}}&{\cos {\alpha _1}\cos {\alpha _2}} \end{array}} \right] \left[ {\begin{array}{*{20}{c}} {{{\dot q}_1}} \\ {{{\dot q}_2}} \\ {{{\dot q}_3}} \end{array}} \right] =\\ & {\boldsymbol{T}} \left[ {\begin{array}{*{20}{c}} {{{\dot q}_1}} \\ {{{\dot q}_2}} \\ {{{\dot q}_3}} \end{array}} \right] \\[-15pt] \end{split} $$ (43) 同理, 可得RVE模型宏观名义应力与虚拟节点反力的转换方程
$$ A \left[ {\begin{array}{*{20}{c}} {{P_{11}}} \\ {{P_{12}}} \\ {{P_{22}}} \end{array}} \right] = {\boldsymbol{T}} \left[ {\begin{array}{*{20}{c}} {{p_1}} \\ {{p_2}} \\ {{p_3}} \end{array}} \right] $$ (44) 式中, q1, q2和p1, p2分别是虚拟节点ghostnode1在水平方向和竖直方向的位移和反力, q3和p3分别是虚拟节点ghostnode2在水平方向的位移和反力.
以RVE模型在竖直方向的宏观名义应力
$ {P_{22}} $ 作为参考应力, 则可以用两个应力比参数来描述作用在RVE模型上3个宏观应力之间的关系$$ \left.\begin{array}{c}\begin{array}{cc}{\rho }_{1} = {P}_{11}/{P}_{22}& (双轴应力比)\end{array}\\ \begin{array}{cc}{\rho }_{2} = {P}_{12}/{P}_{22}& (剪切应力比)\end{array}\end{array} \right\}$$ (45) 将式(45)代入式(34)和式(39), 与
$ {\alpha _1} $ 和$ {\alpha _2} $ 相关的三角函数可以表示为$$ \left.\begin{array}{l} \mathrm{cos}{\alpha }_{1} = \dfrac{1}{\sqrt{1 + {\rho }_{1}^{2}}}\\ \mathrm{sin}{\alpha }_{1} = \dfrac{{\rho }_{1}}{\sqrt{1 + {\rho }_{1}^{2}}} \\ \mathrm{cos}{\alpha }_{2} = \dfrac{\sqrt{1 + {\rho }_{1}^{2}}}{\sqrt{1 + {\rho }_{1}^{2} + {\rho }_{2}^{2}}}\\ \mathrm{sin}{\alpha }_{2} = \dfrac{{\rho }_{2}}{\sqrt{1 + {\rho }_{1}^{2} + {\rho }_{2}^{2}}} \end{array}\right\} $$ (46) 进而, 式(43)和式(44)中的转换矩阵
$ {\boldsymbol{T}} $ 也可以由宏观应力比参数$ {\rho _1} $ 和$ {\rho _2} $ 来表示. 在整个加载历程中, 应力比保持恒定, 则转换矩阵$ {\boldsymbol{T}} $ 不变, 即可实现对二维RVE模型的比例应力加载.由式(43)可以得到RVE模型在右边界(
$ {x_1} = {l_1} $ )的节点位移方程为$$ \left. {\begin{array}{*{20}{l}} {u_{1,{x_1} = {l_1}}} = {u_{1,{x_1} = - {l_1}}} + 2{l_1}(\cos {\alpha _1} \cdot {q_1} -\\ \qquad\sin {\alpha _1}\sin {\alpha _2} \cdot {q_2} + \sin {\alpha _1}\cos {\alpha _2} \cdot {q_3}) \\ \;\;\; {{u_{2,{x_1} = {l_1}}} = {u_{2,{x_1} = - {l_1}}}} \end{array}} \right\} $$ (47) 同理, RVE模型在上边界(
$ {x_2} = {l_2} $ )的节点位移方程可以表达为$$ \left. {\begin{array}{*{20}{l}} {{u_{1,{x_2} = {l_2}}} = {u_{1,{x_2} = - {l_2}}} + 2{l_2}(\cos {\alpha _2} \cdot {q_2} + \sin {\alpha _2} \cdot {q_3})} \\ {u_{2,{x_2} = {l_2}}} = {u_{2,{x_2} = - {l_2}}} + 2{l_2}( - \sin {\alpha _1} \cdot {q_1} - \cos {\alpha _1}\sin {\alpha _2} \cdot {q_2} +\\ \qquad \cos {\alpha _1}\cos {\alpha _2} \cdot {q_3} \end{array}} \right\} $$ (48) 本文的数值仿真工作基于ABAQUS通用有限元分析软件平台(版本: 6.14)开展, 比例加载的周期性边界条件通过编写ABAQUS的用户子程序(user subroutine)实现. 在数值模拟中, 只在虚拟节点ghostnode2的水平(
$ {x_1} $ )方向施加位移载荷$ {q_3} $ , 则该方向的节点反力$ {p_3} $ 为非零值. 当多轴加载的应力比参数已确定, 通过式(49)即可求解得到RVE模型在竖直($ {x_2} $ )方向的宏观名义应力$$ {P_{22}} = \frac{{{p_3}}}{{A\sqrt {1 + \rho _1^2 + \rho _2^2} }} $$ (49) 然后, 通过宏观应力比关系就可得到加载过程中作用在RVE模型上的其他宏观名义应力.
3. 多轴比例加载下的周期性多孔结果断裂行为研究
3.1 周期性多孔结构的RVE模型
针对周期性多孔结构, 本文建立其二维RVE模型, 其中含孔单胞的数量为
$ 2 \times 2 $ , 尺寸为$5\;\rm mm \times 5\;\rm mm$ , 孔洞半径R为$ 1.0\;\rm mm $ , RVE模型的几何构型如图5所示. 其有限元模型采用四边形平面应变单元(CPE4)划分网格, 单元数为115200, 节点数为116637, 模型中对称边界上的节点数均为241个.![]()
图 5 周期性多孔结构RVE模型的几何构型(单位: mm)Figure 5. Geometric configuration of the RVE model for periodic porous structures (unit: mm)RVE模型的裂纹弥散宽度
$ {l_0} = 2 h $ , 其中$ h $ 为模型网格的最小尺寸. 对本文建立的模型, h = 13 μm,$ {l_0} = 26 \;{\text{μm}} $ .本文计算所用的材料拉梅弹性参数分别为
$ \lambda = 121.15\;\rm kN/m{m^2} $ ,$\; \mu = 80.77\;\rm kN/m{m^2} $ , 临界能量释放率$ {G_c} = 2.7 \;\rm N/mm $ [35]. 计算过程中加载步长固定为$ \Delta u = 5 \;{\text{μm}} $ .3.2 与
$ 8 \times 8 $ 常规边界条件模型的计算结果比较为了验证采用比例加载周期性边界条件后所表征的断裂行为的正确性, 本文建立了一个施加常规边界条件的
$ 8 \times 8 $ 周期性多孔结构模型, 模型下端约束竖直方向的位移, 上端施加位移载荷. 并与$ 2 \times 2 $ 的RVE模型在单轴拉伸条件下断裂行为的数值模拟结果相比较.$ 8 \times 8 $ 周期性多孔结构的边界条件与几何构型如图6所示, 模型下端约束竖直方向的位移, 上端施加位移载荷; 其孔洞半径R为$ 1.0\; \rm mm $ , 有限元模型的单元数为204800, 节点数为210498,$ {l_0} = 80\;{\text{μm}} $ , 采用的材料参数与$ 2 \times 2 $ RVE模型一致.![]()
图 6$ 8 \times 8 $ 多孔模型的单轴拉伸边界条件Figure 6. Boundary conditions for uniaxial stretch of$ 8 \times 8 $ periodic porous model对于2×2 RVE模型, 施加
$ {x_2} $ 方向单轴拉伸载荷时, 应力比$ {\rho _1} = {\rho _2} = 0 $ . 在加载过程中, 作用在虚拟节点ghostnode2上的非零反力p3的峰值为${p_{3 \text{-} \max }} = 86\;791\; \rm N$ , 通过式(49)可以计算得到RVE模型在$ {x_2} $ 方向宏观名义应力的极值P22-max为867.91 MPa, 这即是周期性多孔结构在单轴拉伸条件下的极限强度.为了验证
$ 2 \times 2 $ RVE模型网格独立性, 本文在原有115200单元数模型的基础上, 分别建立了单元数为80000和156800的稀疏网格和加密网格模型, 后两者的最小单元尺寸分别为$ h = 15.7\;{\text{μm}} $ 和$ h = 11.2\;{\text{μm}} $ , 裂纹特征宽度参数依然取为$ {l_0} = 2h $ . 3组不同单元数的模型在单轴拉伸载荷下的名义应力-应变曲线如图7所示. 可以发现, 3组模型的曲线吻合度较高, 极限载荷值相差极小. 因此综合考虑模型计算精度和计算成本, 本文选择单元数为115200的$ 2 \times 2 $ RVE模型进行后续研究与讨论.![]()
图 7 单轴拉伸载荷下单元数分别为80000, 115200和156800的周期性多孔结构的2×2 RVE模型名义应力-应变曲线的比较Figure 7. Comparison of nominal stress-strain curves of the RVE model with 80000, 115200 and 156800 elements under tensile loading condition此外, Miehe等[32]在应用断裂相场方法时, 也对模型的网格独立性开展了相似讨论, 同样发现模型网格的疏密程度对于得到的载荷-位移曲线影响程度较小.
如图8所示, 在单轴拉伸载荷作用下, 基于
$ 2 \times 2 $ RVE模型与基于$ 8 \times 8 $ 常规模型得到的名义应力-应变曲线吻合较好. 可以发现, 两条曲线在断裂失效时刻的承载极限${P_{22 \text{-} \max }}$ 和应变值$ {\varepsilon _{22}} $ 都几乎一致.![]()
图 8 单轴拉伸载荷下周期性多孔结构的$ 2 \times 2 $ RVE模型与$ 8 \times 8 $ 模型的名义应力-应变曲线的比较Figure 8. Comparison of nominal stress-strain curves of the$ 2 \times 2 $ RVE model and the$ 8 \times 8 $ periodic porous model under uniaxial tensile loading condition图9分别为周期性多孔结构的
$ 2 \times 2 $ RVE模型与$ 8 \times 8 $ 常规模型在施加竖直($ {x_2} $ )方向单轴拉伸载荷后的渐进断裂相场云图. 在相同的加载条件下, 比较图9(a)和图9(d)、图9(b)和图9(e)、图9(c)和图9(f)的断裂相场云图发现, 周期性多孔结构的$ 2 \times 2 $ RVE模型与$ 8 \times 8 $ 常规模型均是在孔边水平位置萌生裂纹, 并且裂纹均沿水平方向扩展($ \varphi = 1 $ 的路径). 在相同的应变值下, 两种模型的裂纹扩展模式基本一致. 图9(f)中模型的破坏程度不完全一致是由于, 在本文所给的常规单轴拉伸载荷作用下, 受边界条件的影响,$ 8 \times 8 $ 模型内部孔洞的应力分布与模型边缘孔洞的应力分布有一定的差别, 模型内部孔洞的应力状态更接近RVE模型的应力状态.$ 8 \times 8 $ 模型边缘孔洞的应力集中稍大, 因此孔边裂纹的扩展会稍快于中间位置的孔洞. 从宏观上来看, 可以认为周期性多孔结构的$ 8 \times 8 $ 常规模型与$ 2 \times 2 $ RVE模型的裂纹扩展模式基本一致, 且本文图8中两组模型的名义应力-应变曲线吻合较好. 因此, 采用RVE模型研究周期性多孔结构在复杂加载下的力学响应及断裂行为是有效可行的.![]()
图 9 单轴拉伸载荷下的周期性多孔结构的$ 2 \times 2 $ RVE模型与$ 8 \times 8 $ 模型的断裂相场云图Figure 9. Phase-field fracture contours of the$ 2 \times 2 $ RVE model and$ 8 \times 8 $ model for periodic porous structure under uniaxial tensile loading综上所述, 本文采用基于多轴比例加载的周期性边界条件的RVE模型, 可以正确有效模拟周期性结构在复杂载荷作用下的力学响应, 从而评估结构的宏观力学性能. 通过典型算例验证, 表明该多尺度分析方法与相场断裂方法相结合用以研究周期性多孔结构在多轴加载下的损伤断裂行为具有可行性和可靠性.
3.3 双轴比例加载下周期性多孔结构的断裂行为
首先考虑周期性多孔结构在双轴载荷作用下的断裂行为, 即令剪切应力比
$ {\rho _2} = 0 $ . 以RVE模型在竖直方向的拉伸载荷为主应力, 研究水平方向载荷的引入对多孔结构的裂纹萌生、扩展以及承载能力的影响. 这里选取双轴应力比$ {\rho _1} = - 1, $ $ - 0.5, $ $ 0 $ ,$ 0.5 $ ,$ 1 $ ; 当$ \;{\rho _1} < 0 $ 时, 表明在水平方向施加的是压缩载荷.表1给出了当
$ {\rho _2} = 0 $ 时, 周期性多孔结构的$ 2 \times 2 $ RVE模型在不同双轴应力比$ {\rho _1} $ 下的竖直方向的拉伸承载极限${P_{22 \text{-} \max }}$ 和断裂相场云图. 可以发现, RVE模型在承受一个方向的拉伸载荷时, 另一个方向的拉伸载荷会提高结构的承载能力, 而另一个方向的压缩载荷则会促进裂纹的扩展, 降低结构的承载能力.表 1 双轴比例加载下的周期性多孔结构的承载极限${P_{22 \text{-} \max }}$ 与断裂相场云图Table 1. Extreme load${P_{22 \text{-} \max }}$ and phase-field fracture contours for the RVE model under biaxial proportional loading condition$ {\rho _1} $ −1 −0.5 0 0.5 1 P22-max/MPa 702.26 782.25 867.91 947.08 971.39 Phase-
field fracture contours
值得注意的是, 当
$ \;{\rho _1} \lt 1 $ 时, 孔边裂纹扩展路径保持水平. 当$ \;{\rho _1} = 1.0 $ 时, 模型表现出十字正交型裂纹扩展模式. 这是由于多孔结构的RVE模型具有中心对称性, 在双轴等值拉伸载荷作用下, 模型在两个方向产生等量的变形, 水平和竖直方向的孔边应力、应变状态相同. 因此, 当双轴应力比$ \;{\rho _1} = 1.0 $ 时, 孔边裂纹在水平和竖直方向同时萌生并同步扩展.3.4 多轴比例加载下周期性多孔结构的断裂行为
在3.3节的基础上, 引入面内剪切应力, 在周期性边界条件模型中改变剪切应力比
$ {\rho _2} $ 的取值, 讨论周期性多孔结构在同时承受双轴和剪切的复合加载下的断裂行为. 由于在面内剪切方向上, 顺时针和逆时针的剪切载荷具有对称性, 为了减少计算成本, 本文的$ {P_{12}} $ 只选取为顺时针方向, 并且仍以竖直方向的拉伸载荷$ {P_{22}} $ 为主应力. 这里,$\; {\rho _2} $ 的考察范围为$ \left[0, 1\right] $ , 取值分别为$ \;{\rho _2} = 0.25 $ ,$ 0.5 $ ,$ 0.75 $ ,$ 1 .0$ .表2给出了周期性多孔结构的2×2 RVE模型在不同多轴比例加载工况下的断裂相场云图. 观察表2后可以发现, 在不同的双轴应力比
$ {\rho _1} $ 和剪切应力比$ {\rho _2} $ 组合范围内, RVE模型的裂纹萌生和扩展模式存在一定的规律性. 为此, 归纳5种不同的裂纹扩展模式的表征符号和渐进破坏的相场演化云图, 如表3所示.表 2 多轴比例加载下的周期性多孔结构的断裂相场云图Table 2. Phase-field fracture contours of periodic porous structure model under multiaxial proportional loading$ {\rho _1} $, $ {\rho _2} $ −1 −0.5 0 0.5 1.0 0 
0.25 
0.5 
0.75 
1.0 
表 3 5种不同裂纹扩展模式的表征符号和渐进破坏的断裂相场云图Table 3. Representative symbols and progressive phase-field fracture contours for the five distinguished crack propagation modesRepresentative symbol 
crack propagation modes horizontal
typeS-type double arc-type 45°-type orthogonal cross-type progressive phase field
fracture contours
从多孔结构的断裂失效模式来看, 当
$ {\rho _1} < 0 $ 且$ {\rho _2} $ 较小时, 裂纹从孔两边对称萌生, 初始扩展路径与水平轴成一定的夹角, 在扩展过程中与相邻孔边裂纹相互靠近并最终汇合成为S型单裂纹; 当$ {\rho _1} > 0 $ 且$ {\rho _2} $ 较大时, 孔边裂纹萌生点、初始扩展方向与水平轴的夹角进一步增大并趋向于$ {45^ \circ } $ , 裂纹扩展方向随载荷的增加逐渐转向水平, 最终扩展至水平相邻孔边形成双弧线型裂纹.Sarac等[36]对错位排列的金属玻璃异质周期性多孔结构施加单轴拉伸载荷, 并认为孔洞附近同时存在拉应力和剪应力, 实验结果如图10所示, 水平相邻的孔洞之间呈现双弧线型的裂纹, 并逐渐扩展到相邻的孔洞位置. 虽材料不同, 该实验结果与本文多孔结构在相似载荷工况下的模拟结果较一致, 可以认为采用本文所提出方法得到的模拟结果是有效可信的.
![]()
对于S型和双弧线型裂纹模式来说, 随着
$ {\rho _1} $ 和$ {\rho _2} $ 的增大, 裂纹萌生、扩展与水平方向的夹角均逐渐增大. 当双轴应力比$\; {\rho _1} = 1.0 $ 时, 施加面内剪切载荷后, 多孔结构呈现出$ {45^ \circ } $ 斜裂纹扩展模式, 而且剪切应力的增大对$ {45^ \circ } $ 斜裂纹的扩展路径并无显著影响.图11给出了在不同多轴加载比例下, 周期性多孔结构在竖直方向的拉伸承载极限
${P_{22 \text{-} \max }}$ 随应力比ρ1和ρ2变化的曲线. 曲线上的不同标记点符号对应于表3中相应的裂纹扩展模式, 并以黄色“$ \times $ ”符号表示裂纹扩展模式的转变点.![]()
图 11 周期性多孔结构在竖直方向的承载极限${P_{22 \text{-} \max }}$ 随应力比$ {\rho _1} $ 和$ {\rho _2} $ 变化的曲线: (a)保持应力比$ {\rho _2} $ 不变, (b)保持应力比$ {\rho _1} $ 不变, (c)局部放大图Figure 11. Changing curves of 2-directional extreme load${P_{22 \text{-} \max }}$ with$ {\rho _1} $ and$ {\rho _2} $ for periodic porous structures: (a) keep$ {\rho _2} $ unchanged, (b) keep$ {\rho _1} $ unchanged and (c) partial enlarged graph从图11中可以发现, 当双轴应力比
$ {\rho _1} $ 固定时, 随着剪切应力比$ {\rho _2} $ 的增加, 多孔结构在竖直方向的拉伸承载极限显著降低; 而且,$ {\rho _1} $ 越大, 承载极限下降的幅度越明显, 说明本文所研究的周期性多孔结构的抗剪切性能较弱.当
$ 0 \lt {\rho _2} \leqslant 0.5 $ 时, 随着$ {\rho _1} $ 由负向正增加, 即横向载荷由压缩向拉伸状态转变, 多孔结构在竖直方向的拉伸承载能力呈先增加后减少的趋势; 而当$ 0.5 \lt {\rho _2} \leqslant 1.0$ 时, 结构的承载极限随$ {\rho _1} $ 的增加而降低; 说明: 随着剪切载荷的增加, 其对多孔结构极限强度的影响逐渐起到主导作用.此外, 当
$ \;{\rho _1} = 1.0 $ ,$ \; {\rho _2} = 0 $ , 即多孔结构承受双轴等量拉伸载荷时, 其极限承载强度达到最大值.图12总结了周期性多孔结构在复杂多轴比例加载范围内的断裂相图. 根据表3中所列出的5种裂纹扩展模式, 断裂相图总共被划分为5个区域, 并用黄色虚线表示断裂模式转变的相边界, 图12(b)、图12(c)和图12(d)分别为相图局部区域的放大图.
![]()
图 12 周期性多孔结构在不同多轴加载条件下的裂纹扩展模式相图Figure 12. Phase diagram of crack propagation mode for periodic porous structure under different multiaxial proportional loading在双轴和剪切载荷同时作用下, 周期性多孔结构的裂纹扩展模式主要呈现为S型与双弧线型, 且这两种扩展模式的分界线(即相边界)呈近似线性. 当剪切应力比
$ {\rho _2} \leqslant 0.03 $ 时, 裂纹扩展模式由S型向水平裂纹过渡; 当双轴应力比$ {\rho _1} \geqslant 0.94 $ 时, 裂纹扩展模式由双弧线型向$ {45^ \circ } $ 斜裂纹过渡; 在$ {\rho _1} \geqslant 0.95 $ 且$ {\rho _2} \leqslant 0.04 $ 的很小区域内, 裂纹扩展模式呈现十字正交型.观察图12可以发现, 当剪切应力比
$ {\rho _2} $ 趋近于0时, 面内剪切应力$ {P_{12}} $ 相对于竖直方向的拉伸应力$ {P_{22}} $ 已经很小, 可以认为RVE模型退化为承受竖直和水平方向的双轴比例加载, 多孔结构表现出水平裂纹扩展的失效模式; 同理, 当双轴应力比$ {\rho _1} $ 趋近于1时, 模型承受水平和竖直方向的等量加载并耦合面内剪切加载, 多孔结构展示出$ {45^ \circ } $ 斜裂纹扩展的失效模式; 当作用在模型上的多轴应力趋近于等值状态时, 结构表现出十字正交裂纹扩展的失效模式.4. 结论
本文建立了周期性多孔结构的RVE模型, 采用相场断裂方法结合周期性边界条件, 通过引入两个多轴应力比控制参数, 实现模型在复杂应力状态下的恒定比例加载, 并探究其在多轴比例加载工况下的裂纹萌生、扩展行为和承载极限等力学问题. 数值仿真结果得到如下结论.
(1)采用细观力学方法, 建立RVE模型, 施加比例加载的周期性边界条件, 可以正确表征周期性结构的宏观力学性能. 在此基础上, 与相场断裂方法相结合, 可以有效模拟周期性多孔结构在复杂加载条件下的断裂行为.
(2)周期性多孔结构在双轴加载状态下, 剪切应力比
$ {\rho _2} = 0 $ . 结构在承受单轴拉伸载荷基础上, 另一个方向的拉伸载荷会提高结构的拉伸承载能力, 压缩载荷则会降低结构的承载能力, 且双轴加载下裂纹扩展路径保持水平. 当双轴应力比$ \;{\rho _1} = 1.0 $ 时, 多孔结构表征出十字正交裂纹的断裂模式, 这是由结构的几何、载荷和边界条件的对称性导致的.(3)周期性多孔结构在同时承受复杂多轴比例加载条件下, 剪切应力比
$\; {\rho _2} \gt 0 $ . 从承载极限来看, 随着面内剪应力的增加, 结构的极限拉伸强度显著降低, 而且双轴应力比$\; {\rho _1} $ 越大, 强度极限下降幅度越大; 当$ 0 \lt {\rho _2} \leqslant 0.5 $ 时, 随着$\; {\rho _1} $ 增加, 多孔结构的抗拉强度先提高再降低; 当$ 0.5 \lt {\rho _2} \leqslant 1.0 $ 时, 结构的抗拉强度随着$\; {\rho _1} $ 的增加而持续下降, 但下降幅度不大. 从裂纹萌生与扩展行为来看, 周期性多孔结构展现出S型、双弧线型和$ {45^ \circ } $ 斜裂纹型3种模式; 对于前两者,$\; {\rho _1} $ 和$\; {\rho _2} $ 的增大时, 裂纹扩展后与水平方向的夹角逐渐增大; 对于后者,$\; {\rho _2} $ 的增大对其裂纹扩展路径并无显著影响. -
![]()
图 3 弥散裂纹与弹性体受力情况, 左侧为尖裂纹, 右侧为弥散裂纹
Figure 3. Diffused crack and elastomeric stresses, with sharp cracks on the left and diffused crack on the right
![]()
图 5 周期性多孔结构RVE模型的几何构型(单位: mm)
Figure 5. Geometric configuration of the RVE model for periodic porous structures (unit: mm)
![]()
图 6
$ 8 \times 8 $ 多孔模型的单轴拉伸边界条件Figure 6. Boundary conditions for uniaxial stretch of
$ 8 \times 8 $ periodic porous model![]()
图 7 单轴拉伸载荷下单元数分别为80000, 115200和156800的周期性多孔结构的2×2 RVE模型名义应力-应变曲线的比较
Figure 7. Comparison of nominal stress-strain curves of the RVE model with 80000, 115200 and 156800 elements under tensile loading condition
![]()
图 8 单轴拉伸载荷下周期性多孔结构的
$ 2 \times 2 $ RVE模型与$ 8 \times 8 $ 模型的名义应力-应变曲线的比较Figure 8. Comparison of nominal stress-strain curves of the
$ 2 \times 2 $ RVE model and the$ 8 \times 8 $ periodic porous model under uniaxial tensile loading condition![]()
图 9 单轴拉伸载荷下的周期性多孔结构的
$ 2 \times 2 $ RVE模型与$ 8 \times 8 $ 模型的断裂相场云图Figure 9. Phase-field fracture contours of the
$ 2 \times 2 $ RVE model and$ 8 \times 8 $ model for periodic porous structure under uniaxial tensile loading![]()
![]()
图 11 周期性多孔结构在竖直方向的承载极限
${P_{22 \text{-} \max }}$ 随应力比$ {\rho _1} $ 和$ {\rho _2} $ 变化的曲线: (a)保持应力比$ {\rho _2} $ 不变, (b)保持应力比$ {\rho _1} $ 不变, (c)局部放大图Figure 11. Changing curves of 2-directional extreme load
${P_{22 \text{-} \max }}$ with$ {\rho _1} $ and$ {\rho _2} $ for periodic porous structures: (a) keep$ {\rho _2} $ unchanged, (b) keep$ {\rho _1} $ unchanged and (c) partial enlarged graph![]()
图 12 周期性多孔结构在不同多轴加载条件下的裂纹扩展模式相图
Figure 12. Phase diagram of crack propagation mode for periodic porous structure under different multiaxial proportional loading
表 1 双轴比例加载下的周期性多孔结构的承载极限
${P_{22 \text{-} \max }}$ 与断裂相场云图Table 1 Extreme load
${P_{22 \text{-} \max }}$ and phase-field fracture contours for the RVE model under biaxial proportional loading condition$ {\rho _1} $ −1 −0.5 0 0.5 1 P22-max/MPa 702.26 782.25 867.91 947.08 971.39 Phase-
field fracture contours
下载: 导出CSV
表 2 多轴比例加载下的周期性多孔结构的断裂相场云图
Table 2 Phase-field fracture contours of periodic porous structure model under multiaxial proportional loading
$ {\rho _1} $, $ {\rho _2} $ −1 −0.5 0 0.5 1.0 0 0.25 0.5 0.75 1.0
下载: 导出CSV
表 3 5种不同裂纹扩展模式的表征符号和渐进破坏的断裂相场云图
Table 3 Representative symbols and progressive phase-field fracture contours for the five distinguished crack propagation modes
Representative symbol crack propagation modes horizontal
typeS-type double arc-type 45°-type orthogonal cross-type progressive phase field
fracture contours
下载: 导出CSV
-
[1] 张洪武, 王鲲鹏. 材料非线性微-宏观分析的多尺度方法研究. 力学学报, 2004, 36(3): 359-363 (Zhang Hongwu, Wang Kunpeng. Multiscale methods for nonlinear analysis of composite materials. Chinese Journal of Theoretical and Applied Mechanics, 2004, 36(3): 359-363 (in Chinese) [2] Liu QH, Liu XW, Zhang CZ, et al. A novel multiscale porous composite structure for sound absorption enhancement. Composite Structures, 2021, 276(22): 114456
[3] 樊宣青, 刑誉锋. 周期多孔材料力学性能分析中的多尺度均匀化方法的精度// 第十二届全国振动理论及应用学术会议论文集. 中国振动工程学会, 2017 Fan Xuanqing, Xing Yufeng. The accuracy of multi-scale homogenization method in the analysis of mechanical properties of porous materials//Proceedings of the 12th National Conference on Vibration Theory and Applications. Chinese Society for Vibration Engineering, 2017 (in Chinese)
[4] Hayes AM, Wang A, Dempsey BM, et al. Mechanics of linear cellular alloys. Mechanics of Materials, 2004, 36(8): 691-713 doi: 10.1016/j.mechmat.2003.06.001
[5] 杜映洪. 层级多孔结构的力学及断裂性能研究. [硕士论文]. 天津: 天津大学, 2016 Du Yinghong. On the mechanical and fracture properties of the hierarchical cellular structure. [Master Thesis]. Tianjin: Tianjin University, 2016 (in Chinese)
[6] Jelitto H, Schneider GA. Extended cubic fracture model for porous materials and the dependence of the fracture toughness on the pore size. Materialia, 2020, 12(4): 100761
[7] 吴建营. 固体结构损伤破坏统一相场理论, 算法和应用. 力学学报, 2021, 53(2): 301-329 (Wu Jianying. On the unified phase-field theory for damage and failure in solids and structures: theoretical and numerical aspects. Chinese Journal of Theoretical and Applied Mechanics, 2021, 53(2): 301-329 (in Chinese) [8] 卢广达, 陈建兵. 基于一类非局部宏-微观损伤模型的裂纹模拟. 力学学报, 2020, 52(3): 749-762 (Lu Guangda, Chen Jianbing. Cracking simulation based on a nonlocal macro-meso-scale damage model. Chinese Journal of Theoretical and Applied Mechanics, 2020, 52(3): 749-762 (in Chinese) [9] Fan R, Fish J. The rs-method for material failure simulations. International Journal for Numerical Methods in Engineering, 2008, 73(11): 1607-1623 doi: 10.1002/nme.2134
[10] Xu XP, Needleman A. Numerical simulations of fast crack growth in brittle solids. Journal of the Mechanics and Physics of Solids, 1994, 42(9): 1397-1434 doi: 10.1016/0022-5096(94)90003-5
[11] 李录贤, 王铁军. 扩展有限元法(XFEM)及其应用. 力学进展, 2005, 35(1): 5-21 (Li Luxian, Wang Tiejun. The extended finite method and its application. Advances in Mechanics, 2005, 35(1): 5-21 (in Chinese) doi: 10.3321/j.issn:1000-0992.2005.01.002 [12] 王理想, 文龙飞, 肖桂仲等. 非连续问题中单元分割的模板方法. 力学学报, 2021, 53(3): 823-836 (Wang Lixiang, Wen Longfei, Xiao Guizhong, et al. A templated method for partitioning of solid elements in discontinuous problems. Chinese Journal of Theoretical and Applied Mechanics, 2021, 53(3): 823-836 (in Chinese) [13] Wu JY, Qiu JF, Nguyen VP, et al. Computational modeling of localized failure in solids: XFEM vs PF-CZM. Computer Methods in Applied Mechanics and Engineering, 2019, 345(3): 618-643
[14] Griffith AA. The phenomena of rupture and flow in solids. Philosophical Transactions of The Royal Society A: Mathematical Physical and Engineering Sciences, 1920, A221(4): 163-198
[15] Francfort GA, Marigo JJ. Revisiting brittle fracture as an energy minimization problem. Journal of the Mechanics and Physics of Solids, 1998, 46(8): 1319-1342 doi: 10.1016/S0022-5096(98)00034-9
[16] Bourdin B, Francfort GA, Marigo JJ. Numerical experiments in revisited brittle fracture. Journal of the Mechanics and Physics of Solids, 2000, 48(4): 797-826 doi: 10.1016/S0022-5096(99)00028-9
[17] Molnár G, Gravouil A. 2D and 3D Abaqus implementation of a robust staggered phase-field solution for modeling brittle fracture. Finite Elements in Analysis and Design, 2017, 130(5): 27-38
[18] Ambati M, Gerasimov T, Delorenzis L. Phase-field modeling of ductile fracture. Computational Mechanics, 2015, 55(5): 1017-1040 doi: 10.1007/s00466-015-1151-4
[19] Fang JG, Wu CQ, Li J, et al. Phase field fracture in elasto-plastic solids: Variational formulation for multi-surface plasticity and effects of plastic yield surfaces and hardening. International Journal of Mechanical Sciences, 2019, 156(6): 382-396
[20] Geelen JM, Liu Y, Hu T, et al. A phase-field formulation for dynamic cohesive fracture. Computer Methods in Applied Mechanics and Engineering, 2019, 348(6): 680-711
[21] Lo YS, Borden MJ, Ravi CK, et al. A phase-field model for fatigue crack growth. Journal of the Mechanics and Physics of Solids, 2019, 132(11): 103684
[22] Clayton JD, Knap J. Phase field modeling of directional fracture in anisotropic polycrystals. Computational Materials Science, 2015, 98(3): 158-169
[23] 吴建营, 陈万昕, 黄羽立. 基于统一相场理论的早龄期混凝土化-热-力多场耦合裂缝模拟与抗裂性能预测. 力学学报, 2021, 53(5): 1367-1382 (Wu Jianying, Chen Wanxin, Huang Yuli. Computational modeling of shrinkage induced cracking in early-age concrete based on the unified phase-field theory. Chinese Journal of Theoretical and Applied Mechanics, 2021, 53(5): 1367-1382 (in Chinese) [24] Noll T, Kuhn C, Olesch D, et al. 3D phase field simulations of ductile fracture. GAMM-Mitteilungen, 2020, 43(2): 1-16
[25] 王博臣, 侯玉亮, 夏凉等. 基于子结构法与损伤识别的周期性结构脆性断裂相场模拟. 航空学报, 2022, 43(3): 293-304 (Wang Bochen, Hou Yuliang, Xia Liang, et al. Phase field modeling of the brittle fracture of periodic structures based on substructuring and damage identification. Journal of Aeronautics, 2022, 43(3): 293-304 (in Chinese) [26] 郑晓霞, 郑锡涛, 缑林虎. 多尺度方法在复合材料力学分析中的研究进展. 力学进展, 2010, 40(1): 42-52 (Zheng Xiaoxia, Zheng Xitao, Gou Linhu. The Research progress on multiscale method foe the mechanical analysis of composites. Advance in Mechanics, 2010, 40(1): 42-52 (in Chinese) doi: 10.6052/1000-0992-2010-1-J2008-104 [27] Tyrus JM, Gosz M, Desantiago E. A local finite element implementation for imposing periodic boundary conditions on composite micromechanical models. International Journal of Solids & Structures, 2007, 44(9): 2972-2989
[28] 卫宇璇. 周期性复合材料结构力学性能的多尺度分析. [硕士论文]. 吉林: 吉林大学, 2018 Wei Yuxuan. Multi-scale analysis of mechanical properties of periodic composite structures. [Master Thesis]. Jilin: Jilin University, 2018 (in Chinese)
[29] Liang X, Crosby AJ. Uniaxial stretching mechanics of cellular flexible metamaterials. Extreme Mechanics Letters, 2020, 35(2): 100637
[30] Sun Q, Zhou G, Meng Z, et al. Failure criteria of unidirectional carbon fiber reinforced polymer composites informed by a computational micromechanics model. Composites Science and Technology, 2019, 172(4): 81-95
[31] Msekh MA, Sargado JM, Jamshidian M, et al. Abaqus implementation of phase-field model for brittle fracture. Computational Materials Science, 2015, 96(1B): 472-484
[32] Miehe C, Welschinger F, Hofacker M. Thermodynamically consistent phase‐field models of fracture: Variational principles and multi‐field FE implementations. International Journal for Numerical Methods in Engineering, 2010, 83(10): 1274-1311
[33] Fang JG, Wu CQ, Rabczuk T, et al. Phase field fracture in elasto-plastic solids: Abaqus implementation and case studies. Theoretical and Applied Fracture Mechanics, 2019, 103(5): 102252
[34] Guo YD, Huang W, Ma YE. Buckling pattern transition of periodic porous elastomers induced by proportional loading conditions. International Journal of Applied Mechanics, 2021, 13(6): 2150067 doi: 10.1142/S1758825121500678
[35] 刘国威, 李庆斌, 左正. 相场断裂模型分步算法在ABAQUS中的实现. 岩石力学与工程学报, 2016, 35(5): 1020-1029 (Liu Guowei, Li Qingbin, Zuo Zheng. Implementation of phase field fracture model step algorithm in ABAQUS. Journal of Rock Mechanics and Engineering, 2016, 35(5): 1020-1029 (in Chinese) doi: 10.13722/j.cnki.jrme.2015.1264 [36] Sarac B, Schroers J. Designing tensile ductility in metallic glasses. Nature Communications, 2013, 4: 2158 doi: 10.1038/ncomms3158
-
期刊类型引用(4)
1. 许爽，邓庆田，李旺斐，李新波，宋学力，温金鹏. 层合多孔结构的裂纹扩展规律研究. 机械强度. 2025(02): 28-36 . 
百度学术
2. 辛雨柯，邓庆田，宋学力，李新波，温金鹏. 含预制缺陷PLA蜂窝加筋结构承载力分析. 装备环境工程. 2024(01): 26-34 . 百度学术
3. 华军，黄磊，杨亚东，邢小茹，朱正洪. CPCs泡沫微结构仿生构筑及其压缩力学性能研究. 固体力学学报. 2024(05): 665-678 . 百度学术
4. 魏浩林，吴振，任晓辉. 基于各向异性断裂相场模型的自适应多尺度有限元方法. 力学学报. 2024(10): 2902-2912 . 本站查看
其他类型引用(4)
计量
- 文章访问数: 1024
- HTML全文浏览量: 329
- PDF下载量: 184
- 被引次数: 8
