LARGE DEFORMATION AND FRACTURE ANALYSIS OF THIN PLATE BENDING BASED ON PERIDYNAMICS
-
摘要: 薄板结构仅在较小的荷载下就能产生较大的位移、旋转, 甚至引发结构产生裂纹并扩展, 进而发生结构整体断裂, 因此, 建立薄板结构在大变形过程中的裂纹扩展及断裂仿真模型, 具有重要的工程实际意义. 文章建立了用于薄板结构几何大变形和断裂分析的近场动力学(PD)和连续介质力学(CCM)耦合模型. 首先基于冯·卡门假设, 采用更新的拉格朗日法得到薄板在几何大变形增量步下的虚应变能密度增量公式, 并利用虚功原理和均质化假设求出几何大变形微梁键的本构模型参数; 接着分别建立几何大变形薄板PD模型与CCM模型的虚应变能密度增量, 并建立了薄板几何大变形PD-CCM耦合模型; 最后模拟了薄板结构在横向变形作用下的渐进断裂过程, 得到与实验结果高度一致的仿真结果, 验证了所提出的几何非线性PD-CCM耦合模型的精度. 结果表明: 本文所提出的薄板PD-CCM耦合模型具有简单高效, 无需考虑材料参数限制和边界效应的特点, 可以很好地用于预测薄板结构在几何大变形过程中的局部损伤和结构断裂, 有利于薄板结构的断裂安全评价和理论发展.Abstract: Thin plate structures are widely used in the fields of automobiles, ships, and aerospace because of their excellent load-bearing performance, light weight and easy processing. However, in practical applications, thin plate structures often produce large displacement, rotation and even cause crack initiation and growth under small loads, and then the overall structure fractures. Therefore, it is of great engineering practical significance to establish a crack growth and fracture simulation model of thin plate structures in the process of large deformation. In this paper, a peridynamic (PD) and classical continuum mechanics (CCM) coupling model for geometrically nonlinear deformation and fracture analysis of thin plate structures is established. First of all, the updated Lagrangian formula is used to obtain the expression of virtual strain energy density increment of thin plates at each increment step in large deformation analysis under von Karman's hypothesis. Then, the PD constitutive parameters of geometrically nonlinear micro-beam bond are obtained by using the virtual work principle and homogenization hypothesis. After that, the virtual strain energy density increments of the PD model and CCM model for geometrically large deformation thin plate were respectively established, and the geometrically large deformation PD-CCM coupling model of the thin plate was established. Finally, the progressive fracture process of the thin plate structure under the action of lateral deformation is simulated, and the simulation results are highly consistent with the experimental results, which verifies the accuracy of the proposed geometrically nonlinear PD-CCM coupling model. It is shown that the proposed geometrically nonlinear PD-CCM coupling model is simple and efficient without restriction on material parameters and consideration of boundary effects, and can be well used to predict local damage and structural fracture of thin plate structures during geometrically large deformations. It is beneficial to the fracture safety evaluation and theoretical development of thin plate structures.
-
Keywords:
- peridynamics /
- thin plate bending /
- large deformation /
- coupling /
- crack growth
-
引 言
板壳结构以其优良的承载性能、轻量化和易于加工等特性, 被广泛应用于汽车、船舶以及航空航天等领域[1-2]. 而在实际应用中经常可以看到薄板结构的大位移、大转动和小弹性应变行为. 因此, 针对薄板结构因大变形产生失效的现象, 建立可靠的薄板弯曲大变形和断裂分析的数值模拟方法, 具有重要的工程实际意义.
在经典连续介质力学CCM的非线性有限元分析中, 运动方程是通过虚位移原理的线性化得到的[3], 然而由于使用的偏微分方程在不连续条件下会变得无解, 故不适用于求解断裂问题. Silling[4]为此提出了近场动力学PD理论, 该理论用空间积分方程代替偏微分方程, 使得对于不连续模型, PD理论无需任何处理即可参与求解计算[5-6]. 迄今为止, PD理论已经应用于许多领域, 如非弹性材料[7-8]、各向异性材料[9-10]的损伤与失效, 多晶断裂[11], 疲劳裂纹[12-13]和热力学响应[14-15]等.
对于薄板结构, Silling等[16]引入了第一个简化的二维PD模型, 后来由Taylor等[17]开发了适用于薄板的PD公式. O’Grady等[18]引入非常规态基的PD模型, 基于Kirchhoff-Love理论捕捉板的弯曲变形. Yang等[19-20]还为Kirchhoff板开发新的态基模型, 提高了计算效率, 并实现对功能梯度材料下PD模型的开发. Shen等[21-22]通过插值方法建立基于微梁键的PD梁和壳模型. Diyaroglu等[23]引入考虑横向剪切变形的Timoshenko梁和Mindlin板公式, 建立PD梁和板弯曲模型, 此外, Hu等[24]提出的新的态型近场动力学模型解决了考虑横向剪切变形的复合材料层合板的建模问题. 然而上述模型都只考虑小的弹性变形, 不能分析大变形和大转动的板壳结构. Nguyen等[25]最近提出的常规态基几何非线性板的PD模型基于虛位移原理, 利用全拉格朗日法, 建立了PD板的非线性应变能密度和非线性运动方程, 并给出数值计算方法, 验证了模型的准确性, 同时给出损伤预测的数值模型. 但在实际应用中, 更新的拉格朗日法的几何非线性部分表达式比全拉格朗日法更简单, 并且利用共旋理论可以更容易地计算出应变和应力, 因此更新的拉格朗日法更具有实用性. 与此同时, 基于变形的耦合策略[26]还提供一种直接有效的方法极大的节省计算成本, 并且通过耦合的方法可以解决原本无法方便地对PD模型施加边界条件这一问题[27].
据此, 本文采用更新的拉格朗日法, 提出一种新的用于薄板弯曲几何非线性变形和断裂分析的PD和CCM耦合模型. 利用虚功原理和均质化假设推导几何非线性微梁键[22,28]的非线性PD参数, 随后采用不连续有限元方法进行数值求解, 通过建立虚应变能密度增量的平衡, 提出PD和CCM薄板耦合模型, 最后通过模拟在大横向变形下薄板结构的渐进断裂过程来验证该模型的准确性和有效性.
1. 几何非线性板的更新的拉格朗日法
1.1 广义应变和应力
如图1所示, 一个具有4个节点的矩形板单元, 局部坐标系记为
$OXYZ$ . 本文以基于Kirchhoff-Love假设下的薄板结构作为分析对象, 忽略板单元沿厚度方向的剪切变形, 同时引入冯·卡门假设[29], 只计算板单元截面中心线产生的大挠度、小应变变形, 于是Green应变可表示为$$ \begin{split} &{{\boldsymbol{E}}_{{\rm{Green}}}} = \left\{ {\begin{array}{*{20}{c}} {{E_{11}}} \\ {{E_{22}}} \\ {2{E_{12}}} \end{array}} \right\} = \\ &\qquad \left\{ {\begin{array}{*{20}{c}} {\dfrac{{\partial {u_1}}}{{\partial X}} + \dfrac{1}{2}{{\left( {\dfrac{{\partial {u_3}}}{{\partial X}}} \right)}^2} + Z\dfrac{{\partial {\theta _2}}}{{\partial X}}} \\ {\dfrac{{\partial {u_2}}}{{\partial Y}} + \dfrac{1}{2}{{\left( {\dfrac{{\partial {u_3}}}{{\partial Y}}} \right)}^2} - Z\dfrac{{\partial {\theta _1}}}{{\partial Y}}} \\ {\dfrac{{\partial {u_1}}}{{\partial Y}} + \dfrac{{\partial {u_2}}}{{\partial X}} + \dfrac{{\partial {u_3}}}{{\partial X}}\dfrac{{\partial {u_3}}}{{\partial Y}} + Z\left( {\dfrac{{\partial {\theta _2}}}{{\partial Y}} - \dfrac{{\partial {\theta _1}}}{{\partial X}}} \right)} \end{array}} \right\}\end{split} $$ (1) 其中
${u_1},{u_2},{u_3}$ 分别表示未变形板的中面位移,${\theta _1},{\theta _2}$ 分别表示中面法线绕$X$ 轴和$Y$ 轴的旋转角, 由于忽略剪切变形, 这些法线在变形后垂直于中面.上述Green应变可以分解为单元的薄膜应变
$ {\boldsymbol{E}}_L^p $ 、弯曲应变$ {\boldsymbol{E}}_L^b $ 和薄膜应变的非线性分量$ {\boldsymbol{E}}_N^p $ $$ \left.\begin{split} &{{\boldsymbol{E}}}_{L}^{p} = \left\{\begin{array}{c}\dfrac{\partial {u}_{1}}{\partial X}\\ \dfrac{\partial {u}_{2}}{\partial Y}\\ \dfrac{\partial {u}_{1}}{\partial Y} + \dfrac{\partial {u}_{2}}{\partial X}\end{array}\right\}\text{, }{{\boldsymbol{E}}}_{L}^{b} = \left\{\begin{array}{c}\dfrac{\partial {\theta }_{2}}{\partial X}\\ -\dfrac{\partial {\theta }_{1}}{\partial Y}\\ \dfrac{\partial {\theta }_{2}}{\partial Y}-\dfrac{\partial {\theta }_{1}}{\partial X}\end{array}\right\}\\ &{{\boldsymbol{E}}}_{N}^{p} = \left\{\begin{array}{c}\dfrac{1}{2}{\left(\dfrac{\partial {u}_{3}}{\partial X}\right)}^{2}\\ \dfrac{1}{2}{\left(\dfrac{\partial {u}_{3}}{\partial Y}\right)}^{2}\\ \dfrac{\partial {u}_{3}}{\partial X}\dfrac{\partial {u}_{3}}{\partial Y}\end{array}\right\}\end{split}\right\} $$ (2) 对于线弹性各向同性、均质材料, 在平面应力条件下, 应力在
$Z$ 方向为0, 与式(1)中Green应变相关的第二Piola应力矢量为$$ {{\boldsymbol{S}}_{{\rm{Piola}}}} = {{\hat {\boldsymbol{D}}}}{{\boldsymbol{E}}_{{\rm{Green}}}} = {\left\{ {\begin{array}{*{20}{c}} {{S_{11}}}&{{S_{22}}}&{{S_{12}}} \end{array}} \right\}^{\text{T}}} $$ (3) 其中
${{\hat {\boldsymbol{D}}}}$ 为材料参数矩阵, 假设单元杨氏模量为$E$ , 泊松比为$v$ , 有$$ {{\hat {\boldsymbol{D}}}} = \frac{E}{{(1 - {v^2})}}\left[ {\begin{array}{*{20}{c}} 1&v&0 \\ v&1&0 \\ 0&0&{\dfrac{{1 - v}}{2}} \end{array}} \right] $$ (4) 假设板单元厚度为
$h$ , 上述本构方程中应力对厚度方向积分可得到由单位宽度上的薄膜力$ {{\boldsymbol{S}}^p} $ 和弯曲力矩$ {{\boldsymbol{S}}^b} $ 表示的广义应力矢量${\boldsymbol{S}}$ $$ {\boldsymbol{S}} = \left\{ {\begin{array}{*{20}{c}} {{{\boldsymbol{S}}^p}} \\ {{{\boldsymbol{S}}^b}} \end{array}} \right\} $$ (5) $$ \left. \begin{aligned} & {{\boldsymbol{S}}^p} = {\left\{ {\begin{array}{*{20}{c}} {{N_{11}}}&{{N_{22}}}&{{N_{12}}} \end{array}} \right\}^{{\rm{T}}}} = \\ &\qquad{\left[ {\int_{ - h/2}^{h/2} {\left( {\begin{array}{*{20}{c}} {{S_{11}}}&{{S_{22}}}&{{S_{12}}} \end{array}} \right)} {{\rm{d}}}Z} \right]^{{\rm{T}}}} = h{{\hat {\boldsymbol{D}}}}{{\boldsymbol{E}}^p} \\ & {{\boldsymbol{S}}^b} = {\left\{ {\begin{array}{*{20}{c}} {{M_{11}}}&{{M_{22}}}&{{M_{12}}} \end{array}} \right\}^{{\rm{T}}}} = \\ &\qquad{\left[ {\int_{ - h/2}^{h/2} {\left( {\begin{array}{*{20}{c}} {{S_{11}}}&{{S_{22}}}&{{S_{12}}} \end{array}} \right)Z} {{\rm{d}}}Z} \right]^{{\rm{T}}}} = \frac{{{h^3}}}{{12}}{{\hat {\boldsymbol{D}}}}{{\boldsymbol{E}}^b} \end{aligned} \right\} $$ (6) 与上述广义应力矢量对应的广义应变矢量
${\boldsymbol{E}}$ 可由薄膜应变和弯曲应变表示$$ {\boldsymbol{E}} = \left\{ {\begin{array}{*{20}{c}} {{{\boldsymbol{E}}^p}} \\ {{{\boldsymbol{E}}^b}} \end{array}} \right\} = \left\{ {\begin{array}{*{20}{c}} {{\boldsymbol{E}}_L^p + {\boldsymbol{E}}_N^p} \\ {{\boldsymbol{E}}_L^b} \end{array}} \right\} $$ (7) 由此可以得到广义第二Piola应力和广义Green应变表示的本构方程
$$ {\boldsymbol{S}} = {\boldsymbol{DE}} $$ (8) 其中
$$ {\boldsymbol{D}} = \left[ {\begin{array}{*{20}{c}} {h{{\hat {\boldsymbol{D}}}}}&{\boldsymbol{0}} \\ {\boldsymbol{0}}&{\dfrac{{{h^3}}}{{12}}{{\hat {\boldsymbol{D}}}}} \end{array}} \right],\,\,\,{\boldsymbol{0}} = \left[ {\begin{array}{*{20}{c}} 0&0&0 \\ 0&0&0 \\ 0&0&0 \end{array}} \right] $$ (9) 1.2 广义应变增量分解
板单元上任意一点在
$t + \Delta t$ 时刻的位移分量可定义为$$ \left.\begin{gathered} {{\bar u}_1} = {u_1} + \Delta {u_1},\;\;\;{{\bar u}_2} = {u_2} + \Delta {u_2},\;\;\;{{\bar u}_3} = {u_3} + \Delta {u_3}\; \\ {{\bar \theta }_1} = {\theta _1} + \Delta {\theta _1},\;\;\;{{\bar \theta }_2} = {\theta _2} + \Delta {\theta _2} \\ \end{gathered}\right\} $$ (10) 其中
${u_1},{u_2},{u_3},{\theta _1},{\theta _2}$ 和$\Delta {u_1},\Delta {u_2},\Delta {u_3},\Delta {\theta _1},\Delta {\theta _2}$ 分别表示$t$ 时刻的位移分量以及$t$ 时刻到$t + \Delta t$ 时刻的位移分量增量,$t + \Delta t$ 时刻的广义应变可以写成$$ {{\bar {\boldsymbol{E}}}} = {\boldsymbol{E}} + \Delta {\boldsymbol{E}} $$ (11) 由式(7)可知广义应变增量
$ \Delta {\boldsymbol{E}} $ 可表示为$$ \Delta {\boldsymbol{E}} = \left\{ {\begin{array}{*{20}{c}} {\Delta {{\boldsymbol{E}}^p}} \\ {\Delta {{\boldsymbol{E}}^b}} \end{array}} \right\} = \left\{ {\begin{array}{*{20}{c}} {\Delta {\boldsymbol{E}}_L^p + \Delta {\boldsymbol{E}}_N^p} \\ {\Delta {\boldsymbol{E}}_L^b} \end{array}} \right\} $$ (12) 式中广义应变增量可被分为线性增量应变
$\Delta {{\boldsymbol{E}}_L}$ 和非线性增量应变$\Delta {{\boldsymbol{E}}_N}$ , 即$$ \Delta {\boldsymbol{E}} = \Delta {{\boldsymbol{E}}_L} + \Delta {{\boldsymbol{E}}_N} = \left\{ {\begin{array}{*{20}{c}} {\Delta {\boldsymbol{E}}_L^p} \\ {\Delta {\boldsymbol{E}}_L^b} \end{array}} \right\} + \left\{ {\begin{array}{*{20}{c}} {\Delta {\boldsymbol{E}}_N^p} \\ {\boldsymbol{0}} \end{array}} \right\} $$ (13) 1.3 虚功增量
$t + \Delta t$ 时刻的增量虚功方程为$$ \int_{{a_0}} {{{ {\delta {{\bar {\boldsymbol{E}}}}} }^{\text{T}}}{{\bar {\boldsymbol{S}}}}{\text{d}}a} = \int_{{v_0}} {\delta {{{{\bar {\boldsymbol{u}}}}}^{\text{T}}}{{\boldsymbol{f}}_v}{\text{d}}v} + \int_{{a_0}} {\delta {{{{\bar {\boldsymbol{u}}}}}^{\text{T}}}{{\boldsymbol{f}}_a}{\text{d}}a} $$ (14) 其中,
${{\bar {\boldsymbol{S}}}}$ 和${{\bar {\boldsymbol{u}}}}$ 分别表示$t + \Delta t$ 时刻的广义应力矢量和位移矢量,$ {{\boldsymbol{f}}_v} $ 和$ {{\boldsymbol{f}}_a} $ 分别表示板单元体力和面力矢量, 对应积分域${v_0}$ 和${a_0}$ 分别表示板单元的体积和表面积. 根据式(13), 式(14)左侧可表示为$$ \begin{split} & \int_{{a_0}} {{{ {\delta {{\bar {\boldsymbol{E}}}}} }^{\text{T}}}{{\bar {\boldsymbol{S}}}}{\text{d}}a} = \int_{{a_0}} \delta \left( {\Delta {{\boldsymbol{E}}^{{\rm{T}}}}} \right){{\bar {\boldsymbol{S}}}}{\text{d}}a =\\ &\qquad \int_{{a_0}} \delta \left( {\Delta {{\boldsymbol{E}}_L}^{{\rm{T}}} + \Delta {{\boldsymbol{E}}_N}^{{\rm{T}}}} \right)({\boldsymbol{S}} + \Delta {\boldsymbol{S}}){\text{d}}a = \\ &\qquad \int_{{a_0}} \delta \left( {\Delta {{\boldsymbol{E}}_L}^{{\rm{T}}}} \right){\boldsymbol{S}}{\text{d}}a + \int_{{a_0}} \delta \left( {\Delta {{\boldsymbol{E}}_N}^{{\rm{T}}}} \right){\boldsymbol{S}}{\text{d}}a + \\ &\qquad\int_{{a_0}} \delta \left( {\Delta {{\boldsymbol{E}}_L}^{{\rm{T}}} + \Delta {{\boldsymbol{E}}_N}^{{\rm{T}}}} \right)(\Delta {\boldsymbol{S}}){\text{d}}a \end{split} $$ (15) 第1项为结构反力引起的虚应变能增量
$$ \Delta \omega _{CCM}^R = \int_{{a_0}} \delta \left( {\Delta {{\boldsymbol{E}}_L}^{{\rm{T}}}} \right){\boldsymbol{S}}{\text{d}}a $$ (16) 第2项为初始应力引起的虚应变能增量
$$ \Delta \omega _{CCM}^S = \int_{{a_0}} \delta \left( {\Delta {{\boldsymbol{E}}_N}^{{\rm{T}}}} \right){\boldsymbol{S}}{\text{d}}a $$ (17) 第3项可以表示为
$$ \begin{split} & \int_{{a_0}} \delta \left( {\Delta {{\boldsymbol{E}}_L}^{{\rm{T}}} + \Delta {{\boldsymbol{E}}_N}^{{\rm{T}}}} \right)(\Delta {\boldsymbol{S}}){\text{d}}a = \\ &\qquad\int_{{a_0}} \delta \left( {\Delta {{\boldsymbol{E}}_L}^{{\rm{T}}} + \Delta {{\boldsymbol{E}}_N}^{{\rm{T}}}} \right){\boldsymbol{D}}(\Delta {{\boldsymbol{E}}_L} + \Delta {{\boldsymbol{E}}_N}){\text{d}}a = \\ &\qquad \Delta \omega _{CCM}^L + \Delta \omega _{CCM}^N \end{split} $$ (18) 其中
$ \Delta \omega _{CCM}^L $ 和$ \Delta \omega _{CCM}^N $ 分别为线性应变和非线性应变引起的虚应变能增量$$ \left.\begin{split} & \Delta \omega _{CCM}^L = \int_{{a_0}} \delta \left( {\Delta {{\boldsymbol{E}}_L}^{{\rm{T}}}} \right){\boldsymbol{D}}(\Delta {{\boldsymbol{E}}_L}){\text{d}}a \\ & \Delta \omega _{CCM}^N = \int_{{a_0}} \left[ \delta \left( {\Delta {{\boldsymbol{E}}_L}^{{\rm{T}}}} \right){\boldsymbol{D}}(\Delta {{\boldsymbol{E}}_N}) + \delta \left( {\Delta {{\boldsymbol{E}}_N}^{{\rm{T}}}} \right) \right.\\ &\qquad \left.{\boldsymbol{D}}(\Delta {{\boldsymbol{E}}_L}) +\delta \left( {\Delta {{\boldsymbol{E}}_N}^{{\rm{T}}}} \right){\boldsymbol{D}}(\Delta {{\boldsymbol{E}}_N}) \right] {\text{d}}a \end{split}\right\} $$ (19) 假设在每一个增量步下每个点的应变增量和初始应力都是均匀的,
$t$ 时刻的广义初始应力可由共旋理论进行直接累加得到$$ {\boldsymbol{S}} = \sum\limits_{k = 0}^n {{{(\Delta {\boldsymbol{S}})}_k}} = \sum\limits_{k = 0}^n {{\boldsymbol{D}}{{(\Delta {\boldsymbol{E}})}_k}} $$ (20) 其中
$k$ 为更新的拉格朗日法下$t$ 时刻之前的增量步数,$n$ 为$t$ 时刻之前的总增量步,$ {(\Delta {\boldsymbol{S}})_k} $ 和$ {(\Delta {\boldsymbol{E}})_k} $ 分别为第$k$ 增量步时的应力增量和应变增量.因此在CCM理论的框架内, 在一个增量步内的虚应变能密度增量可以表示为
$$\left. \begin{split} & \Delta W_{CCM}^L\left( {{{\boldsymbol{x}}_i}} \right) = \frac{1}{h}\delta \left( {\Delta {\boldsymbol{E}}_L^{\text{T}}} \right){\boldsymbol{D}}\left( {\Delta {{\boldsymbol{E}}_L}} \right) \\ & \Delta W_{CCM}^R\left( {{{\boldsymbol{x}}_i}} \right) = \frac{1}{h}\delta \left( {\Delta {\boldsymbol{E}}_L^{\text{T}}} \right){\boldsymbol{S}} = \frac{1}{h}\delta \left( {\Delta {\boldsymbol{E}}_L^{\text{T}}} \right)\sum\limits_{k = 0}^n {\boldsymbol{D}} {(\Delta {\boldsymbol{E}})_k} \\ & \Delta W_{CCM}^S\left( {{{\boldsymbol{x}}_i}} \right) = \frac{1}{h}\delta \left( {\Delta {\boldsymbol{E}}_N^{\text{T}}} \right){\boldsymbol{S}} = \frac{1}{h}\delta \left( {\Delta {\boldsymbol{E}}_N^{\text{T}}} \right)\sum\limits_{k = 0}^n {\boldsymbol{D}} {(\Delta {\boldsymbol{E}})_k} \\ & \Delta W_{CCM}^N\left( {{{\boldsymbol{x}}_i}} \right) = \frac{1}{h}\left[ \delta \left( {\Delta {\boldsymbol{E}}_L^{\text{T}}} \right){\boldsymbol{D}}\left( {\Delta {{\boldsymbol{E}}_N}} \right) + \delta \left( {\Delta {\boldsymbol{E}}_N^{\text{T}}} \right){\boldsymbol{D}}\left( {\Delta {{\boldsymbol{E}}_L}} \right) + \right.\\ &\qquad \left.\delta \left( {\Delta {\boldsymbol{E}}_N^{\text{T}}} \right){\boldsymbol{D}}\left( {\Delta {{\boldsymbol{E}}_N}} \right) \right] \end{split} \right\}$$ (21) 2. 基于微梁键的非线性近场动力学模型
2.1 经典键基近场动力学理论
在经典键基PD[4]模型中,
$t$ 时刻物质点${{\boldsymbol{x}}_i}$ 的运动方程可以表示为$$ \rho {\boldsymbol{\ddot u}}({{\boldsymbol{x}}_i},t) = \int_{{H_{{{\boldsymbol{x}}_i}}}} {{\boldsymbol{f}}\left( {{\boldsymbol{u}}({{\boldsymbol{x}}_j},t} ) - {\boldsymbol{u}}({{\boldsymbol{x}}_i},t),{{\boldsymbol{x}}_j} - {{\boldsymbol{x}}_i}\right){{\rm{d}}}{V_j}} + {\boldsymbol{b}}({{\boldsymbol{x}}_i},t) $$ (22) 其中
$ \rho $ 为参考构形中材料的密度,$ {\boldsymbol{u}} $ 为位移场矢量,${\boldsymbol{f}}$ 是物质点${{\boldsymbol{x}}_j}$ 施加在物质点${{\boldsymbol{x}}_i}$ 上的力矢量,${H_{{{\boldsymbol{x}}_i}}}$ 是物质点${{\boldsymbol{x}}_i}$ 的近场域,${\boldsymbol{b}}({{\boldsymbol{x}}_i},t)$ 表示施加在物质点${{\boldsymbol{x}}_i}$ 上的体积力函数. 为了便于描述, 两个物质点的相对位置和位移表示为${\boldsymbol{\xi }}$ 和${\boldsymbol{\eta }}$ $$ {\boldsymbol{\xi }} = {{\boldsymbol{x}}_j} - {{\boldsymbol{x}}_i},\,\,\,\,\,{\boldsymbol{\eta }} = {\boldsymbol{u}}({{\boldsymbol{x}}_i},t) - {\boldsymbol{u}}({{\boldsymbol{x}}_i},t) $$ (23) 根据文献[4]定义的微弹性材料模型, 成对力由标量微势能函数
$\omega ({\boldsymbol{\eta }},{\boldsymbol{\xi }})$ 导出$$ {\boldsymbol{f}}({\boldsymbol{\eta }},{\boldsymbol{\xi }}) = \frac{{\partial \omega ({\boldsymbol{\eta }},{\boldsymbol{\xi }})}}{{\partial {\boldsymbol{\eta }}}} $$ (24) 在线弹性理论下, 成对力的大小与键的相对伸长率之间呈线性关系, 微势能函数表示为
$$ \omega ({\boldsymbol{\eta }},{\boldsymbol{\xi }}) = \frac{1}{2}c({\boldsymbol{\eta }},{\boldsymbol{\xi }}){s^2}\xi $$ (25) 其中
$c$ 是作为本构参数的微模量常数,$\xi $ 是原始键长,$s$ 是键的相对伸长率, 其定义如下$$ s = \frac{{\left| {{\boldsymbol{\xi }} + {\boldsymbol{\eta }}} \right| - \left| {\boldsymbol{\xi }} \right|}}{{\left| {\boldsymbol{\xi }} \right|}} $$ (26) 给定物质点
${{\boldsymbol{x}}_i}$ 的应变能密度可以通过对近场域积分获得$$ {W_{PD}}({{\boldsymbol{x}}_i}) = \frac{1}{2}\int_{{H_{{{\boldsymbol{x}}_i}}}} {\omega ({\boldsymbol{\eta }},{\boldsymbol{\xi }})} {{\rm{d}}}{V_{{{\boldsymbol{x}}_j}}} $$ (27) 其中系数
${1}/{2}$ 是由于键中的每个物质点具有总键能的一半. 通过均匀变形条件下的经典线弹性均质体的应变能密度${W_{CCM}}({{\boldsymbol{x}}_i})$ 与${W_{PD}}({{\boldsymbol{x}}_i})$ 相等可以求得微模量参数$c$ .2.2 几何非线性近场动力学板模型的微模量
对于大变形条件下的几何非线性PD模型, 可以采用虚应变能密度增量来建立, 假设在每一个增量步中, PD模型下的每个材料点的增量应变和初始应力都是均匀的, 因此, 几何非线性PD板的虚应变能增量可以很容易地从均质化假设得到.
假设几何非线性PD板模型中的键是一种类似于CCM理论框架下梁单元的微梁, 该微梁键也具有几何非线性特性[28]. 如图2所示, 微梁键的局部坐标系记为
$oxyz$ , 节点${{\boldsymbol{x}}_i}$ 和其相邻点$ {{\boldsymbol{x}}_j} $ 位于半径为$\delta $ 的近场域${H_{{{\boldsymbol{x}}_i}}}$ 中, 键的中心线上任意一点的位移分别定义为${u'_1}(x)$ ,${u'_2}(x)$ 和${u'_3}(x)$ .对于该微梁键上的任意点在
$t + \Delta t$ 时刻的广义应变可以写成$$ {\boldsymbol{\bar \varepsilon }} = {\boldsymbol{\varepsilon }} + \Delta {\boldsymbol{\varepsilon }} $$ (28) 其中
${\boldsymbol{\varepsilon }}$ 为$t$ 时刻微梁键的广义应变,$\Delta {\boldsymbol{\varepsilon }}$ 为微梁键$t$ 到$t + \Delta t$ 时刻的广义应变增量, 可以定义为$$ \Delta {\boldsymbol{\varepsilon }} = {\left[ {\begin{array}{*{20}{c}} {\Delta {\varepsilon _x}}&{\Delta {\kappa _y}}&{\Delta {\kappa _z}}&{\Delta {\phi _x}} \end{array}} \right]^{\text{T}}} $$ (29) 其中
$\Delta {\varepsilon _x}$ ,$\Delta {\kappa _y}$ ,$\Delta {\kappa _z}$ 和$\Delta {\phi _x}$ 分别表示微梁键的轴向应变、$y$ 和$z$ 方向的曲率以及单位长度的扭转角增量. 按照冯卡门大挠度理论, 假设微梁键也产生大挠度、小应变变形, 即$$ \left.\begin{split} & \Delta {\varepsilon _x} = \frac{{\partial \Delta {{u'}_1}}}{{\partial x}} + \frac{1}{2}{\left( {\frac{{\partial \Delta {{u'}_2}}}{{\partial x}}} \right)^2} + \frac{1}{2}{\left( {\frac{{\partial \Delta {{u'}_3}}}{{\partial x}}} \right)^2} \\ & \Delta {\kappa _y} = \frac{{{\partial ^2}\Delta {{u'}_3}}}{{\partial {x^2}}},\;\;\;\Delta {\kappa _z} = \frac{{{\partial ^2}\Delta {{u'}_2}}}{{\partial {x^2}}},\;\;\;\Delta {\phi _x} = \frac{{\partial \Delta {{\theta '}_1}}}{{\partial x}} \end{split}\right\} $$ (30) 参考式(13), 式(30)中的广义应变增量亦可分为线性应变增量和非线性应变增量, 即
$$ \left.\begin{split} & \Delta {\boldsymbol{\varepsilon }} = \Delta {{\boldsymbol{\varepsilon }}_L} + \Delta {{\boldsymbol{\varepsilon }}_N} \\ & \Delta {{\boldsymbol{\varepsilon }}_L} = {\left\{ {\begin{array}{*{20}{c}} {\dfrac{{\partial \Delta {{u'}_1}}}{{\partial x}}}&{\dfrac{{{\partial ^2}\Delta {{u'}_3}}}{{\partial {x^2}}}}&{\dfrac{{{\partial ^2}\Delta {{u'}_2}}}{{\partial {x^2}}}}&{\dfrac{{\partial \Delta {{\theta '}_1}}}{{\partial x}}} \end{array}} \right\}^{{\rm{T}}}} \\ & \Delta {{\boldsymbol{\varepsilon }}_N} = {\left\{ {\begin{array}{*{20}{c}} {\dfrac{1}{2}{{\left( {\dfrac{{\partial \Delta {{u'}_2}}}{{\partial x}}} \right)}^2} + \dfrac{1}{2}{{\left( {\dfrac{{\partial \Delta {{u'}_3}}}{{\partial x}}} \right)}^2}}&0&0&0 \end{array}} \right\}^{{\rm{T}}}} \end{split}\right\} $$ (31) 定义PD微模量参数矩阵
$$ {{\boldsymbol{D}}_{{\rm{PD}}}} = \left[ {\begin{array}{*{20}{c}} {{c_{ax}}}&0&0&0 \\ 0&{{c_{by}}}&0&0 \\ 0&0&{{c_{bz}}}&0 \\ 0&0&0&{{c_{tor}}} \end{array}} \right] $$ (32) 与微梁键的广义应变增量相对应的广义应力增量可由非线性微梁键的本构关系表示为
$$ \Delta {{\boldsymbol{S}}_{{\rm{PD}}}} = {{\boldsymbol{D}}_{{\rm{PD}}}}\Delta {\boldsymbol{\varepsilon }} $$ (33) 其中广义应力增量
$\Delta {{\boldsymbol{S}}_{{\rm{PD}}}}$ 可由微梁键的轴向应力、$y$ 和$z$ 方向的弯矩以及单位长度的扭矩增量表示$$ \Delta {{\boldsymbol{S}}_{{\rm{PD}}}} = {\left\{ {\begin{array}{*{20}{c}} {\Delta F_x^{{\rm{PD}}}}&{\Delta M_y^{{\rm{PD}}}}&{\Delta M_z^{{\rm{PD}}}}&{\Delta T_x^{{\rm{PD}}}} \end{array}} \right\}^{{\rm{T}}}} $$ (34) 参考式(21), 在PD理论框架下, 同一增量步下的虚应变能密度的增量可以表示为
$$\left. \begin{split} &\Delta {W}_{{\rm{PD}}}^{L}\left({{\boldsymbol{x}}}_{i}\right) = {\displaystyle {\int }_{\hslash \left({{\boldsymbol{x}}}_{i}\right)}\dfrac{1}{2}{\displaystyle {\int }_{{x}_{i}}^{{x}_{j}}\delta }\left(\Delta {{\boldsymbol{\varepsilon}} }_{L}^{\text{T}}\right){{\boldsymbol{D}}}_{{\rm{PD}}}\left(\Delta {{\boldsymbol{\varepsilon}} }_{L}\right)\text{d}x}\text{d}{V}_{{{\boldsymbol{x}}}_{j}}\\ &\Delta {W}_{{\rm{PD}}}^{R}\left({{\boldsymbol{x}}}_{i}\right) = {\displaystyle {\int }_{\hslash \left({{\boldsymbol{x}}}_{i}\right)}\dfrac{1}{2}{\displaystyle {\int }_{{x}_{i}}^{{x}_{j}}\delta }\left(\Delta {{\boldsymbol{\varepsilon}} }_{L}^{\text{T}}\right){{\boldsymbol{S}}}_{{\rm{PD}}} \text{d}x}\text{d}{V}_{{{\boldsymbol{x}}}_{j}}\\ &\Delta {W}_{{\rm{PD}}}^{S}\left({{\boldsymbol{x}}}_{i}\right) = {\displaystyle {\int }_{\hslash \left({{\boldsymbol{x}}}_{i}\right)}\dfrac{1}{2}{\displaystyle {\int }_{{x}_{i}}^{{x}_{j}}\delta }\left(\Delta {{\boldsymbol{\varepsilon}} }_{N}^{\text{T}}\right){{\boldsymbol{S}}}_{{\rm{PD}}} \text{d}x}\text{d}{V}_{{{\boldsymbol{x}}}_{j}}\\ &\Delta {W}_{{\rm{PD}}}^{N}\left({{\boldsymbol{x}}}_{i}\right) = \displaystyle {\int }_{\hslash \left({{\boldsymbol{x}}}_{i}\right)}\dfrac{1}{2}\displaystyle {\int }_{{x}_{i}}^{{x}_{j}}\left\{\delta \left(\Delta {{\boldsymbol{\varepsilon}} }_{L}^{\text{T}}\right){{\boldsymbol{D}}}_{{\rm{PD}}}\left(\Delta {{\boldsymbol{\varepsilon}} }_{N}\right) + \right. \\ &\qquad \left.\delta \left(\Delta {{\boldsymbol{\varepsilon}} }_{N}^{\text{T}}\right){{\boldsymbol{D}}}_{{\rm{PD}}}\left(\Delta {{\boldsymbol{\varepsilon}}}_{L}\right) +\delta \left(\Delta {{\boldsymbol{\varepsilon}}}_{N}^{\text{T}}\right){{\boldsymbol{D}}}_{{\rm{PD}}}\left(\Delta {{\boldsymbol{\varepsilon}} }_{N}\right)\right\}\text{d}x\text{d}{V}_{{{\boldsymbol{x}}}_{j}}\end{split} \right\}$$ (35) 根据文献[30], 存在微梁键与非线性板单元的广义应变增量转换矩阵
${{\boldsymbol{T}}_\xi }$ 使得$$ \Delta {\boldsymbol{\varepsilon }} = {{\boldsymbol{T}}_\xi }\Delta {\boldsymbol{E}} $$ (36) 其中
$$\left.\begin{split} & {{\boldsymbol{T}}_\xi } = \left[ {\begin{array}{*{20}{c}} {{{\boldsymbol{T}}_{p1}}}&{\boldsymbol{0}} \\ {\boldsymbol{0}}&{{{\boldsymbol{T}}_{b1}}} \\ {{{\boldsymbol{T}}_{p2}}}&{\boldsymbol{0}} \\ {\boldsymbol{0}}&{{{\boldsymbol{T}}_{b2}}} \end{array}} \right] \\ & {{\boldsymbol{T}}_{p1}} = \left[ {\begin{array}{*{20}{c}} {{{\cos }^2}\theta }&{{{\sin }^2}\theta }&{\sin \theta \cos \theta } \end{array}} \right] \\ & {{\boldsymbol{T}}_{b1}} = \left[ {\begin{array}{*{20}{c}} {{{\cos }^2}\theta }&{{{\sin }^2}\theta }&{\sin \theta \cos \theta } \end{array}} \right] \\ & {{\boldsymbol{T}}_{p2}} =\dfrac{{6(2\chi - 1)}}{\xi } \left[ \sin \theta \cos \theta \;\;\;-\sin \theta \cos \theta\;\;\; -\left( {{\cos }^{2}}\theta- \frac{1}{2} \right) \right] \\ & {{\boldsymbol{T}}_{b2}} = \left[ {\begin{array}{*{20}{c}} { - \sin \theta \cos \theta }&{\sin \theta \cos \theta }&{{{\cos }^2}\theta - \dfrac{1}{2}} \end{array}} \right] \end{split}\right\} $$ (37) 其中
$\chi $ 为形状参数,$\theta $ 表示微梁键坐标系与单元坐标系的夹角.根据式(35)中PD理论下的虚应变能密度增量与式(21)中CCM理论下的虚应变能密度增量相等, 可求得PD微模量参数为
$$ \left.\begin{split} & {c_{ax}} = \frac{{6E}}{{\text{π} h{\delta ^3}\left( {1 - v} \right)}},\,\,{c_{by}} = \frac{{Eh}}{{2\text{π} {\delta ^3}(1 - v)}} \\ & {c_{bz}} = \frac{{E\left( {1 - 3v} \right)}}{{6\text{π} h\delta \left( {1 - {v^2}} \right)}},\,\,{c_{tor}} = \frac{{Eh\left( {1 - 3v} \right)}}{{2\text{π} {\delta ^3}\left( {1 - {v^2}} \right)}} \end{split} \right\}$$ (38) 2.3 断裂准则
本文采用基于应变能思想的断裂准则来模拟板的损伤和断裂[31]. 在这种条件下, 其临界能量释放率取决于微梁键的能量失效准则, 如图3所示, 断裂能量可以通过断开单位断裂表面上的所有键释放的能量来确定. 本文所提出的非线性PD板模型属于二维问题, 假设每个增量步中每根键的临界微势能为
$\Delta {\omega _c}$ , 表示为$$ \Delta {\omega _c} = \frac{1}{2}\int_{{x_i}}^{{x_j}} {\Delta {{\boldsymbol{\varepsilon }}^{\text{T}}}} {{\boldsymbol{D}}_{PD}}\Delta {\boldsymbol{\varepsilon }}{\text{d}}x + \int_{{x_i}}^{{x_j}} {\Delta {{\boldsymbol{\varepsilon }}^{\text{T}}}} {{\boldsymbol{S}}_{PD}}{\text{d}}x $$ (39) 在二维条件下断裂能
${G_c}$ 的积分方程可定义为$$ {G_c} = 2h\sum\limits_{k = 1}^n {\int_0^\delta {\int_z^\delta {\int_0^{{{\cos }^{ - 1}}(z/\xi) } {\Delta {\omega _c}} } } \xi {\text{d}}\varphi {\text{d}}\xi {\text{d}}z} $$ (40) 假设在点
${{\boldsymbol{x}}_i}$ 的近场域内产生均匀变形, 即近场域内板中性层上的所有广义应变都相等, 则式(39)可以重写为$$ \Delta {\omega _c} = \frac{1}{2}\Delta {{\boldsymbol{\varepsilon }}^{\text{T}}}{{\boldsymbol{D}}_{{\rm{PD}}}}\Delta {\boldsymbol{\varepsilon }}\xi + \Delta {{\boldsymbol{\varepsilon }}^{\text{T}}}{{\boldsymbol{S}}_{{\rm{PD}}}}\xi $$ (41) 将式(41)代入断裂能积分方程(40)可以得到
$$ {G_c} = \frac{{h{\delta ^4}}}{2}\sum\limits_{k = 1}^n {\left( {\frac{1}{2}\Delta {{\boldsymbol{\varepsilon }}^{\text{T}}}{{\boldsymbol{D}}_{{\rm{PD}}}}\Delta {\boldsymbol{\varepsilon }} + \Delta {{\boldsymbol{\varepsilon }}^{\text{T}}}{{\boldsymbol{S}}_{{\rm{PD}}}}} \right)} $$ (42) 最后, 微梁键的能量断裂准则可以表示为
$$ {\omega _{{\rm{bond}}}} = \sum\limits_{k = 1}^n {\left( {\frac{1}{2}\Delta {{\boldsymbol{\varepsilon }}^{\text{T}}}{{\boldsymbol{D}}_{{\rm{PD}}}}\Delta {\boldsymbol{\varepsilon }} + \Delta {{\boldsymbol{\varepsilon }}^{\text{T}}}{{\boldsymbol{S}}_{{\rm{PD}}}}} \right)} = \frac{{2{G_c}}}{{h{\delta ^4}}} $$ (43) 3. 几何非线性PD-CCM板耦合模型
Lubineau等[26]提出了一种结合PD和CCM优点的耦合方法. 如图4所示, 一个完整的模型区域
$\varOmega $ 由耦合方法内的3个子区域组成: 纯CCM模型区域${\varOmega _1}$ 、纯PD模型区域${\varOmega _2}$ 和耦合区域${\varOmega _m}$ .在PD和CCM耦合区域
${\varOmega _m}$ 中, 利用耦合函数$\alpha ({\boldsymbol{x}})$ [26]定义了CCM模型和PD模型之间的关系, 板中性层上一点${{\boldsymbol{x}}_i}$ 在增量步下的虚应变能密度增量可以表示为$$\left. \begin{split} & \Delta W_{Hb}^L({{\boldsymbol{x}}_i}) = \frac{1}{h}\delta \left( {\Delta {\boldsymbol{E}}_L^{\text{T}}} \right){{\boldsymbol{D}}^{Hb}}\left( {\Delta {{\boldsymbol{E}}_L}} \right) + \\ &\qquad \int_{\hbar \left( {{{\boldsymbol{x}}_i}} \right)} {\frac{1}{2}\int_{{x_i}}^{{x_j}} {\alpha ({\boldsymbol{x}})\delta } \left( {\Delta {\boldsymbol{\varepsilon }}_L^{\text{T}}} \right){{\boldsymbol{D}}_{PD}}\left( {\Delta {{\boldsymbol{\varepsilon }}_L}} \right){\text{d}}x} {\text{d}}{V_{{{\boldsymbol{x}}_j}}} \\ &\Delta W_{Hb}^R({{\boldsymbol{x}}_i}) = \frac{1}{h}\delta \left( {\Delta {\boldsymbol{E}}_L^{\text{T}}} \right)\sum\limits_{k = 0}^n {{{\boldsymbol{D}}^{Hb}}} {(\Delta {\boldsymbol{E}})_k} + \\ &\qquad\int_{\hbar \left( {{{\boldsymbol{x}}_i}} \right)} {\frac{1}{2}\int_{{x_i}}^{{x_j}} {\alpha ({\boldsymbol{x}})\delta } \left( {\Delta {\boldsymbol{\varepsilon }}_L^{\text{T}}} \right){{\boldsymbol{D}}_{PD}}\sum\limits_{k = 0}^n {{{(\Delta {\boldsymbol{\varepsilon }})}_k}} {\text{d}}x} {\text{d}}{V_{{{\boldsymbol{x}}_j}}} \\ & \Delta W_{Hb}^S({{\boldsymbol{x}}_i}) = \frac{1}{h}\delta \left( {\Delta {\boldsymbol{E}}_N^{\text{T}}} \right)\sum\limits_{k = 0}^n {{{\boldsymbol{D}}^{Hb}}} {(\Delta {\boldsymbol{E}})_k} + \\ &\qquad\int_{\hbar \left( {{{\boldsymbol{x}}_i}} \right)} {\frac{1}{2}\int_{{x_i}}^{{x_j}} {\alpha ({\boldsymbol{x}})\delta } \left( {\Delta {\boldsymbol{\varepsilon }}_N^{\text{T}}} \right){{\boldsymbol{D}}_{PD}}\sum\limits_{k = 0}^n {{{(\Delta {\boldsymbol{\varepsilon }})}_k}} {\text{d}}x} {\text{d}}{V_{{{\boldsymbol{x}}_j}}} \\ & \Delta W_{Hb}^N({{\boldsymbol{x}}_i}) = \frac{1}{h}\Bigr[ \delta \left( {\Delta {\boldsymbol{E}}_L^{\text{T}}} \right){{\boldsymbol{D}}^{Hb}}\left( {\Delta {{\boldsymbol{E}}_N}} \right) + \delta \left( {\Delta {\boldsymbol{E}}_N^{\text{T}}} \right){{\boldsymbol{D}}^{Hb}}\left( {\Delta {{\boldsymbol{E}}_L}} \right) + \\ &\qquad \delta \left( {\Delta {\boldsymbol{E}}_N^{\text{T}}} \right){{\boldsymbol{D}}^{Hb}}\left( {\Delta {{\boldsymbol{E}}_N}} \right) \Bigr] + \int_{\hbar \left( {{{\boldsymbol{x}}_i}} \right)} \frac{1}{2}\int_{{x_i}}^{{x_j}} \alpha ({\boldsymbol{x}})\Bigr[ \delta \left( {\Delta {\boldsymbol{\varepsilon }}_L^{\text{T}}} \right){{\boldsymbol{D}}_{PD}} \\ &\qquad \left( {\Delta {{\boldsymbol{\varepsilon }}_N}} \right) + {\delta \left( {\Delta {\boldsymbol{\varepsilon }}_N^{\text{T}}} \right){{\boldsymbol{D}}_{PD}}\left( {\Delta {{\boldsymbol{\varepsilon }}_L}} \right) + } \\ &\qquad {\delta \left( {\Delta {\boldsymbol{\varepsilon }}_N^{\text{T}}} \right){{\boldsymbol{D}}_{PD}}\left( {\Delta {{\boldsymbol{\varepsilon }}_N}} \right)} \Bigr] {\text{d}}x {\text{d}}{V_{{{\boldsymbol{x}}_j}}} \end{split}\right\} $$ (44) 由均质化假设可知耦合区域下的虚应变能密度增量与CCM理论下的虚应变能密度增量也相同, 即式(44)与式(21)相等, 由此可得
$$ {{\boldsymbol{D}}^{Hb}} = {\boldsymbol{D}} - \frac{h}{2}\int_{\hbar \left( {{{\boldsymbol{x}}_i}} \right)} {\int_{{x_i}}^{{x_j}} {\alpha ({\boldsymbol{x}}){{\boldsymbol{T}}_\xi }^{\text{T}}{{\boldsymbol{D}}_{PD}}{{\boldsymbol{T}}_\xi }} {\text{d}}x} {\text{d}}{V_{{{\boldsymbol{x}}_j}}} $$ (45) 考虑耦合函数
$\alpha ({\boldsymbol{x}})$ 是线性耦合的标量函数[26], 式(45)可被简化为$$ {{\boldsymbol{D}}^{Hb}} = {\boldsymbol{D}} - \frac{h}{2}\int_{\hbar \left( {{{\boldsymbol{x}}_i}} \right)} {\xi \frac{{\alpha ({{\boldsymbol{x}}_i}) + \alpha ({{\boldsymbol{x}}_j})}}{2}{{\boldsymbol{T}}_\xi }^{\text{T}}{{\boldsymbol{D}}_{PD}}{{\boldsymbol{T}}_\xi }} {\text{d}}{V_{{{\boldsymbol{x}}_j}}} $$ (46) 根据不同区域的材料参数, 建立不同区域的刚度矩阵和结构反力矢量, 进而得到整个结构的刚度矩阵和结构反力矢量, 最后可以建立各增量步的全局平衡方程.
4. 数值算例
4.1 预置裂纹钢板的平面外拉伸撕裂
为了验证所提出的几何非线性PD模型模拟裂纹扩展的能力, 在本节中, 研究了具有预置裂纹的正方形低碳钢板在横向载荷作用下的渐进损伤过程. 这个问题的实验过程可以参考文献[32]. 如图5所示, 钢板的长度和宽度均为
$L = 0.203\;{{\rm{m}}}$ , 厚度为$h = 0.8\;{{\rm{mm}}}$ . 钢板在$X = L/2$ 处有预置初始裂纹, 裂纹长度为${a_0} = 3\;{{\rm{cm}}}$ , 低碳钢板的弹性模量为$E = 210\;{{\rm{GPa}}}$ , 泊松比为$v = 0.3$ , 密度为$7.85 \;{\rm{t}}/{{{{\rm{m}}}^{3}}}$ , 临界能量释放率为$Gc = 255\;{{\rm{kN}}/{\rm{m}}}$ [33]. 钢板左右两端固定约束, 在位于裂纹处两端施加增量载荷$F$ .针对裂纹扩展区域, 采用网格大小为
$\Delta x = 1\;{{\rm{mm}}}$ 的网格离散, 近场域尺寸大小设置为$\delta = 3\Delta x$ , 采用显式中心差分法计算载荷逐步增加的裂纹扩展过程, 时间步长为$\Delta t = 0.1 \,\text{μ}{{\rm{s}}}$ , 这足以满足数值积分的稳定性条件. 经过50000个增量步, 载荷逐步增加到$1800\;{{\rm{N}}}$ , 得到了裂纹扩展过程的准静态结果.如图6所示, 展示了钢板上的损伤演变过程. 从图中可以看出, 当垂直载荷持续增加时, 裂纹路径沿着正
$Y$ 方向扩展. 当施加的载荷逐步增加达到最大时, 裂纹扩展到上边缘的位置, 此外, 由于较大的变形, 观察到较大的裂纹开口. 同时图7给出了载荷-位移曲线对比验证结果, 可以看出该数值结果与文献[32]中的实验结果取得了很好的一致性.![]()
图 7 钢板撕裂的载荷-位移对比曲线Figure 7. Load-displacement comparison curve of steel plate tearing4.2 预置裂纹纸张的平面外拉伸撕裂
在本节中, 通过非线性PD模型计算了纸张的平面外拉伸撕裂过程. 如图8所示, 纸张的长度为
$L = 0.12\;{{\rm{m}}}$ , 宽度为$W = 6 \;{{\rm{cm}}}$ , 厚度为$h = 0.5\;{{\rm{mm}}}$ . 纸张的3条边是固定的, 一条边是自由的, 在自由边两侧$W/3$ 处有对称的初始裂纹, 初始裂纹长度为${a_0} = 3 \;{{\rm{cm}}}$ . 纸张材料的弹性模量为$E = 5.69 \;{{\rm{GPa}}}$ , 泊松比为$v = 0.2$ , 密度为$1.2 \;{{\rm{t}}/}{{{\rm{m}}}^{3}}$ , 临界能量释放率为$Gc = 8.8 \;{{\rm{kN}}/{\rm{m}}}$ [34]. 在两条裂纹之间施加竖直向上的端部剪力增量载荷F.针对裂纹扩展区域, 采用网格大小为
$\Delta x = 0.5\;{{\rm{mm}}}$ 的网格离散, 近场域尺寸大小设置为$\delta = 3\Delta x$ , 采用显式中心差分法计算载荷逐步增加的裂纹扩展过程, 时间步长为$\Delta t = 0.1\;\text{μ}{{\rm{s}}}$ , 经过30000个增量步, 载荷逐步增加到20${{\rm{N}}}$ , 得到了裂纹扩展过程的准静态结果.如图9所示, 展示了纸张上的损伤演变过程, 随着施加的载荷不断增加, 裂纹之间的部分呈现非线性变化弯曲, 并且在达到一定程度时候裂纹开始扩展, 两侧裂纹对称分布呈现对角扩展, 在达到最大载荷时两条裂纹相交, 纸张被剥离成为两个部分. 根据观察, 本研究所展示的非线性PD模型数值计算的损伤结果与Silling等[16]获得的数值结果高度相似.
5. 结 论
本文提出了一种新的PD-CCM薄板弯曲耦合模型, 通过引入几何非线性微梁键, 使得该模型可以考虑几何非线性板弯曲变形的影响. 采用更新的拉格朗日法, 得到了每个增量步PD模型的虚应变能密度增量. 基于均质化和虚功原理, 通过建立每个增量步下虚应变能密度增量的平衡, 获得了几何非线性薄板弯曲模型的PD材料参数. 最后, 将耦合方法应用于几何非线性薄板弯曲PD模型和CCM模型. 通过几个算例验证了该模型的有效性和准确性, 并得出以下结论:
(1)相对于全拉格朗日法, 本文采用更新拉格朗日法表示应变的几何非线性部分更为简捷, 并且使用共旋理论可以更容易地计算应变和应力;
(2)相对于传统键基PD模型, 本文所采用的几何非线性微梁键考虑了键的弯曲和扭转微势能, 模型中的材料参数不受限制;
(3)PD模型, 特别是用于分析几何非线性问题的PD模型, 与CCM模型相比, 具有相当大的计算量, 并且工程师无法方便地将边界条件应用于PD模型, 而本文采用耦合的策略提供了一种有效的方法来节省计算成本并方便地满足边界条件.
-
![]()
图 7 钢板撕裂的载荷-位移对比曲线
Figure 7. Load-displacement comparison curve of steel plate tearing
-
[1] 高云凯, 杨欣, 金哲峰. 轿车车身刚度优化方法研究. 同济大学学报 (自然科学版), 2005, 33(8): 1095-1097 Gao Yunkai, Yang Xin, Jin Zhefeng. Study on method for optimizing car body stiffness. Journal of Tongji University (Natural Science), 2005, 33(8): 1095-1097 (in Chinese)
[2] Kim CS, Shin JG, Kim EK, et al. A study on classification algorithm of rectangle curved hull plates for plate fabrication. Journal of Ship Production and Design, 2016, 32(3): 166-173 doi: 10.5957/jspd.2016.32.3.166
[3] Sussman T, Bathe KJ. 3D-shell elements for structures in large strains. Computers & Structures, 2013, 122: 2-12
[4] Silling SA. Reformulation of elasticity theory for discontinuities and long-range forces. Journal of the Mechanics and Physics of Solids, 2000, 48(1): 175-209 doi: 10.1016/S0022-5096(99)00029-0
[5] Agwai A, Guven I, Madenci E. Predicting crack propagation with peridynamics: a comparative study. International Journal of Fracture, 2011, 171(1): 65-78 doi: 10.1007/s10704-011-9628-4
[6] Kilic B, Madenci E. Prediction of crack paths in a quenched glass plate by using peridynamic theory. International Journal of Fracture, 2009, 156(2): 165-177 doi: 10.1007/s10704-009-9355-2
[7] Nguyen CT, Oterkus S. Ordinary state-based peridynamic model for geometrically nonlinear analysis. Engineering Fracture Mechanics, 2020, 224: 106750 doi: 10.1016/j.engfracmech.2019.106750
[8] Madenci E, Oterkus S. Ordinary state-based peridynamics for plastic deformation according to von Mises yield criteria with isotropic hardening. Journal of the Mechanics and Physics of Solids, 2016, 86: 192-219 doi: 10.1016/j.jmps.2015.09.016
[9] Hu YL, Madenci E. Bond-based peridynamic modeling of composite laminates with arbitrary fiber orientation and stacking sequence. Composite Structures, 2016, 153: 139-175 doi: 10.1016/j.compstruct.2016.05.063
[10] Madenci E, Yaghoobi A, Barut A, et al. Peridynamic modeling of compression after impact damage in composite laminates. Journal of Peridynamics and Nonlocal Modeling, 2021, 3(4): 327-347 doi: 10.1007/s42102-021-00054-1
[11] Zhu N, De Meo D, Oterkus E. Modelling of granular fracture in polycrystalline materials using ordinary state-based peridynamics. Materials, 2016, 9(12): 977 doi: 10.3390/ma9120977
[12] Zhang G, Le Q, Loghin A, et al. Validation of a peridynamic model for fatigue cracking. Engineering Fracture Mechanics, 2016, 162: 76-94 doi: 10.1016/j.engfracmech.2016.05.008
[13] Jung J, Seok J. Mixed-mode fatigue crack growth analysis using peridynamic approach. International Journal of Fatigue, 2017, 103: 591-603 doi: 10.1016/j.ijfatigue.2017.06.008
[14] Oterkus S, Madenci E, Agwai A. Peridynamic thermal diffusion. Journal of Computational Physics, 2014, 265: 71-96 doi: 10.1016/j.jcp.2014.01.027
[15] Nguyen CT, Oterkus S. Peridynamics for the thermomechanical behavior of shell structures. Engineering Fracture Mechanics, 2019, 219: 106623 doi: 10.1016/j.engfracmech.2019.106623
[16] Silling SA, Bobaru F. Peridynamic modeling of membranes and fibers. International Journal of Non-Linear Mechanics, 2005, 40(2-3): 395-409 doi: 10.1016/j.ijnonlinmec.2004.08.004
[17] Taylor M, Steigmann DJ. A two-dimensional peridynamic model for thin plates. Mathematics and Mechanics of Solids, 2015, 20(8): 998-1010 doi: 10.1177/1081286513512925
[18] O’Grady J, Foster J. Peridynamic plates and flat shells: A non-ordinary, state-based model. International Journal of Solids and Structures, 2014, 51(25-26): 4572-4579 doi: 10.1016/j.ijsolstr.2014.09.003
[19] Yang Z, Vazic B, Diyaroglu C, et al. A Kirchhoff plate formulation in a state-based peridynamic framework. Mathematics and Mechanics of Solids, 2020, 25(3): 727-738 doi: 10.1177/1081286519887523
[20] Yang Z, Oterkus E, Oterkus S. A state-based peridynamic formulation for functionally graded Kirchhoff plates. Mathematics and Mechanics of Solids, 2021, 26(4): 530-551 doi: 10.1177/1081286520963383
[21] Shen G, Xia Y, Li W, et al. Modeling of peridynamic beams and shells with transverse shear effect via interpolation method. Computer Methods in Applied Mechanics and Engineering, 2021, 378: 113716 doi: 10.1016/j.cma.2021.113716
[22] Shen G, Xia Y, Hu P, et al. Construction of peridynamic beam and shell models on the basis of the micro-beam bond obtained via interpolation method. European Journal of Mechanics-A/Solids, 2021, 86: 104174 doi: 10.1016/j.euromechsol.2020.104174
[23] Diyaroglu C, Oterkus E, Oterkus S, et al. Peridynamics for bending of beams and plates with transverse shear deformation. International Journal of Solids and Structures, 2015, 69: 152-168
[24] Hu YL, Yu Y, Madenci E. Peridynamic modeling of composite laminates with material coupling and transverse shear deformation. Composite Structures, 2020, 253: 112760 doi: 10.1016/j.compstruct.2020.112760
[25] Nguyen CT, Oterkus S. Ordinary state-based peridynamics for geometrically nonlinear analysis of plates. Theoretical and Applied Fracture Mechanics, 2021, 112: 102877 doi: 10.1016/j.tafmec.2020.102877
[26] Lubineau G, Azdoud Y, Han F, et al. A morphing strategy to couple non-local to local continuum mechanics. Journal of the Mechanics and Physics of Solids, 2012, 60(6): 1088-1102 doi: 10.1016/j.jmps.2012.02.009
[27] Han F, Lubineau G, Azdoud Y, et al. A morphing approach to couple state-based peridynamics with classical continuum mechanics. Computer Methods in Applied Mechanics and Engineering, 2016, 301: 336-358 doi: 10.1016/j.cma.2015.12.024
[28] Zheng G, Li L, Han F, et al. Coupled peridynamic model for geometrically nonlinear deformation and fracture analysis of slender beam structures. International Journal for Numerical Methods in Engineering, 2022, 123(16): 3658-3680 doi: 10.1002/nme.6984
[29] Fung YC. Foundations of Solid Mechanics. Prentice Hall, 1965
[30] Zheng G, Yan Z, Xia Y, et al. Peridynamic shell model based on micro-beam bond. CMES-Computer Modeling in Engineering & Sciences, 2023, 3: 1975-1995
[31] Foster JT, Silling SA, Chen W. An energy based failure criterion for use with peridynamic states. International Journal for Multiscale Computational Engineering, 2011, 9(6): 675-687 doi: 10.1615/IntJMultCompEng.2011002407
[32] Muscat-Fenech CM, Atkins AG. Out-of-plane stretching and tearing fracture in ductile sheet materials. International Journal of Fracture, 1997, 84(4): 297-306 doi: 10.1023/A:1007325719337
[33] Areias P, Rabczuk T, Msekh MA. Phase-field analysis of finite-strain plates and shells including element subdivision. Computer Methods in Applied Mechanics and Engineering, 2016, 312: 322-350 doi: 10.1016/j.cma.2016.01.020
[34] Mai YW, He H, Leung R, et al. Fracture Mechanics: 26th Volume. Philadelphia: ASTM, 1995: 587-599
下载:
计量
- 文章访问数: 810
- HTML全文浏览量: 273
- PDF下载量: 230
