CN108763841B

CN108763841B - 一种基于对偶边界元和应变能优化分析的弹性断裂仿真方法

Info

Publication number: CN108763841B
Application number: CN201810815571.9A
Authority: CN
Inventors: 李帅; 崔柳; 杨晨; 于洋
Original assignee: Qingdao Research Institute Of Beihang University; Beihang University
Current assignee: Qingdao Research Institute Of Beihang University; Beihang University
Priority date: 2018-07-24
Filing date: 2018-07-24
Publication date: 2022-04-05
Anticipated expiration: 2038-07-24
Also published as: CN108763841A

Abstract

本发明提出一种基于对偶边界元和应变能优化分析的弹性断裂仿真方法，对有关三维线性弹性断裂的两个基本问题进行了处理。通过使用对偶边界元方法来表示含有重叠裂口表面的模型，并采用半解析方法来求解超奇异边界积分方程。通过引入接触力模型来表示来自接触载荷和裂纹前端应变能释放的连续属性。通过在断裂持续时间段内对裂纹的延展距离计算断裂因子来对裂口端的延展进行精确控制。

Description

一种基于对偶边界元和应变能优化分析的弹性断裂仿真方法

技术领域

本发明属于物理仿真技术领域，具体涉及基于对偶边界元和应变能优化分析的弹性断裂仿真方法。

背景技术

边界元方法作为一种强有力的数值计算工具已经在工程应用的诸多领域得到广泛应用。研究者通过将面向不同应用的基础解带入边界积分方程(Boundary IntegralEquations,BIEs)，使得边界元方法具有处理各种应用(特别是在无限场域和开放表面等问题上)。

在断裂仿真方面，对偶边界元方法(Dual Boundary Element Method， Dual-BEM)在裂口处一侧施加传统位移边界积分方程(Displacement BIEs)的基础上，针对裂口处的另一侧施加独立的受力边界积分方程(Traction BIEs) 来解决重叠裂口表面的问题。Dual-BEM只对模型的表面进行离散，并针对模型裂口处重叠表面两侧设定不同的边界积分方程来解决系统退化问题，在线性弹性断裂仿真应用中具有非常明显的优势，但是也必须要面对求解各种奇异边界积分的情况。

对模型断裂过程的建模包含几何方法和物理方法。刚体断裂的模拟可以通过预定义的方法或者动态计算应力来计算如何断裂。在模型上预定义的方法可以事先设计出一个断裂模式，这个模式可以是从某一断裂点开始，或者裂纹没有集中到某一断裂点，在断裂时会用预定义的断裂模型来代替，这样速度快但缺少真实性。而从物体的内应力动态计算断裂的方法，可以产生很真实的从某一点开始的很真实的断裂模式，但需要很大的计算量，在实时应用中不太适用。脆性断裂通常指硬体之间在裂纹张开和扩展时，只有很小的弹性变形。

三维断裂仿真的主要挑战在于裂纹延展过程中的几何方面的更新。一般地，断裂过程中的裂纹增长幅度一般会小于网格分辨率两个数量级，此外，考虑到裂口延展的自由性，对几何网格的更新带来了巨大的困难。更重要的是，网格更新所带来的复杂性会导致在仿真过程中频繁出现不连续性和奇异性，这对断裂仿真的数值求解方法提出了新的要求。如果采用体网格细分方法来表示裂口，受网格最小颗粒度和数值连续性的限制，不可能准确地符合实际的裂口，而这种病态的网格反映到物理系统中会导致系统的退化。有限元方法作为一种在较为常用的断裂建模方法，对此类问题的支持是十分无力的。对于裂口处的细分网格操作所导致的大规模网格，有限元方法的效率会出现大幅度的下降，低质量的网格甚至会影响有限元方法的收敛率。不仅如此，传统有限元方法的形函数是不能表示奇异性结果的，如何表示奇异的应力是有限元方法面临的一大挑战。

发明内容

为解决现有技术存在的问题，本发明提出一种基于对偶边界元和应变能优化分析的弹性断裂仿真方法，该方法主要利用对偶边界元方法来表示含有重叠裂口表面的模型，并采用半解析方法来求解超奇异边界积分方程。在断裂的物理仿真过程中，通过引入接触力模型来表示来自接触载荷和裂纹前端应变能释放的连续属性，并且在断裂持续时间段内对裂纹的延展距离计算断裂因子来对裂口端的延展进行精确控制。

本发明是采用以下的技术方案实现的：一种基于对偶边界元和应变能优化分析的弹性断裂仿真方法，包括如下步骤：

步骤(1)、将目标模型及其裂纹进行对偶边界积分方程建模；

步骤(2)、基于分段光滑表面的边界积分方程离散；

步骤(3)、对边界积分方程中奇异的被积函数采用奇异性消减技术进行统一处理；

步骤(4)、边界积分方程中自由项的求解；

步骤(5)、在仿真过程中，使用连续接触模型对断裂产生过程进行准确建模；

采用弹性球体接触的Hertz模型建立连续接触时间段，在接触时间段内，通过定义随时间变化的外力来表示接触过程中受力的变化，通过用户指定的仿真时间步以及得到的连续接触时间段，将断裂置于一个连续的接触时间段内，来表示外力作用在模型表面的动态过程；

步骤(6)、在仿真过程中，使用应力强度因子来描述裂纹附近的应力场；

在模型受连续外力作用下进行形变仿真，得到裂口表面边界元素的位移，从而得到裂口的张开位移，在裂纹的前沿处对延展点建立局部坐标系，将张开位移投影到局部坐标基向量上，并计算应力强度因子，当延展点被多个裂纹共享时，通过对多个裂纹中的应力强度因子求平均值得到最终的应力强度因子；

步骤(7)、在产生断裂的过程中，根据表面应力分析进行断裂的初始化，并且在裂纹上基于应力强度因子进行延展；

针对边界元素的表面应力分析，首先获得当前元素的形变梯度，并得到三角面片的First Piola-Kirchhoff应力张量，当最大应力超过了材料强度，在当前三角面片上生成新的裂口；对裂纹上的点计算其应力强度因子，当有效应力强度大于材料的断裂韧性时进行裂纹的延展，在裂纹的延展过程中，计算延展的速度和延展的角度。

进一步地，步骤(1)包括：采用Betfi-Somigliana等式来对目标模型的区域Ω建立边界积分方程；源点处的边界积分方程通过源点和域边界的场点之间的距离所构成的Kelvin基础解来表示；将域中的裂纹建模成一个张开型表面，使用对偶边界元将三组表面建模成独立的位移边界积分方程，并将重叠表面上的位移边界积分方法对源点求导，得到受力边界积分方程。

进一步地，步骤(2)包括：对在计算机图形学领域内广泛存在的分段光滑表面模型进行线性离散；边界积分方程在三角面片集合上的求解采用正常的数值积分方法；对于具有数值奇异性面片集合的边界积分方程求解，需在源点处建立剔除域。

进一步地，步骤(3)包括：采用奇异性消减技术对边界积分方程中奇异的被积函数进行统一的处理；将极坐标系下的奇异的被积函数进行Taylor展开，使得被积函数的奇异项显式地表达，通过把这些奇异项从被积函数减去，得到一个不含奇异项的正规的被积函数；对被移除的奇异项求得它们的解析解，并将结果带回到边界积分方程。

进一步地，步骤(4)包括：采用球形剔除域，对于源点周围的边界元素所在的切平面集合与球形剔除域进行相交操作，得到位于球面上的顶点集合和轮廓弧线集合；定义如下变量集合：切平面集合，切平面集合的外法线集合，两两相交切平面的夹角弧度集合；得到源点所在局部表面的集合结构相关的系数矩阵。

与现有技术相比，本发明的优点和积极效果在于：

1、提出了一种采用对偶边界元方法的仿真断裂，该方法对重叠断裂表面的两侧分别使用位移边界积分方程和受力边界积分方程来表示，对两个重叠表面的采样点进行独立地建模。并对随之产生的各种奇异积分方程，采用半解析方法和鲁棒的积分变换方法来剔除奇异积分方程对数值求解的影响。

2、引入连续的接触力模型，能够将断裂仿真中施加优势负载的过程显式地表示成一个连续变化的时间区间，并能获得连续变化的负载，对断裂生成过程进行时间序上的连续控制，扩展断裂的应用场合。

3、提出了一种针对断裂延展过程中对断裂前端延伸和应变能量的显著性特征分析方法。一方面简化了裂纹生成过程中参数的计算，提高计算效能，另一方面，使得裂纹生成能够体现对其先验数据的影响。

附图说明

图1为断裂仿真流程图；

图2为边界元积分方程域示意图；

图3为分段光滑表面网格模型示意图；

图4为非光滑表面局部坐标投影示意图；

图5为裂纹延展的局部坐标系；

图6为裂口负载模式；

图7为全局坐标系、局部坐标系、保角坐标系变换示意图；

图8为计算由五个切平面构成的自由项的几何结构示意图；

图9为人类股骨(髀骨)断裂仿真结果图；

图10为人体第五腰椎拥有不同断裂韧性断裂仿真结果图。

具体实施方式

本发明提出一种基于对偶边界元和应变能优化分析的弹性断裂仿真方法，下面结合具体方法步骤，介绍本发明原理。

1)分段光滑表面的对偶边界元方法建模

步骤(1)、采用Betfi-Somigliana等式来对目标模型的区域Ω建立边界积分方程。源点处的边界积分方程通过源点和域边界的场点之间的距离所构成的 Kelvin基础解U_ij和T_ij来表示，其中Kelvin基础解表示了在源点处施加了一个单位负载(Unit Load)后，经过材质影响后，反应到场点处的位移。当场点趋近源点时，基础解U_ij具有弱奇异性(WeakSingularity),此类弱奇异项可以通过极坐标变换(Polar Coordinate Transformation,PCT)的方法将U_ij转换成可积的被积函数(Integrable Integrand)。T_ij具有强奇异性(Strong Singularity), 因此需要引入柯西主值积分(Cauchy Principal Value，CPV)，来进行处理。当目标模型的域Ω中含有裂纹时，将域中的裂纹建模成一个张开型表面，域Ω边界Γ变成了三部分(Γ＝Γ_r+Γ_c++Γ_c-)其中Γ_c+,Γ_c-是重叠的表面，它们的位置相同，方向相反，在这里我们假设在Γ_c+,Γ_c-上没有受力(Traction-Free)。使用对偶边界元(Dual-BEM)将三组表面建模成独立的位移边界积分方程(DBIEs)，并将重叠表面上的位移边界积分方法对源点求导，得到受力边界积分方程 (TBIEs)。

步骤(2)、本发明所用到的边界元方法主要针对的是计算机图形学领域，需要对在计算机图形学领域内广泛存在的分段光滑表面模型进行线性离散 (LinearDiscretization)。由于使用分段光滑表面作为模型的边界，导致了源点s是一个不连续(Non-smooth)的点，被多个三角面片共享(如图3(a))。由于在边界积分方程中的被积函数F对场点x(积分点)到源点s的距离具有奇异性，将分段连续表面Γ分成了两个部分：奇异面片集合Γ_s和三角面片集合Γ_r，其中奇异面片集合Γ_s是由共享源点s的面片组成的(如图3(a)中的绿色部分)，余下的灰色部分构成了三角面片集合Γ_r。边界积分方程在三角面片集合上的求解采用正常的数值积分方法，并能获得良好的计算精度。对于奇异面片集合的边界积分方程求解，方法需要在源点s处建立剔除域(Excluded Domain)，该区域是一个以源点s为圆心的半径为ε的球体(如图3(b))，方法采用球型剔除域是为了利用它在对自由项c_ij计算时的优点。分段光滑表面的对偶边界元方法建模往往同样条函数(NURBS)一起使用，样条函数对于模型网格的表示要求较高，虚拟手术模型网格数据(或者说图形学领域广泛使用的网格模型)中，一般使用分段光滑网格(三角网格)，本发明在边界元方法与图形学离散方法之间建立了链接，使得边界元方法可以应用在图形学领域，这种链接不是简单的尝试，而是提出了在三角网格上有效使用边界元方法需要进行的数值优化方法，否则会导致运算不收敛的情况。

步骤(3)、经过步骤(2)的处理得到了一个在分段光滑表面上定义的边界积分方程。方程中含有各种奇异的被积函数。为了能够对这些奇异的被积函数进行统一的处理，采用奇异性消减技术(Singularity Subtraction Technique， SST)。将极坐标系下的奇异的被积函数进行Taylor展开，使得被积函数的奇异项可以显式地(Explicitly)表达，通过把这些奇异项从被积函数减去，得到一个不含奇异项的正规的(Regular)被积函数。对被移除的奇异项求得它们的解析解(Analytic Solution)，并将结果带回到边界积分方程。通过使用奇异性消减方法可以避免公式边界积分方程中的极限操作。

2)边界元方程的数值求解

步骤(4)、在边界元建模过程中，Somigliana恒等式一般包含了两部分，一部分是通过边界积分方程对具有奇异属性的kelvin基础解进行求解，另一部分是与源点所在局部表面的几何结构相关的系数矩阵C。对于本发明所涉及的在分段光滑表面上建模三维弹性断裂问题，采用显式计算技术，该技术支持在分段光滑表面上显式地计算被任意多个切平面所共享的边界顶点处的系数矩阵 C。本步骤采用了专为多边形网格开发的自由项计算技术，相比于传统的应用场景，如CAD设计，该领域的集合模型是基本固定的，都是由各种基本几何形状进行拼接，所以，它们计算自由项的方法是，事先(穷举)计算好(使用积分等运算量大的方法)，使用的时候，去引用(查表)，本发明采用的方法能够对于任意多边形几何结构继续快速计算，很好地适用于断裂过程中出现的复杂的几何情况。

3)基于连续接触模型的断裂仿真方法

步骤(5)、为了更好地对断裂的产生过程进行准确地建模，引入连续接触模型。采用弹性球体接触的Hertz模型来建立连续接触时间段。通过定义随时间变化的外力来表示接触过程中受力的变化，由于Hertz接触模型是一种非粘着接触模型(Non-adhesiveContact Model)，将接触外力视为一个正弦函数 (Sinusoidal Function)。通过用户指定的仿真时间步以及得到的连续接触时间段，将断裂置于一个连续的接触时间段内，来表示外力作用在模型表面的动态过程。

步骤(6)、对于裂纹的扩展，因为应力在裂纹尖端具有奇异性，所以使用应力强度因子(Stress Intensity Factors,SIFs)来描述裂纹附近的应力场。在模型受连续外力作用下进行形变仿真，得到裂口表面边界元素的位移，从而得到裂口的张开位移(CrackOpening Displacement,COD)。计算应力强度因子需要在在裂纹的前沿处对延展点建立局部坐标系，对于延展点c，根据其所在裂纹及裂纹所属边界元素的法线以及平行于裂纹的切线，建立相互垂直的局部基向量，将张开位移投影到局部坐标基向量上，并计算应力强度因子；当延展点被多个裂纹共享时，通过对多个裂纹中的应力强度因子求平均值得到最终的应力强度因子。

步骤(7)、对于物体断裂过程的建模需要处理裂纹的初始生成(CrackInitiation)和延展(Crack Propagation)。断裂的生成是基于表面应力的分析，对模型拥有的所有边界元素逐一计算其主应力和主应力方向，首先获得当前元素的形变梯度(Deformation Gradient)，并得到三角面片的First Piola-Kirchhoff应力张量，当最大应力超过了材料强度，在当前三角面片的质心，以最大应力的方向为裂口表面的法线，插入新的裂口。对裂纹上的点计算其应力强度因子，当有效应力强度大于材料的断裂韧性时进行裂纹的延展，在裂纹的延展过程中，计算延展的速度和延展的角度。

下面结合附图和具体实施方式对本发明做进一步详细地说明。

图1给出了基于对偶边界元和应变能优化分析的弹性断裂仿真过程的总体处理流程，本发明实施例一种基于对偶边界元和应变能优化分析的弹性断裂仿真方法，其步骤包括：

步骤(1)对于三维的线性弹性问题，本方法采用Betfi-Somigliana等式来对目标模型的区域Ω(如图2所示)建立边界积分方程(见公式1)。

在公式1中，源点(Source Point)s处的边界积分方程可以通过s和域Ω边界

的场点(Field Point)x之间的距离r(r＝|s-x|)所构成的Kelvin基础解U_ij和

来表示。其中，Kelvin基础解U_ij表示了在源点 s处施加了一个单位负载(Unit Load)后，经过材质影响后，反应到场点x处的位移。在本方法中，当场点x趋近源点s时(r→0)，公式1中的U_ij具有弱奇异性(Weak Singularity),U_ij＝O(r^-1),此类弱奇异项可以通过极坐标变换(Polar Coordinate Transformation,PCT)的方法将转换成可积的被积函数 (IntegrableIntegrand)。不幸的是，当r→0时，T_ij具有强奇异性(Strong Singularity),T_ij＝O(r^-2),因此需要引入柯西主值积分(Cauchy Principal Value，CPV，用符号

表示)，来进行处理。公式1中的u和t分别表示场点处的位移和受力，下标i,j表示三维空间中的方向(i,j＝1,2,3)。

当目标模型的域Ω中含有裂纹时(如图2(b)所示)，本方法将裂纹建模成为一个张开型表面(Open Surface)。对于图2(b)中的情况，域Ω边界Γ变成了三部分(Γ＝Γ_r+Γ_c++Γ_c-)其中Γ_c+,Γ_c-是重叠的表面，它们的位置相同，方向相反，在这里我们假设在Γ_c+,Γ_c-上没有受力(Traction-Free)。对偶边界元(Dual-BEM)通过将三组表面建模成独立的位移边界积分方程(DBIEs)，并将重叠表面Γ_c-上的位移边界积分方法对源点s求导，得到受力边界积分方程 (TBIEs)来解决了这个问题(如公式2，公式3所示)。可以看到，受力边界积分方程中的被积函数

和

由于求导操作使得它们在源点s处具有了强奇异性,V_ij＝O(r^-2)和超奇异性(Hyper Singularity)， W_ij＝O(r^-3)，分别用柯西主值积分符号

和Hadamard有限部分(Hadamard Finite Part，HFP)积分符号

来表示。

公式2中Γ_c+的下标“+”和公式3中Γ_c-的下标“-”是为了区别来自裂口处两侧的重叠表面。公式2和公式3中左侧出现的被称为自由项(Free Terms)，该项是由Somigliana等式的推导产生的，c_ij反应了源点s处的几何结构。

步骤(2)由于所用到的边界元方法主要针对的是计算机图形学领域，需要对在计算机图形学领域内广泛存在的分段光滑表面模型(如图3(a)所示)进行线性离散(LinearDiscretization)。假定所使用的模型是由C¹连续的三角面片(Patches)组成的分段光滑表面，并对其进行线性断裂分析。由于使用分段光滑表面作为模型的边界，导致了源点s是一个不连续(Non-smooth)的点，被多个三角面片共享(如图3(a)所示)。不仅如此，由于在边界积分方程中的被积函数F(在这里，F可以是U,T,W,V中的任何一个)对场点x(积分点)到源点s的距离具有奇异性，很自然地将分段连续表面Γ分成了两个部分：奇异面片集合Γ_s和三角面片集合Γ_r，其中奇异面片集合Γ_s是由共享源点s的面片组成的(如图3(a)中的绿色部分)，余下的灰色部分构成了三角面片集合Γ_r。边界积分方程在三角面片集合上进行求解时可以采用正常的数值积分方法(如 Gauss Quadrature Rule)，并能获得良好的计算精度。对于奇异面片集合的边界积分方程求解，方法需要在源点处建立剔除域(Excluded Domain)，该区域是一个以源点s为圆心的半径为ε的球体(如图3(b)所示)，极坐标下的形式是ρ＝ε。采用球型剔除域是为了利用它在对自由项c_ij计算时的优点。经过上面的处理，在奇异面片集合Γ_s上公式1变成了公式4的形式。

公式4是Somigliana等式的理论形式(Theoretical Form)，为了将公式 4转换成在数值计算上可用的形式需要进行两处修改。第一是将每一个满足C¹连续的任意三角面片从全局空间映射(Mapping)到内蕴几何空间(Intrinsic Geometric Space)中的规范域(本发明使用直角三角形)，以便对其进行数值积分。为了处理奇异性被积函数，我们采用极坐标(Polar Coordinates)来表示参量空间。在图4中显示了位于非光滑边界处的顶点被多个面片共享的源点s 和s的剔除域T_ε在坐标变换中的状态。在图4(b)中，可以看到在极坐标系下，是以源点s在内蕴几何空间中的映像(Image)为原点的。需要注意的是，源点s 处的球型剔除域(Spherical Excluded Domain)经过坐标变换之后并不一定在内蕴几何空间中保持球型(如图4(b)所示)。经过极坐标坐标变换，公式4 的被积函数U_ij(s,x)u_j(x)和T_ij(s,x)t_j(x)变成了

和

其中N 代表场点x在内蕴几何空间中的关于源点s的形函数，J表示从全局空间到内蕴几何空间的雅克比变换(Jacobian Transformation)，ρ代表极半径(Polar Radius),θ代表极角(Polar Angle)。第二是将公式4变成在多个面片共享源点s情况下的边界积分方程(见公式5)，在转变的过程中要保持映射形式在相邻面片之间的连续性。在公式5，N表示的是共享面片的数量。

步骤(3)经过了步骤(2)的处理，得到了一个在分段光滑表面上定义的边界积分方程，公式5，方程中含有各种奇异的被积函数。为了能够对这些奇异的被积函数进行统一的处理，我们采用奇异性消减技术(Singularity Subtraction Technique，SST)。奇异性消减方法首先将极坐标系下的奇异的被积函数进行Taylor展开，使得被积函数的奇异项可以显式地(Explicitly)表达，通过把这些奇异项从被积函数减去，得到一个不含奇异项的正规的 (Regular)被积函数。最后对被移除的奇异项求得它们的解析解(AnalyticSolution)，并将结果带回到边界积分方程。通过使用奇异性消减方法可以避免公式5中的极限操作。为了方便描述，本方法只对超奇异(Hyper-singular) 被积函数

进行奇异性消减操作，见公式6。为了使用奇异性消减方法，本方法假定位移向量(DisplacementVector)u在场点x趋近于源点s时满足

连续条件，记做u∈C^0,α(s)。在位移向量u满足

连续条件的基础上，方法能够对公式6中的奇异性积分项

进行级数展开操作，从而得到公式7。

在公式(5-7)中

和

是两个只与极角(Polar Angle)θ相关的函数。执行奇异性消减操作的过程可以表示成公式8，公式8中右侧的双重积分 (DoubleIntegral)方程已经消减了(Subtraction)奇异项

其外侧的极限操作与积分方程实现了解耦，可以被当做正常积分方程进行数值积分处理。

公式8右侧的单层积分(Single Integral)是需要被进行解析求解的奇异项部分。最终步骤(2)得到的公式5可以被处理成只含有常规积分方程的形式，见公式9。对于强奇异项(Strong Singular Term)，弱奇异项(Weak Singular Term)而言，它们的处理方式和超奇异项是相同的，而且处理上会更简单一些。

步骤(4)针对本方法所面临的分段光滑表面上三维弹性断裂边界元建模问题的求解提供稳定高效的数值求解方法。从数值积分和自由项求解两方面对现有方法提供进一步的数值优化。经过步骤(3)的处理，初始的含有超奇异性被积函数的边界积分方程，公式1，已经通过奇异性消减方法(SST)转变成为规则的数值积分方程，见公式9。在公式9中，可以看到在对含有奇异性被积函数的边界积分方程进行奇异性消减处理中，使用了极坐标来对参照空间进行表示 (见步骤(2))以利用极坐标在处理奇异性被积函数方面的优势。但是也需要注意到，虽然极坐标表示的规则的数值积分方法同样可以采用标准的高斯数值积分方法(Standard Gaussian Quadrature Rules)来对其进行求解，但是需要通过加大对积分点的采样密度来满足数值计算精度。不仅如此，对于边界元素(本发明使用三角面片，Triangle Element)具有较大纵横比(Aspect Ratio) 的情况，元素所对应的被积函数仍然会具有奇异属性。出现这种情况的原因在于通过公式7所得的对被积函数的Laurent级数展开形式中，函数

虽然已经与极轴变量ρ进行了解耦，但是由于其仍然是极角θ的函数，在对极角θ进行积分操作过程中仍然会显现出来自极角θ奇异性。针对上述问题，通过将公式7 得到的极角θ的函数投影变换到稳定的保角空间中，来提高现有边界积分方程的数值求解稳定性。

对于具有较大纵横比的边界元素来说，公式7中所得的函数具有如下表达形式：

其中分子处的V_i(θ)是θ的函数，p是由下标i确定的常数。函数A(θ)是一个反映了边界元素几何形状信息的函数，其展开式如公式11所示，可以看到公式11中的

表示了边界元素的全局坐标到局部坐标(ξ₁,ξ₂)上映射的基向量。我们通过定义两个空间基向量长度之比λ和两个空间基向量的夹角 γ的余弦cos(γ)来对投影空间进行度量，详见公式13。由公式13可以得到两个辅助系数

和μ，见公式14。对公式11使用倍角公式并将公式13带入，可以得到新的函数A(θ)的表达式，公式12。可以看到，当μ趋近于1时，系数θ可以通过使

导致函数A(θ)→0，进而导致函数

拥有奇异属性。从几何的角度看，有两个因素导致了上述情况(μ→1)的出现：(1)由于投影空间基向量u₁和u₂之间趋近于重叠(夹角γ→0)或者平行(夹角γ→π) 导致；(2)由于|u₁|和|u₂|的比值趋于0或者无穷，即纵横比过大，导致 (λ+λ^-1)²→∞。

解决此类问题的方法就是在现有全局坐标系(图7)向局部坐标系(图7b) 投影的基础上再进行一次保角变换(Conformal Transformation)，如图7c所示，通过保角变换使得新的投影空间基向量

满足如下条件：

当将新的基向量带入到公式12中，函数A(θ)由之前的关于极角θ的函数变成了一个常数项

由此带来的好处在于：(1)由于A(θ)是一个常数项，公式10中的

对于极角θ的奇异性已经被解耦了；(2)一个简单的公式10可以降低沿着极角θ进行数值积分的复杂度。

为了实现这一目标，需要计算两个辅助系数

和

(见公式16)，用来表示新的投影空间中点的坐标，如图7(c)所示，其中λ和cos(γ)由公式(5-23) 定义。于是可以得到新的投影矩阵T(见公式17)，对现有的局部坐标系(见图7b)使用投影矩阵T进行投影，得到了新的局部坐标(见图7c)，可以看出，新的保角空间的基向量

和

满足公式15。由于进行了新的保角空间投影，所以需要对公式9进行修改，得到公式18。公式中的|T^-1|由公式17所定义，极角 θ变成了新的投影空间中的极角。

对于本方法所涉及的在分段光滑表面上建模三维断裂问题，采用自由项显式计算技术，支持在分段光滑表面上显式地计算被任意多个切平面所共享的边界顶点处的自由项系数矩阵。

为了描述方便，下面通过一个示例来说明对系数矩阵C的计算过程，该过程可以被运用到分段光滑表面的通用情况中，包括重叠断裂面的情况。如图8 所示，在示例中假设当前源点s被五个边界元素(满足C¹连续条件的任意三角面片)所共享。由于采用了球型剔除域(见步骤(2))，所以5个边界元素所在的切平面集合π与球型剔除域进行相交操作，得到了位于球面上的顶点集合为 v＝{v₁,v₂,v₃,v₄,v₅}和轮廓弧线集合为γ＝{γ_1,2,γ_2,3,γ_3,4,γ_4,5,γ_5,1}。在这里，将顶点和轮廓弧线进行顺时针排序，如图8所示。为了方面描述，定义如下变量集合： (1)定义切平面集合π＝{π_1,2,π_2,3,π_3,4,π_4,5,π_5,1}；(2)定义平面集合π的外法线集合为n＝{n_1,2,n_2,3,n_3,4,n_4,5,n_5,1}；(3)定义两两相交切平面的夹角弧度集合为 α＝{α₁,α₂,α₃,α₄,α₅}，对于α_i的计算参见公式19：

α_i＝π+sgn((n_i-1,i×n_i,i+1)·r_i)arccos(n_i-1,i·n_i,i+1) (19)

在公式19中，r_i＝r(v_i)＝v_i-s，s表示源点。sgn(x)是符号函数。

至此，我们可以得到一个由n个切平面相交所构成的局部几何结构所对应的系数矩阵C的通用公式20。使用公式20得到的是一个针对源点s的3×3矩阵 C＝{c_ij}，其中i,j＝1,2,3。

步骤(5)为了更好地对断裂的产生过程进行准确地建模，本方法引入连续接触模型，将断裂置于一个连续的接触时间段内，来表示外力作用在模型表面的动态过程。在接触时间段内，通过定义随时间变化的外力来表示接触过程中受力的变化。为了建立连续接触时间段，本发明采用弹性球体接触的Hertz模型，见公式21。

其中，m代表质量，E代表弹性模量，R代表接触球的半径，v代表模型与外界接触时的速度，c是Hertz模型的常系数(本发明所使用的参数，c＝2.87)。通过用户指定的仿真时间步Δt以及公式21得到的连续接触时间段可以将断裂过程视为一个时间段内的仿真操作。

由于Hertz接触模型是一种非粘着接触模型(Non-adhesive Contact Model)，本方法将接触外力视为一个正弦函数(Sinusoidal Function)，见公式22。

其中，f_max代表接触过程中模型表面所受到的最大外力。

步骤(6)针对裂纹的扩展，因为应力在裂纹尖端具有奇异性，所以使用应力强度因子(Stress Intensity Factors,SIFs)来描述裂纹附近的应力场。在模型受连续外力作用下进行形变仿真，得到裂口表面边界元素的位移，从而得到裂口的张开位移(CrackOpening Displacement,COD)，并在裂纹处根据张开位移计算应力强度因子(SIF)。由此可见，裂纹延展是以裂纹尖端为基础，根据应力强度因子和模型的材料属性进行扩展的过程。

计算应力强度因子需要在裂纹的前沿处对延展点建立局部坐标系。对于延展点c所处的局部坐标系来说，首先要确定延展点c所在的裂纹ct以及裂纹所属的边界元素elt_ct，得到边界元素elt_ct的法线n₁以及平行于裂纹ct的切线n₃，通过n₁与n₃的向量积(CrossProduct)得到垂直于裂纹ct的法线n₂，如图5所示。在图5中，n₁，n₂，n₃是三个相互垂直的局部基向量，其中n₁垂直于边界元素elt_ct，对应于张开方向的受力(见图6中ModeI)，n₂垂直于裂纹ct，对应于剪切方向的受力(见图6中ModeII)，n₃平行于裂纹ct，对应于滑动方向的力(见图6 中ModeIII)。

对于Δu可以通过边界元素elt_ct中不属于裂纹ct的顶点p来得到。得到张开位移Δu之后需要将其投影到局部坐标基向量上，得到Δu_I，Δu_II和Δu_III，然后根据公式13计算应力强度因子K_I，K_II和K_III，公式23中的μ代表剪切模量 (Shear Modulus)，μ＝E/(2+2v)，E代表杨氏模量，v代表泊松比，r代表顶点p到裂纹ct的距离。当延展点被多个裂纹共享时，可以通过对多个裂纹中的应力强度因子求平均值得到最终的应力强度因子。

步骤(7)对于物体断裂过程的建模需要处理裂纹的初始生成(Crack Initiation)和延展(Crack Propagation)。对于边界元方法来说，断裂的生成是基于表面应力的分析，计算边界元素(三角面片)的最大主应力(Maximal Principal Stress)，当最大主应力超过材料强度时就会导致断裂现象。对模型拥有的所有边界元素逐一计算其主应力和主应力方向，然后检查当前边界元素是否需要产生裂纹，当满足产生裂纹条件时，以当前主应力方向作为裂口表面的法线初始化裂口，并以这些裂纹为基础可以进行裂纹的延展。当没有新的裂纹面产生时裂纹生成结束。

在裂纹初始化阶段，为了执行边界元素的表面应力分析，首先要获得当前元素的形变梯度(Deformation Gradient)，考虑到本发明所采用的的边界元素为空间三角形，它的三个顶点在形变之前和形变之后的状态不能完全确定空间三角形的仿射变换(这是因为缺少了垂直于空间三角形的维度)。为了解决这个问题，本发明通过为空间三角形增加其垂直方向的点，见公式24。假设空间三角形的三个变形前的顶点为v_i,i＝1,2,3，三个变形后的顶点为

公式 24中的v₄就是新增的垂直于空间三角面片的点，使用公式24中同样的方法可以得到变形后的

得到垂直方向的点后，就可以构建空间三角面片的边矢量矩阵V＝[v₂-v₁ v₃-v₁ v₄-v₁]和

使用公式24得到形变梯度F_3×3。通过对形变梯度F执行奇异值分解(Singular Value Decomposition，SVD)得到F_3×3＝U∑V^T，将奇异值矩阵∑带入公式26计算三角面片的First Piola-Kirchhoff应力张量P，其中P代表单位面积上的应力，应力的方向对应于V的列向量。

P＝2μ(∑-I)+λTr(∑-I)I (26)

当应力张量P中的最大应力超过了材料强度，需要在当前三角面片上生成新的裂口。首先定位当前三角面片的质心，以最大应力的方向(来自V的列向量)为裂口表面的法线，插入新的裂口。

当裂纹初始化之后，就可以在仿真过程中对裂纹进行延展操作。首先，通过对裂纹上的点使用公式23计算其应力强度因子K_I，K_II和K_III。然后，将K_I， K_II和K_III带入公式27得到有效应力强度K_eff(Effective Stress Intensity)。通过将有效应力强度K_eff与材料的断裂韧性K_c进行对比，当K_eff大于K_c时进行裂纹的延展。

在裂纹的延展过程中，需要计算延展的速度v_cp和延展的角度θ_cp，可以通过公式28和公式30得到。在公式28中c_R代表了Rayleigh波速，

ρ是材料的密度。为了计算延展角度θ_cp，需要先将K_I和K_III通过公式29进行合并得到，这是由于K_I和K_III主导了裂纹向前的延展，K_II导致裂纹侧向旋转延展。将K_I,III和K_II带入公式30得到了θ_cp，最后需要将延展角度θ_cp转换成全局坐标下的方向。

本方法原型系统的实现使用的软件平台为Microsoft visual studio 2010 与OpenGL，硬件平台为有Intel Core i7 3.60GHz CPU、NVIDIA GeForce GTX 770GPU、4G内存。方法效果图如图9，图10所示。

以上所述，仅是本发明的较佳实施例而已，并非是对本发明作其它形式的限制，任何熟悉本专业的技术人员可能利用上述揭示的技术内容加以变更或改型为等同变化的等效实施例应用于其它领域，但是凡是未脱离本发明技术方案内容，依据本发明的技术实质对以上实施例所作的任何简单修改、等同变化与改型，仍属于本发明技术方案的保护范围。

Claims

1.一种基于对偶边界元和应变能优化分析的弹性断裂仿真方法，其特征在于：包括如下步骤：

步骤(1)、将目标模型及其裂纹进行对偶边界积分方程建模；

步骤(2)、基于分段光滑表面的边界积分方程离散；

步骤(4)、边界积分方程中自由项的求解；

2.根据权利要求1所述的基于对偶边界元和应变能优化分析的弹性断裂仿真方法，其特征在于步骤(1)包括：采用Betfi-Somigliana等式来对目标模型的区域Ω建立边界积分方程；源点处的边界积分方程通过源点和域边界的场点之间的距离所构成的Kelvin基础解来表示；将域中的裂纹建模成一个张开型表面，使用对偶边界元将三组表面建模成独立的位移边界积分方程，并将重叠表面上的位移边界积分方法对源点求导，得到受力边界积分方程。

3.根据权利要求1所述的基于对偶边界元和应变能优化分析的弹性断裂仿真方法，其特征在于步骤(2)包括：对在计算机图形学领域内广泛存在的分段光滑表面模型进行线性离散；边界积分方程在三角面片集合上的求解采用正常的数值积分方法；对于具有数值奇异性面片集合的边界积分方程求解，需在源点处建立剔除域。

4.根据权利要求1所述的基于对偶边界元和应变能优化分析的弹性断裂仿真方法，其特征在于步骤(3)包括：采用奇异性消减技术对边界积分方程中奇异的被积函数进行统一的处理；将极坐标系下的奇异的被积函数进行Taylor展开，使得被积函数的奇异项显式地表达，通过把这些奇异项从被积函数减去，得到一个不含奇异项的正规的被积函数；对被移除的奇异项求得它们的解析解，并将结果带回到边界积分方程。

5.根据权利要求1所述的基于对偶边界元和应变能优化分析的弹性断裂仿真方法，其特征在于步骤(4)包括：采用球形剔除域，对于源点周围的边界元素所在的切平面集合与球形剔除域进行相交操作，得到位于球面上的顶点集合和轮廓弧线集合；定义如下变量集合：切平面集合，切平面集合的外法线集合，两两相交切平面的夹角弧度集合；得到源点所在局部表面的集合结构相关的系数矩阵。