RESEARCH ON DEFORMATION OF HYPERELASTIC MATERIALS BASED ON A NEW VIRTUAL ELEMENT METHOD
-
摘要: 虚单元法是一种先进的求解固体力学问题的数值方法, 在过去10年间, 该算法在线弹性问题中得到了较为广泛地开发和应用. 文章尝试给出一种通用的高阶虚单元法格式, 可用于计算超弹性问题以及更普遍的非线性问题. 与传统求解力学问题的虚单元法的思想不同, 其主要思想是对泊松方程求解映射算子, 并将该映射算子直接用于位移场的近似, 从而可求解众多非线性力学问题. 由于采用了标量场的映射算子来近似矢量场, 因此该算法格式简单, 并且可以轻易扩展到高阶格式或者三维问题求解. 将从泊松方程出发, 介绍虚单元法中椭圆映射算子的计算方法, 在此基础上, 详细推导虚单元法在求解超弹性问题时的具体格式, 并给出切线刚度矩阵的计算方法. 最后, 给出了几个典型的超弹性数值算例, 从而证明该虚单元法格式的有效性.Abstract: The virtual element method is an advanced numerical method for solving solid mechanics problems. In the past ten years, the numerical method has been widely developed and applied in linear elasticity problems. This work attempts to give a general high-order virtual element method format that can be used to calculate hyperelastic problems and more general nonlinear problems. Different from the traditional virtual element method for solving mechanical problems, its main idea is to solve the projection operator for the Poisson equation and use the projection operator directly for the approximation of the displacement field, so that it can solve many nonlinear mechanical problems. Since the projection operator of the scalar field is used to approximate the vector field, the method has a simple format and can be easily extended to high-order formats or three-dimensional problems. This work will start from the Poisson equation and introduce the calculation method of the elliptical projection operator in the virtual element method. On this basis, the specific format of the virtual element method in solving hyperelastic problems will be derived in detail, and the calculation of the tangent stiffness matrix will be given. Finally, this paper gives several typical hyperelastic numerical examples to prove the effectiveness of the virtual element method format.
-
引 言
与传统有限元法相比, 多边形有限元法是一种更灵活的求解技术, 可以更好地处理拥有较复杂计算域的问题. 现有的多边形算法包括但不限于多边形有限元法(polygonal finite element method, PFEM)[1-2]、间断伽辽金法(discontinuous Galerkin methods, DG)[3-4]、比例边界有限元法(scaled boundary finite element method, SBFEM)[5]、扩展有限元法(extended FEMs, XFEM)[6]和广义有限元法(generalised FEMs, GFEM)[7]. 除了以上提出的算法之外, 意大利学者于2013年提出的虚单元法(virtual element method, VEM)[8-10]获得了越来越多人的关注. 虚单元法来源于MFD (mimetic finite difference methods)[11-13], 并可视为一种一般形式的有限元法. 和传统有限元法相比, 该算法所使用的单元可以是任意形状的多边形和多面体, 因此在处理众多问题时, 具有一定的灵活性和优势. 到目前为止, 该算法已经在线弹性力学[14-18]、超弹性问题[19-22]、接触问题[23-26]、动力学问题[27-29]、相场断裂问题[30-31]和特征值问题[32-33]等领域得到了应用. 到目前为止, 虚单元法在国内的研究尚不广泛, 力学相关中文论文可见文献[34-36]. 此外, 无稳定项虚单元法(self-stabilization或者stabilization-free)也得到了学者的广泛研究, 可见文献[37-40].
虚单元法有两个重要特点, 一是可以使用任意多边形(单元的边数大于3的单元为多边形单元)或者多面体单元, 因此放松了有限元对单元形状的限制; 二是在计算单元刚度矩阵时, 需要构造额外的稳定项. 和有限元法不同, 虚单元法无需计算单元内的插值函数$ \phi $, 只要求多边形边界上为连续多项式. 由于插值函数没有显式表达式, 需要将基函数$ \phi $由其在多项式空间中的投影$\prod \phi $代替, 而这里的投影算子需要用到单元边界上的插值函数进行计算. 对于多边形单元, 这将造成单元刚度矩阵的秩不满足要求($ {n_u} - {n_\varepsilon } > 3 $,${n_u}$和${n_\varepsilon }$分别为位移插值函数和应变插值函数的维数[41]), 故出现零能模式. 因此, 需要在刚度矩阵之外, 额外增加稳定项.
在传统的虚单元法中, 需要根据指定的方程, 构造正交条件, 求解指定方程的映射算子. 一般来说, 直接根据非线性问题的控制方程计算映射算子是几乎不可能的. 因此在进行非线性力学问题计算时, 通常首先计算线弹性方程的映射算子, 再进一步对应变进行近似, 从而得到非线性问题的离散形式. 由于线弹性方程是一个矢量方程, 构造高阶格式比较麻烦(主要集中在内部矩量的定义和映射限制条件的施加). 因此现有文献多使用的格式为低阶格式($k = 1$时不需要内部矩量). 除此之外, 非线性问题的稳定项不仅依赖于映射算子, 还依赖于当前状态(当前材料或者应力状态), 无论从推导到计算都比较繁琐[23]. 当前涉及到非线性问题的文献, 其切线刚度矩阵的格式、稳定项的选择和内力项的计算过程都没有很清晰的描述. 本文将以超弹性问题为例, 对以上问题详细展开说明.
和传统弹性力学问题虚单元法中直接计算位移的映射算子不同, 本文利用泊松方程的映射算子, 对位移场的分量分别进行近似, 从而直接构造虚单元法求解非线性力学问题的求解格式. 通常来说, 针对泊松方程, 虚单元法的映射算子计算更简单, 也更容易扩展到高阶格式, 甚至任意阶格式. 利用计算所得的高阶映射算子, 可进一步构造对位移场的映射算子, 从而直接代入到线性化之后的平衡方程中, 计算刚度矩阵. 对于稳定项, 由于方程最终进行了线性化, 因此可直接通过自由度相乘的形式计算稳定刚度矩阵[8,22]. 一般来说, 以上提到的计算格式具有很强的可扩展性. 根据泊松方程的映射算子, 便可直接求解大部分非线性力学问题.
该文分为以下几个部分: 第1章给出超弹性问题的基本理论和选用的超弹性材料本构; 第2章主要解释如何计算虚单元法中的不同映射算子; 第3章介绍如何利用映射算子, 进一步推导出超弹性问题的虚单元法求解格式; 在第4章中, 给出一些典型算例, 并进行分析. 最后给出总结.
1. 基本理论
1.1 超弹性问题的控制方程及其弱形式
在本节中考虑二维超弹性的静力及动力问题. 对于给定二维几何体$\varOmega \in {\mathbb{R}^2}$, $ \varGamma = \partial \varOmega = {\varGamma _D} \cup {\varGamma _t} $为模型的边界, ${\varGamma _D}$为Dirichlet边界, 在其上指定位移; ${\varGamma _t}$为Neumann边界, 在其上指定面力.
如图1所示, 假定${\boldsymbol{X}}$为任意物质点在初始构型下的坐标, $ {\boldsymbol{x}} $为该物质点在当前构型下的坐标, 则其关系可描述为
$$ {\boldsymbol{x}} = {\boldsymbol{\varphi}} \left( {{\boldsymbol{X}},t} \right) ={\boldsymbol{ X}} + {\boldsymbol{u}}\left( {{\boldsymbol{X}},t} \right) $$ (1) 其中$ {\boldsymbol{u}} $为位移. 位移便于在参考坐标系下对物体变形进行描述, 可定义变形梯度张量${\boldsymbol{F}}$为
$$ {\boldsymbol{F}} = {\text{Grad}}{\boldsymbol{\varphi}} = \frac{{{\text{d}}{\boldsymbol{x}}}}{{{\text{d}}{\boldsymbol{X}}}} = {\boldsymbol{I}} + \nabla {\boldsymbol{u}} $$ (2) 其中${\boldsymbol{I}}$为单位矩阵, $ \nabla (\cdot) = \dfrac{\partial }{\partial {X}_{j}}{(\cdot)}_{i} $为矢量的梯度算子. 注意这里的梯度算子参考初始构型. 引入右柯西-格林(right Cauchy-Green)变形张量${\boldsymbol{C}} = {{\boldsymbol{F}}^{\mathrm{T}}} \cdot {\boldsymbol{F}}$, 格林-拉格朗日(Green-Lagrange)应变${\boldsymbol{E}}$有以下形式
$$ {\boldsymbol{E}} = \frac{1}{2}\left( {{\boldsymbol{C}} - {\boldsymbol{I}}} \right) $$ (3) 参考构型下变形体$\varOmega $平衡方程为
$$ {\text{Div}}{\boldsymbol{P}} + {\boldsymbol{f}} = {\boldsymbol{0}} $$ (4) 其中${\boldsymbol{P}}$为第一P-K应力(first Piola-Kirchhoff stress), ${\boldsymbol{f}}$为体积力. 由于应力张量${\boldsymbol{P}}$是非对称的, 在数值计算中往往采用对称的第二P-K (second Piola-Kirchhoff) 应力张量${\boldsymbol{S}},$ ${\boldsymbol{S}}$与${\boldsymbol{P}}$存在如下转换关系
$$ {\boldsymbol{S}} = {{\boldsymbol{F}}^{ - 1}} \cdot {\boldsymbol{P}} $$ (5) 并且, ${\boldsymbol{S}}$和${\boldsymbol{E}}$是功共轭的, 即
$$ {\boldsymbol{S}} = \frac{{\partial \varPsi }}{{\partial {\boldsymbol{E}}}} $$ (6) 其中$\varPsi $为超弹性应变能密度. 超弹性问题的弱形式为
$$ \delta W = \mathop \smallint \nolimits_{{V_0}} \delta {\boldsymbol{E}}:{\boldsymbol{S}}{\text{d}}{V_0} $$ (7) 式中, $\delta {\boldsymbol{E}}$为格林-拉格朗日应变的变分, 并且有
$$ \delta {\boldsymbol{E}} = \frac{1}{2}\left( {\delta {{\boldsymbol{F}}^{\mathrm{T}}} \cdot {\boldsymbol{F}} + {{\boldsymbol{F}}^{\mathrm{T}}} \cdot \delta {\boldsymbol{F}}} \right) $$ (8) 平衡方程式(4)的数值求解可以转化成对式(7)中弱形式的数值求解. 通过伽辽金法可以将式(7)转化成代数方程组. 由于式(7)是非线性的, 在采用牛顿迭代法对其所对应的代数方程组进行数值求解时, 需对其进行线性化, 以得到刚度矩阵. 式(7)线性化后, 有
$$ \Delta \left( {\delta \mathcal{W}} \right) = \int_{{V_0}} {\delta {\boldsymbol{E}}:\Delta {\boldsymbol{S}}\, + \Delta \left( {\delta {\boldsymbol{E}}} \right):{\boldsymbol{S}}} {\text{d}}{V_0} $$ (9) 第二P-K应力${\boldsymbol{S}}$的线性增量可以进一步表示为
$$ \Delta {\boldsymbol{S}} = \frac{{\partial {\boldsymbol{S}}}}{{\partial {\boldsymbol{E}}}}:\Delta {\boldsymbol{E}} = {{\boldsymbol{D}}}:\Delta {\boldsymbol{E}} $$ (10) 其中$ \boldsymbol{D} $为4阶本构张量. 对式(8)进行线性化, 并结合变形梯度F的定义, 可以得到
$$ \begin{split} &\Delta \left( {\delta {\boldsymbol{E}}} \right) = \frac{1}{2}\left( {\delta {{\boldsymbol{F}}^{\mathrm{T}}} \cdot \Delta {\boldsymbol{F}} + \Delta {{\boldsymbol{F}}^{\mathrm{T}}} \cdot \delta {\boldsymbol{F}}} \right)= \\&\qquad \frac{1}{2}\left[ {\delta {{\left( {\nabla {\boldsymbol{u}}} \right)}^{\mathrm{T}}} \cdot \Delta \left( {\nabla {\boldsymbol{u}}} \right) + \Delta {{\left( {\nabla {\boldsymbol{u}}} \right)}^{\mathrm{T}}} \cdot \delta \nabla {\boldsymbol{u}}} \right] \end{split} $$ (11) 将式(10)和式(11)代入式(9)中, 并根据S和E的对称性, 可以得到
$$ \begin{split} &\Delta \delta \mathcal{W} = \displaystyle\int \nolimits_{{V_0}} \delta {\boldsymbol{E}}:{{\boldsymbol{D}}}:\Delta {\boldsymbol{E}}{\text{d}}V+ \\&\qquad \displaystyle\int \nolimits_{{V_0}} {\boldsymbol{S}}:\left[ {\delta {{\left( {\nabla {\boldsymbol{u}}} \right)}^{\mathrm{T}}} \cdot \Delta \left( {\nabla {\boldsymbol{u}}} \right)} \right]{\text{d}}{V_0} \end{split} $$ (12) 材料模型选用可压缩Neo-Hookean超弹性模型, 其应变能$\varPsi $表达式为
$$ \varPsi = \frac{\mu }{2}\left( {{I_C} - 3} \right) - \mu {\text{ln}}J + \frac{\lambda }{2}{\left( {{\text{ln}}J} \right)^2} $$ (13) 其中$ J = \det {\boldsymbol{F}} $, $\lambda $和$\mu $为拉梅常数. 由式(6)可以得到Neo-Hookean超弹性模型的第二P-K应力表达式
$$ {\boldsymbol{S}} = 2\frac{{\partial \varPsi }}{{\partial {\boldsymbol{C}}}} = \mu \left( {{\boldsymbol{I}} - {{\boldsymbol{C}}^{ - 1}}} \right) + \lambda \left( {{\text{ln}}J} \right){{\boldsymbol{C}}^{ - 1}} $$ (14) 对第二P-K应力进行求导, 进而得到对应的本构张量
$$ \left.\begin{split} &{{\boldsymbol{D}}} = 2\frac{{\partial {\boldsymbol{S}}}}{{\partial {\boldsymbol{C}}}} = \lambda {{\boldsymbol{C}}^{ - 1}} \otimes {{\boldsymbol{C}}^{ - 1}} + 2\left( {\mu - \lambda {\text{ln}}J} \right){\boldsymbol{L}} \\& {L_{IJKL}} = - \frac{{\partial {C^{ - 1}}_{IJ}}}{{\partial {C_{KL}}}} \end{split}\right\} $$ (15) 根据以上弱形式及其线性化后的形式, 可以建立虚单元法非线性求解超弹性静力问题的数值格式. 虚单元法相比有限元, 可以使用任意不规则单元, 因此在处理超弹性问题上具有一定优势. 传统虚单元法在求解此类问题时, 需要构造位移场或者应变的投影, 计算映射算子的过程性比较复杂. 本文将讨论如何使用泊松方程的映射算子, 构造超弹性问题的虚单元法格式.
2. 基于泊松方程的投影算子计算
根据现有文献, 在使用虚单元法处理力学问题时, 需要引入矢量多项式空间${\mathcal{P}_k}$. 为了简化计算, 分别将位移场的分量近似, 因此只需要构造标量场的映射算子即可. 本文从泊松方程出发, 计算标量场的映射算子, 并构造位移场的映射, 从而计算多边形单元内的位移梯度, 进而过渡到能够处理超弹性问题.
2.1 模型离散与虚单元空间
虚单元法需要将模型$\varOmega $离散成一系列互不重叠的任意多边形单元$ E \in {\mathcal{J}_h} $. 单元E的边数为n, 半径为${h_E}$ (特征尺寸), 面积为$|E|$. 对于任意给定阶数k, 为了表达自由度, 定义$ {{{\boldsymbol{M}}}_k}\left( E \right) $为单元E上的尺度单项式集
$$ {{{\boldsymbol{M}}}_k}\left( E \right): = \left\{ {{{\left( {\frac{{{\boldsymbol{X}} - {{\boldsymbol{X}}_E}}}{{{h_E}}}} \right)}^s},\left| {\boldsymbol{s}} \right| \leqslant k} \right\} $$ (16) 其维数为$\left( {k + 1} \right)\left( {k + 2} \right)/2$,$|s| = {s_1} + {s_2}$, ${{\boldsymbol{X}}_E}$为单元的E形心, 此外
$$ {{\boldsymbol{X}}^s}: = X_1^{{s_1}}X_2^{{s_2}} $$ (17) 由于虚单元法没有显式的形函数表达式, 因此需要构造多项式空间的局部投影. 对于泊松方程, 借助以下虚单元空间来构造投影算子
$$ \begin{split} & {\mathcal{V}_k}\left( E \right): = \Big\{ {{u_h} \in {\mathcal{H}^1}\left( E \right):\Delta u \in {\mathcal{P}_{k - 2}}\left( E \right){\text{ in }}E,} \\& \qquad {u{|_{\partial E}} = {\mathcal{B}_k}\left( {\partial E} \right)} \Big\} \end{split} $$ (18) 其中$ {\mathcal{P}_k} $为最高阶次不超过k的标量多项式空间, ${\mathcal{B}_k}\left( {\partial E} \right)$为边界元空间
$$ {\mathcal{B}_k}\left( {\partial E} \right): = \left\{ {{u_h} \in C\left( {\partial E} \right):{u_e} \in {\mathcal{P}_k}(e),\;\;e \in \partial E} \right\} $$ (19) 接下来需要根据$ {\mathcal{V}_k}\left( E \right) $设定自由度. 对于多边形的任意条边, 需要k + 1个点的值以确定k阶多项式, 因此需要$ n + \left( {k - 1} \right)n = nk $个自由度. 此外, 由于$\Delta u = {\mathcal{P}_{k - 2}}$, 因此需要在单元内部额外增加${{k(k - 1)} \mathord{\left/ {\vphantom {{k(k - 1)} 2}} \right. } 2}$个自由度. 对于$k = 2$, 自由度的设置如图2所示. 最终, 函数空间$ {\mathcal{V}_k}\left( E \right) $的维数为
$$ {N_E}: = {\text{dim}}{\mathcal{V}_k}\left( E \right) = nk + \frac{{k\left( {k - 1} \right)}}{2} $$ (20) 对于单元角点及边界上的节点, 其自由度可设置为节点值, 对于内部点, 其自由度可定义为
$$ \frac{1}{{\left| E \right|}}\mathop \smallint \nolimits_E {u_h}{m_\alpha }{\text{d}}\varOmega ,\quad \forall {m_\alpha } \in {\mathcal{M}_{k - 2}}\left( E \right) $$ (21) 2.2 椭圆投影算子
对应虚单元空间$ {\mathcal{V}_k}\left( E \right) $, 定义以下椭圆投影算子
$$ \Pi _k^\nabla \left( E \right):{\mathcal{V}_k}\left( E \right) \to {\mathcal{P}_k}\left( E \right),u \mapsto \Pi _k^\nabla u $$ (22) 其中$\Pi _k^\nabla $表示椭圆投影算子. $ E \in {\mathcal{J}_h} $, $\Pi _k^\nabla $可通过以下正交条件计算
$$ \int \nolimits_E \nabla \left( {\Pi _k^\nabla {u_h} - {u_h}} \right) \cdot \nabla p\,{\text{d}}\varOmega = 0 $$ (23) 因此可进一步得到
$$ \int \nolimits_E \nabla \Pi _k^\nabla {u_h} \cdot \nabla p\,{\text{d}}\varOmega = \int \nolimits_E \nabla {u_h} \cdot \nabla p\,{\text{d}}\varOmega $$ (24) 式(23)和式(24)中, 均有$\forall {u_h} \in {\mathcal{H}^1}\left( E \right)$, $\forall p \in {\mathcal{P}_k}\left( E \right)$. 由于式(24)左侧矩阵的第一行全为0, 为了求解映射算子$\Pi _k^\nabla \left( E \right)$, 需要补充以下限制条件
$$ {P_0}:{\mathcal{H}^1}\left( E \right) \to {\mathcal{P}_0}\left( E \right),{\text{ }}{P_0}\left( {\Pi _k^\nabla {u_h} - {u_h}} \right) = 0 $$ (25) 对于不同阶次算法, 限制条件算符${P_0}$有以下不同的选择
$$ {P_0}\left( {{u_h}} \right): = \left\{ {\begin{array}{*{20}{l}} {\dfrac{1}{n} \displaystyle\sum \limits_{i = 1}^n {u_h},\quad {\text{ for }}k = 1} \\ {\dfrac{1}{{\left| E \right|}}\int \nolimits_E {u_h}{\text{d}}\varOmega, \quad {\text{ for }}k \geqslant 2} \end{array}} \right. $$ (26) 在实施过程中, 例如对于式(24)的右侧, 其矩阵第一行需要变换成式(26)的形式(对于一阶VEM, 矩阵第一行的值全为$1/{N_E}$, 对于二阶VEM, 需要用到内部自由度的定义, 见式(21)).
对式(24)的右侧进行分部积分, 并借助高斯公式, 有
$$ \begin{split} &\int \nolimits_E \nabla \Pi _k^\nabla {u_h} \cdot \nabla p\,{\text{d}}\varOmega = \\&\qquad - \int \nolimits_E {u_h} \cdot {\nabla ^2}p\,{\text{d}}\varOmega + \int \nolimits_{\partial E} {u_h} \cdot \frac{{\partial p}}{{\partial {\boldsymbol{n}}}}{\text{d}}\varGamma \\ \end{split} $$ (27) 由于节点自由度已知, 因此可计算椭圆投影基函数. 假设虚单元空间$ {\mathcal{V}_k}\left( E \right) $的基函数(节点基)为${\phi _k}\left( E \right)$, 多项式空间${\mathcal{P}_k}\left( E \right)$的基函数为$ {\boldsymbol{m}} \in {{{\boldsymbol{M}}}_k}\left( E \right) $, 则节点基函数的投影为
$$ {\Pi _k^\nabla {\phi _i}}{ = \displaystyle\sum \limits_{\alpha = 1}^{{N_P}} {a_{\alpha ,i}}{m_\alpha } = \displaystyle\sum \limits_{j = 1}^{{N_E}} {s_{j,i}}{\phi _j}} $$ (28) 其中${N_P}: = {\text{dim}}{\mathcal{P}_k}\left( E \right)$. 上式可写成以下矩阵向量形式
$$ \begin{split} &\left[ {\begin{array}{*{20}{l}} {\Pi _k^\nabla {\phi _1},}&{\Pi _k^\nabla {\phi _2},}& \cdots &{\Pi _k^\nabla {\phi _{{N_K}}}} \end{array}} \right]= \\&\qquad {\boldsymbol{\Pi}} _k^\nabla {{\boldsymbol{\phi}} ^{\mathrm{T}}} = {{\boldsymbol{m}}^{\mathrm{T}}}{\boldsymbol{\Pi}} _{k*}^\nabla = {{\boldsymbol{\phi}} ^{\mathrm{T}}}{\boldsymbol{\Pi }}_k^\nabla \\ \end{split} $$ (29) 其中$ {\boldsymbol{\Pi}} _k^\nabla $为椭圆投影算子在节点基函数下的矩阵表示, $ {\boldsymbol{\Pi}} _{k*}^\nabla $为椭圆投影算子在多项式基下的矩阵表示. 根据两种基函数的转换
$$ {{\boldsymbol{m}}^{\mathrm{T}}} = {{\boldsymbol{\phi}} ^{\mathrm{T}}}{\boldsymbol{D}} $$ (30) 可得到以下关系
$$ {\boldsymbol{\Pi}} _k^\nabla = {\boldsymbol{D\Pi}} _{k{\text{*}}}^\nabla $$ (31) 其中${\boldsymbol{D}}$为过渡矩阵(尺寸为${N_E} \times {N_P}$), 其表达式为
$$ {{\boldsymbol{D}}_{j\alpha }} = {\text{do}}{{\text{f}}_j}\left( {{m_\alpha }} \right),{\text{ }}j = 1,2, \cdots ,{N_E};\alpha = 1,2, \cdots ,{N_P} $$ (32) 为了计算投影算子矩阵, 在式(27)中考虑基函数, 则有
$$ \begin{split} & \int \nolimits_E \nabla {\boldsymbol{m}} \cdot \nabla \Pi _k^\nabla {{{\boldsymbol{\phi}} }^{\mathrm{T}}}\,{\mathrm{d}}\varOmega = \\&\qquad - \int \nolimits_E {\nabla ^2}{\boldsymbol{m}} \cdot {{{\boldsymbol{\phi}} }^{\mathrm{T}}}\,{\mathrm{d}}\varOmega + \int \nolimits_{\partial E} \left( {\nabla {\boldsymbol{m}} \cdot {\boldsymbol{n}}} \right){{{\boldsymbol{\phi}} }^{\mathrm{T}}}{\mathrm{d}}\varGamma \\ \end{split} $$ (33) 其中
$$ \nabla {{\boldsymbol{\phi}} ^{\mathrm{T}}} = \left[ {\begin{array}{*{20}{c}} {\dfrac{{\partial {{\boldsymbol{\phi}} ^{\mathrm{T}}}}}{{\partial X}}} \\ {\dfrac{{\partial {{\boldsymbol{\phi}} ^{\mathrm{T}}}}}{{\partial Y}}} \end{array}} \right],\nabla {{\boldsymbol{m}}^{\mathrm{T}}} = \left[ {\begin{array}{*{20}{c}} {\dfrac{{\partial {{\boldsymbol{m}}^{\mathrm{T}}}}}{{\partial X}}} \\ {\dfrac{{\partial {{\boldsymbol{m}}^{\mathrm{T}}}}}{{\partial Y}}} \end{array}} \right] $$ (34) 将式(29)代入式(33), 可得
$$ \begin{split} &\int \nolimits_E \nabla {\boldsymbol{m}}\nabla {{\boldsymbol{m}}^{\mathrm{T}}}{\mathrm{d}}\varOmega {\boldsymbol{\Pi}} _{k*}^\nabla = \\&\qquad {\text{ }} - \int \nolimits_E {\nabla ^2}{\boldsymbol{m}}{{\boldsymbol{\phi}} ^{\mathrm{T}}}\,{\mathrm{d}}\varOmega + \int \nolimits_{\partial E} \left( {\nabla {\boldsymbol{mn}}} \right){{\boldsymbol{\phi}} ^{\mathrm{T}}}{\mathrm{d}}\varGamma \end{split} $$ (35) 其中${{\boldsymbol{n}}}$为单元边界外法线. 对于一阶VEM, 式(35)右端第一项为0; 对于二阶VEM, 需要通过内部自由度计算. 因此, 式(35)可进一步写为
$$ {\boldsymbol{G\Pi}} _{k*}^\nabla = {\boldsymbol{B}} $$ (36) 其中$ {\boldsymbol{G}} $是尺寸为 ${N_P} \times {N_p}$ 的矩阵, ${\boldsymbol{B}}$是尺寸为${N_E} \times {N_p}$的矩阵. 根据自由度的假设, 式(35)的右端是可计算的, 因此矩阵${\boldsymbol{B}}$是可计算的. 值得一提的是, 由于没有考虑限制条件, 矩阵${\boldsymbol{G}}$的第一行元素都为0. 将矩阵元素为0的方程用限制条件(式(26))代替, 可得
$$ \tilde {\boldsymbol{G}}{\boldsymbol{\Pi}} _{k*}^\nabla = \tilde {\boldsymbol{B}} $$ (37) 则椭圆投影算子矩阵可通过下式计算
$$ {\boldsymbol{\Pi}} _{k*}^\nabla = {\tilde {\boldsymbol{G}}^{ - 1}}\tilde {\boldsymbol{B}} $$ (38) 针对泊松方程问题, 我们得到了如式(38)及式(31)所示的投影算子矩阵的计算方法. 通过该投影算子, 对于任意标量$u$的梯度$\nabla u$可通过如下近似
$$ \begin{gathered} \nabla \left( {\Pi _k^\nabla u} \right) = \nabla \left( {\Pi _k^\nabla {{\boldsymbol{\phi}} ^{\mathrm{T}}}} \right){\boldsymbol{u}} = \nabla {{\boldsymbol{\phi}} ^{\mathrm{T}}}{\boldsymbol{\Pi}} _k^\nabla {\boldsymbol{u}} = \nabla {{\boldsymbol{m}}^{\mathrm{T}}}{\boldsymbol{\Pi }}_{k*}^\nabla {\boldsymbol{u}} \end{gathered} $$ (39) 虚单元法的基本思想是将变量${\boldsymbol{u}}$分为映射分量$\Pi _k^\nabla u$和余量, 即
$$ {u_h} = \Pi _k^\nabla {u_h} + \left( {{u_h} - \Pi _k^\nabla {u_h}} \right) $$ (40) 因此, 在构造双线性格式的时候, 会得到额外的稳定项.
对于弹性力学问题, 和传统的VEM求解格式不同, 本文不需要针对矢量场构造投影算子, 只需要根据式(38)及式(31)计算所得投影算子进行变换, 从而计算出多边形单元内的位移梯度、应变等, 从而可适用于求解非线性力学问题. 接下来以超弹性问题为例, 利用式(38)及式(31)计算所得投影算子, 详细推导VEM求解超弹性问题的计算过程.
3. 超弹性问题的虚单元法求解格式
3.1 超弹性问题的切线刚度矩阵
对于超弹性材料的静力学问题, 线性化之后的内力虚功为
$$ \begin{split} &\Delta \delta \mathcal{W} = \int \nolimits_{{V_0}} \delta {\boldsymbol{E}}:{{\boldsymbol{D}}}:\Delta {\boldsymbol{E}}{\text{d}}V+ \\& \qquad \int \nolimits_{{V_0}} {\boldsymbol{S}}:\left[ {\delta {{\left( {\nabla {\boldsymbol{u}}} \right)}^{\mathrm{T}}} \cdot \Delta \left( {\nabla {\boldsymbol{u}}} \right)} \right]{\text{d}}{V_0} \end{split} $$ (41) 对于二维问题, 为了便于数值计算, 将格林-拉格朗日应变张量${\boldsymbol{E}}$和第二P-K应力分别表示成如下向量形式
$$ \hat {\boldsymbol{E}} = {\left[ {\begin{array}{*{20}{c}} {{E_{XX}}}&{{E_{YY}}}&{2{E_{XY}}} \end{array}} \right]^{\mathrm{T}}},\;\;\hat {\boldsymbol{S}} = {\left[ {\begin{array}{*{20}{c}} {{S_{XX}}}&{{S_{YY}}}&{{S_{XY}}} \end{array}} \right]^{\mathrm{T}}} $$ (42) 式中$\hat {\boldsymbol{S}}$和$\hat {\boldsymbol{E}}$仍然是功共轭的. 此外, 对于式(41)中等号右端第二项, 有
$$ {\boldsymbol{S}}:\left[ {\delta {{\left( {\nabla {\boldsymbol{u}}} \right)}^{\mathrm{T}}} \cdot \Delta \left( {\nabla {\boldsymbol{u}}} \right)} \right] = \delta {{\boldsymbol{\theta}} ^{\mathrm{T}}} \cdot {{\boldsymbol{I}}} \cdot \Delta {\boldsymbol{\theta}} $$ (43) 其中
$$ {\boldsymbol{\theta}} = {\left[ {\begin{array}{*{20}{c}} {\dfrac{{\partial u}}{{\partial X}}}&{\dfrac{{\partial u}}{{\partial Y}}}&{\dfrac{{\partial v}}{{\partial X}}}&{\dfrac{{\partial v}}{{\partial Y}}} \end{array}} \right]^{\mathrm{T}}},\;\;{{\boldsymbol{I}}} = \left[ {\begin{array}{*{20}{c}} {\boldsymbol{S}}&{\boldsymbol{0}} \\ {\boldsymbol{0}}&{\boldsymbol{S}} \end{array}} \right] $$ (44) 考虑
$$ \delta \hat {\boldsymbol{E}} = {\boldsymbol{A}} \cdot \delta {\boldsymbol{\theta}} ,\quad \Delta \hat {\boldsymbol{E}} = {\boldsymbol{A}} \cdot \Delta {\boldsymbol{\theta}} $$ (45) 其中
$$ {\boldsymbol{A}} = \left[ {\begin{array}{*{20}{c}} {{F_{11}}}&0&{{F_{21}}}&0 \\ 0&{{F_{12}}}&0&{{F_{22}}} \\ {{F_{12}}}&{{F_{11}}}&{{F_{22}}}&{{F_{21}}} \end{array}} \right] $$ (46) 则可得到以下关系
$$ \begin{aligned} \delta \boldsymbol{E}^{\mathrm{T}} : {{\boldsymbol{D}}}: \Delta \boldsymbol{E} = \delta \hat{\boldsymbol{E}}^{\mathrm{T}} \cdot \hat{{{\boldsymbol{D}}}} \cdot \Delta \hat{\boldsymbol{E}} = \delta \boldsymbol{\theta}^{\mathrm{T}} \cdot \boldsymbol{A}^{\mathrm{T}} \cdot \hat{{{\boldsymbol{D}}}} \cdot \boldsymbol{A} \cdot \Delta \boldsymbol{\theta}\end{aligned} $$ (47) 将式(45)和式(47)代入式(41), 最终得到
$$ \begin{split}& \Delta \delta \mathcal{W} = \int_{V_0} \delta \hat{\boldsymbol{E}}^{\mathrm{T}} \cdot \hat{{{\boldsymbol{D}}}} \cdot \Delta \hat{\boldsymbol{E}} \mathrm{d} V_0 + \int_{V_0} \delta \boldsymbol{\theta}^{\mathrm{T}} \cdot {{\boldsymbol{I}}} \cdot \Delta \boldsymbol{\theta} \mathrm{d} V_0= \\ & \qquad\int_{V_0} \delta \boldsymbol{\theta}^{\mathrm{T}} \cdot\left(\boldsymbol{A}^{\mathrm{T}} \cdot \hat{{{\boldsymbol{D}}}} \cdot \boldsymbol{A} + { { {\boldsymbol{I}} }}\right) \cdot \Delta \boldsymbol{\theta} \mathrm{d} V_0 \\[-4pt]\end{split} $$ (48) 和有限元法不同, 虚单元法中并不能直接给出单元内插值形函数的具体表达式, 但是根据之前构造的形函数向多项式空间的投影算子, 可近似计算变量梯度. 由于式(39)中的映射算子是基于泊松方程求得, 如果将其直接应用到超弹性问题, 便可大幅度减小问题分析的复杂度. 这也是本文一直所推荐的格式. 考虑式(44), 有
$$ \begin{split} &\Delta {{\boldsymbol{\theta}} _h} \approx \left[ {\begin{array}{*{20}{c}} {\nabla {{\boldsymbol{\phi}} ^{\mathrm{T}}}{\boldsymbol{\Pi}} _k^\nabla }&{{{\boldsymbol{0}}^{\mathrm{T}}}} \\ {{{\boldsymbol{0}}^{\mathrm{T}}}}&{\nabla {{\boldsymbol{\phi}} ^{\mathrm{T}}}{\boldsymbol{\Pi}} _k^\nabla } \end{array}} \right]\left\{ {\begin{array}{*{20}{c}} {\Delta {{\boldsymbol{u}}_h}} \\ {\Delta {{\boldsymbol{v}}_h}} \end{array}} \right\}= \\ &\qquad \left[ {\begin{array}{*{20}{c}} {\nabla {{\boldsymbol{m}}^{\mathrm{T}}}{\boldsymbol{\Pi}} _{k*}^\nabla }&{{{\boldsymbol{0}}^{\mathrm{T}}}} \\ {{{\boldsymbol{0}}^{\mathrm{T}}}}&{\nabla {{\boldsymbol{m}}^{\mathrm{T}}}{\boldsymbol{\Pi}} _{k*}^\nabla } \end{array}} \right]\left\{ {\begin{array}{*{20}{c}} {\Delta {{\boldsymbol{u}}_h}} \\ {\Delta {{\boldsymbol{v}}_h}} \end{array}} \right\} =\\ &\qquad \nabla {\boldsymbol{\phi}} _d^{\mathrm{T}}{{\boldsymbol{\Pi}} _k}\left\{ {\begin{array}{*{20}{c}} {\Delta {{\boldsymbol{u}}_h}} \\ {\Delta {{\boldsymbol{v}}_h}} \end{array}} \right\} = \nabla {\boldsymbol{m}}_d^{\mathrm{T}}{{\boldsymbol{\Pi}} _{k{\text{*}}}}\left\{ {\begin{array}{*{20}{c}} {\Delta {{\boldsymbol{u}}_h}} \\ {\Delta {{\boldsymbol{v}}_h}} \end{array}} \right\} \end{split} $$ (49) 其中$\nabla {{\boldsymbol{m}}^{\mathrm{T}}}$和$\nabla {{\boldsymbol{\phi}} ^{\mathrm{T}}}$在式(34)中给出. 此外, 式(49)中的其他项为
$$ \begin{split} \boldsymbol{\Pi}_k=\left[\begin{array}{cc}\boldsymbol{\Pi}_k^{\nabla} & \boldsymbol{0} \\ \boldsymbol{0} & \boldsymbol{\Pi}_k^{\nabla}\end{array}\right], \quad \boldsymbol{\Pi}_{k^*}=\left[\begin{array}{cc}\boldsymbol{\Pi}_{k^*}^{\nabla} & \boldsymbol{0} \\ \boldsymbol{0} & \boldsymbol{\Pi}_{k^*}^{\nabla}\end{array}\right]\qquad\qquad\;\; \end{split} $$ (50) $$ \nabla \boldsymbol{m}_d^{\mathrm{T}}=\left[\begin{array}{cc}\nabla \boldsymbol{m}^{\mathrm{T}} & \boldsymbol{0}^{\mathrm{T}} \\ \boldsymbol{0}^{\mathrm{T}} & \nabla \boldsymbol{m}^{\mathrm{T}}\end{array}\right], \quad \nabla \boldsymbol{\phi}_d^{\mathrm{T}}=\left[\begin{array}{cc}\nabla \boldsymbol{\phi}^{\mathrm{T}} & \boldsymbol{0}^{\mathrm{T}} \\ \boldsymbol{0}^{\mathrm{T}} & \nabla \boldsymbol{\phi}^{\mathrm{T}}\end{array}\right] $$ (51) 如式(40)中所描述, 式(49)中的左右两侧并不是严格相等. 若考虑余项, 定义$ {{\boldsymbol{U}}_h} = {\left[ {\begin{array}{*{20}{c}} {{\boldsymbol{u}}_h^{\mathrm{T}}}&{{\boldsymbol{v}}_h^{\mathrm{T}}} \end{array}} \right]^{\mathrm{T}}} $, 则有
$$ \begin{split} &\Delta \boldsymbol{\theta}_h =\nabla \boldsymbol{\phi}^{\mathrm{T}} \cdot {\boldsymbol{\Pi}}_k \cdot \nabla \boldsymbol{U}_h+\left(\Delta \boldsymbol{\theta}-\nabla \boldsymbol{\phi}^{\mathrm{T}} \cdot {\boldsymbol{\Pi}}_k \cdot \nabla \boldsymbol{U}_h\right) =\\ & \qquad \nabla \boldsymbol{m}^{\mathrm{T}} \cdot {\boldsymbol{\Pi}}_{k^*} \cdot \nabla \boldsymbol{U}_h+\nabla \boldsymbol{\phi}^{\mathrm{T}} \cdot\left(\boldsymbol{I}_k-{\boldsymbol{\Pi}}_k\right) \cdot \nabla \boldsymbol{U}_h \end{split} $$ (52) 其中${{\boldsymbol{I}}_k}$为尺寸和$ {{\boldsymbol{\Pi}} _k} $相同的单元矩阵.
为了得到刚度矩阵, 将式(52)代入式(48), 可得
$$ {\text{ }}\Delta \delta \mathcal{W} = \left\{ {\begin{array}{*{20}{c}} {\delta {\boldsymbol{u}}_h^{\mathrm{T}}}&{\delta {\boldsymbol{v}}_h^{\mathrm{T}}} \end{array}} \right\} \cdot {\boldsymbol{\varLambda }} \cdot \left\{ {\begin{array}{*{20}{c}} {\Delta {{\boldsymbol{u}}_h}} \\ {\Delta {{\boldsymbol{v}}_h}} \end{array}} \right\} $$ (53) 其中
$$ \begin{aligned} & \boldsymbol{\varLambda} = \boldsymbol{\Pi}_{k^*}^{\mathrm{T}} \cdot \int_{V_0} \nabla \boldsymbol{m}_d\left(\boldsymbol{A}^{\mathrm{T}} \cdot \hat{{{\boldsymbol{D}}}} \cdot \boldsymbol{A} + {{\boldsymbol{I}}}\right) \nabla \boldsymbol{m}_d^{\mathrm{T}} \mathrm{ ~ d} V_0 \cdot \boldsymbol{\Pi}_{k^*} + \\ &\quad \left(\boldsymbol{I}_k-\boldsymbol{\Pi}_k\right)^{\mathrm{T}} \cdot \int_{V_0} \nabla \boldsymbol{\phi}_d\left(\boldsymbol{A}^{\mathrm{T}} \cdot \hat{{{\boldsymbol{D}}}} \cdot \boldsymbol{A} + {{\boldsymbol{I}}}\right) \nabla \boldsymbol{\phi}_d^{\mathrm{T}} \mathrm{ ~ d} V_0 \cdot\left(\boldsymbol{I}_k-\boldsymbol{\Pi}_k\right)\end{aligned} $$ (54) 则其切线刚度矩阵需要分为两部分
$$ {\boldsymbol{K}}_{K2}^t = {\boldsymbol{K}}_{K2}^{tc} + {\boldsymbol{K}}_{K2}^{ts} $$ (55) 其中${\boldsymbol{K}}_{K2}^{tc}$为协调刚度矩阵, ${\boldsymbol{K}}_{K2}^{ts}$为稳定刚度矩阵.
由式(54)等号右端第一、二项可分别推导出协调刚度矩阵和稳定刚度矩阵. 协调刚度矩阵为
$$ {\boldsymbol{K}}_{K2}^{tc} = {\boldsymbol{\Pi}} _{k{\text{*}}}^{\mathrm{T}} \cdot {{\boldsymbol{G}}_0} \cdot {{\boldsymbol{\Pi}} _{k{\text{*}}}} $$ (56) 其中
$$ \boldsymbol{G}_0 = \int_{V_0} \nabla \boldsymbol{m}_d\left(\boldsymbol{A}^{\mathrm{T}} \cdot \hat{{{\boldsymbol{D}}}} \cdot \boldsymbol{A} + {{\boldsymbol{I}}}\right) \nabla \boldsymbol{m}_d^{\mathrm{T}} \mathrm{ ~ d} V_0 $$ (57) 式中的积分, 对于$k = 1$时可直接使用单点积分进行计算. 对于$k > 1$, 可通过将多边形单元$E$划分成三角形单元计算, 如图3所示. 式中的核函数为多项式函数, 因此也可以构造特殊格式的高斯公式, 将面积分转化成形单元边界上的线积分.
此外, 稳定刚度矩阵为
$$ {\boldsymbol{K}}_{K2}^{ts} = {\left( {{{\boldsymbol{I}}_k} - {{\boldsymbol{\Pi}} _k}} \right)^{\mathrm{T}}} \cdot {\boldsymbol{G}}_0^s \cdot \left( {{{\boldsymbol{I}}_k} - {{\boldsymbol{\Pi}} _k}} \right) $$ (58) 其中
$$ \boldsymbol{G}_0^s = \int_{V_0} \nabla \boldsymbol{\phi}_d\left(\boldsymbol{A}^{\mathrm{T}} \cdot \hat{{{\boldsymbol{D}}}} \cdot \boldsymbol{A} + {{\boldsymbol{I}}}\right) \nabla \boldsymbol{\phi}_d^{\mathrm{T}} \mathrm{ ~ d} V_0 $$ (59) 虚单元法没有形函数${\boldsymbol{\phi}} $的具体表达式, 但是在合理的${h_E}$的假设下, 可以保证
$$ \int \nolimits_{{V_0}} \nabla {\boldsymbol{\phi}} \cdot \nabla {{\boldsymbol{\phi}} ^{\mathrm{T}}}{\text{d}}{V_0} = O\left( 1 \right) $$ (60) 因此可根据式(58)直接得到稳定刚度矩阵, 其表达式如下
$$ {\boldsymbol{K}}_{K2}^{ts} = \alpha {\left( {{{\boldsymbol{I}}_k} - {{\boldsymbol{\Pi}} _k}} \right)^{\mathrm{T}}}\left( {{{\boldsymbol{I}}_k} - {{\boldsymbol{\Pi}} _k}} \right) $$ (61) 其中$\alpha $为保证${\boldsymbol{K}}_{K2}^{tc}$和${\boldsymbol{K}}_{K2}^{ts}$有相同量级的参数, 具体表示为
$$ \alpha = \frac{{{\alpha _0}}}{{{N_{en}}}}\sqrt {\mathop \sum \limits_{i = 1}^{2{N_{en}}} {{\left( {{\boldsymbol{K}}_{ii}^{tc}} \right)}^2}} $$ (62) 其中${N_{en}}$为矩阵${{\boldsymbol{\Pi}} _k}$的维数, ${\alpha _0}$为稳定化参数. 对于非线性问题, 稳定项的选择对结果的影响较大, 尤其是对于一阶VEM. 一般来说, 稳定化参数的选择过小会有较高的精度, 但是会导致迭代的不收敛. 稳定化参数的选择过大, 非线性迭代的收敛效果较好, 但是精度会有所降低. 通过数值实验, 本文设置$2 \leqslant {\alpha _0} \leqslant 4$. 对于高阶虚单元法, 稳定项的选择对结果的影响较小.
3.2 超弹性问题中内力的计算
和有限元法相同, 在虚单元法中, 内力项的计算需要额外的稳定项. 为了导出内力向量${{\boldsymbol{F}}_{{\text{int}}}}$的表达式, 将式(45)代入式(7)中, 可得
$$ \begin{split} &\left\{ {\begin{array}{*{20}{c}} {\delta {{\boldsymbol{u}}^{\mathrm{T}}}}&{\delta {{\boldsymbol{v}}^{\mathrm{T}}}} \end{array}} \right\} \cdot {{\boldsymbol{F}}_{{\text{int}}}} = \int \limits_{{V_0}} \delta {{\hat {\boldsymbol{E}}}^{\mathrm{T}}} \cdot \hat {\boldsymbol{S}}{\text{d}}{V_0}= \\ &\qquad \int \nolimits_{{V_0}} {\left( {{\boldsymbol{A}} \cdot \delta {\boldsymbol{\theta}} } \right)^{\mathrm{T}}} \cdot \hat {\boldsymbol{S}}{\text{d}}{V_0} \end{split} $$ (63) 进一步考虑式(49)和式(52), 可得
$$ \begin{split} &\left\{ {\begin{array}{*{20}{l}} { \delta {{\boldsymbol{u}}^{\mathrm{T}}}}&{\delta {{\boldsymbol{v}}^{\mathrm{T}}} } \end{array}} \right\} \cdot {{\boldsymbol{F}}_{{\text{int}}}} = \left\{ { \begin{array}{*{20}{c}} {\delta {{\boldsymbol{u}}^{\mathrm{T}}}}&{\delta {{\boldsymbol{v}}^{\mathrm{T}}}} \end{array} } \right\} \cdot \\&\qquad \left( {\int \nolimits_{{V_0}} \nabla {{\boldsymbol{\Pi}} ^{\mathrm{T}}} \cdot {{\boldsymbol{A}}^{\mathrm{T}}} \cdot \hat {\boldsymbol{S}} + {\boldsymbol{F}}_{{{\mathrm{int}}} }^s} \right) \end{split} $$ (64) 则内力向量为
$$ {{\boldsymbol{F}}_{{\text{int}}}} = \int \nolimits_{{V_0}} \nabla {{\boldsymbol{\Pi}} ^{\mathrm{T}}} \cdot {{\boldsymbol{A}}^{\mathrm{T}}} \cdot \hat {\boldsymbol{S}} + {\boldsymbol{F}}_{{{\mathrm{int}}} }^s $$ (65) 其中${\boldsymbol{F}}_{{\text{int}}}^s$为内力向量的稳定项.
和刚度矩阵的稳定项推导方法类似, 可以得到${\boldsymbol{F}}_{{\text{int}}}^s$的具体表达式
$$ {{\boldsymbol{F}}_{{\text{int}}}^s}{ = {{\left( {{{\boldsymbol{I}}_k} - {{\boldsymbol{\Pi}} _k}} \right)}^{\mathrm{T}}} \cdot {\boldsymbol{G}}_0^{Fs} \cdot \left( {{{\boldsymbol{I}}_k} - {{\boldsymbol{\Pi}} _k}} \right) \cdot \left\{ {\begin{array}{*{20}{c}} {{{\boldsymbol{u}}_h}} \\ {{{\boldsymbol{v}}_h}} \end{array}} \right\}} $$ (66) 其中
$$ \boldsymbol G_0^{Fs} = \int_{V_0} \nabla \boldsymbol{\phi}_d\left(\boldsymbol{A}^{\mathrm{T}} \cdot \hat{{{\boldsymbol{D}}}} \cdot \boldsymbol{A}\right) \nabla \boldsymbol{\phi}_d^{\mathrm{T}} \mathrm{d} V_0 $$ (67) 和式(59)相似, 上式也是不可计算的. 但是在合理的${h_E}$的假设下, 可参考刚度矩阵的稳定项, 假设${\boldsymbol{G}}_0^s \approx {\boldsymbol{G}}_0^{Fs}$, 有
$$ {\boldsymbol{F}}_{{\text{int}}}^s \approx {\boldsymbol{K}}_{K2}^{ts} \cdot \left\{ {\begin{array}{*{20}{c}} {{{\boldsymbol{u}}_h}} \\ {{{\boldsymbol{v}}_h}} \end{array}} \right\} $$ (68) 将式(68)代入式(65), 可得到最终的内力向量.
综上所述, 本文所采用的新型虚单元法求解超弹性问题的基本流程如表1所示.
表 1 算法流程框架Table 1. Algorithmic process frameworkAlgorithm flow framework: Novel virtual element method for solving hyperelastic problems 1. Read mesh, define boundary conditions, set material parameters. 2. Calculate the element stiffness matrix and internal force vector:
Calculate the projection operator ${\boldsymbol{\Pi}} _{k*}^\nabla $ based on Eqs. (31) and (38);
Calculate the second P-K stress and constitutive tensor based on Eqs. (14) and (15);
Calculate coordination stiffness matrix based on Eqs. (38), (50), (51), (56), (57) and ${\boldsymbol{F}}$, $ {{\boldsymbol{D}}} $;
Calculate stabilization term based on Eq. (61);
Calculate the tangent stiffness matrix ${\boldsymbol{K}}_{K2}^t$ based on Eq. (55).
3. Use Newton iteration method to solve the displacement and output results for each loading step.3.3 应变及应力的计算
在通过VEM计算出节点的位移之后, 可将位移代入式(49), 则单元内的应变为
$$ \hat {\boldsymbol{E}} = {\boldsymbol{A}} \cdot {\boldsymbol{\theta}} = {\boldsymbol{A}} \cdot \nabla {\boldsymbol{m}}_d^{\mathrm{T}} \cdot {{\boldsymbol{\Pi}} _{k{\text{*}}}}\left\{ {\begin{array}{*{20}{c}} {{{\boldsymbol{u}}_h}} \\ {{{\boldsymbol{v}}_h}} \end{array}} \right\} $$ (69) 然后第二P-K应力${\boldsymbol{S}}$可通过式(14)计算, 柯西应力${\boldsymbol{\sigma}} $可通过第二P-K应力${\boldsymbol{S}}$计算
$$ {\boldsymbol{\sigma}} = \frac{1}{J}{\boldsymbol{F}} \cdot {\boldsymbol{S}} \cdot {\boldsymbol{F}} $$ (70) 其中涉及的变形梯度${\boldsymbol{F}}$可由下式计算
$$ {\boldsymbol{F}}{\text{ = }}{\left( {\nabla {{\boldsymbol{m}}^{\mathrm{T}}}{\boldsymbol{\Pi}} _{k{\text{*}}}^\nabla \left[ {\begin{array}{*{20}{c}} {{{\boldsymbol{u}}_h}}&{{{\boldsymbol{v}}_h}} \end{array}} \right]} \right)^{\mathrm{T}}} + {\boldsymbol{I}} $$ (71) 4. 数值算例
基于上述的求解超弹性材料静力和动力学问题的数值格式, 编写相应的Matlab程序, 采用该程序模拟了超弹性材料的变形行为, 并与有限元的数值结果进行了对比, 比较了两种算法的收敛性和精确性. 所有算例的虚拟单元法的多边形网格都是基于开源软件PolyMesher[42]生成.
4.1 Cook膜问题
本小节研究变截面悬臂梁的Cook问题, 如图4所示. 悬臂梁右端承受均布4 N/mm的剪切载荷, 其几何尺寸分别为: L = 48 mm, ${H_1}$ = 44 mm, ${H_2}$ = 16 mm. 超弹性本构的材料参数$ \lambda $ 和 $ \mu $分别为100 MPa和40 MPa. 在本算例中, 稳定项参数选${\alpha _0} = 2$. 同时, 本算例还与另一种VEM格式进行对比(无稳定项虚单元法, stabilization-free virtual element method, SFVEM [40]).
![]()
图 4 Cook膜问题的边界和加载条件Figure 4. Boundary and loading conditions of Cook membrane problem对于本算例, 分别采用四边形网格和多边形网格进行计算. 为了展示所提出的虚单元法数值格式的性能, 对该几何模型进行均匀网格加密, $x$方向和$y$方向的网格划分份数为${2^N}$. 为了比较虚单元法和有限元法的收敛性, 本文将模型右上角(图4中的点A)的纵向位移作为参考. 对于SFVEM的不同阶数的问题, 应变模式的阶数选择为$l$. 不同算法所得结果如图5所示.
![]()
图 5 虚单元和有限元收敛性的比较Figure 5. Comparison of convergence between virtual element method and finite element method图5中横坐标表示网格单元的细分数, 纵坐标为Cook膜自由端顶点(图4中的点A)沿y方向的最大位移. 从图中可以看出, 对于不同网格, 一阶虚单元法所得结果要比一阶有限元法更好. 对于二阶虚单元法, 其在较低密度网格下便可得到较优秀的结果. 对于四边形网格, SFVEM所得结果和有限元法很接近, 但是由于采用了高阶的应变模式, 因此计算量较大. 二阶VEM求解得到的应力分量($ {\sigma _{xy}} $)云图如图6所示.
![]()
图 6 Cook膜问题在当前构型下的应力$ {\sigma }_{xy} $云图Figure 6. Contour plot of stress $ {\sigma }_{xy} $ for Cook membrane problem at current configuration4.2 冲压问题
本算例研究冲压问题的虚单元法表现. 如图7所示, 方形板顶端和左端水平方向的平动自由度被约束; 底端竖直方向的平动自由度被约束, 均布载荷大小为800 N/mm. 几何尺寸设置为: L = H = 1 mm. 超弹性本构的材料参数$ \lambda $和$ \mu $分别为400.75 MPa和92.5 MPa. 本算例采用四边形和多边形网格进行计算, 稳定项参数选${\alpha _0} = 4$.
为了验证该算法的收敛性, 对该几何模型进行均匀网格加密, $x$方向和$y$方向的网格划分份数为$2 \times {2^N}$和${2^N}$. 通过比较不同网格下模型左上角点(图7中点A)的纵向位移${u_y}$, 比较了虚单元法和有限元的收敛性, 如图8所示.
![]()
图 8 虚单元和有限元收敛性的比较Figure 8. Comparison of convergence between virtual element method and finite element method从图8中可以看出, 对于该问题, 一阶虚单元法在网格密度较低的情况下所得结果精度较差, 但是随着网格的加密, 所得位移结果和一阶有限元法类似. 二阶虚单元法所得结果精度较高. 对于$ N = 3 $ 和 $N = 4$, 不同网格的位移 ${u_y}$ 云图如图9所示.
![]()
图 9 不同网格下冲压问题的位移uy云图Figure 9. Contour plot of displacement ${u_y}$of punch problem for different meshes4.3 半圆拱受压问题
最后, 考虑半圆拱受压的不稳定问题. 如图10所示, 模型的外半径为10 mm, 内半径为9 mm. 模型的右端底部固定, 左端底部受单点固支, 此外其顶部承受竖直向下的集中力$F$, 初始载荷$F = - 0.1\; {\mathrm{N}}$. 超弹性材料参数$ \lambda $和$ \mu $分别为5.56 MPa和1.85 MPa. 使用多边形网格对该半圆拱模型进行离散. 本算例采用一阶VEM进行计算, 稳定项参数选${\alpha _0} = 2$.
![]()
图 10 半圆拱的边界条件和几何尺寸Figure 10. Boundary conditions and geometric size of semicircular arch容易看出在此类边界条件下, 该结构为非稳定结构, 因此采用牛顿法求解容易发生“突跳”和“回弹”现象. 弧长法相比牛顿法可以跟踪位移载荷曲线极值点以后的位移-载荷路径, 因此可以捕捉结构的失稳现象. 因此在采用虚单元法进行模拟时, 需使用弧长法对位移和加载因子进行求解. 弧长法中的初始载荷因子和对应的弧长都为0.025. 在本算例中, 采用一阶虚单元法进行计算, 因此可直接使用点积分计算式(57), 计算效率会高于传统FEM.
模拟得到的半圆拱顶部的载荷-位移曲线如图11所示. 载荷-位移曲线中特殊位置的Mises应力云图如图12(a) ~ 图12(c)所示.
![]()
图 12 半圆拱的载荷-位移曲线及特殊点的应力云图Figure 12. Loading-displacement curve and stress contours of semicircular arch5. 结 论
本文给出了一种求解超弹性问题的一般格式的高阶虚单元法. 由于虚单元法不需要获取单元内部的插值基函数, 因此可以使用任意多边形单元进行几何离散. 和传统的虚单元法不同的是, 通过求解泊松方程的椭圆映射算子, 对矢量场的每个分量分别进行近似, 从而可直接求解非线性弹性力学问题. 由于不需要针对位移或者应变构造椭圆投影算子, 因此该格式更简洁, 且扩展性更强. 通过本文给出的例子可以看出, 虚单元法对前处理更具灵活性和多样化, 其数值结果和传统有限元法很接近, 因此, 该算法具有一定的有效性. 本文给出的求解格式不仅适用于超弹性问题, 还适用于其他非线性力学问题, 如弹塑性问题等.
-
![]()
图 4 Cook膜问题的边界和加载条件
Figure 4. Boundary and loading conditions of Cook membrane problem
![]()
图 5 虚单元和有限元收敛性的比较
Figure 5. Comparison of convergence between virtual element method and finite element method
![]()
图 6 Cook膜问题在当前构型下的应力$ {\sigma }_{xy} $云图
Figure 6. Contour plot of stress $ {\sigma }_{xy} $ for Cook membrane problem at current configuration
![]()
图 8 虚单元和有限元收敛性的比较
Figure 8. Comparison of convergence between virtual element method and finite element method
![]()
图 9 不同网格下冲压问题的位移uy云图
Figure 9. Contour plot of displacement ${u_y}$of punch problem for different meshes
![]()
图 10 半圆拱的边界条件和几何尺寸
Figure 10. Boundary conditions and geometric size of semicircular arch
![]()
图 12 半圆拱的载荷-位移曲线及特殊点的应力云图
Figure 12. Loading-displacement curve and stress contours of semicircular arch
表 1 算法流程框架
Table 1 Algorithmic process framework
Algorithm flow framework: Novel virtual element method for solving hyperelastic problems 1. Read mesh, define boundary conditions, set material parameters. 2. Calculate the element stiffness matrix and internal force vector:
Calculate the projection operator ${\boldsymbol{\Pi}} _{k*}^\nabla $ based on Eqs. (31) and (38);
Calculate the second P-K stress and constitutive tensor based on Eqs. (14) and (15);
Calculate coordination stiffness matrix based on Eqs. (38), (50), (51), (56), (57) and ${\boldsymbol{F}}$, $ {{\boldsymbol{D}}} $;
Calculate stabilization term based on Eq. (61);
Calculate the tangent stiffness matrix ${\boldsymbol{K}}_{K2}^t$ based on Eq. (55).
3. Use Newton iteration method to solve the displacement and output results for each loading step.
下载: 导出CSV
-
[1] Sukumar N, Tabarraei A. Conformal polygonal finite elements. International Journal for Numerical Methods in Engineering, 2004, 61: 2045-2066 doi: 10.1002/nme.1141
[2] Nguyen-Xuan H. A polygonal finite element method for plate analysis. Computers & Structures, 2017, 188: 45-62
[3] Pietro DAD, Eru A. Mathematical Aspects of Discontinuous Galerkin Methods. Springer Berlin, Heidelberg, 2012
[4] Hesthaven J, Warburton T. Nodal Discontinuous Galerkin Methods: Algorithms, Analysis, and Applications. Springer Science & Business Media, 2007
[5] Song C. The Scaled Boundary Finite Element Method. John Wiley & Sons Ltd, 2018
[6] Moes N, Dolbow J, Belytschko T. A finite element method for crack growth without remeshing. International Journal for Numerical Methods in Engineering, 1999, 46: 131-150
[7] Strouboulis T, Babuška I, Copps K. The design and analysis of the generalized finite element method. Computer Methods in Applied Mechanics and Engineering, 2000, 181: 43-69
[8] Veiga L, Brezzi F, Cangiani A, et al. Basic principles of virtual element methods. Mathematical Models and Methods in Applied Sciences, 2012, 23(1): 199-214
[9] Veiga L, Brezzi F, Marini L, et al. Mixed virtual element methods for general second order elliptic problems on polygonal meshes. ESAIM : Mathematical Modelling and Numerical Analysis, 2014, 503: 721-747
[10] Ahmad B, Alsaedi A, Brezzi F, et al. Equivalent projectors for virtual element methods. Computers & Mathematics with Applications, 2013, 66: 376-391
[11] Veiga L, Lipnikov K, Manzini G. The Mimetic Finite Difference Method for Elliptic Problems. Springer Cham, 2014
[12] Brezzi F, Lipnikov K, Simoncini V. A family of mimetic finite difference methods on polygonal and polyhedral meshes. Mathematical Models and Methods in Applied Sciences, 2005, 15: 1533-1551, 10
[13] Tishkin V, Samarskii A, Favorskii A, et al. Operational finite-difference schemes. Differential Equations, 1981, 17: 854-862
[14] Veiga L, Brezzi F. Virtual elements for linear elasticity problems. SIAM Journal on Numerical Analysis, 2013, 52: 794-812
[15] Gain A, Talischi C, Paulino C. On the virtual element method for three-dimensional elasticity problems on arbitrary polyhedral meshes. Computer Methods in Applied Mechanics and Engineering, 2013, 282: 132-160
[16] Artioli E, Veiga L, Lovadina C, et al. Arbitrary order 2D virtual elements for polygonal meshes: Part I, elastic problem. Computational Mechanics, 2017, 60: 643-657
[17] Dassi F, Lovadina C, Visinoni M. A three-dimensional Hellinger-Reissner virtual element method for linear elasticity problem. Computer Methods in Applied Mechanics and Engineering, 2020, 364: 112910
[18] Mengolini M, Benedetto M, Aragón A. An engineering perspective to the virtual element method and its interplay with the standard finite element method. Computer Methods in Applied Mechanics and Engineering, 2019, 350: 995-1023
[19] Chi H, Veiga L, Paulino G. Some basic formulations of the virtual element method (VEM) for finite deformations. Computer Methods in Applied Mechanics and Engineering, 2017, 318: 148-192
[20] Wriggers P, Reddy B, Rust W, et al. Efficient virtual element formulations for compressible and incompressible finite deformations. Computational Mechanics, 2017, 60: 253-268
[21] van Huyssteen D, Reddy B. A virtual element method for isotropic hyperelasticity. Computer Methods in Applied Mechanics and Engineering, 2020 367: 113134
[22] De Bellis M, Wriggers P, Hudobivnik B. Serendipity virtual element formulation for nonlinear elasticity. Computers & Structures, 2019, 223: 106094
[23] Wriggers P, Rust W, Reddy B. A virtual element method for contact. Computational Mechanics, 2016, 58: 1039-1050
[24] Aldakheel F, Hudobivnik B, Artioli E, et al. Curvilinear virtual elements for contact mechanics. Computer Methods in Applied Mechanics and Engineering, 2020, 372: 113394
[25] Shen W, Ohsaki M, Zhang J. A 2-dimentional contact analysis using second-order virtual element method. Computational Mechanics, 2022, 70: 225-245
[26] Cihan M, Hudobivnik B, Korelc J, et al. A virtual element method for 3D contact problems with non-conforming meshes. Computer Methods in Applied Mechanics and Engineering, 2022, 402: 115385 doi: 10.1016/j.cma.2022.115385
[27] Park K, Chi H, Paulino G. On nonconvex meshes for elastodynamics using virtual element methods with explicit time integration. Computer Methods in Applied Mechanics and Engineering, 2019, 356: 669-684
[28] Park K, Chi H, Paulino G. Numerical recipes for elastodynamic virtual element methods with explicit time integration. International Journal for Numerical Methods in Engineering, 2019, 121(1): 1-31
[29] Cihan M, Hudobivnik B, Aldakheel F, et al. Virtual element formulation for finite strain elastodynamics. Computer Modeling in Engineering and Sciences, 2021, 129: 1154-1180
[30] Liu T, Aldakheel F, Aliabadi M. Virtual element method for phase field modeling of dynamic fracture. Computer Methods in Applied Mechanics and Engineering, 2023, 411: 116050
[31] Aldakheel F, Hudobivnik B, Wriggers P. Virtual element formulation for phase-field modeling of ductile fracture. International Journal for Multiscale Computational Engineering, 2019, 17: 181-200
[32] Čertík O, Gardini F, Manzini G, et al. The p- and hp-versions of the virtual element method for elliptic eigenvalue problems. Computers & Mathematics with Applications, 2019, 79: 2035-2056
[33] Meng J, Wang G, Mei L. Mixed virtual element method for the helmholtz transmission eigenvalue problem on polytopal meshes. IMA Journal of Numerical Analysis, 2023, 43(3): 1685-1717
[34] 刘传奇, 许广涛, 魏宇杰. 虚单元计算方法的最新理论与应用进展. 力学进展, 2022, 52(4): 874-913 (Liu Chuanqi, Xu Guangtao, Wei Yujie. Virtual element method: Theory and applications. Advances in Mechanics, 2022, 52(4): 874-913 (in Chinese) Liu Chuanqi, Xu Guangtao, Wei Yujie. Virtual element method: Theory and applications. Advances in Mechanics, 2022, 52(4): 874-913 (in Chinese)
[35] 林姗, 杨永涛, 孙冠华等. 弹塑性力学问题的虚单元法. 固体力学学报, 2020, 41(1): 30-40 (Lin Shan, Yang Yongtao, Sun Guanhua, et al. Elastoplastic mechanical analysis based on the virtual element method. Chinese Journal of Solid Mechanics, 2020, 41(1): 30-40 (in Chinese) Lin Shan, Yang Yongtao, Sun Guanhua, et al. Elastoplastic mechanical analysis based on the virtual element method. Chinese Journal of Solid Mechanics, 2020, 41(1): 30-40 (in Chinese)
[36] 江 巍, 徐建城, 王乐华等. 基于虚单元法的非连续变形分析方法新格式. 岩石力学与工程学报, 2022, 41(1): 106-119 (Jiang Wei, Xu Jiancheng, Wang Lehua, et al. A novel formulation of discontinuous deformation analysis enlightened by virtual element method. Chinese Journal of Rock Mechanics and Engineering, 2022, 41(1): 106-119 (in Chinese) Jiang Wei, Xu Jiancheng, Wang Lehua, et al. A novel formulation of discontinuous deformation analysis enlightened by virtual element method. Chinese Journal of Rock Mechanics and Engineering, 2022, 41(1): 106-119 (in Chinese)
[37] Berrone S, Borio A, Marcon F. Lowest order stabilization free virtual element method for the Poisson equation. arXiv Preprint: 2103.16896
[38] Berrone S, Borio A, Marcon F. A stabilization-free virtual element method based on divergence-free projections. Computer Methods in Applied Mechanics and Engineering, 2024, 424: 116885
[39] Xu BB, Wriggers P. 3D stabilization-free virtual element method for linear elastic analysis. Computer Methods in Applied Mechanics and Engineering, 2024, 421: 116826
[40] Xu BB, Peng F, Wriggers P. Stabilization-free virtual element method for finite strain applications. Computer Methods in Applied Mechanics and Engineering, 2023, 417: 116555
[41] Lamperti A, Cremonesi M, Perego U. A Hu-Washizu variational approach to self-stabilized virtual elements: 2D linear elastostatics. Computational Mechanics, 2023, 71: 1-21
[42] Talischi C, Paulino G, Pereira A, et al. PolyMesher: A general-purpose mesh generator for polygonal elements written in MATLAB. Structural and Multidisciplinary Optimization. 2012, 45: 309-328
计量
- 文章访问数: 174
- HTML全文浏览量: 59
- PDF下载量: 69
