LS-DYNA IGA同几何分析介绍（上）

将同几何分析(IGA)引入到有限元分析(FEA)的框架中，目的是使数值分析模型与计算机辅助设计(CAD)的几何模型相同。与标准的、低阶的有限元单元相比，许多研究论文已经证明了使用更高阶和更高连续性的基函数是有益和优越的分析特性。B-样条曲线(B-splines)和非均匀有理B样条曲线 (NURBS)是CAD中使用最广泛的几何描述，近几年基于NURBS的有限元技术快速发展并被应用到LS-DYNA中。同时，LS-DYNA也开发了一些适合分析的非结构样条，以满足更复杂的工业几何形状和精细化分析的需要。本文将主要介绍：

• LS-DYNA中IGA方法概述

• IGA的基本概念

• 如何使用IGA

• IGA在LS-DYNA中的一些功能（下一期分享）

• 演示IGA在碰撞模拟、钣金成形到频域分析中的应用等（下一期分享）

基本概念

IGA的提出基于一项现实研究，美国德克萨斯大学计算与应用数学系主任、航空航天工程与工程力学教授 T. J. R. Hughes 在论文中引用了桑迪亚国家实验室在研究中发现，工程计算从模型建模、划分网格再到计算结束整个过程中，划分网格工作花费了大量时间，占据了整个流程近 70-80% 时间，而如果能找到一种解决办法减少网格划分的工作量，计算效率将大大提高，这是最初提出IGA想法的原因。

传统有限元最初由CAD几何模型进行网格划分，然后进行计算，这里需要大量的时间在网格划分上。对于曲面复杂的模型有限元分析虽然细化了网格，但与初始CAD模型只是近似，并不完全一样，在计算过程中可能出现结果不准确的问题。同几何分析是利用CAD相同的几何， CAD几何模型与用于计算上的几何模型是同一个几何模型。

传统有限元中用的是同参数法（ISOPARAMETRIC），即用来计算的几何模型和形变几何模型用相同的参数近似，都是用低阶拉格朗日多项式（比如LS-DYNA线性单元一样）。而同几何法（ISOGEOMERTRIC）的区别在于，用于计算的几何模型和CAD模型一致，好处是计算完成之后根据计算结果，如果需要反馈回去修改原来的几何模型，可以直接在计算模型上修改，而有限元方法可能需要用计算模型近似或mapping回CAD 模型等的操作。总结来说，IGA是用FEA的方式去计算，使用准确的几何模型，以简化建模的过程，提高计算效率，同时计算用的几何与CAD设计用的是同一个几何，这样设计和分析、模型修改的数据是互通，使得双方反馈更流畅，流程更简易。

IGA中最常用的样条曲线为NURBS（Non-Uniform Rational B-splines）是一种B-样条曲线，此外还有T-splines，LR-splines，LRB-splines，HR-splines, THR-splines等曲线，后四种可以做局部的网格细化。本文将主要介绍NURBS，是用得比较多的，也是在LS-DYNA中最早采用的样条曲线。

基础知识

随着越来越多的专家学者开始研究IGA，产生了各种各样不同的样条曲线，但并非所有的样条曲线都适合分析。图中展示了部分适合用来进行工程分析的样条曲线，这些曲线未来也将会添加到LS-DYNA中，里面包括：

• Unstructured Splines，可以模拟复杂的形状；

• Hierarchical splines可进行局部网格细化；

• B-splines是一种比较规则的曲线；

• Bernstein-Bézier forms等

B-样条曲线（B-splines）和非均匀有理B样条曲线 (NURBS)是CAD中使用最广泛的几何描述，下文将就此进行介绍。

图中的公式为B-splines的基函数，i代表基函数序号，p为阶数，ξ为节点向量（是一系列参数坐标），B-样条曲线的特征包括：

• 不减的

• 递归方式计算

• 数值大于等于零的，没有负值

• 基于节点向量和阶数

• 单位分解

B样条曲线与普通有限元参数坐标形式对比。

B-样条曲线的表达类似传统有限元，选取一些控制点，前面乘上基函数，然后相加。

如何定义B样条曲线单元？包含两个概念控制点Control point（红色圆点）以及节点向量Knot vector（蓝色菱形点）。每一个单元是一段非零的节点向量，如上图，一共有5个单元。

与有限元的网格细分功能一样，B样条曲线也可以做网格细分，包含h-refinement以及p-refinement， h-refinement是指element数量会增加，但曲线的连续性不变，通过插入knot value（如0与1之间插入0.5）增加单元数来实现。同时还可以保证curve形状不变。h-refinement的Control point会增加。

p-refinement代表阶数continuity增加，但element数量不变。实现方法：通过加入相同的knot({0,0,0}变成{0,0,0,0}，此时element数量不变但阶数提高了1，由于knot vector变化它的Control point数量和原来相比也相应增加（位置可能也不一样），这是为了保证curve在refine之后形状保持不变。

与有限元网格细化不同的是，有限元中没有k-refinement概念（k-refinement为同时进行h-refinement和p-refinement），有限元中没有对应的细分网格方法。IGA可进行k-refinement，在增加element数量的同时，又能提高阶数。

以上图右上角案例为例，原来的knot vector是{0,0,1,1}，阶数p为1。第一种方案先插入knot，使得参数坐标变成{0,0,1/3,2/3,1,1}，3个element，阶数p仍然为1。然后提高阶数，通过插入重复的knot，参数坐标变成{0,0,0,1/3,1/3,2/3,2/3,1,1,1}，此时阶数p为2。

若反过来采用先提高阶数的第二种方案，{0,0,1,1}变为{0,0,0,1,1,1}，随后再插入knot，最后得到结果完全不一样的参数坐标{0,0,0,1/3,2/3,1,1,1}，最后也是得到3个单元，但他们的基函数是不一样的。因此从这里可以得出h-refinement和p-refinement顺序不能随意改变，现实中通常采用第二种方案（先提高阶数，再增加knot）。k-refinement结合了h-refinement和p-refinement以确保C（p-1）的光滑性以及连续性。第二种方案单元和单元之间可以提高阶数（p等于2），但是左边的第一种方案单元和单元之间的连接处，p仍然是1并没有变化(光滑性没有变化)。因此，若需要将整个区域进行全部的Pure k-refinement，除非只有一个单元。

为什么我们需要NURBS呢？如果是一个圆，我们用B-样条曲线就可以。如果是椭圆形或抛物形，双曲线形，B-样条曲线就无法准确的来描述这类图形，这里就需要用到NURBS。从名字可以看出为非均匀的，每个Control point贡献的权重可能不一样。

所以在NURBS里面会引进一个权重参数w，NURBS curve形式与B-样条曲线类似，只是多考虑了权重参数的线性组合。右边展示了权重不同导致某个曲线变化的案例，若改变第9个Control point权重，黄色曲线代表 w9 等于0的情况，紫色曲线代表 w9 等于10的情况，也说明Control point的权重不一样，那么画出来的曲线是不一样的。

前面是介绍曲线，这里再介绍曲面和体。NURBS surfaces也是线性组合（linear combination），不同点在于它是二维的，正如上图所示。公式中二位的基函数是两个B-Spline曲线考虑权重后基函数的积：

以右图为例，它展示了R方向与 S方向分别有各自B-Spline的基函数，将他们相乘最后可以得到曲面的基函数。整个曲面呈现向量积的形状，也比较容易得到所有的Control points数量。假设一个方向Control points数量为n个，另一个方向是m个，总数量就是n*m。总共有两个互相独立的Knot vector，以及两个互相独立的基函数。

NURBS Solids概念相同，区别在于它是三维（R、S、T三个方向都分别有各自的基函数N、M、L，然后利用向量积得到NURBS Solids的基函数）。NURBS Solids特性与NURBS surfaces基本一致。

上面介绍的B-Spline和NURBS是CAD中最常用的样条曲线。由于我们在unstructured splines时会通过Bézier extraction转化成B-样条曲线，所以这里介绍一下Bernstein-Bézier forms。图中展示了Bernstein-Bézier forms的基函数表达式，与B-样条曲线类似，也是通过递归方式计算得出基函数。右侧图表分别展示了阶数为1、2、3、4时基函数的表现形式。基函数的个数是它的阶数加上1。同样地，Bernstein-Bézier forms基函数与B-样条曲线的基函数有着类似的特性，均为递归方式得到，均为非负值，满足单位分解（partition of unity和总为1），线性独立，continuity连续性p的值可以无穷大，p值越大曲线越平滑。

Bézier曲线是P+1个 Bernstein多项式基函数的线性组合，右图展示的曲线包含4个Control point（p为3），起始点P1和P4在曲线之上，而P2和P3在曲线之外，P1 P2为曲线的切线，P3 P4为曲线的切线（element boundary）。

这里对比B-样条曲线与Bézier曲线，研究它们是否可以相互转化？上图分别展示了B-样条曲线和Bézier曲线的形状，同样为二阶函数，knot vector相同均为{0,0,0,1,1,2,2,3,3,3}，可以看到B-样条曲线一共有5个基函数，5个Control point，而Bézier曲线却有7个基函数。尽管B-样条曲线与Bézier曲线的基函数和Control point数量不同，但却可以画出形状相同的曲线（右侧红色曲线），可以说明它们之间存在转化关系。

B-样条曲线横轴第一个element通过P1 P2 P3 三个ctronl point来控制，对应的基函数为N1、N2、N3。Bézier曲线第一个element由B1、B2、B3构成，B-样条曲线与Bézier曲线之间的转换通过C^e完成，C^e为Bézier extraction operator。图中最后一行为C^e的关系式。第二第三个等都有各自不同的C^e关系式，可以将Bézier unstructured splines转化成B-样条曲线，方便后续计算。

Unstructured splines。右侧案例可以看到，大部分的地方还是tensor product形式，其中红点Extraordinary point（EP）和T-junctions的地方不是tensor product，导致整个Patch并不是tensor product，但是每个局部的element还是可以写成tensor product形式，利用C^e将Bézier extraction形式转换成B-样条曲线，方便后续计算。对于Unstructure splines，三维的案例也是目前比较活跃的研究领域，因为二维案例中想要改变任何一个形状，可以对它进行trimming（LS-DYNA中比较成熟的方法之一），但是对于solid，如果用同样的算法的来trimming会相当复杂，这里并不推荐，但如果用Unstructure splines，在三维中也可以模拟任何形状，所以在三维中，使用Unstructure splines是比较常见的。

LS-DYNA中的同几何分析

IGA的概念由Dave Benson于2008年开始引进LS-DYNA，最初称为GENERALIZED EMEMENTS，直到2010年通过 *ELEMENT_SHELL_NURBS_PATCH 开发NURBS Shells，2015年李利萍博士开始通过*ELEMENT_SOLID_NURBS_PATCH开发NURBS Solids，2019年引入新的关键字命名*IGA（从R13版本开始）。这个keyword可以和CAD design的data Structure联系起来，如几何结构信息，拓扑结构信息等，而且也比较容易加边界条件，所以引入这个新的关键字*IGA。

LS-DYNA中的IGA经过多年发展，目前已可以支持SMP/MPP，Contact以及各种boundary condition，trimmed shells等。

NURBS Patch，也有element概念（与有限元对应更方便理解）。例如图中的蓝色部分为NURBS element从方便理解上，可以对应有限元中的element。假设有限元中某个壳单元有4个节点来定义，蓝色区域代表一个NURBS element，而这里的节点就对应NURBS中的Control point（图中蓝色点），Control point数量多少由阶数决定，图中案例阶数为二，说明一个方向有3个Control point，总共就有3*3，9个Control points贡献到这个element上。右上图是这个NURBS element的基函数。

其他的element也是由9个control point，单元之间几个共同的control point，使得element之间有一种overlap重叠(使其具有higher continuity高连续性)，重叠的大小取决于阶数。

LS-DYNA *IGA keyword可以与标准CAD data Structure兼容，会包括几何信息和拓扑信息等。几何信息包括定义NURBS Patch， Control point，knot vector等；拓扑信息包括Face，edge，point等（更易于施加边界）。*IGA 关键字目前包括壳单元，固体单元，将来也会包括trimmed solid，同时也包括Bézier format 的壳单元和固体单元。

有限元关键词与IGA关键词信息对比，以Shell单元为例。

有限元关键词与IGA关键词信息对比，以Solid单元为例。

（NISR/NISS定义Number of interpolation elements per NURBS-Element(r-/s-direction)，Shell例子（Solid会多一个NIST）在R direction或S direction，蓝色区域为NURBS Element，绿色区域为interpolation element。这里定义NISR和NISS均为2，代表R方向与S方向分别使用2个element，来interpolate 这个IGA element。使用interpolation element的原因在于，post processing目前无法直接在IGA上看结果，但是可以将IGA上的结果map 到interpolation element上，且因为interpolation element是Finite element mesh, 也可以用有限元的方式来读取结果。

此外，interpolation element还可以应用于接触分析， IGA Surface可以在接触中用来检测计算两个接触的物体还有多远的距离开始接触，interpolation element也同样可以用来检索，且方法更简便。两个不同的部件之间连接也能应用interpolation element，具体的用法下文有一个例子介绍。）

*IGA_INCLUDE关键字可以引入unstructured splines，视频中案例展示了几个solid的模型，由于solid模型有一些extraordinary point和孔洞，这几个模型无法用NURBS或只用一个Patch模拟出完整的形状。而IGA中的unstructured splines却可以模拟这些复杂的形状。右侧为左边模型的模态计算结果，以及另外两个案例的模态分析结果。

---

更多内容将在下一期继续分享，欢迎关注我们！

文章来源：2022 LS-DYNA网络研讨会，作者：李利萍博士，Ansys高级研发工程师

视频链接：LS-DYNA IGA等几何分析介绍

技术校对：王强， Ansys高级应用工程师；整理编辑：俞琴