第1章 绪论
　　1.1 引言
　　数值分析作为当今公认的四大科学研究范式(理论、试验、计算、数据)之一，起到连接其他三类科学研究范式的桥梁作用，在科学探究之路中持续扮演着重要的角色。等几何分析方法(isogeometric analysis， IGA)作为数值分析领域的前沿方向之一，自?Hughes及其合作者于2005年[1]在计算力学领域顶级期刊Computer Methods in Applied Mechanics and Engineering上发表以来，因其克服了数值计算中的几何离散误差(其他误差来源主要有模型近似误差、数值计算误差等)，受到众多学者的持续关注和跟进研究，其知识版图和应用领域不断地拓展，近二十年来其发展的研究谱系如图1.1所示。
　　图1.1 等几何分析方法研究谱系
　　前期的研究主要集中在将等几何分析方法的思想引入单*一类数值方法中，充分利用等几何模型的精确几何表达、形函数的非负高阶次和高连续性、CAD与CAE之间的无缝连接等优点，对传统数值分析方法进行改进、加强，或者重塑其分析框架，如等几何有限元法(isogeometric analysis finite element method， IGAFEM)、等几何边界元法(isogeometric analysis boundary element method， IGABEM)、等几何物质点法(isogeometric analysis material point method， IGAMPM)等。然而，正所谓“尺有所短，寸有所长”，每类等几何数值方法都有其发展技术瓶颈和应用局限性。开展各类等几何先进数值方法的耦合算法(如本书中的等几何有限元-边界元耦合方法)研究，可以充分发挥每种方法的优点，达到互补的目的。
　　本章*先介绍本书中提到的等几何有限元法、等几何边界元法及等几何有限元-边界元耦合方法的相关研究进展和应用，然后介绍本书的章节内容安排。
　　1.2 等几何有限元法
　　传统数值方法在分析问题时，*先需要在计算机辅助设计(CAD)软件中建立几何模型，然后将模型导入计算机辅助工程(CAE)软件中进行网格划分。CAD采用非均匀有理B样条(non-uniform rational B-splines， NURBS)函数，CAE软件采用拉格朗日形函数。CAE数值计算后不能直接对几何模型进行修改，必须再次在CAD软件中修改后重新导入。*终的计算结果可视化展示采用的是双线性插值函数。上述不同几何表达之间频繁的数据传输和交换既麻烦又耗时，极大地增加了分析成本。如图1.2所示，在汽车、航空航天和造船工业中，大约80%的分析时间用于CAE网格生成。2005年，Hughes等[1]提出了等几何分析(IGA)的概念。IGA的核心思想是将CAD软件几何建模的NURBS函数直接用作数值分析的形函数。IGA完全消除了网格划分的概念，既能精确地描述几何模型，消除几何离散误差和获得更好的数值解，又大大节省了数值计算成本。
　　图1.2 美国桑迪亚(Sandia)国家实验室模型生成和分析过程中每个组成部分的相对时间成本[1]
　　如图1.3所示，IGAFEM消除了现有数值计算方法中存在的几何离散误差，实现了CAE和CAD系统的无缝连接和统一，大幅度提高了数值分析和设计的效率[2]。IGA是将CAD中普遍采用的NURBS基函数作为形函数用于数值分析，避免了烦琐的网格划分步骤和不同系统间频繁的数据交换，大大缩短了产品的开发周期。等几何概念的提出主要是实现了工程设计与计算模拟之间数学模型表达方式的统一。2017年，来自挪威科技大学的Stahl与其合作者将这一概念拓展到了结果可视化的领域[2]，即相同的数学表达既用于几何描述和数值分析，也用于后处理过程中计算结果的可视化。他们采用Bézier分解将结构性样条和非结构性样条单元均转化为Bézier单元，进行可视化操作，这样做的好处是利于并行计算和处理局部细分样条下的可视化问题。IGA采用CAD系统中的样条函数当作形函数，同时用于模型几何表示和未知场变量近似，其主要优点如下所述。
　　图1.3 传统有限元法和等几何有限元法的区别[2]
　　(a) 传统有限元法中，设计、分析和可视化被分离，并使用不同的数学表达；(b) 等几何分析中，设计、分析和可视化过程使用相同的数学表达
　　(1) 消除了传统分析中由拉格朗日(Lagrange)多项式插值引起的几何离散误差，同时也消除了烦琐的网格划分步骤。单元细化过程可以方便地通过插入节点(h细分)和升阶(p细分)来实现，而不再需要与CAD系统进行额外的数据通信。
　　(2) 采用很少的单元就能精确描述不同尺度比例结构的几何形状，而且其描述的几何形状不随离散单元的多少而发生改变，并且可以通过其*有的“k细分”实现单元之间的高阶连续性。
　　(3) 精确的空间离散使得与*面几何相关的运动学量(*率、切线和法向量等)可以在每个计算点处精确确定，而不会产生任何误差。
　　(4) NURBS形函数的高阶次和高阶连续性能自然满足各种力学理论中对基本变量C1或更高插值连续性的要求，如薄板壳理论、Cahn-Hilliard方程[3]等。
　　(5) 进行k细分的高阶NURBS基函数在对光滑或者非光滑函数进行插值近似时不会出现不稳定的振荡现象[4]，如Runge现象。高阶NURBS形函数由于其非负性和凸包性的特点，在处理大位移和大滑动下的相互作用面问题时具有明显的优势[5]。
　　(6) 高阶连续性的样条形函数使得应变、应力等二阶张量的计算和可视化不再需要通过在传统有限元分析中广泛应用的节点平均方法来获得光滑的节点变化值[2]。
　　然而，在应用基于NURBS的IGAFEM时仍存在一些不足。
　　(1) NURBS基函数由于其张量拓扑形式的特点而缺乏有效的局部细分能力，在分析具有局部特征问题时会在非局部区域产生多余的控制点和引入冗余的自由度，因此基于NURBS基函数的IGA不易于直接分析含缺陷的结构。一些具有固有局部细分特性的样条函数被提出，尽管它们还没有与现有的CAD环境完全集成，如T样条[6，7]、PHT(polynomial splines over hierarchical T-meshes)样条[8]、LR-B样条(locally refined-B splines)[9]和细分*面(subdivision surfaces， SS)[10]等。上述样条不仅克服了NURBS在分析局部特征问题时的不足，而且还继承了NURBS基函数非负性、单位分解性和局部支撑性等利于数值分析的优点。
　　(2) 复杂几何体，尤其是那些包含局部特征(如裂纹和孔洞)的几何体，通常需要由多个NURBS片组成，这些片之间通常具有不相适应的离散网格。虽然IGA可以容易地实现单个片中更高阶的单元连续性，但是实施片与片之间的高阶连续性一直是一个开放且具有挑战性的问题。许多学者提出了一系列的解决方法，如静态凝聚法[11]、罚函数法[12]、Mortar法[13]和Nitsche法[14]等。
　　(3) IGAFEM在数值计算时需要将CAD系统中几何模型的边界表征向适合于分析的内部区域实体表征进行转换，这种转换的存在使得IGAFEM在处理复杂模型时并未能真正实现CAD与CAE系统间的无缝连接，而且这种转换的精度和有效性目前仍是一个开放性的问题[15-1617]。鉴于此，一些学者提出了基于边界样条表示的等几何分析方法(isogeometric B-rep analysis， IBRA)以实现从设计到分析过程的统一[18，19]，允许直接对复杂的CAD模型进行数值分析，无须重新建模和进行网格转换。
　　1.2.1 在壳体分析中的应用
　　壳体结构分析的两个主要且广泛使用的公式分别来自Kirchhoff-Love[20]和Reissner-Mindlin[21]壳体理论。区分薄壳和厚壳的常用标准是壳体*率半径R和厚度h的比值。按照这个标准，壳体结构通常分为厚壳(R/h<20)和薄壳(R/h?20)。Kirchhoff-Love壳体理论忽略了横向剪切变形，仅需要位移自由度就能描述壳体的变形行为，其仅适用于变形行为主要由弯*和薄膜应变决定的薄壳结构[22]。该理论的两个关键问题是：①变分公式中涉及形函数的二阶导数，这就要求壳单元之间至少具有C1连续性，这是传统FEM难以实现的[23，24]，这也是仅需要单元间C0连续性的Reissner-Mindlin壳单元在传统FEM中广泛使用的重要原因；②涉及旋转和力矩的边界条件的施加变得更加复杂，因为它涉及位移自由度的导数，而系统方程中只包含位移自由度[25]。与无旋转自由度的Kirchhoff-Love理论相比，Reissner-Mindlin壳理论考虑了横向剪切应变，可以有效地模拟薄壳和中厚壳。在其壳体的运动学描述中，使用转角作为附加未知量，有助于施加力矩边界条件[26]。此外，平衡方程中只涉及未知场变量的一阶导数。因此，单元间变形满足C0连续性要求就足够了。当然，也有学者利用上述两种壳理论优缺点之间的互补关系提出了混合壳体理论(blended shell theory)，并应用于模拟金属成形过程[27]，即低*率区域可以使用无旋转的薄壳理论进行有效分析，而高*率区域可以使用Reissner-Mindlin公式进行分析。
　　近些年来，许多学者结合IGA对上述两种壳体理论进行了大量的研究，并从多个方面说明了IGA在壳体领域相对于传统方法具有的令人信服的优势，具体可参考等几何Reissner-Mindlin壳[21，26，28-2930]和等几何Kirchhoff-Love壳的文章[20，22，31，32]。一方面，可以在每个积分点处精确计算与壳体基础几何相关的运动学量(*率、切线和法向量等)。另一方面，NURBS的高阶光滑性使得壳单元的高阶公式可以直接实现，特别是要求C1连续性的Kirchhoff-Love薄壳公式。等几何的出现使得Kirchhoff-Love壳体理论经历了一次复苏。尽管如此，与传统方法相比，NURBS更高的连续性和精确的几何描述仍然为Reissner-Mindlin壳单元的实施提供了显著的优势[33]。与传统有限元Reissner-Mindlin壳单元相比，等几何Reissner-Mindlin壳单元出现剪切自锁问题的可能性要小得多[34]，只需提高NURBS基函数的阶次即可简单和有效地减小自锁问题产生的影响。此外，在商业有限元软件LS-DYNA中已成功地实现了基于等几何Reissner-Mindlin壳单元的数值分析[35]。
　　1.2.2 在含缺陷问题中的应用
　　传统FEM在处理裂纹问题时主要存在如下几个方面的挑战[36]：①需要相适应的网格拓扑来表示裂纹，裂纹扩展路径必须事先给定且只能沿着单元边界；②基于多项式的插值形函数无法描述裂尖应力奇异场，需要构建特殊类型的单元，如四分之一节点单元[37]，并且在裂纹尖端区域要划分密集的网格，这无疑增加了计算成本；③FEM网格存在几何离散误差，因此依赖于网格投影的计算误差会不断累积。为了克服传统方法中裂纹求解对网格的严重依赖性和描述裂尖附近应力场的困难，Mo?s等[38]利用形函数的单位分解性，在位移插值近似表达式中引入包含裂纹附近解的某些已知特性的富集函数，提出了扩展有限元法(extended finite element method， XFEM)。在该方法中，通过引入阶跃函数和裂尖渐近函数而扩充了原来的形函数，分别对裂纹面的不连续和裂纹尖端的奇异性进行描述，允许离散裂纹的表达*立于底层的求解网格，且富集函数的添加只涉及裂纹附近的局部单元，这大大减轻甚至消除了重新划分网格的负担，同时也提高了求解精度。