刘传奇 许广涛 ,2 魏宇杰 ,2,3,*

1 中国科学院力学研究所非线性力学国家重点实验室,北京 100190

2 中国科学院大学未来技术学院,北京 100049

3 中国科学院大学工程科学学院,北京 100049

1 发展历程

虚单元方法(virtual element method,VEM)的基本驱动来源于对任意形状多边形的处理.传统意义的有限元或有限差分需要将具备显著几何特征的物理实体离散化而求解,这在一定程度上丧失了对实体几何信息的“宏观”描述.而实际的工程需求中,越来越多的计算涉及到处理特定几何体结构,如非凸性多边形的变形以及复杂结构的接触等问题.2013年,由意大利 Milano-Bicocca 大学的 L.Beirão da Veiga 教授组首先提出可适用于任意形状多边形的虚单元数值方法(Beirão Da Veiga et al.2013a,2013b;Beirão da Veiga et al.2014).由于该方法和有限元方法的高度兼容性,同时在处理悬挂节点(适用于含任意数目节点的单元)、接触(均可转化为点-点接触)、晶体变形(对形状、晶格匹配等无特殊要求)等问题方面具有特定优势,近年来得到计算力学领域的广泛关注与发展.国际上不同的团队将该方法从理论上进行了拓展,并应用到不同的力学问题,如弹性问题 (Gain et al.2014,Artioli et al.2017a)、接触问题 (Wriggers et al.2016,Wriggers and Rust 2019)、高阶单元 (Beirão da Veiga et al.2017a)、大变形 (Beirão da Veiga et al.2015,Wriggers and Rust 2019)、断裂问题 (Hussein et al.2018,Aldakheel et al.2019a),等工程科学/力学问题的计算.

对比有限元方法,虚单元法也需要对几何空间进行离散,通过形成并求解线性方程组对实际问题进行近似.不同之处在于: (1)有限元中所采用的近似函数为显式多项式函数;虚单元法中,近似函数除了多项式函数外,还有在单元边界为连续多项式且单元内部满足某些条件(如进行拉普拉斯运算后为多项式)的函数.这些函数在单元域内无显式表达,这也是虚单元方法中“虚”字的由来.(2)有限元方法中自由度均是近似函数的取值,虚单元方法通过定义合理的自由度,在形成刚度矩阵时避免计算单元内部近似函数的取值;(3)有限元的刚度矩阵中仅含有一项,虚单元法的刚度矩阵包含了协调矩阵和稳定矩阵,以此来保证计算的收敛.

为了实现摘要中所提出的目标,本文将采取有别于传统综述论文的行文方式,大幅度地增加了关于虚单元法基本原理和最新理论的详细介绍.这样处理的目的是双重性的: 一方面,读者能够了解该方法在一些典型问题中所体现的有别于有限元方法的能力;与此同时,对该方法感兴趣的读者将可以通过对本文的通读和理解,全面掌握虚单元法的理论基础和实现途径,为利用该方法来发展适用于自身面临的科学与工程计算问题提供一份全面的参考资料.基于这一设想,文章的基本框架为: 第 2 章简略介绍有限元法与虚单元法的区别与联系,以便具有力学背景的同行对虚单元法有个直观了解;第 3 章详细介绍虚单元法基本原理方面的发展,以熟知的泊松方程为对象,介绍了虚单元法的核心理论;第 4 章针对线弹性问题的求解,全面综述了求解过程中涉及自由度、单元刚度矩阵构建、相应计算流程的具体方法和步骤,并给出了典型的计算实例来考察这一计算方法;第 5 章详细介绍了当前虚单元在非线性问题、断裂问题、接触问题、多场耦合等涉及复杂几何边界处理的典型应用;第 6 章总结该方法的利弊并展望虚单元法的发展前景.

为便于介绍,需要在此做一些符号及其运算的基本设定.以二维问题为例,下文中黑斜体表示张量（张量函数）或矩阵,比如位移场u,刚度矩阵K;斜体表示标量函数,比如泊松方程的解u;若带有下标则表示矢量的分量形式,比如u=uiei(i=1,2),下文为描述简单省去单位基矢量ei.∇()表示空间梯度运算,比如∇u=du/dxi为矢量;∇·() 表示散度运算,比如∇·u=dui/dxi为标量,这里重复下标采用了爱因斯坦求和约定;Δ()=∇·∇() 表示拉普拉斯运算,比如为标量.

2 虚单元法与有限元法的异同概述

在详细阐述虚单元法的基本理论前,需要对比有限元法并指出两种方法的异同,以便工程力学背景的同行对该方法有个直观了解.需要指出,该章意为对比基本概念,相关定义并不严格,更严格推导与解释详见下两章.

考虑Ω域内定义的标量场u(x),x ∈Ω,为得到该场的近似解uh(x),两种方法均需进行空间离散Ω ≈∪iEi,其中Ei为单元网格.

网格要求的不同: 有限元需要网格为特定节点数目的凸单元,如三角形单元、四边形单元等;虚单元法中单元可为任意节点数目,且对单元凸/凹性没有限制.其原因为: 有限元中需要进行单元映射,单元为凸才能保证映射矩阵的雅可比行列式为正,即单元体积恒正;而虚单元法无需进行单元映射,可直接采用空间离散单元的几何信息进行近似.

在单元E上,近似函数uh(x) 均表示成形函数与待解自由度的乘积,即

其中,nd为待解自由度的数目,ϕi(x) 为节点i的形函数,dofi(uh) 为节点i处的自由度.根据变分原理,刚度矩阵分量的一般形式为

式中推导采用了高斯公式,f,g,l为根据具体问题确定的函数.有限元确定刚度矩阵主要依赖式(2)中第二个等号左边直接进行体积积分,而虚单元主要依赖式(2)中第二个等号右边,需要在单元边界上进行线积分以及单元内部进行体积积分(通过自由度设定可避免).因而,形函数与自由度的设定均有不同.

形函数ϕi(x)形式的不同: 有限元中单元形函数为显式的多项式函数;虚单元法形函数除了多项式函数也可以为在单元边界上∂E为多项式、单元内部满足拉普拉斯运算后 Δϕi为多项式的函数,该函数可没有显式表达.如此,虚单元方法将形函数的形式进行了扩充.

自由度设定的不同: 有限元中待解量为函数在节点处的取值,即:dofi(uh)dofi(uh)=uh(xi);虚单元法中,在单元边界上的节点处的待解量为函数取值,但在单元内部的待解量为形函数的拉普拉斯运算与多项式的乘积.如此,可避免对没有显式表达的形函数进行体积分.

具体实施起来,有限元法由于形函数形式确定,可直接对式(2)中第二个等号左边进行体积积分,刚度矩阵因而仅有一项.需要指出,在特定单元与积分格式会造成“沙漏”,需要进行积分格式的调整或者增加稳定项.而对于虚单元方法需进行线积分以及体积分,由于形函数在边界上为多项式,因而线积分容易求得;通过自由度设定可避免进行单个形函数的体积分,但对于牵扯两个形函数ϕi,ϕj的积分便需要引入映射操作符 Π 来使得 Πϕi为多项式,如此

式中,为协调刚度矩阵分量,为稳定刚度矩阵分量,其具体形式可根据收敛准则给出.

刚度矩阵的不同: 有限元的刚度矩阵可直接通过体积积分获得,仅含有一项;虚单元法的刚度矩阵包含协调刚度矩阵和稳定刚度矩阵,其中稳定刚度矩阵通过收敛准则确定.

至此,通过对比有限元中的概念,简述了虚单元法的基本理论,第三章从数学角度探讨虚单元法的完备性;第四章从应用角度介绍虚单元方法在弹性问题中的应用流程,若对数学完备不感兴趣的读者,可直接跳至第四章.

3 虚单元法的基本原理与其在泊松方程中的应用

为解释虚单元法的基本理论,特别是数学上的完备性,需要补充部分基本理论,涉及到函数空间、离散方程推导的基本流程、解的存在唯一性与收敛条件、以及多项式函数空间与虚单元法特定的近似函数形式等的讨论.为便于理解且从信息的全面性考虑,这里作一些基本介绍.

3.1 函数空间

实数的集合称为实数集,类似的,函数的集合称为函数空间.比如函数u(x) 和v(x) 隶属于同一个函数空间(为表述更简洁,下文中省去自变量x),若在作用域Ω内定义了两者的内积

则称这样的函数空间为内积空间.对于内积空间的任意函数u,进一步定义一个实数‖u‖Lp(Ω)表示其大小,称为函数的模

其中,p为整数.若‖u‖Lp(Ω)＜∞,则称 Lebesgue 空间(下标L的来源).在本文讨论中仅涉及L2空间(即p=2,并略去了作用域Ω的符号).进一步,若函数及其导数满足

其中,k为整数,则称为 Hilbert 空间(符号H的来源).在本文讨论中仅涉及H1空间(即k=1 ).因为考虑了函数的导数,需要将函数的模进一步细化为半模(单竖线表示)与模(双竖线表示),分别定义为

其中,下标表示最高求导阶次.

3.2 离散化方程组(泊松方程及其离散化)

无论何种数值方法最终都将变为求解线性方程组,因而本小节将以简单的泊松方程为例,简述推导离散方程组的基本流程,以便后续的讨论.

记空间域Ω的边界为∂Ω,考虑

将包含真实物理解的函数空间记为V,则方程(8)对应的弱形式可通过分部积分并考虑边界条件推得,即寻找u ∈V,使得

其中,a(u,v),(f,v) 分别称为双线性格式与线性格式,与式(4)相对应,分别定义为

将包含有限个元素的近似函数空间记为Vh ⊂V,待解问题变为: 寻找uh ∈Vh,使得

需要指出,由于Vh ⊂V,对于∀uh,vh ∈Vh,除了在Vh空间内定义双线性格式ah(uh,vh),还可以在V空间内定义相应的a(uh,vh).为避免引入过多的数学定义,而不能专注于问题本身,在下文不引起混淆的情况下,把ah中的下标h略去①h 一般代表空间离散尺寸,考虑到单元离散后必对应近似函数空间,为避免引入过多符号,直接采用h符号表示近似空间,且在定义变量时是上标、在定义运算时是下标.,但在考虑解的收敛性时需要加以区分.

进一步,将Vh中线性无关的基函数记为ϕi(i=1,2,···,n) (考虑到Vh的定义,n为有限值),则对于函数uh,vh ∈Vh,将有如下线性组合形式

其中,Uj,Vi为实数.将式(12)代入到式(11)中,并考虑到vh的任意性,则有

由于双线性格式的对称性,记Kij=a(ϕi,ϕj),Fi=(f,ϕi),若采用矩阵写法K=(Kij),U=(Uj),F=(Fi),则式(13)可写为

再进一步,将空间域Ω离散为有限个空间单元,Ω ≈Ωh ≡∪iEi,此时,作用于整个空间域Ω上的函数空间Vh将由作用在单元E上的函数空间Vh|E集合构成,对应着整体刚度矩阵与右端项将由单元内积分(亦可以理解为全域积分,但在单元外部的函数取值为零)组装得到,即

其中,∑为组装操作,由于不易产生歧义,直接采用了相加符号.

至此,上述理论是有限元法与虚单元法共同的理论基础.

3.3 存在性、唯一性和收敛条件

对于数值求解,自然地需要考虑式(9)、式(11)是否存在解且解是否唯一;当将计算域进行空间离散后(记空间离散特征尺寸为h),还需要考虑收敛性,即h→0 时,‖u-uh‖→0 是否成立.分别讨论如下:

(1)解的存在性与唯一性

对于式(9),真实解所在的空间为在边界取值为零的一阶 Hilbert 空间,即V=H01,下标 0 表示函数在边界处取值为零.如果存在C ＞0,使得

称双线性格式a(·,·) 为连续的,类似地,连续性需要|(f,v)|≤||f||·||v||m.进一步如果存在α ＞0,使得

称a(·,·) 为强制的(coercive) 或者Hm-椭圆的(elliptic).可以证明强制的双线性格式对应的方程存在且具有唯一解,即著名的 Lax-Milgram 定理.需要指出,式(16)、式(17)中不等式右端项采用模度量或半模度量是等价的.对于式(9),可以证明

因而,式(9)存在唯一解.上述理论与相关证明,可参阅 (Braess 2007,Brenner and Scott 2008).对于式(11),其解的存在性亦需要验证ah(uh,vh) 是否是强制的,即: 是否满足式(16)和式(17),相关讨论将在下一节展开.

(2)解的收敛性

判断解是否收敛,需要验证协调性(consistency)和稳定性(stability)(Ralston and Rabinowitz 2001).协调性指当h→0 时,式(14)的解满足微分方程与边界条件,即趋向真实解;稳定性指式(14)有解,即:K非奇异.因而可以通过给定一些简单的真实解,来判断方程的解是否满足协调性与稳定性,即进行分片测试(patch test) (Taylor et al.1986).

3.4 多项式函数空间与虚单元法特定的近似函数形式

3.4.1 多项式函数空间

函数空间中,最直观的为多项式函数构成的空间.记Pk为最高阶次不超过k的多项式函数构成的空间,Pk为对应的矢量函数空间.

(1) 在一维问题中,记ξ为无量纲坐标,则有

Pk共有k+1 个线性无关函数,此外,记P-1=0;

(2) 在二维几何空间的标量问题中,记ξ,η为无量纲坐标,则有

Pk共有 (k+2)(k+1)/2 个线性无关的函数;

(3) 在二维几何空间的矢量问题中,则有

共有 (k+2)(k+1) 个线性无关的矢量函数.对于小变形下的弹性问题,将P1变为

如此,并不改变函数空间内线性无关多项式矢量函数的数目,但利用前三项对位移矢量u进行近似时,可对应着应变为零的情况,即发生刚体位移.比如,若u=C(-η,ξ)T,其中,C为常数,对应的应变为ϵ=[∇u+(∇u)T]/2=0.如此调整,在生成刚度矩阵时,仅需少量操作,便可避免刚性位移,具体讨论将在下一章给出.

3.4.2 虚单元法特定的近似函数形式

首先考虑如下问题,在二维多边形E域内,标量函数u ∈H1满足

即函数u进行拉普拉斯运算后得到最高阶次不超过k-2 的多项式,并且在多边形边界上u为连续的且不超过k次的多项式,那么函数u一定是多项式么？

通过反证法,显然答案是否定的.下文将介绍虚单元法的函数空间由满足上述条件的函数构成(多项式函数以及非多项式函数),因为所对应的非多项式没有显式形式(或者不易求出),而被称为虚单元法.与之对应,有限元法的函数空间仅为多项式空间,因而,不同方法的核心区别是近似函数空间Vh的差异.

3.5 虚单元法在泊松方程中的应用

由于泊松方程待解量为标量函数,相对简单,本节将以泊松方程为例,介绍虚单元法的基本思路与流程.相比于下一节介绍的弹性方程,本节略侧重于数学理论,以保证方法的数学完备.本节按照逆序的思路,首先给出式(11)存在唯一解以及收敛的条件,其次引出虚单元方法来实现上述条件,最后给出具体的实施细节.

3.5.1 目标

空间离散后,对于任意的Vh ⊂V,假定

即整体的双线性格式可由单元内的双线性格式组装得到,如此仅需对单元内的双线性格式进行讨论.对于任意的多边形单元E,相应的空间离散特征尺度h,以及定义在单元上的近似空间Vh|E,如果存在一个正整数k≥1,使得Pk(E)⊂Vh|E,且有

(1) 对于∀p ∈Pk(E) 以及∀vh ∈Vh|E,满足

(2) 存在两个不依赖于E以及h的正数α*和α*,对于∀vh ∈Vh|E,满足

此时,可以证明式(11)存在唯一解,且有足够的精度,保证

条件(1)保证真实解为最高阶次不超过k的多项式时,近似解可逼近真实解,即协调性条件;条件(2)保证ah的度量与a相同,此时能保证生成的刚度矩阵是满秩的,称为稳定性条件.上述证明可参考 (Beirão Da Veiga et al.2013a).

为满足条件(1),首先定义函数空间VE,k为式(19)的解所构成的空间.其次,定义映射操作符 ΠEk:VE,k→Pk(E)⊂VE,k,满足: 对于∀v ∈VE,k有

至此,对于∀u,v ∈Vh|E,令aEh(u,v)=aE(ΠEk u,ΠEk v),则能满足条件(1),即式(21).

为满足条件(2),假定一个对称正定的双线性格式SE(u,v),对于∀v ∈VE,k满足 ΠEk v=0 的前提下,满足

其中,正数c0,c1与单元和单元尺度无关.根据式(25),对于∀u ∈VE,k,则有:令v=u-ΠEk u,代入式(26),则对于∀u ∈VE,k,有

至此,令

此时,考虑式(25),易证明式(28)的定义满足条件(1),即式(21);考虑到式(26)以及aE(ΠEk v,ΠEk v)+aE(v-ΠEk v,v-ΠEk v)=aE(v,v),亦可证明式(28)满足条件(2),即式(22),其中α*=max{1,c1},α*=min{1,c0}.

基于上述讨论,虚单元法的总体思路可以简述为: 在满足式(19)的VE,k空间中,定义满足式(24)的映射操作符 ΠEk:VE,k→Pk(E)⊂VE,k,并选用满足式(27)的SE来定义双线性格式(28).

3.5.2 实施方案

虚单元法的具体实施中,包括三个核心细节: 根据VE,k设定自由度;构造SE;以及计算线性格式(离散方程组的右端项).现分别予以讨论.

(1)根据 VE,k设定自由度

考虑多边形单元E,记其形心为xE,边数为n,特征尺寸为hE.

对于式(19)的第一个条件,对于多边形的任意边,需该边上k+1 个点的函数取值则可确定最高阶次为k的多项式函数,这些点可由该边的 2 个顶点,以及边内k-1 个平均分布的点(或位于边内线积分点的 位置)构成(当k=1 时,每个多项式函数均为线性函数,则无需边内节点).由于多边形的顶点被 2 条边共用,因而,共需要n+(k-1)n=nk个自由度即可确定∂E上连续的多项式函数.

对于式(19)的第二个条件,对于任意给定的f,则式(19)的解是唯一的.因而若满足f ∈Pk-2(E),由于Pk-2共有k(k-1)/2 个线性无关的多项式函数(预备知识 3.4.1 中已予以讨论),需要增加k(k-1)/2个自由度来保证式(19)的解的函数空间的维度,这些自由度设定在单元内部点上.

因而,共需Ndof=dimVE,k=nk+k(k-1)/2 个自由度.图1表示了五边形(n=5 )单元在k=2时,自由度的设定位置.

图1 k=2的单元自由度设定位置(E:单元； e:单元边界)

在确定完自由度位置后,现在确定自由度的取值.将节点与单元边界上的设定自由度的点的集合记为nf,单元内部设定自由度的点的集合记为nb,并记xi处自由度取值为 dofi(vh).

若xi ∈nf,则该点处自由度取值为近似函数的取值,即: dofi(vh)=vh(xi);

若xi ∈nb,则该点处自由度取值需要考虑双线性格式的定义,即,对于∀uh,vh ∈VE,k

上式推导中采用了分部积分以及高斯定理,n为单元的外法向,且有:n·∇vh=∂vh/∂n.由于VE,k为满足式(19)的解构成的空间,因而,在∂E上vh,∂uh/∂n为多项式,可由xi ∈nf上设定的自由度确定(可设定边内自由度的点位于线积分点的位置,从而可以精确计算上式右端第一项的单元边界积分);单元内部积分仅包含右端第二项,而考虑到式(19)中 Δuh ∈Pk-2(E),因而把xi ∈nb处的自由度取值设定为

其中,AE为单元面积.需要指出,一个xi ∈nb对应一个p,Pk-2共有k(k-1)/2 项,则共有k(k-1)/2个单元内部点设定了自由度.如此,可以避免体积分来确定刚度矩阵,且避免了确定uh在单元内的取值.因而,虚单元的一个优势为: 通过设定合理的自由度取值,在规避求解形函数在单元内部的取值的同时,避免进行体积分来确定刚度矩阵.

对自由度的取值进一步讨论如下.

① 式(30)中p为无量纲的多项式,因而式(30)的量纲与vh的量纲相同.特别的,当p=1时,该自由度对应着vh在单元内部的平均值.在单元E内的无量纲度量可以为

② 与有限元类似,每个自由度对应着一个基函数ϕi(x),使得

此时,亦有:dofi(ϕj)=δij.由于vh(x)无显式表达式,则ϕi(x)亦无显式表达式.

③ 考虑到式(28),与上述自由度对应的局部刚度矩阵为

(2)稳定项SE的构造

式(33)右端第一项生成的刚度矩阵并不满秩,因而需要增加稳定项.考虑到稳定性需要满足式(27),最直观的设定为

需要指出,在上述推导中,ϕi和 ΠEk ϕi被表示成基函数ϕr的线性组合,即:对于 ΠEk ϕi同理.在合理的hE的设定下,可以保证aE(ϕr,ϕr)=O(1),因而可以设定

来满足式(27).需要指出稳定项的设定有多种形式,更多的讨论可参考 (Beirão da Veiga et al.2017b).

(3)离散方程组右端项的构造

考虑到单元离散,有

根据k的取值不同,fh的定义与 (fh,ϕi)E的构造可分为三种情况

① 对于k=1,fh为f在单元E内的平均值,即

此时

上述推导中,考虑了ϕi(xv)=δiv.

② 对于k=2,类似式(24),需要定义线性格式的映射操作符 Π0k:L2(E)→Pk(E),对于∀w ∈L2(E),满足

此时,fh=Π0kf,类似可定义 Π0kϕi.如此,反复利用式(39),可构造单元右端项为

其中,Π0k可以通过单元内部自由度式(30)求出,更多实施细节可参考 4.2.4 节.

③ 对于斜体>k＞2,类似的,定义映射操作符Π0k-2:L2(E)→Pk-2(E)fh=Π0k-2f,以及Π0k-2ϕi.构造单元右端项为

其中,Π0k-2的计算可参考 4.2.4 节.

如此,可证明(Beirão Da Veiga et al.2013a)

至此,虚单元求解泊松问题时单元内部的刚度矩阵为式(33),右端项根据k的取值在式 (38)、式(40)、式(41)三者中选择合适的公式计算,进行组装后,便可进行近似求解.

4 虚单元法处理线弹性问题

上一章以泊松方程为例,介绍了虚单元法的核心理论,侧重解释为什么引入一些处理,比如自由度的设定、多项式映射等,并强调数学完备.本章以弹性问题为例,通过对比有限元,具体介绍虚单元法的应用过程.为独立于上一章节,本章直接从应用角度理解虚单元,部分概念略有重复,但更侧重物理意义.

4.1 弹性问题描述

以二维连续弹性体为例,如图2所示,物体Ω在准静态、小变形条件下的平衡条件为

图2 二维弹性边值问题

其中,σ为柯西应力张量,f是为单位体积的体力.物体边界记为∂Ω,边界条件为

其中,n为边界法向,为指定应力,且有∂Ω=∂Ωt ∪∂Ωu以及∂Ωt ∩∂Ωu=Ø.本构方程为σ=Cϵ,其中C为四阶材料弹性常数,ϵ为应变张量,定义为位移梯度的对称部分,即ϵ=[∇u+(∇u)T]/2.

控制方程的弱形式为: 寻找u ∈V使得①计算数学的微元符号常采用 dx(如上章所述),本章更侧重于力学解释,因而体积微元为 dΩ,界面微元为d∂Ω

其中,v为试矢量,V为包含真实解的矢量函数空间,其中每个分量隶属于H1(Ω).通过分部积分,考虑应力的对称性以及代入边界条件,式(45)可写成如下形式

其中,a(u,v) 与L(v) 分别为弹性问题的双线性与线性形式,定义为

同样,记含有限个基函数的近似空间为Vh,相应地,近似的双线性格式ah(uh,vh) 和线性格式Lh(vh)形式与式(47)、式(48)类似,这里不再重复.

4.2 离散方程

如前所述,虚单元法与有限元法类似,也需要将计算域近似为有限个单元构成的区域,即:Ω ≈Ωh ≡∪iEi.不同于有限元中单元需要为凸多边形来保证雅可比矩阵的行列式为正,在虚单元中单元Ei可以为任意形状多边形(凸、凹).本节将重复虚单元法的映射操作符与自由度设定,并给出刚度矩阵与右端项的具体矩阵表达.

4.2.1 映射操作符

对于矢量函数空间,同样定义映射操作符 ΠEk,将定义在单元E上的近似函数空间Vh(E) 映射到最高阶次不超过k的多项式空间Pk(E),ΠEk:Vh(E)→Pk(E),下文中在不产生歧义的情况下,Π≡ΠEk.式(24)亦称为正交条件,即:

由于双线性形式在某种程度上表示了单元内部的能量,因而正交条件表示了以能量为度量时,uh与 Πuh的误差并不能通过不超过k次的多项式空间进行度量.

如此,利用式(49),定义在单元上的双线性格式为:uh,vh ∈Vh

右端第一项为映射操作符对uh,vh操作后在单元E内的积分,第二项为进行映射操作后没能考虑的能量,又被称为稳定项,需要能够确保aE(uh,vh) 为强制的,保证存在唯一解.

4.2.2 自由度设定

在弹性问题中,近似函数空间中的分量亦需满足式(19).对于nv边形单元,且边界上近似函数为最高阶次不超过k的多项式时,需在

●nv个节点布置函数取值的自由度;

● 每条边界上存在k-1 个边界内部节点,在其上布置函数取值的自由度;

● 在单元内部布置vh与最高阶次不超过k-2 的多项式内积的自由度(矩形式)

如第三章关于自由度的讨论所述,对于标量问题共需nd=nvk+k(k-1)/2 个自由度,因而对于二维矢量问题,共需 2nd个自由度.

在一个维度的函数空间内定义操作符 dofi(vh),i=1,2,···,nd,表示vh在第i个自由度的取值,ϕi(x) 为第i个自由度对应的基函数,如前所述,此时有

当设定自由度的点处均有两个自由度(二维矢量)时,ϕ表示基函数的向量表达,有ϕ1=[ϕ1,0]T,ϕ2=[0,ϕ1]T,···,ϕ2nd-1=[ϕnd,0]T,ϕ2nd=[0,ϕnd]T,且有

需要指出,尽管虚单元法的空间离散格式(53)与有限元法相同,但由于节点形函数并无显示表达,单元内部任意位置处的位移并不能直接获取.但网格节点的位移是待解量,可通过节点位移后处理插值获取位移场、应力场等信息.

4.2.3 单元刚度矩阵

首先回顾传统有限元中单元刚度矩阵的构造.二维弹性问题的自由度为节点位移

其中,下标 1,2 分别表示两个维度,上标 1,2,···,nd为单元节点编号.形函数矩阵为

应变为:ϵ=BuE=∂NuE,其中∂在二维 Voigt 表示下,为

从式(47)出发,并考虑vh的任意性,可以推得有限单元法中单元刚度矩阵为

其中,C是材料弹性常数矩阵.

虚单元法中的刚度矩阵构造过程与有限元法类似,但如式(50)所示,需包含两部分

其中,KcE为协调刚度矩阵(与映射操作符的协调性相关),KsE为稳定刚度矩阵(与稳定项相关).

(1)协调刚度矩阵

最高阶次不超过k的多项式基矢量集为:{pα}α=1,2,···,nk,且pα ∈Pk(E),如 3.4.1 中所述,nk=(k+2)(k+1).映射操作符 Π 对式(53)中的矢量基函数ϕi的操作定义为多项式基的线性组合,即

式中,si,β为常数.按照正交条件式(49),需要满足

令α=1,2,···,nk,则可以得到如下矩阵形式

其中

以及

代入所有的矢量基函数ϕi,i=1,2,···,2nd,bi组装成矩阵,映射系数si组装成矩阵Π˜,即

此时,由式(61),得

该条件可视为vh和 Πvh在pα(α=1,2,3) 度量下的平均值相等.如此,矩阵和矩阵G需要做如下修正 (Mengolini et al.2019)

对于矩阵G,前三行外的其他元素均可通过数值积分直接求得,从而可完全确定矩阵G的各个元素.

其中,为单元边界∂E的外法向量.

对于式(70)右端第一项,如前所述,pα ∈Pk(E),则:∇·(C ϵ(pα))∈Pk-2(E),进一步考虑线性组合

其中,dα,β为有量纲的系数.考虑到内部自由度的定义式(51),可以看出,式(70)右端第一项与单元内部自由度相关,且有 dof2nvk+β(ϕi)=δ2nvk+β,i.

对于式(70)右端第二项,单元边界积分为各条边上线积分的贡献之和,每条边上的线积分按如下数值积分计算

其中,ej为单元边界∂E的第j条边,是该边外法线矢量.可将边界与顶点设定的自由度置于每条边界上的数值积分点ξr(其相对应的权重为wr),如此,可提高计算精度.

作为总结,矩阵G的前三行按照式(69)求得,剩余行可直接数值积分得到;矩阵的前三行按照式(68)求得,剩余行通过式(70)～式(72)求得;矩阵按照式(66)求得.进一步,将矩阵G作一变换,令其前三行所有元素设为 0,其他行所有元素不变,记此矩阵为.

如此,协调刚度矩阵的分量形式可表示为

相应地,协调刚度矩阵为

(2)稳定刚度矩阵

与有限元法不同,虚单元法中需要增加稳定刚度矩阵.

首先,定义矩阵D2nd×nk为多项式基pα在自由度设定位置上的自由度取值,即

对应的矩阵形式为

考虑到自由度取值可为函数本身取值,亦可为函数与多项式乘积的积分(矩形式).由此,可以确定矩阵D的各个元素:

① 对于边界自由度: 1 ≤i≤2nvk,dofi(pα) 定义为基矢量pα在第i个自由度所属节点的矢量分量值.

②对于单元内部自由度: 2nvk+1 ≤i≤2nd,考虑到单元内部自由度的定义式(51)

按照数值积分即可求得矩阵D的剩余元素.

进一步,定义

与式(35)相对应,二维弹性问题的稳定项可以取

将式(78)中Πij的定义,以及单位矩阵Iij=dofi(ϕj) 的定义代入式(79),则稳定项的矩阵形式为

可以证明(Beirão da Veiga et al.2014,2013b),若存在一个因子σE ＞0,满足

其中,常数σ*,σ*均不依赖于单元特征尺寸h,则稳定项式(80)乘以该因子σE依然能够保证计算结果收敛.通常取σE=τh·tr(kcE)＞0,从而得到稳定刚度矩阵

式中,τh为稳定系数,对于弹性问题,τh=0.5,更多选择可参考 (Gain et al.2014);tr(·) 为矩阵的迹,引入此项可保证单元内部的能量相对于单元尺寸和材料弹性常数为相对正确的比例(Artioli et al.2017a).

最后,需要指出矩阵G恰好可以表示成矩阵和矩阵D的乘积:G=.因此,如果在计算过程中首先求出了矩阵和矩阵D,就无需再通过式(63)进行求解.当然,按式(63)求出矩阵G,可帮助检查两种方式的求解结果是否相同.

4.2.4 右端项

如式(48),单元内部的线性格式为

由于近似函数在单元边界上为多项式,式中第二项无需特别处理,但对于第一项,近似函数在单元内部无显式表达,因此需要进行近似处理.在不产生歧义的情况下,本节下述右端项均指与体积力f对应的右端项.

与泊松方程类似,对于弹性问题,右端项也需要进行合理近似

根据k的取值不同,fh的定义与 (fh,ϕi)E的构造可分为三种情况

(1) 对于k=1,fh为f在单元E内的平均值,即

(2) 对于k=2,类似式(49)定义的双线性格式的正交条件,首先定义线性格式的映射操作符Π0k:VE,k→Pk(E),对于∀ϕi ∈VE,k,满足

其中

式中,si,β为常数.下面将讨论如何计算si,β.

将式(88)代入式(87),可得方程组

将其写成矩阵形式

其中

以及

代入所有矢量基函数ϕi,i=1,2,···,2nd,将qi0组装成矩阵Q0,将映射系数si0组装成矩阵Πk0,即

方程组(90)可整理为

至此,系数矩阵Π0k可由式(95)求得.下面讨论H0,Q0的求解.对于矩阵可通过已知的基矢量pα(α=1,2,···,nk) 完全确定.对于Q0的前nk-2行

其中

考虑到单元内部自由度的定义式(51),可将基函数在单元内部的自由度写成

因此,矩阵Q中的元素均可以确定

对于Q0的第nk-2+1 行至第nk行(需要特别注意nk-2和nk分别是Pk-2和Pk空间的线性无关基函数的个数,两者并非相差 2),需要计算基函数ϕi和基矢量pα在单元内部的积分,然而ϕi没有显式表达,考虑到式(60)与式(89),由于已按照式(66)求得,可直接采用的第nk-2+1行至第nk行进行替换,即

如此,矩阵Q0nk×2nd的计算公式如下

通过式(95)求出Π0k,反复利用式(87),可构造右端项为

(3) 对于k ＞2,类似的,定义映射操作符 Π0k-2:VE,k→Pk-2(E),对于∀ϕi ∈VE,k,满足

其中

式中,ri,β为常数.将式(104)代入式(103),则有方程组

其矩阵形式为

代入所有矢量基函数ϕi,i=1,2,···,2nd,将qi组装成矩阵Q,恰好为式(97),将映射系数ri组装成矩阵Π0k-2,即

方程组(106)可以整理为

其中,矩阵Hnk-2×nk-2可通过已知的基矢量pα(α=1,2,···,nk-2) 完全确定;矩阵Qnk-2×2nd可通过式(97)、式(99)确定.

此时,通过式(109)求出操作符 Π0k-2的矩阵表达式Π0k-2,可构造右端项

如上所述,可以在三种不同情况下求得基函数ϕi所对应的离散方程右端项.如此,单元内与体积力f对应的右端项为

4.3 计算流程

本小节给出单元矩阵与右端项的计算流程,如统一输入量为: 内插阶数k,多边形单元E,单元节点坐标X=[x1,x2,···,xnv],积分点坐标为ξ,权重为γ,材料弹性常数矩阵C,体积力f;输出量为: 单元矩阵kE,单元右端项rE(包含与体积力f对应的右端项rfE).

Π ←D ˜ΠVh(E)(10) //映射到 空间的操作符式(78)(11) //计算刚度矩阵kc E ←˜ΠT ˜G ˜Π(12) //协调刚度矩阵式(74)ks E ←τhtr(kc E)(I -Π)T(I -Π)(13) //稳定刚度矩阵式(82)kE ←kc E +ks (14)(15) // 计算右端项rf E,rE ←02nd×1(16) //初始化k=1 E(17) if: then Xfrf ErE(18) 根据 和 计算,//右端项式(86)、式(83)(19) end(20) else if then H0αβ ←∫E pα·pβ dEH0nk×nk(21) // 矩阵式(92)Qαi ←∫E pα·ϕi dEQnk-2×2nd k=2(22) // 矩阵式(97)Q0 ←Q,H0,G, ˜BQ0(23) // 矩阵式(101)Π0k ←(H0)-1 Q0nk×2ndΠ0k(24) // 的 矩阵式(95)XfΠ0krf ErE(25) 根据,,计算,//右端项式(102)、式(83)(26) end(27) else Hαβ ←∫E pα·pβ dEHnk-2×nk-2 nk×2nd(28) // 矩阵式(105)Qαi ←∫E pα·ϕi dEQnk-2×2nd (29) // 矩阵式(97)Π0 k-2 ←H-1Qnk-2×2ndΠ0k-2(30) // 的 矩阵式(109)k-2rf ErE(31) 根据,,计算,//右端项式(110)、式(83)Xf Π0(32) end

4.4 应用实例

如前所述,虚单元法可处理任意节点数目的多边形且对单元凸凹性没有限制.为做出直观展示,本小节以含有“力”字单元的悬臂梁模型为例,讨论收敛性.本小节结果为阶次为k=1 的情况.

考虑到悬臂梁的位移场存在理论解(Timoshenko and Goodier 1951).如图3所示,在平面应变状态下单位厚度的悬臂梁左端,根据理论解施加位移边界条件,右端施加沿 y 轴分布的剪应力:

图3 悬臂梁模型

此时,水平向与垂向的位移为

为展示虚单元的优势,采用含“力”字单元的网格,分别计算四种网格密度,如图4(a) 所示.图4(b) 为计算得到的水平向位移场,可以看出位移场连续且光滑.

图4 网格与水平向位移场.(a)网格,(b)水平向位移

图5 收敛曲线

图6展示了不同纵向剖面(如图3所示x=2.0 m,x=4.0 m,x=7.2 m)的水平向位移与理论解的比较,可以看出虚单元方法在不均匀的多边形网格下,随着网格的加密能够趋于收敛.

图6 不同纵向剖面的水平向位移与理论解比较

5 虚单元方法的深入应用

上一章系统梳理了理解虚单元基本理论的预备知识;侧重数学推导,阐述了虚单元在泊松方程的应用;侧重物理解释,介绍了虚单元在二维弹性问题的应用过程.如前言所述,虚单元方法的相关研究已从计算数学领域逐渐拓展到工程科学/力学领域.本章对虚单元方法在计算力学领域的研究进展进行归纳与总结,包括材料与几何非线性、接触与界面问题、断裂问题等.为避免内容过度重合,每部分仅引用近 5年的代表性工作.

5.1 材料与几何非线性

本节首先梳理小变形的相关工作,侧重于单元曲直、阶数的研究进展,其次按照准静态问题与动力学问题进行区分,讨论材料与几何非线性问题的研究进展.需要指出,晶体塑性的相关工作在下一节接触与界面问题中予以讨论.

5.1.1 小变形问题

对于低阶直线单元,Artioli et al.(2017a) 提出虚单元方法求解线弹性问题的标准化代码实现流程,详细介绍了协调项、稳定项和载荷右端项的构造过程,结果表明,同阶的虚单元方法与有限元方法具有相近的准确性和相同的收敛阶数,且对网格畸变几乎不敏感.此后,Artioli et al.(2017b) 将其工作扩展到非线性材料中,研究了黏弹性模型、含各向同性与运动硬化的 Mises 塑性模型、以及记忆合金本构模型.Beirão da Veiga et al.(2015) 对小变形条件下,虚单元的非线性弹性本构与 Mises 屈服准则的弹塑性本构的应用过程进行了总结.Gain et al.(2014) 将虚单元法应用到三维的小变形弹性问题中,并对比了不同类型网格的收敛性.

针对单元阶数,Beirão da Veiga et al.(2017a) 提出了高阶单元方法求解对流反应方程,Beirão da Veiga et al.(2016)类比高阶 Serendipity 有限元方法,提出了高阶 Serendipity 虚单元方法求解弹性问题.与高阶有限元方法类似,高阶虚单元法对网格畸变问题具有更高的鲁棒性,并且由于增加了自由度,计算精度更高.为避免离散曲线边界时需要“以直代曲”,针对单元曲直,Beirão da Veiga et al.(2019) 将单元边界的描述转化为自然坐标,提出了曲边单元,并对椭圆方程进行了求解.Artioli et al.(2020) 将其方法应用到二维线弹性问题中,但曲线形状不能任意.Wriggers et al.(2020) 提出可将参考构型的单元(直线边)映射到初始构型(曲线边)中,实现任意形状的单元离散,对于映射操作可采用等参映射或者 NURBS 映射.而后,Wriggers et al.(2021b) 将该方法扩展到高阶单元中.

5.1.2 准静态问题的几何与材料非线性

对于弹性条件下的几何非线性问题,Wriggers et al.(2017) 通过对应变能的讨论,提出了有限变形的稳定项构建方法,通过可压缩和不可压缩材料的模拟,验证了其方法的收敛性和鲁棒性,如图7所示.同一时期,Chi et al.(2017) 通过构造局部位移空间准确求得体积变化,根据切线模量的迹构建稳定项,讨论了混合单元(位移、压力)格式与纯位移格式的有限变形模型.图8是其模型对弹性多孔材料的模拟结果.针对不可压缩材料,Taylor-Hood 是混合单元格式,采用不同的阶数插值位移与压力,以保证 LBB 稳定条件,Wriggers et al.(2021a) 将 Taylor-Hood 单元引入到虚单元法中,构造了不可压缩材料的有限变形格式.

图7 有限变形下的冲压问题.(a)模型,(b)可压缩材料的冲压变形,(c)不可压缩材料冲压变形(修改自 (Wriggers et al.2017))

图8 弹性多孔材料的有限变形.(a)模型,(b)剪应变为 1.345 下的最大应变分布(修改自(Chi et al.2017))

若进一步考虑材料非线性,Wriggers and Hudobivnik (2017) 采用低阶单元从最小能量增量出发,构造了二维弹塑性本构下有限变形的虚单元格式,并对精度和收敛性进行了讨论.Hudobivnik et al.(2019) 将该方法扩展到三维问题中,并对比了不同收敛算法对结果的影响.图9是对冲压以及扭转过程的模拟结果.

图9 准静态大变形弹塑性问题的等效塑性应变分布.(a)冲压,(b)扭转 (修改自(Hudobivnik et al.2019))

需要指出,大部分虚单元方法在处理非线性问题时,均采用低阶单元,De Bellis et al.(2019)成功将高阶单元应用到非线性弹性问题中.

5.1.3 动力学问题的几何与材料非线性

对于动力学问题,Cihan et al.(2021b) 指出不同于刚度矩阵需要稳定项,虚单元方法的质量矩阵只需要进行映射操作即可,其中时间积分可按照隐式 Newmark 格式进行.图10是不考虑阻尼条件下,悬臂梁的振动模拟结果.

图10 悬臂梁振动问题.(a)模型,(b)垂向位移时间演化(修改自 (Cihan et al.2021b))

对于弹塑性的动力学问题,Cihan et al.(2021a) 按照 Hu-Washizu 泛函求极值的基本思路,提出了混合单元格式,以避免弹塑性不可压条件引入的体积自锁问题,其中的时间积分亦采用隐式 Newmark 格式.图11为泰勒杆与冲压过程中的塑性应变.若采用显式时间积分,Park et al.(2020) 对弹性不可压缩材料、Park et al.(2019) 对弹塑性材料分别进行了讨论,其中临界时间步长都是通过矩阵最大的特征值确定.

图11 弹塑性动力学问题的塑性应变.(a)泰勒杆,(b)冲压 (修改自 (Cihan et al.2021a))

上述讨论的都是各向同性材料,对于横观各向同性材料（比如纤维增强）,虚单元法也被逐步应用到小变形 (Wriggers et al.2018b,Reddy and van Huyssteen 2019) 与有限变形模拟中(Wriggers et al.2018a,van Huyssteen and Reddy 2021).

5.2 接触与界面问题

由于虚单元法对网格形状没有要求,因而在网格配对处理上可得到极大简化,在处理接触与界面问题上具有先天优势.本节分别从接触问题、非均匀材料的界面问题以及衍生的晶体塑性三个方面总结虚单元的研究进展.

5.2.1 接触问题

采用虚单元法处理接触问题,是该方法一个非常重要的应用.处理接触问题需要判断两个物体间的相对位置,对于有限元网格(含节点、网格边界线),如图12(a)～图12(c)所示,常用的方法有点-点、点-线、Mortar 型,其中点-点型最容易实现,但需要网格匹配.由于虚单元对单元节点数目没有要求,可任意增加节点,如图12(d)所示,可以将潜在接触面的节点向目标面做映射生成新的节点,采用点-点格式进行接触判断,如此可避免复杂的网格匹配、助于各体间独立建模.需要指出,虚单元法中增加了单元节点后,增加节点的单元形式会发生改变,待解自由度增加,而其他单元没有变化,这是与有限元节点插值最大的区别.

图12 判断两物体间相对位置的常用方法.(a)点-点,(b)点-线,(c)Mortar 型,(d)虚单元法插节点

对于两体接触,Wriggers et al.(2016) 对虚单元法在接触问题中的应用进行了梳理,分别采用拉格朗日乘子法和罚函数法施加约束,其结果表明虚单元法利用插点格式处理接触问题,精度高、操作简单.图13是采用不同方法对接触问题进行的分片测试,可以看出相对于有限元中的点-线格式,虚单元法求得的应力均匀,精度更高.Wriggers and Rust (2019) 将该方法扩展到大变形的摩擦接触问题中,图14为其模拟的赫兹接触过程.对于接触体具有曲线边界的情况,Aldakheel et al.(2020) 引入了曲边单元处理接触,其中插入的节点需要按照曲边单元进行处理.

图13 接触问题的分片测试.(a)虚单元法,(b)有限元点-线格式(修改自(Wriggers et al.2016))

图14 大变形条件下的赫兹接触问题(修改自 (Wriggers and Rust,2019))

颗粒材料中颗粒间的相互作用涉及多体接触问题,常采用离散元进行描述,其核心假定为颗粒为刚体.虚单元法在网格上的任意性以及接触处理的优势,使得快捷准确地捕捉变形颗粒间的相互作用变为可能.Gay Neto et al.(2021) 首次做出了尝试,采用完全的隐式时间积分、并且研究了刚度与质量矩阵参数的敏感性.图15模拟多面体颗粒集合的沙漏过程.该方法为颗粒材料研究提供一个新的思路与方向.

图15 颗粒材料沙漏过程(修改自 (Gay Neto et al.2021))

5.2.2 不考虑断裂的非均匀材料界面问题

对于非均匀材料,考虑两方面的研究: 材料界面与均质化.

关于材料界面问题,Lo Cascio et al.(2021) 采用虚单元法描述基材、边界元法描述内嵌材料,边界元获取的应力-位移方程转化为力-位移方程,组装到虚单元法的方程组中,实现两者的耦合.该方法有效地描述了非均匀材料的变形及损伤演化,图16为典型算例的模拟结果.

图16 虚单元法与边界元法模拟非均匀材料.(a)网格,(b)单轴拉伸条件下的应力分布(修改自 (Lo Cascio et al.2021))

多晶材料的均质化是进行多尺度建模的核心.对于有限元,增加晶粒的各向异性需增加大量的自由度,而虚单元法将任意形状的晶粒视为一个网格,自由度数目自然大大降低.Marino et al.(2019)对线性和非线性的均质化格式开展研究.若考虑电-磁-力耦合的多晶材料,Böhm et al.(2021b) 计算了不同晶格结构和不同程度的各向异性材料,发现虚单元方法在低自由度条件下计算的均质化性能都表现出很好的精度,如图17所示.

图17 多晶结构均质化.(a)虚单元法网格,(b)粗糙的四面体网格,(c)精细的四面体网格,(d) BaNiO3，中度各向异性，六方晶系晶胞(修改自 (Böhm et al.2021b))

5.2.3 晶体塑性

晶体塑性采用虚单元方法进行模拟具有巨大潜力.虚单元法中单元的任意性可以完美地契合晶粒结构,无需进行精细的网格剖分,但该部分工作尚属起步阶段.Böhm et al.(2021a) 采用虚单元模拟了 FCC 晶格结构受单向拉伸的滑移过程,其中映射操作在晶粒的表面进行以产生协调矩阵,稳定项需要将晶格内部划分为多个四面体.图18为单元分解过程与单滑的模拟结果.工业应用方面,Behrens et al.(2020)将基于虚单元的晶体塑性模型应用到轴承衬套的生产设计中,以描述多种材料连接处多晶材料的损伤与疲劳性质,图19为模拟结果.

图18 虚单元法模拟晶体塑性.(a)单个晶格网格分解为界面网格与内部四面体网格,(b)单轴拉伸条件下 FCC 晶格结构的剪应变(修改自 (Böhm et al.2021a))

图19 虚单元方法分析轴承衬套中不同材料连接区域的多晶塑性模型.(a)轴承衬套示意,(b)三种材料连接区域的代表性体积单元,(c) FCC 晶格在 (111) 滑移面上剪切应变率 的分布,(d)等效von Mises 应力分布(修改自 (Behrens et al.2020))

5.3 断裂问题

虚单元法中网格的任意性简化了空间离散的难度、并能随意增加/减少节点.虚单元法处理断裂问题可以从直接单元切割(强间断)以及与其他方法(相场、内聚力、扩展有限元等)的结合两个角度进行总结.

5.3.1 直接单元切割

针对直接单元切割,Hussein et al.(2019)提出了沿裂纹扩展方向增加新的单元节点并进行单元分割的裂纹扩展技术,对比不同形状的网格,验证其方法的有效性.图20为单元切割算法与单轴拉伸下的裂纹扩展.直接切割算法的优势是可直接确定裂纹面的位移间隔,在采用位移差求解应力强度因子时大大简化了程序设计,并且由于可任意增加节点,虚单元法在裂纹交叉、分叉、三维曲面裂纹的扩展模拟中具有潜力.

图20 虚单元法在断裂问题中的应用.(a)单元切割,(b)裂纹扩展(修改自 (Hussein et al.2019))

5.3.2 与其他方法结合模拟断裂

断裂问题一直是计算力学的热点,在虚单元法提出以前,适用于不同场景的描述断裂的方法已被广泛发展.

与相场法耦合.考虑到裂纹扩展方向可以容易从相场方法获得,Hussein et al.(2020) 将相场法与单元切割融合,采用虚单元法中单元切割表征裂纹扩展,未开裂部分的网格采用自适应方案,来求解相场.该方法为裂纹扩展模拟提供了高效、鲁棒性强的模拟方案.图21展示了该方法模拟的裂纹扩展过程.

图21 单元切割与自适性相场耦合模拟裂纹扩展.(a)自适应网格求解相场,(b)单元切割模拟裂纹扩展(修改自 (Hussein et al.2020))

与内聚力模型结合.Marfia et al.(2022)将虚单元法与内聚力模型相结合,按照单元域内平均的最大主应力确定裂纹方向,通过内聚力模型以及单元切割,模拟裂纹的起裂与扩展过程,如图22所示.进一步考虑非均匀材料的界面引起的断裂问题,Benedetto et al.(2018) 应用零厚度的界面单元模拟应力-裂纹张开过程,结果如图23所示.

图22 虚单元法与内聚力模型模拟 L 形状板的裂纹扩展(修改自 (Marfia et al.2022))

图23 虚单元法与内聚力模型模拟耦合.(a)接触脱离模拟,(b)非均匀材料三点弯曲梁模拟(修改自(Benedetto et al.2018))

与扩展有限元法结合.扩展有限元法是将近似函数空间进行扩充,采用间断函数以及奇异函数描述物理场的间断.Benvenuti et al.(2019) 将扩展有限元法与虚单元进行结合,对拉普拉斯方程的裂纹尖端进行了描述,其中满足拉普拉斯方程的间断函数以及奇异函数被扩充到基函数空间中,并推导了扩充映射操作符,以便将扩充的基函数空间映射到线性多项式与扩充函数上.图24所示为含裂纹的薄膜受 III 型载荷作用的模拟结果.Benvenuti 等 (2022)将扩展有限元法与虚单元结合的方法进一步应用到线弹性断裂问题中,着重讨论了应力强度因子的求解.

图24 虚单元法与扩展有限元结合模拟含裂纹的薄膜在 III 型载荷作用下变形.(a)100 个网格,(b)1600 个网格(修改自 (Benvenuti et al.2019))

5.4 多场耦合、拓扑优化与地下裂隙网络流体流动

本节总结了虚单元法在其他领域的应用,包括热-力耦合作用、拓扑优化、以及地下裂隙网络流体流动模拟等.

5.4.1 热-力耦合作用

对于热-力耦合问题,Dhanush and Natarajan (2019) 采用 Matlab 生成多边形网格以及Abaqus 的输入文件,将虚单元法成功应用到 Abaqus 中来处理小变形下的热-力耦合问题,并给出了详细的数据格式.进一步,Aldakheel et al.(2019b) 将虚单元法中有限变形下的弹塑性模型引入到热-力耦合问题中,如此,自由能包含弹性能、弹-热耦合部分、纯热部分、以及塑性势能.图25为三维钢螺栓成型过程模拟.

图25 三维钢螺栓成型的热力耦合过程.(a)～(f)等效塑性应变,(g)～(l)绝对温度场 T (Aldakheel et al.2019b)

5.4.2 拓扑优化

拓扑优化问题需要求解各种几何形状的边值问题,由于虚单元法放松了对网格形状的要求,简化了复杂形状的网格剖分难度,特别适用于优化问题.Gain et al.(2015) 对比了虚单元法与有限元法的优化结果,指出虚单元法只需要少量的单元,便可得到足够精度.图26展示了对悬臂梁的优化结果,六面体网格需要 50 000 多个网格,而非规则多面体仅需 10 000 个网格.对于复合材料的大变形优化问题,Zhang et al.(2020) 引入虚单元法提高计算效率,其中材料性能的确定是通过插值能量函数确定,而非插值材料参数,且该函数可以解决粗糙网格面临的单元畸变问题.进一步,对于纤维增强的软物质的拓扑优化问题,Zhang 等 (2021) 将基材与纤维参数作为两组设计参数,其中纤维的方向是预先设定的一组离散方向,并采用虚单元法求解有限变形的边值问题.图27为三点弯曲梁的优化结果.

图26 拓扑优化问题.(a)悬臂梁模型,(b)六面体网格,(c)非规则多面体网格(修改自 (Gain et al.2015))

图27 拓扑优化问题.(a)目标承载体,(b)优化结果(修改自 (Zhang et al.2021))

5.4.3 地下裂隙网络流体流动模拟

地下裂隙网络中流体的流动是地质力学中的重要问题.相对于裂隙流动,孔隙渗流过程一般忽略不计,但三维空间中二维裂隙平面相互交叉、独立,若采用传统有限元协调网格,无法对每个裂隙单独建模,必需整体考虑,如此网格生成变得十分困难.考虑到虚单元法对网格节点数目没有要求,因而处理连接处的悬挂节点等十分自然,可以大大简化网格生成的难度.Benedetto et al.(2016) 采用虚单元法对多种裂隙网络下的流动进行了尝试,如图28所示.此后,Benedetto et al.(2017) 将该方法扩展到混合单元与高阶单元.

图28 地下裂隙流体流动问题.(a)网格生成,(b)裂隙网络中水头分布(修改自 (Benedetto et al.2016))

6 结论与展望

正如在该综述的开头所强调的,通过有别于传统综述论文的行文方式,详述了虚单元方法的基本原理和最新理论进展,同时也展示了目前该方法在一些典型问题中所体现的有别于有限元方法的潜力: 包括它对单元凸凹性的处理,以及因此延伸到如悬挂节点、接触、裂纹扩展、多晶体变形等问题方面的应用.基于处理非凸性单元任意多边形单元的考虑,虚单元方法引入了和有限元方法迥异的处理方法,最为显著的是(1)虚单元法中的近似函数既包含了多项式函数,同时也引入了在单元域内无法显式表达的函数,这也是虚单元方法中“虚”字的由来;(2) 自由度取值在节点以及单元边界上为近似函数的取值,在单元内部为矩形式(近似函数与多项函数的乘积在单元内的积分),避免计算非多项式在单元内部的取值;(3)为保证计算的收敛性,虚单元法的刚度矩阵包含了协调矩阵和稳定矩阵两部分,对应于多项式的精确解和多项式未被考虑的部分,比有限元方法更为复杂.

虚单元法尚属于发展阶段,在几何与材料非线性问题、接触与界面问题、断裂问题等诸多方面具有发展潜力,尤其需要关注它在以下几个方面的应用.

(1) 裂纹扩展: 断裂问题一直是计算力学的热点,在虚单元法提出以前,适用于不同场景的描述断裂的方法已被广泛发展.当前虚单元法在这一类问题方面展现了非常好的兼容性.将相场法与单元切割融合,采用虚单元法中单元切割表征裂纹扩展 (Hussein et al.2020),可以为裂纹扩展模拟提供了高效、鲁棒性强的模拟方案.将虚单元法与内聚力模型相结合,通过内聚力模型以及虚单元法中的单元切割 (Marfia et al.2022),可高效模拟裂纹的起裂与扩展过程.将扩展有限元法与虚单元进行结合 (Benvenuti et al.2019,2022),可模拟诸多的线弹性断裂问题.

(2) 接触问题: 虚单元法对单元形状没有要求,可任意增加接触面节点,因而在处理接触与界面问题时,可避免网格配对带来的挑战,简化计算与单元划分问题.(Wriggers et al.2016) 对虚单元法在接触问题中的应用进行了梳理,分别采用拉格朗日乘子法和罚函数法施加约束,其结果表明虚单元法利用插点格式处理接触问题,精度高、操作简单.对涉及大变形的摩擦接触问题Wriggers and Rust (2019)和具有曲线接触边界的情况 (Aldakheel et al.2020),虚单元法也具有高效的处理方案.同样对于涉及多体接触的问题,虚单元法在网格上的任意性以及接触处理方面,也具超越离散元的优势,可快捷准确捕捉变形颗粒间的相互作用 (Gay Neto et al.2021).工程中涉及到类似的接触问题不胜枚举.

(3) 多晶体变形: 利用虚单元法来模拟晶体塑性具有巨大潜力.这是因为虚单元法中对多边形单元的包容将可以和晶粒结构完美地契合,一方面无需进行精细地网格剖分,可以确保具有晶体取向的晶粒的整体性,同时又可以保持晶粒内部不同区域之间的非局域作用,使得变形梯度、位错塞积等晶体变形特征能够得到妥善处理.尽管目前一些研究组已经采用虚单元方法开展了晶格塑性模拟方面的应用 (Böhm et al.2021a,Behrens et al.2020),该方面的工作还亟待开发,以充分利用虚单元法对非规则单元处理方面的优势.

以上是笔者从自身的科研背景和需求出发所展望的虚单元计算方法的发展潜力,因此它不可能是全面的.同时需要强调的是,虚单元方法在这些特定问题方面,具有超越传统有限元方法的某些优势,因此为饱受这些问题困扰的科研工作者提供了一种新的解决方案,但这并不意味着它能全面地替代有限元方法的成熟性、便利性以及历经长时间发展带来的高可靠性.它作为一个特定模块,结合相应的有限元方法,将丰富和提升计算力学的处理能力和便利性.通过对两者的融会贯通,促进发展出适应我国力学计算需求的新型算法与高性能软件.

致 谢作者们感谢国家自然科学基金委基础科学中心项目“非线性力学的多尺度问题研究”(资助号 11988102)支持和面上项目支持(资助号 12172368,刘传奇),刘传奇同时感谢中科院百人计划项目支持;作者们对意大利 Milano-Bicocca 大学 L.Beirão da Veiga 教授和美国加州大学N.Sukumar 教授关于论文算法富有成效的讨论深表谢意.