虚拟单元方法中的构型力及其在裂纹问题中的应用

《Computer Methods in Applied Mechanics and Engineering》：Configurational forces in the virtual element method with applications to cracks

　　凯文·施密茨 | 安德烈亚斯·里科尔
卡塞尔大学，力学研究所，Moenchebergstr. 7，卡塞尔，34125，德国

**摘要**
虚拟元素方法（VEM）是有限元方法（FEM）的推广，允许使用任意形状的多边形元素进行空间离散化，甚至包括非凸多边形。这对于涉及移动不连续性、局部应力集中或具有明显不同尺寸的多相问题具有显著优势。构型力（CF）理论为在热力学框架内量化虚拟缺陷位移的驱动力提供了非常灵活的基础，同时也能指示数值计算中的局部网格质量。因此，VEM与CF理论结合在一起，为具有强烈不连续性和尖锐裂纹尖端的裂纹生长数值模拟提供了有效的基础。面临的挑战包括对虚拟元素内形状函数知识的缺乏、在任意多边形元素上的数值积分问题，以及需要自动生成多边形网格的需求。能量-动量平衡的弱形式是通过变分方法推导出来的，该方法将作用在自由表面上的节点CF与驱动缺陷的CF区分开来。此外，高阶虚拟元素固有的内部自由度必须在CF的背景下进行解释。非物理的、所谓的虚假CF被用于改进裂纹尖端的网格设计，最终生成准确的裂纹尖端载荷和裂纹路径模拟结果。

**1. 引言**
在经典力学中，牛顿力无论其物理起源如何，都是通过加速度或反作用力的形式来体验的，后者是由于几何或运动学约束引起的。在公理化背景下数学描述时，它们与材料属性无关。相比之下，构型力（CF）与物体或连续体的配置有关[1]，[2]。它们的作用归因于不均匀性的不可逆虚拟位移，从而导致物体总势能的减少。因此，CF在非平衡热力学的背景下被解释和数学处理。它们构成了作用于材料或结构不连续性（如点缺陷或线缺陷、不同介质之间的界面或裂纹）迁移的能量驱动力。

在断裂力学中，裂纹尖端的CF与J积分[3]相关，J积分是一个能量驱动的载荷量，常用于数值评估裂纹[4]，[5]，[6]，[7]，[8]。其大小与临界能量释放率的阈值相比，决定了裂纹生长的潜力，而其方向构成了裂纹偏转判据的基础。然而，J积分和CF的计算在数值上存在根本差异。CF在数值断裂力学中的优势包括：与节点相关的计算效率高，与其他类型缺陷的相互作用在同一理论基础上处理，以及CF可作为自动局部网格优化的指标[9]，[10]，[11]。

虚拟元素方法（VEM）[12]，[13]是一种稳定的Galerkin方法，它是有限元方法（FEM）到多面体离散化的有效推广。与经典FEM不同，VEM不依赖于基函数的显式构造，而是使用适当的投影将其映射到多项式子空间。这使得可以使用任意阶数的凸和非凸多边形进行离散化，这些多边形与标准有限元兼容，并且即使在大网格畸变的情况下也能产生准确和稳健的结果。因此，VEM已被应用于工程领域的各种问题，包括从线性弹性[14]到非线性和非弹性问题[15]，[16]，以及弹性动力学问题[17]。此外，由于VEM能够直接包含悬挂节点和非凸元素，因此特别适用于模拟具有移动不连续性的问题，例如断裂力学中离散建模的裂纹。

在这方面，VEM首次被Nguyen-Thanh等人[18]应用于线性弹性断裂力学（LEFM）中的裂纹载荷分析，随后Hussein等人[19]将其用于脆性裂纹生长的高效模拟。基于粘聚区模型的裂纹生长实现后来由Artioli等人[20]、Marfia等人[21]和Choi等人[22]提出。在最近的工作中，Chen等人[23]使用基于Hu-Washizu的VEM公式，结合了类似于[19]中提出的元素分割技术。具体来说，作者使用了最初自稳定的虚拟元素网格，并通过选择性沙漏稳定化在裂纹传播过程中可能生成的不稳定元素，从而形成了一个用于模拟二维（2D）脆性裂纹生长的稳健高效算法框架。实际上，这些作者中的几位在VEM中使用了J积分作为等效域积分[24]，而最近的工作[25]还使用了路径独立的轮廓积分来计算矢量J积分。然而，据我们所知，CF理论至今尚未在VEM框架内得到详细阐述。

具体而言，CF被用作指导h-细化优化策略的指标，并指导裂纹生长模拟。与标准方法（特别是FEM）相比，VEM能够通过消除连接粗略区域和细化区域所需的过渡元素或约束，直接实现局部网格细化。此外，VEM即使在存在畸变元素的情况下也表现出稳健性，这些畸变元素可能在局部细化区域自然出现。此外，多项研究证明了VE适用于离散裂纹建模，这得益于用简化的裂纹尖端附近局部网格修改代替频繁的域重网格化[19]，[25]。所提出框架的潜在缺点主要与VEM形式本身有关。例如，由于虚拟元素内的形状函数不是显式提供的，可能会出现困难。此外，本工作中使用的标准VEM公式依赖于稳定化程序，这些程序引入了额外的计算工作和必须仔细选择的数值参数。此外，还需要在一般多边形元素上进行数值积分，这与标准有限元计算在参考元素上的计算不同。

本文的目标是将VEM和CF理论的优势结合起来，创建一种在线性弹性框架内进行高效和准确模拟的数值方法。考虑了线性和二次近似阶数，并在CF的背景下进行了详细阐述。理论框架中的难点在于评估CF需要形状函数的导数，而这些导数在虚拟元素公式中并不是显式提供的。此外，必须将自由边界处的物理合理的CF与驱动缺陷的CF区分开来。基于断裂力学中由强烈不连续性和奇异应力控制的裂纹的经典定义，进行了多项基准测试。研究了非物理的、所谓的虚假CF，包括与内部自由度相关的CF，以评估它们在裂纹尖端载荷中的作用及其在自动局部网格细化判据中的可行性。最后，VEM-CF框架被用于在自适应细化网格上的裂纹生长数值模拟，其中多边形元素被裂纹面切割，裂纹尖端被允许穿透元素。

**2. 线性弹性的虚拟元素方法**
**2.1. 连续介质设置**
考虑一个占据域Ω?R2的线性弹性体。在假设小变形的情况下，其机械行为受经典连续介质力学原理的支配。在没有惯性的情况下，平衡条件由线性动量守恒给出：
(1)
div(σ) + b = 0，
其中σ表示柯西应力张量，b表示单位体积的体力。材料响应被建模为线性弹性，本构关系遵循胡克定律：
(2)
?σ = C·?(u)，
其中C是四阶弹性张量，?(u)是从位移场u导出的无穷小应变张量。应变张量定义为位移梯度的对称部分：
(3)
??(u) = 1/2 · grad(u) + gradT(u)，
其中gradT(·)表示转置梯度。这一设置为后续的变分公式和离散化方案提供了基础。

**2.2. 离散问题**
动量平衡方程(1)的弱形式为：
(4)
{
find u ∈ [H1(Ω)]2，
satisfying a(u, v) = L(v) ? v ∈ [H1(Ω)]2，
}
其中H1(Ω)表示定义在Ω上的第一阶Sobolev空间，v是[H1(Ω)]2中的向量测试函数。选择虚拟位移作为测试函数，a(u, v)是虚拟弹性能量，而L(v)表示施加的体力和外力的虚拟功。连续域Ω?R2由NE个不重叠的多边形Ωi组成的划分Ωh近似，每个多边形的域为Ωi，见图1：
(5)
Ω ≈ Ωh = ??=1?Ωi，
其中Ωi∩Ωj = ? ?i, j ∈ {1, ?,?}, i ≠ j。
与分解Ωh相关的虚拟元素网格表示为Th，E∈Th表示各个元素。为了构建有限维近似，引入了全局虚拟元素空间Vh(Ωh)，它是H1(Ωh)的一个符合子空间，具有C?-连续的分段定义：
(6)
Vh(Ωh) = {
v ∈ H1(Ωh) ∩ C?(Ωh) : v | E ∈ Vh(E) ? E ∈ Th，
}
其中Vh(E)表示每个多边形上待定义的1D局部虚拟元素空间。按照Bubnov-Galerkin方法，所考虑问题的离散弱形式为：
(7)
{
find uh ∈ [Vh(Ωh)]2，
satisfying ah(uh, vh) = Lh(vh) ? vh ∈ [Vh(Ωh)]2，
}
其中[Vh(Ωh)]2表示向量值位移解uh和测试函数vh的2D近似空间。此外，方程(4)和(7)的组成部分分别按元素分解为：
(8)
a(h)(·, ·) = ∑?∈Th aE(h)(·, ·)
和
L(h)(·) = ∑?∈Th LE(h)(·)。
元素级的对称双线性形式aEh(·, ·)表示离散化的虚拟弹性能量：
(9)
?aEh(uh, vh) = ∫?E σ(uh) : ?(vh) dE，
其中σ(uh)是由离散位移场uh产生的应力张量，?(vh)是相应的测试应变，均在线性弹性假设下。线性形式LEh(·)表示施加的体力和外力的离散虚拟功：
(10)
LEh(vh) = ∫?E bh · vhdE + ∫?Et th · vh dE，
其中bh是元素边界?E上的近似体力，th是边界上的近似牵引向量，其中规定了非均匀Neumann边界条件，见图1。

**2.3. 线性弹性的虚拟元素方法**
**2.3.1. 连续介质设置**
考虑一个占据域Ω?R2的线性弹性体。在假设小变形的情况下，其机械行为受经典连续介质力学原理的支配。在没有惯性的情况下，平衡条件由线性动量守恒给出：
(1)
div(σ) + b = 0，
其中σ表示柯西应力张量，b表示单位体积的体力。材料响应被建模为线性弹性，本构关系遵循胡克定律：
(2)
?σ = C·?(u)，
其中C是四阶弹性张量，?(u)是从位移场u导出的无穷小应变张量。应变张量定义为位移梯度的对称部分：
(3)
??(u) = 1/2 · grad(u) + gradT(u)，
其中gradT(·)表示转置梯度。这一设置为后续的变分公式和离散化方案提供了基础。

**2.4. 离散问题**
动量平衡方程(1)的弱形式为：
(4)
{
find u ∈ [H1(Ω)]2，
satisfying a(u, v) = L(v) ? v ∈ [H1(Ω)]2，
}
其中H1(Ω)表示定义在Ω上的第一阶Sobolev空间，v是[H1(Ω)]2中的向量测试函数。选择虚拟位移作为测试函数，a(u, v)是虚拟弹性能量，而L(v)表示施加的体力和外力的虚拟功。连续域Ω?R2由NE个不重叠的多边形Ωi组成的划分Ωh近似，见图1：
(5)
Ω ≈ Ωh = ??=1?Ωi，
其中Ωi∩Ωj = ? ?i, j ∈ {1, ?,?}, i ≠ j。
与分解Ωh相关的虚拟元素网格表示为Th，E∈Th表示各个元素。为了构建有限维近似，引入了全局虚拟元素空间Vh(Ωh)，它是H1(Ωh)的一个符合子空间，具有C?-连续的分段定义：
(6)
Vh(Ωh) = {
v ∈ H1(Ωh) ∩ C?(Ωh) : v | E ∈ Vh(E) ? E ∈ Th，
}
其中Vh(E)表示每个多边形上待定义的1D局部虚拟元素空间。按照Bubnov-Galerkin方法，所考虑问题的离散弱形式为：
(7)
{
find uh ∈ [Vh(Ωh)]2，
satisfying ah(uh, vh) = Lh(vh) ? vh ∈ [Vh(Ωh)]2，
}
其中[Vh(Ωh)]2表示向量值位移解uh和测试函数vh的2D近似空间。此外，方程(4)和(7)的组成部分分别按元素分解为：
(8)
a(h)(·, ·) = ∑?∈Th aE(h)(·, ·)
和
L(h)(·) = ∑?∈Th LE(h)(·)。
元素级的对称双线性形式aEh(·, ·)表示离散化的虚拟弹性能量：
(9)
?aEh(uh, vh) = ∫?E σ(uh) : ?(vh) dE，
其中σ(uh)是由离散位移场uh产生的应力张量，?(vh)是相应的测试应变，均在线性弹性假设下。线性形式LEh(·)表示施加的体力和外力的离散虚拟功：
(10)
LEh(vh) = ∫?E bh · vhdE + ∫?Et th · vh dE，
其中bh是元素边界?Et上的近似体力，th是边界上的近似牵引向量，其中规定了非均匀Neumann边界条件，见图1。

**2.3.2. 本地近似空间**
第2.3节至2.7节的理论框架主要基于Beir?o da Veiga等人[12]，[13]的工作，以及Artioli等人[14]、Mengolini等人[26]和Dassi[27]的工作。设Pk(E)和Pk?2(E)分别是定义在多边形分解的通用元素E上的最高为k和k?2的多项式向量空间。此外，定义了约定P?1(E)={0}。然后，对于选定的精度阶数k∈N≥1，引入定义在E上的有限维空间Vh?V：
(11)
Vh(E) = {
v ∈ H1(E) : Δv ∈ Pk?2(E), v | e ∈ Pk(e) ? e ∈ ?E，
}
其中e表示E的边，形成其边界?E。这一选择，结合方程(6)的定义，对于Vh中的每个通用标量函数vh，有：
(12)
• vh是E上的第一阶Sobolev空间H?的函数，并且在Ωh中连续；
• Δvh是元素内部k?2阶的多项式；
• vh是元素边界?E上每条边e的k阶多项式。

**2.4. 自由度**
在平面弹性框架中，需要2ndE个值来表示[Vh(E)]2中的向量函数vh。虽然有多种选择，但我们使用以下自由度集：
对于多项式阶数k ≥ 1：
• 在多边形顶点节点处有2nv个vh值。
此外，对于k ≥ 2：
• 在多边形边上(k+1)点Gauss-Lobatto求积规则的积分点处有2(k?1)nv个vh值；
• 由|E|?1 ∫E vh · mdE定义的k(k?1)个内部矩，?m ∈ [Pk?2(E)]2。
图2(b)展示了一个示例五边形上的二次插值自由度集。这种特定选择是为了确保与经典有限元方法的兼容性，并便于FEM和VEM的直接耦合。在这个背景下，节点自由度（DOFs）沿着元素边界位于相同的位置，而形状函数的多项式阶数与相应的拉格朗日型有限元假设相同。接下来，引入了一组算子(13)dofi:[Vh(Ωh)]2→Vh(Ωh)(14)vh?dofi(vh)，i=1,?,2nd，其中dofi(vh)用于计算vh在第i个自由度处的值，2nd是离散化的总自由度数。利用这种表示法，元素级上虚拟函数的VEM插值可以表示为(15)vh|E=∑i=12ndEdofi(vh)φi，向量形状函数φi∈[Vh(E)]2满足克罗内克δ性质dofj(φi)=δij，确保每个形状函数唯一对应一个自由度。与传统的有限元公式不同，这些形状函数或局部虚拟元素空间的基函数不会被显式定义。实际上，VEM形式主义主要是为了依赖于函数空间的隐式定义而构建的，这将在后续的阐述中指出。

2.5. 投影算子
VEM的核心思想是使用适当的投影算子从隐式定义的虚拟元素基函数中提取多项式分量。具体来说，将使用算子(16)Πk,E:[Vh(E)]2→[Pk(E)]2vh?vΠ:=Πk,E(vh)，它将属于[Vh(E)]2的函数vh映射到属于向量值多项式空间[Pk(E)]2?[Vh(E)]2的多项式vΠ上。因此，在域分解的每个元素上引入了一个局部空间，即(17)Pk(E):={q:E→R|q(x,y)=∑i+j≤kaijxiyj,aij∈R}。相应地，2D多项式空间[Pk(E)]2包含所有其分量属于Pk(E)的向量函数。这些空间的构建在这里不会详细说明，但第2节末尾会提供更多信息。方程(16)中的投影算子在VEM框架中非常重要，因为它可以在不显式知道方程(11)中假设函数的空间变化的情况下获得基函数的基本信息，从而进行精确计算。通常，这些算子是基于适当的条件获得的，这些条件可以在某种范数下最小化形状函数逼近的误差。因此，根据所考虑的问题和所需的精度属性，可以有多种Πk,E的定义，参见[14]、[28]、[29]。有了这样的投影器，[Vh(E)]2中的任意函数，例如形状函数φi，可以重写为多项式部分和余数的和，如下所示(18)φi=φiΠ+φi?φiΠ(19)=φiΠ+φiR，其中φiΠ:=Πk,E(φi)是φi在[Pk(E)]2上的投影，φiR:=φi?Πk,E(φi)表示未被投影捕获的部分，因为它是[Vh(E)]2?[Pk(E)]2中的函数。

2.6. 元素刚度矩阵
根据[12]，使用方程(19)将双线性形式(9)分解为一致性和稳定性部分(20)aEh(uh,vh)=aE(uΠ,vΠ)+SE(uh,vh)，其中aE(uΠ,vΠ)是一致性项，包含了多项式投影uΠ和vΠ。这一项确保了精确的逼近，并保证了VEM的最佳收敛性质。更准确地说，对于精度为k的VEM，与k-1阶多项式应变相关的虚拟弹性能量可以精确计算，满足(21)aEh(u?h,vh)=aE(u?h,vh)?u?h∈[Pk(E)]2。因此，该方法通过了k阶补丁测试，这一特性称为k一致性[12]。方程(20)中的第二项SE(uh, vh)被称为稳定性形式。这一项考虑了uh和vh的贡献，这些贡献没有被多项式投影捕获。通常，在求解原始变量的方程组时，需要确保VEM形式主义的稳定性。然而，在某些情况下，根据元素几何形状和所需的精度，可能不需要稳定性，或者甚至在构建无稳定性元素时有意避免稳定性，参见[30]、[31]。与投影算子类似，选择稳定性形式有多种可能性。这些方法的共同点是，SE的计算避免了显式知道所包含函数uh和vh的需要，而是使用了对aEh(uh,vh)的虚假能量贡献的粗略近似。在线性弹性方面，所选择的技术不会显著影响方法的精度，但使用最适合模拟目的的形式仍然是有利的。有关理论细节，请参考[12]，有关计算细节，请参考[13]、[26]和[27]。

利用方程(15)和(19)的表示法引入元素级插值(22)uΠ|E=∑i=12ndEdofi(uh)φiΠ,vΠ|E=∑i=12ndEdofi(vh)φiΠ，并将方程(22)与重新表述的双线性形式方程(20)一起插入方程(8)和(7)，将弱形式转换为线性组合(23)∑E∈Th∑i=12ndEdofi(uh)︸uih∑j=12ndE(aE(φiΠ,φjΠ)︸[KEC]ij+SE(φih,φjh)︸[KES]ij)dofj(vh)︸vjh=∑E∈Th∑j=12ndELE(φjΠ)︸rjhdofj(vh)︸vjh，涉及确定BVP位移解的元素相关贡献。这里，[KEC]ij和[KES]ij,i,j=1,?,2ndE表示元素一致刚度矩阵和稳定性刚度的笛卡尔分量，构成元素刚度矩阵(24)KE=KEC+KES。符号uih、vjh和rjh表示在元素自由度处收集离散位移和测试函数值的向量分量以及外部载荷。一旦计算出局部量，就像在FEM中一样，使用组装算法构建它们的全局对应物：(25)K=AE∈ThKEandr=AE∈ThrE。最后，使用标准的有限元求解器程序求解组装好的方程组，其形式为(26)K·U=r，其中U表示收集全局自由度处位移的向量。

2.7. 实现细节
如前几节所述，VEM的构建提供了几个自由度。在这种情况下，相关量的具体选择，即多项式空间的基础、定义投影算子和稳定性形式的条件，在下文中进行了概述。有关数值计算的更多细节，请参考[13]、[26]和[27]，它们提供了所涉及程序的全面描述，并为VEM实现提供了良好的起点。

2.7.1. 多项式空间
对于分解中的每个元素E，[Pk(E)]2的基础Mk(E)定义为(27)Mk(E)={mα(ξE,ηE)}α=1,?,nk，由nk=(k+1)(k+2)个线性独立的基向量mα组成，以下称为缩放单项式。在每个元素上，这些基向量根据缩放坐标定义(28)ξE=x?xˉEhE和ηE=y?yˉEhE，其中xˉE、yˉE表示多边形的质心坐标，hE表示E的直径。限制k≤2，定义多项式空间的缩放单项式集合分别为(29)M0(E)={[10],[01]，(30)M1(E)=M0(E)∪{[?ηξ],[ηξ],[ξ0],[0η]，(31)M2(E)=M1(E)∪{[ξ20],[0ξ2],[ξη0],[0ξη],[η20],[0η2]}，为了提高可读性，省略了缩放坐标的上标(·E)。

2.7.2. 投影算子
在目前的工作中，使用了方程(16)中引入的通用投影算子定义，即(32)?Πk,E?={aE(q,vh?Πk,E?(vh))=0，?q∈[Pk(E)]2P0,E(q?,vh?Πk,E?(vh))=0，?q?∈[Pk(E)]2，满足?(q?)=0，这主要是基于aE(·, ·)的Galerkin正交性条件，在推导VEM方程时是最自然的选择。然而，还需要一个额外的条件P0,E来完整投影算子，以提供唯一投影刚体模式集q?∈[Pk(E)]2的有效条件，对于这些模式有aE(q?,vh)=aE(q?,Πk,E?(vh))≡0。因此，定义了算子(33)P0,E(q?,vh):=∫?Eq?·vhde，它产生了涉及的多项式q?的积分值乘以函数vh。在计算投影函数值时，使用算子在基础Mk(E)中的矩阵表示，即(34)Π?E?∈Rnk×2ndE，其分量为[Π?E?]αi=sα，使得Π?(φi)=∑α=1nksαimα?i=1,?,2ndE。此外，还需要在[Vh(E)]2的基础上表示相同的算子，定义为(35)ΠE?∈R2ndE×2ndE，其分量为[ΠE?]ji=dj，使得Π?(φi)=∑j=12ndEdjiφj?i=1,?,2ndE。

2.7.3. 一致性刚度矩阵
使用方程(9)、(23)、(27)和(34)的定义，计算元素刚度的一致性部分如下(36)KEC=(Π?E?)TG?EΠ?E?，其中G?是一个对称矩阵，包含了与缩放单项式(29) - (31)相关的虚拟弹性能量贡献：(37)G?E∈Rnk×nk，其分量为[G?E]αβ=aE(mα,mβ)?α,β=1,?,nk。

2.7.4. 稳定性刚度矩阵
使用所谓的dofi-dofi稳定性方法完成局部刚度矩阵，这是一种通过近似虚假能量贡献SE(uh, vh)来构建可计算稳定性形式的常用方法：(38)KES=τtr(KEC)(I?ΠE?)T(I?ΠE?)。该方法基于方程(35)中引入的算子表示。矩阵分量的数量级通过刚度矩阵的一致性贡献的迹来缩放，并加入了一个用户定义的参数τ。后者表示一个正数，在这里选择τ=1/2，这是文献[14]中推荐用于线性弹性的值。

3. 材料空间中的构型力理论
构型力（CFs）理论提供了一个热力学框架，用于研究作用在材料力学缺陷上的驱动力。在这种情况下，缺陷代表了边界值问题（BVP）的特征，表现为应力和应变的不均匀性和不连续性，例如由于裂纹、位错或相界引起。与牛顿力不同，构型力的作用不是直接可观察的，而是表现为对材料配置变化的反应。在初始域Ω中由于不可逆缺陷位移而耗散的能量，如图3所示，可以通过虚拟功来表示(39)∫ΩδχdV=∫Ωp(X)·δz(X)dV，其中δχ表示在从初始状态Ω变化到新的材料配置Ω?时，在不可逆虚拟缺陷位移δz(X)处所做的特定虚拟耗散功。因此，p(X)表示材料位置X处的体积特定CF。

为了获得计算离散CFs的公式，首先展示了如何严格推导出包含缺陷的机械BVP的能量-动量平衡的弱形式，遵循[32]中介绍的程序。它与[7]、[8]、[33]等中的现有方法不同，因为测试函数是物理上合理的，最终导致在CFs背景下边界积分的明确解释，例如在裂纹面上。此外，构型体力严格与与缺陷相关的CFs分开，因此结果的驱动力的符号表示了预期的缺陷运动方向。研究将再次限制在无限小的变形上，即X≡x，其中x是变形配置中的坐标。假设Ω再次受到其Neumann边界上的机械表面牵引力t和内部的体力b的作用。此外，现在假设Ω内存在任意分布的（点）缺陷，除此之外，Ω由均匀弹性材料组成。这些缺陷被认为是可移动的，因此能够根据任意的平滑虚拟位移场δz(x)改变其在材料配置中的位置。假设（准）静态条件并忽略热效应，BVP的能量平衡由(40)δu=δwrev+δχ给出，其中δu表示系统的内能密度变化，δwrev是特定的虚拟可逆功。其他导致能量耗散的现象可以忽略，但可以直接纳入能量平衡中。因此，结构的总势能Πtot的减少等于沿虚拟缺陷位移耗散的能量，即(41)?δΠtot=∫ΩδχdV=∫Ω(δu?δwrev)dV。由于虚拟缺陷位移导致的内能密度变化因此为(42)δu=grad(u)·δz。同样的程序也适用于可逆虚拟功：(43)??δwrev=σ:δ?=σ:grad(?)·δz=σ:grad(grad(u))·δz。这里，利用了Cauchy应力σ=σT的对称性来用位移梯度grad(u)替换应变ε。使用乘积规则重写方程(43)，考虑平衡条件方程(1)，并将结果关系和方程(42)代入方程(41)，得到全局变分能量平衡(44)?δΠtot=∫Ω(div(uI?gradT(u)·σ)?gradT(u)·b)·δzdV，其中I是单位矩阵。这个方程给出了由于虚拟缺陷位移导致的能量耗散的变分公式，可以从边界值问题的解中计算得出。数学上，方程(44)构成了能量-动量平衡的弱形式，其中δz充当了足够平滑的测试函数。括号中的第一项代表Eshelby的构型应力或能量-动量张量Q[34]，而后者项现在表示为构型体力g的向量，遵循例如[35]或[33]中的概念。因此，全局能量-动量平衡重写为(45)?δΠtot=∫Ω(p·δzdV=∑γ=1nγPγ·δzγ=∫Ω(div(Q)+g)·δzdV，其中CFs Pγ作用在一组位于xγ∈Ω的虚拟位移缺陷γ=1,?,nγ上。考虑单个缺陷时，δzγ≡δz表示其虚拟位移，因此(47)P·δz=∫Ω(div(Q)+g)·δzdV。通过比较方程(47)的两边基本上可以得到CF的方程，但由于Eshelby张量在不连续处的发散，这不可行。作为数值计算的起点，弱形式是不可避免的，如第4节中将详细说明。对于没有缺陷或任何类型不均匀性的域Ω，P消失，得到局部能量-动量平衡(48)div(Q)+g=0。后者作为已建立的CFs弱公式推导的起点[7]、[33]、[35]，在加权残差方法的框架内引入了一个随机测试函数。关于断裂力学，裂纹面和裂纹前沿在CF理论中都构成缺陷。然而，裂纹驱动力仅限于沿着裂纹前沿的贡献，并且与J积分[32]密切相关。尽管J积分自20世纪60年代末以来在断裂力学中已经得到广泛应用[3]，通常通过沿远离裂纹尖端的轮廓积分或在包含裂纹尖端的足够大的域内积分来计算，但CF在理论上是在裂纹尖端直接确定的。此外，CF分布能够提供对底层网格质量的有良好解释的指标，这些指标非常适合指导网格优化算法[9]、[35]。这一点在与VEM结合使用时特别有用，因为即使是以非常非传统的方式操作网格，也由于对网格设计的要求相对较低而可以轻松实现。

4. VEM离散化中的构型力
4.1. 弱能量-动量平衡的离散公式
从方程(46)出发，对计算域Ωh进行部分积分并引入外向单位法向量n，得到
(49)∫Ωhp·δzdV=∫?Ωh(Q·n)·δzdS+∫Ωhg·δzdV?∫ΩhQ:grad(δz)dV，
将未知的构型力密度p与材料空间中的场量联系起来，这些场量是BVP解的结果。正如我们在3D有限元分析中关于CFs的先前工作[32]所介绍的，并且与有限元框架中关于CFs的基础工作[7]、[8]、[9]、[35]一致，采用了一组形状函数φi∈[Vh(Ωh)]2来离散化δz，
(50)δz=∑i=12nddofi(δz)φi，
从而得到一组二阶方程，
(51)dofi(δz)∫Ωhp·φidV=dofi(δz){∫?Ωh(Q·n)·φidS+∫Ωhg·φidV?∫ΩhQ:grad(φi)dV}。
要求方程(51)独立于足够平滑的测试函数dofi(δz)的离散值成立，
(52)∫Ωhp·φidV=∫?Ωh(Q·n)·φidS+∫Ωhg·φidV?∫ΩhQ:grad(φi)dV，i=1,?,2nd，
从而将材料力场的加权积分值与基于相关自由度的形状函数的数值可计算右侧联系起来。因此，方程(52)左侧的积分是作用在i-th自由度上的构型力的离散表示，随后表示为
(53)Pih=dofi(P):=∫Ωhp·φidV，
其中h表示与涉及的虚拟元素的特征网格大小有明显的关系。因此，结合(53)的方程(52)可以写为
(54)Pih=∫?Ωh(Q·n)·φidS+∫Ωhg·φidV?∫ΩhQ:grad(φi)dV。
此外，为了便于实现，引入了运算符
(55)supp(i):={E∈Th:ihassupportonE}，
将自由度i与其支持的元素块相关联。因此，可以根据该元素块上的积分以及随后的组装来重写方程(54)，得到结果力，
(56)Pih=AE∈supp(i){∫?E∈?Ωh(Q·n)·φide+∫Eg·φidE?∫EQ:grad(φi)dE}，
例如，在有限元框架中，这可以基于标准程序轻松实现。然而，在VEM背景下，如何使用方程(56)仍有一些问题，因为
•涉及的形状函数φi在元素内部是未知的，
•需要对方程进行任意多边形的积分来评估积分。
为了解决这些问题，以下小节详细介绍了在VEM背景下实现方程(56)的方法。

4.2. VEM实现
CF的计算作为后处理步骤实现，基于场Qh和gh，这些场可以追溯到BVP位移解的投影。此外，还需要考虑部分未知形状函数的合理表示φ?i。因此，问题是要推导出一个可计算的形式
(57)Pih≈f(grad(uh),φ?i)，
为每个自由度提供良好的近似或最佳情况下的精确CF评估。为了后续分析方程(56)，定义了积分
(58)I1:=∫?E(Qh·n)·φide，
(59)I2:=∫Egh·φidE，
(60)I3:=∫EQh:grad(φi)dE。
进一步分析方程(58) - (60)中各组成部分的性质，可以发现多项式的度数根据
?Qh|E=12I(σΠ:?Π)?gradT(uΠ)·σΠ∈[P2k?2(E)]2×2，
gh|E=?gradT(uΠ)·bh∈[Pk?1+deg(bh)(E)]2，
φi|?E∈[Pk(?E)]2，
φi|E∈[Vh(E)]2，
grad(φi)|E∈[L2(E)]2×2，
其中deg(bh)表示施加在元素上的体积力的多项式度数。因此，积分方程(58)
(61)I1∝∫?Epde,p∈P3k?2(?E)需要沿着元素边界评估已知的多项式，这可以通过使用适当的求积方案来实现。然而，第二项和第三项，方程(59)和(60)，需要更多的关注。以下阐述仅限于bh在每个元素上允许的deg(bh)≤k?1的情况，以使I2和I3中的多项式度数相等，从而得到
(62)I2,I3∝∫EpvhdE,p∈P2k?2(E)和vh∈Vh(E)，
需要将虚拟函数vh（这不是显式已知的）与E上的多项式进行积分。为了精确评估所使用的域积分，因此需要一个运算符Πk,E*，它满足性质
(63)∫EpvhdE=∫EpΠk,E*(vh)dE，?vh∈Vh(E)，?p∈{P0(E)，如果k=1，则P2(E)；如果k=2。
在方程(63)中，Πk,E*是所谓的“L2-正交投影器，在≤k的多项式空间上”[28]，对于k≤2，它与通过方程(32)中引入的投影算子Πk,E?是一致的。然而，这并不适用于一般的k。总之，方程(56)的数值评估如下：
(64)Pih≈AE∈supp(i){∫?E∈?Ωh(Qh·n)·φide+∫Egh·Πk,E?(φi)dE?∫EQh:grad(Πk,E?(φi))dE}，k≤2。
另外，值得注意的是，在最低阶情况k=1时，每个元素中的Qh是常数。因此，可以将格林定理应用于方程(60)，通过边界积分来评估它，得到
Pih≈AE∈supp(i){∫?E∈?Ωh(Qh·n)·φide+∫Egh·Π1,E?(φi)dE?∫?EQh:(n?φi)de}，k=1，其中?表示外向法向量n和形状函数的叉积。最后，需要评论方程(64)中积分的数值实现。考虑到第一个被积函数，可以在元素边界?E上使用高斯求积。然而，仍然需要注意被积函数是一个3k?2次的多项式，如前所述。具体来说，如果VEM近似阶数k>1，可以直接检查每个元素边缘上的k+1个可用积分点（边缘节点，参见图2）是否足以通过高斯-洛巴托求积准确评估积分。为了纠正这种不匹配，限制k不超过2，采用了边界积分贡献的解析预积分。这提供了基于高斯-洛巴托积分点的数值积分公式的权重，这些权重是由所选自由度决定的，详见小节2.4。因此，积分表示为元素边缘贡献的总和：
(65)I1=∫eqh·φide=∑e∈?E∫eqh·φide，
其中qh=Qh·n∈P2k?2(?E)和φi∈Pk(?E)。
在关于k的假设下，有2k?2≤k，因此每个边缘上的qh由其值在2(k+1)个自由度处唯一确定。因此，在每个边缘e上表示为
(66)qh|e=∑j=12(k+1)dofj(qh)φj*，
其中φj*是边缘e上的(2k?2)阶近似的形状函数。将方程(66)代入方程(65)，每个边缘上的积分重写为：
(67)∫eqh·φide=∑j=12(k+1)dofj(qh)∫eφj*·φide
(68)=∑j=12(k+1)dofj(qh)wji。
因此，修改后的积分权重wij在预处理步骤中获得，因为它们仅依赖于k和元素几何形状，而不是BVP的解。然后，通过在后处理阶段应用积分公式(68)来高效评估方程(64)的边界积分。
根据方程(59)和(60)，积分I2和I3的贡献由一定度数的多项式组成，这些多项式必须在对多边形域上进行积分。在任意多边形的情况下，这些积分不能像在有限元框架中处理三角形或四边形那样直接通过将被积函数映射到参考域来计算。尽管如此，例如，可以使用虚拟元素的三角剖分来使用数值积分公式，例如高斯求积。其他可能性包括基于格林定理的解析公式[36]的积分，或者结合矩积分[37]的积分项的泰勒级数展开[37]，或者为多边形积分量身定制的求积方案[38]、[39]。对于当前的实现，采用了后一种方法，更具体地使用了Chin和Sukumar[38]、[40]最初提出的缩放边界求积方案，这些权重和积分点在预处理期间获得，而数值积分在后处理阶段进行。

5. 数值结果
本节展示了所提出框架在线性弹性断裂力学理论中的应用能力。考虑了一系列基准问题，这些问题涉及验证推导出的CF公式、CFs作为h-细化的误差指标的适用性，以及它们作为裂纹生长模拟的矢量驱动力的使用。

5.1. 单轴拉伸下的双材料界面
5.1.1. 模型
使用受单轴拉伸作用的双材料二维域来验证CF公式。问题如图4所示。假设材料行为为线性弹性、各向同性，无量纲材料参数为EI=2，νI=0和EII=1，νII=0，得到的应变分别为：
(69)εxx(x,y)={σxx/EI,0≤x<1；σxx/EII,14 线性Voronoi单元，稳定
V2 2 >4 二次Voronoi单元，稳定
模拟结果如图5(a) - (c)所示。节点CF与离散化以及Eshelby张量的Frobenius范数||Q||F一起展示。与解析值(75)∥Q(x)∥F≈{0.35,0≤x<1；0.71,1
下载：下载高分辨率图像（414kb）
下载：下载全尺寸图像
图5. 双材料界面问题的数值研究。(a) - (c)：不同网格和不同近似阶数得到的eshelby应力的frobenius范数||q||f的等高线图。(d) 网格(a) - (c)界面节点处的非零cf分量值

5.2. 混合模式加载静止裂纹的分析
构型力的一个显著优点是它们能够提供裂纹尖端的矢量驱动力，从而估计裂纹载荷和裂纹生长方向，这与矢量j积分[41]一致。然而，由于裂纹尖端附近存在奇点，计算这些量具有挑战性。特别是当使用多项式近似位移场时，如标准vem中那样，如果不采取适当措施，很难在奇点附近获得准确的cf。通常通过在后处理过程中结合解析知识或应用数值技术（如h-、p-或r-细化及其组合）来改进。这个例子展示了cf框架在处理由于u∝r而具有弱奇点的问题时的特性。cfs既用于指导h-自适应细化，也用于量化裂纹尖端的载荷。

5.2.1. 模型
一个长度为a的通用裂纹受到平面应变下的混合模式加载。该问题通过考虑一个半径为rˉ=1的圆来建模，该圆包含裂纹沿θ=±π,r∈[0,rˉ]的线段，在位于裂纹尖端的极坐标系中。假设材料为线性弹性、均匀且各向同性，无量纲参数为e=1和ν=0.3。为了在整个计算域内强制一个特征位移解，边界条件是根据所需的裂纹尖端载荷反向规定的。为此，考虑了williams系列[42]的前两个元素，包括用应力强度因子(sifs) [43] ki和kii表示的近尖端位移场。此外，还叠加了一个恒定应力txx，通常称为t应力，以将方程(64)中的裂纹表面积分包含进来，否则该积分将消失。相应的位移场表示为[44]
(76)u(r,θ)=12μ(r2π(kifi(1)(θ)+kiifii(θ))+rtxx4fi(2)(θ))，
其中μ表示剪切模量。角度函数定义为
(77)fi(n)(θ)=[[κ+(?1)n+n2]cosn2θ?n2cos(n2?2)θ[κ?(?1)n?n2]sinn2θ+n2sin(n2?2)θ]，
(78)fii(θ)=[[?κ?32]sin12θ+12sin?32θ[κ?32]cos12θ+12cos?32θ]，
其中κ=3?4ν。远端载荷指定为ki=0.1，kii=?0.1和txx=?ki/2π。 下载：下载高分辨率图像（414kb） 下载：下载全尺寸图像 图5. 双材料界面问题的数值研究。(a) - (c)：不同网格和不同近似阶数得到的eshelby应力的frobenius范数||q||f的等高线图。(d) 网格(a) - (c)界面节点处的非零cf分量值 5.2. 混合模式加载静止裂纹的分析 构型力的一个显著优点是它们能够提供裂纹尖端的矢量驱动力，从而估计裂纹载荷和裂纹生长方向，这与矢量j积分[41]一致。然而，由于裂纹尖端附近存在奇点，计算这些量具有挑战性。特别是当使用多项式近似位移场时，如标准vem中那样，如果不采取适当措施，很难在奇点附近获得准确的cf。通常通过在后处理过程中结合解析知识或应用数值技术（如h-、p-或r-细化及其组合）来改进。这个例子展示了cf框架在处理由于u∝r而具有弱奇点的问题时的特性。cfs既用于指导h-自适应细化，也用于量化裂纹尖端的载荷。 5.2.1. 模型 一个长度为a的通用裂纹受到平面应变下的混合模式加载。该问题通过考虑一个半径为rˉ=1的圆来建模，该圆包含裂纹沿θ=±π,r∈[0,rˉ]的线段，在位于裂纹尖端的极坐标系中。假设材料为线性弹性、均匀且各向同性，无量纲参数为E=1和ν=0.3。为了在整个计算域内强制一个特征位移解，边界条件是根据所需的裂纹尖端载荷反向规定的。为此，考虑了Williams系列[42]的前两个元素，包括用应力强度因子(SIFs) [43] ki和kii表示的近尖端位移场。此外，还叠加了一个恒定应力txx，通常称为t应力，以将方程(64)中的裂纹表面积分包含进来，否则该积分将消失。相应的位移场表示为[44] (76)u(r,θ)=12μ(r2π(KIfI(1)(θ)+KIIfII(θ))+rTxx4fI(2)(θ))， 其中μ表示剪切模量。角度函数定义为 (77)fi(n)(θ)=[[κ+(?1)n+n2]cosn2θ?n2cos(n2?2)θ[κ?(?1)n?n2]sinn2θ+n2sin(n2?2)θ]， (78)fii(θ)=[[?κ?32]sin12θ+12sin?32θ[κ?32]cos12θ+12cos?32θ]， 其中κ=>
下载：下载高分辨率图像（414kb）
下载：下载全尺寸图像
图5. 双材料界面问题的数值研究。(a) - (c)：不同网格和不同近似阶数得到的eshelby应力的frobenius范数||q||f的等高线图。(d) 网格(a) - (c)界面节点处的非零cf分量值

5.2. 混合模式加载静止裂纹的分析
构型力的一个显著优点是它们能够提供裂纹尖端的矢量驱动力，从而估计裂纹载荷和裂纹生长方向，这与矢量j积分[41]一致。然而，由于裂纹尖端附近存在奇点，计算这些量具有挑战性。特别是当使用多项式近似位移场时，如标准vem中那样，如果不采取适当措施，很难在奇点附近获得准确的cf。通常通过在后处理过程中结合解析知识或应用数值技术（如h-、p-或r-细化及其组合）来改进。这个例子展示了cf框架在处理由于u∝r而具有弱奇点的问题时的特性。cfs既用于指导h-自适应细化，也用于量化裂纹尖端的载荷。

5.2.1. 模型
一个长度为a的通用裂纹受到平面应变下的混合模式加载。该问题通过考虑一个半径为rˉ=1的圆来建模，该圆包含裂纹沿θ=±π,r∈[0,rˉ]的线段，在位于裂纹尖端的极坐标系中。假设材料为线性弹性、均匀且各向同性，无量纲参数为e=1和ν=0.3。为了在整个计算域内强制一个特征位移解，边界条件是根据所需的裂纹尖端载荷反向规定的。为此，考虑了williams系列[42]的前两个元素，包括用应力强度因子(sifs) [43] ki和kii表示的近尖端位移场。此外，还叠加了一个恒定应力txx，通常称为t应力，以将方程(64)中的裂纹表面积分包含进来，否则该积分将消失。相应的位移场表示为[44]
(76)u(r,θ)=12μ(r2π(kifi(1)(θ)+kiifii(θ))+rtxx4fi(2)(θ))，
其中μ表示剪切模量。角度函数定义为
(77)fi(n)(θ)=[[κ+(?1)n+n2]cosn2θ?n2cos(n2?2)θ[κ?(?1)n?n2]sinn2θ+n2sin(n2?2)θ]，
(78)fii(θ)=[[?κ?32]sin12θ+12sin?32θ[κ?32]cos12θ+12cos?32θ]，
其中κ=3?4ν。远端载荷指定为ki=0.1，kii=?0.1和txx=?ki/2π。>最终，方程（76）被作为狄利克雷边界条件施加在Ωh的外边界上，而裂纹面保持无牵引力状态。下载：下载高分辨率图像（228KB）下载：下载全尺寸图像图6. 一个通用裂纹模型，其长度为a，受到预定义的裂纹尖端载荷（KI, KII）的作用，以及尺寸为rˉ的圆形计算域Ωh。5.2.2. 网格自适应性初始网格如图6所示，由64个Voronoi单元组成，其中包括裂纹尖端的一个非凸元素。采用了k=1（V1单元）和k=2（V2单元）的近似方法。h细化策略基于求解-估计-标记-细化的过程，在每个求解步骤后评估自由度（DOFs）处的CF，并将其用作细化指标。采用基于中点的细化技术，在单元边缘的中点和单元质心处引入新节点。新单元的构建确保每个单元包含四个节点。如果一个待细化的单元E有一个在先前步骤中已经被细化的相邻单元，则使用共同边上的现有中点节点，而不是创建新节点，以防止生成严重变形的单元。图7显示了一个代表性网格，其中基于方程（64）进行后处理并将结果组装成矢量量。从物理角度来看，构型力场应仅包括裂纹尖端力以及作用在两个裂纹面上的力。然而，由于方程（64）本质上已经补偿了裂纹面上的贡献，因此应该只剩下裂纹尖端的贡献。尽管如此，所谓的虚假CFs主要出现在裂纹尖端周围的有限域中。这些CFs缺乏物理意义，已知它们指示应力和应变不连续性，这可能与网格质量不足有关。因此，如果某个单元至少包含一个满足条件（79）∥Pnodeh∥≥η∥Pmaxh∥的节点，则该单元被标记为需要细化，其中∥Pnodeh∥是矢量节点CF的幅度，∥Pmaxh∥表示当前网格迭代的矢量CFs的最大绝对值。参数η用于改变细化的阈值。此外，如果没有任何节点满足方程（79）或最小单元面积与最大单元面积的比率|E|min/|E|max≤0.02，则终止该过程。图8和图9显示了使用k=2和η=0.05获得的一个代表性结果。在中间步骤中，非凸裂纹尖端单元被三角剖分，如图8和图9（b）所示，以实现基于实现的中点细化规则的稳健自动细化，这要求单元质心位于多边形内核内。这对于凸多边形是预先保证的，因此我们仅对凸多边形进行细化，即使非凸多边形仍然可以满足要求。下载：下载高分辨率图像（181KB）下载：下载全尺寸图像图7. 裂纹尖端周围区域的CF分布。下载：下载高分辨率图像（444KB）下载：下载全尺寸图像图8. η=0.05时增量细化过程的各个阶段，以及k=2情况下CF的幅度。下载：下载高分辨率图像（413KB）下载：下载全尺寸图像图9. η=0.05时增量细化过程的各个阶段，以及k=2情况下CF矩（CFMs）的幅度。图8和图9（c）显示，网格在裂纹尖端附近被强烈细化。此外，节点CF和与内部力矩自由度相关的值在域的内部减少，而在裂纹尖端附近集中。图10提供了裂纹尖端区域的详细视图以及断裂力学分析相关的应力。比较初始和最终应力，它们相差大约一个数量级，空间分辨率显著提高。下载：下载高分辨率图像（375KB）下载：下载全尺寸图像图10. 对初始网格和细化网格上施加的混合模式载荷的相关应力进行比较。此外，还研究了位移解。图11展示了用于评估裂纹开口位移（COD）Δu|i,i=1,2的量，在位于裂纹面上两个不同节点对处，其中一个节点对位于r1≈0.06（靠近裂纹尖端），另一个位于距离奇异性r2≈0.31的位置。图12显示了COD幅度的相对误差（80）err(Δuh|i)=∥Δuh|i?Δu∥∥Δu∥，结果针对不同阈值η以及整个网格的均匀细化进行了绘制。图13展示了特定选择η=0.05和k=1及k=2时的两个结果网格。对于k=1（图12（a）和k=2（图12（b）），自适应细化加速了COD误差的减少，特别是在靠近奇异性的节点对处，这一点用实线表示。然而，当η=0.25时，在中间距离r2处的准确性下降，且改进不再一致。这种效应是由于在裂纹尖端附近高度局部化的网格细化造成的，这在较大距离上几乎不影响位移解。然而，较小的阈值在相关位置提供了持续且显著的改进。总体而言，结果表明CFs是驱动与VEM和混合模式裂纹问题相关的局部网格优化的可靠指标。下载：下载高分辨率图像（140KB）下载：下载全尺寸图像图11. 裂纹处的术语和位移跳变Δu|i的定义。下载：下载高分辨率图像（419KB）下载：下载全尺寸图像图12. 不同细化阈值η下裂纹面上位移跳变Δu|i的相对误差。节点对位于r1≈0.06，r2≈0.31。下载：下载高分辨率图像（462KB）下载：下载全尺寸图像图13. η=0.05时的最终网格。5.2.3. 裂纹尖端载荷分析基于先前获得的自适应网格，随后基于CFs评估裂纹尖端载荷。为此采用了两种方法：首先，直接使用从方程（64）在裂纹尖端自由度处获得的CF，这是一种计算效率高的方法来获得裂纹驱动力。其次，引入了对CF公式的轻微调整，并结合了在裂纹尖端周围有限区域内常用的域类型评估方法，例如在[4]和[32]中提出的方法。对于研究的第一部分，裂纹尖端CF被重写为（81）Ptip=Pxtip[1tan(θtip)]，其中θtip=arctan(Pytip/Pxtip)，以分别考虑向量的幅度和角度。从断裂力学的角度来看，Pxtip理论上提供了自相似无限小裂纹扩展的能量释放率，因此对应于矢量J积分的第一分量。此外，θtip表示在能量最小化[45]背景下缺陷迁移的方向，并且通常用于根据J积分向量准则[46]，[47]评估裂纹倾斜角。因此，所考虑问题的分析参考公式为：（82）Pxref=J1=1?ν2E(KI2+KII2)=0.0182，（83）Pyref=?1?ν2E2KIKII=0.0182，（84）θref=arctan(J2J1)=π4。数值结果θ?tip=θtip/θref和Px?tip=Pxtip/Pxref在图14中展示，对于k=1和k=2以及不同细化阈值下的自由度数量进行了绘制。归一化角度θ?tip收敛到令人满意的值，尽管对于k=1来说不是单调的，而Px?tip也显示出收敛趋势，但得到的值比初始估计值有更大的误差。此外，θ?tip和Px?tip的最终值几乎不受η的影响。对于k=1，θ?tip的值约为0.9，对于k=2，最终结果接近1。然而，Px?tip的最终值低于0.5，基本上不受近似阶数和细化阈值的影响。尽管裂纹尖端CF在评估裂纹尖端载荷方面不够准确，但结果对于使用CFs作为裂纹起始角度的指标来说是令人鼓舞的。通过直接将方程（64）应用于裂纹尖端节点，可以以令人满意的准确性和最小的计算工作量获得角度。下载：下载高分辨率图像（471KB）下载：下载全尺寸图像图14. 网格细化过程中的裂纹尖端CF。实线表示归一化的裂纹倾斜角，虚线表示直接从方程（64）计算出的CF向量的归一化第一分量。研究的第二部分旨在基于节点CF准确评估裂纹尖端载荷。该方法主要基于所考虑域的全局CF平衡[4]，[32]。因此，考虑在包含裂纹尖端的有限域Ωγ?Ωh中的虚假CFs，得到（85）P=∑xi∈ΩγPh|xi≡J，其中Ph|xi用于表示在位置xi处组装的CF向量。此外，P等同于J积分向量J。对于k=1，因此总结了作用在Ωγ内所有节点上的矢量CFs。对于k=2，需要纳入额外的力矩自由度，尽管它们代表的是单元级的平均值，参见图9，这些值并不与唯一的空间位置xi相关联。因此，力矩自由度被分配给每个单元的质心（xˉE,yˉE）处的矢量CF。为了方便起见，Ωγ被定义为以裂纹尖端为中心、半径为rγ的圆形域。因此，包括位于质心处的力矩在内的所有与裂纹尖端的欧几里得距离≤rγ的CF向量都包含在求和中。方程（85）已知可以显著提高精度，特别是针对第一分量Px。然而，先前的工作表明，在混合模式情况下，如果不对方程（64）中的边界积分进行仔细处理，Py则不是这样。为此，通过在相应的积分中直接纳入无牵引力裂纹面的条件σnn=σtn=0来修改这一项，其中σnn和σtn是边界上的法向和切向应力。因此，在Ωγ内的所有位于裂纹面上的节点（裂纹尖端节点除外）处设置σxy=σyy=0。图15（a）和（b）显示了相对于评估域半径的Px和Py的结果。计算是基于图6中所示的初始粗网格以及之前获得的η=0.01的细化网格进行的，对于k=1和k=2。结果表明，在rγ=0时最初无用的结果得到了显著改进，并且对于较大的rγ获得了与域无关的结果。在粗网格上，仅获得了中等精度。相比之下，细化网格获得了Px的极好一致性结果和Py的良好结果。值得注意的是，Px的计算当然可以与准确预测的角度θtip结合，以获得Py的非常好结果。这些结论强调，基于CF的评估能够提供裂纹尖端载荷的幅度和裂纹生长方向。最后，应该注意的是，在[32]中开发了一种方法，用于在FEM背景下提高Py的准确性，这也适用于VEM计算。下载：下载高分辨率图像（440KB）下载：下载全尺寸图像图15. 从域评估方法计算的CFs的比较，对于初始（a）和细化（b）网格。5.3. 裂纹生长模拟：带有孔的SENB试样第三个示例展示了VEM-CF框架在脆性裂纹生长模拟中的能力。使用了一个常用的基准问题，即带有孔的单边缺口弯曲（SENB）试样。文献[48]中详细记录了该问题的多种配置的实验获得的裂纹轨迹。实验装置由一块20英寸×8英寸（50.8厘米×20.32厘米）的PMMA板和三个直径为0.5英寸的孔组成，这些孔位于距载荷施加点一定距离的位置。裂纹插入的位置由距载荷施加垂直线的距离b0指定。此外，引入了不同的初始裂纹长度a0，导致不同配置下的裂纹路径明显不同。数值建模基于均匀、弹性且各向同性的材料，其杨氏模量为474000 psi（约3268 MPa），在平面应变条件下。泊松比假设为0.35，这对于PMMA来说是合适的选择。施加的载荷设置为F=1磅（lb）（约4.45 N）。研究了三种不同的初始配置。裂纹路径的数值模拟以准静态、增量方式进行。增量过程包括求解当前的边界值问题（BVP），在后处理中计算裂纹倾斜角，最后将裂纹延长一个有限长度Δa，从而得到下一个增量中需要求解的新BVP。VEM方法在此方面的优势在于能够直接追踪裂纹路径，正如Hussein等人[19]所证明的那样。选择了一种类似的技术来模拟裂纹扩展，通过切割现有网格来实现，同时保留了上一步中的大部分元素信息，包括局部刚度矩阵和投影算子。主要使用边长为0.25英寸的均匀四边形网格，孔洞处除外。为了简化计算，使用了k=1的线性元素。此外，在裂纹周围的区域控制了网格质量，将元素细化到直径不小于hE=0.05英寸。裂纹扩展角θc的评估基于公式(86)θc=arctan(PyPx)进行，其中CF是根据公式(85)计算得出的，适用于半径为rγ=0.4英寸的圆形域。最后，每步的裂纹延长量Δa根据评估的扩展角通过公式(87)Δa(θc)={Δamax?8(Δamax?Δamin)π|θc|,|θc|≤π/8,Δamin,else}来控制，以在较大偏转角的情况下“减缓”裂纹扩展。参数选择为Δamax=0.1英寸和Δamin=0.01英寸。图16显示了a0=1.5英寸和b0=5英寸时的初始BVP及所用网格。图17展示了与此初始设置相对应的不同增量阶段的模拟结果，以及裂纹尖端附近的细节。图17(b)显示了裂纹尖端即将穿透孔洞之前的情况。图18(a)至(c)展示了与[48]中的实验参考路径相比的最终裂纹轨迹。可以看出，在所有三种情况下都获得了极好的一致性。裂纹的初始倾斜角、接近的孔洞以及孔洞的穿透点都被准确预测。需要指出的是，由于根据图14，即使仅从裂纹尖端节点也可以相当准确地预测角度，因此评估域半径rγ可以选择得较小。下载：下载高分辨率图像（489KB）下载：下载全尺寸图像图16. a0=1.5英寸，b0=5英寸时的初始BVP。下载：下载高分辨率图像（406KB）下载：下载全尺寸图像图17. 不同模拟阶段的裂纹路径和非凸形裂纹尖端元素。下载：下载高分辨率图像（239KB）下载：下载全尺寸图像图18. 不同配置(a0, b0)下模拟裂纹路径与实验观察结果的比较[48]。6. 结论在变分框架内推导出了能量-动量平衡的弱形式，并将其转化为数值框架，从而在虚拟元素离散化中得到作为缺陷驱动量的节点CF。通过分离自由表面（如裂纹面）的贡献，可以方便地获得与缺陷相关的CF。虚拟元素中缺乏形状函数的问题通过在CF的离散化公式中引入投影算子来弥补。积分是在基于广义求积方案的任意多边形元素上进行的。在二次近似阶数下，作为内部自由度的矩会产生虚假的CF，这些虚假CF必须纳入精确的裂纹尖端载荷分析中。结合节点虚假CF，它们为应力集中位置的自动网格改进提供了适当的指标。可以证明，混合模式裂纹尖端载荷的计算是准确的，但不可避免地需要实施局部网格调整。增量裂纹偏转几乎不受虚假CF的影响，裂纹路径的模拟与实验结果非常吻合。CRediT作者贡献声明Kevin Schmitz：撰写 – 审稿与编辑，撰写 – 原始草稿，可视化，验证，软件，方法论，研究，形式分析。Andreas Ricoeur：撰写 – 审稿与编辑，撰写 – 原始草稿，监督，项目管理，方法论，概念化。