等几何分析细化收敛性研究

谢宇洋，卜宁远

(上海理工大学机械工程学院，上海 200093)

826647688@qq.com;BNY19971997@163.com

1 引言(Introduction)

现阶段各个领域的物理性能研究大多根据有限元分析(FEM)软件仿真计算，如温度、变形和应力等性能。有限元方法在实现分析优化时，CAD数据与CAE数据的相互转换会造成模型数据特征丢失，同时点的离散性会导致优化模型的边界曲线外轮廓无法精确表达。在产品的研发中没有设计方法可满足从设计到分析在同一个模型数据下进行，而等几何方法的提出则满足了CAD模型和CAE模型的统一表达和计算。

等几何分析(IGA)是HUGHES等人基于等参单元思想提出的分析方法，它提供了计算的可能性，将分析方法集成到计算机辅助设计(CAD)工具中。等几何模型都是由非均匀有理B样条(NURBS)参数化表示的，其几何模型与分析模型不需要转化，避免了模型数据转换误差，其好处是从设计到分析所花费的时间减少，从而节省了工业生产中的时间成本。在等几何分析中，CAD和CAE的紧密耦合需要这两个领域的知识，即计算几何和计算力学，可以解决有限元优化中边界运动所导致的模型边界不规则、网格扭曲等问题，目前已应用于板壳分析、流体力学、接触分析、优化设计、电磁场等。

本文介绍了等几何分析的理论，使用体参数化造型方法生成三维几何模型。最后，通过与有限元方法的变形和应力分布等分析结果对比验证等几何分析算法的正确性，以及计算结果的收敛性差异。

2 几何模型构建(Geometric modeling)

2.1 NURBS曲面

NURBS曲面的定义：

2.2 多片分析模型

同理，NURBS实体表达式定义为：

各参数所代表的含义同NURBS曲面定义的参数。

NURBS实体基函数表达式为：

图1 圆筒几何模型图Fig.1 Cylinder geometry model diagram

3 等几何分析理论(Theory of isogeometric analysis)

3.1 等参思想

在CAD中，实体实际上只是边界的表面，而不是内部的建模，而在FEM中是一个三维实体，因此从CAD实体到FEM实体的过渡需要一个CAD的步骤边界表示法转化为有限元实体表示法。在实际应用中，很多情况下几何模型比较复杂，此时选用曲边单元代替直边单元来进行几何域的逼近。为了积分求解的简便，要在规整单元域中进行，必须在规整单元域和曲边单元域之间建立参数映射，称之为参数域和物理域之间的参数变换，即为等参变换，单元为等参单元。

在等几何分析中，需要用到两级映射，第一是母单元到参数空间映射，第二是参数空间到物理域映射。

母单元到参数空间映射公式：

参数空间到物理域映射公式：

雅可比矩阵行列式公式：

3.2 线弹性力学分析

式中，左端为等效力，为刚度矩阵，为控制点位移。其中，在等几何分析中的单元刚度K和矩阵的表达式为：

3.3 等几何分析流程

(1)针对复杂几何模型无法直接进行整体网格划分的特点，采用剪裁操作或划分子域。

(2)在各个子域内进行体参数化处理。利用NURBS建模技术，通过表面控制点的体参数化得到内部控制点，然后根据网格评判标准对控制点网格进行优化，最后得到可用于分析的控制点网格模型。

(3)选择合适的网格细分策略进行模型网格细分操作，以提高分析的结果精度，可以对关键受力点进行局部细分操作。

(4)进行单元刚度矩阵的列式求解，选用合适的高斯积分点，引入等参概念，将位移、应力等要求解的未知量表示为与几何模型相同的基函数表达形式。

(5)根据FEM中对号入座装配原则，将单元刚度矩阵组装成总体刚度矩阵。

(6)将实际载荷与约束进行等效，转移到控制点上，对总体刚度矩阵形成的方程组进行求解，然后求出控制点位移的值，进而求出其他未知变量的值。

(7)分析完成后，显示各种结果云图，并进行相应后处理。

等几何分析流程图如图2所示。

图2 等几何分析流程图Fig.2 Flow chart of IGA

4 数值结果对比(Comparison of numerical results)

接下来根据前述等几何分析理论，通过编程实现结构的静力学分析，以证明该方法的正确性。梁结构在工程中应用广泛，在机械、土木行业的大规模承载结构和电子元件等微结构设计中发挥着重要作用。梁结构采用有限元方法进行仿真分析时，需要将几何模型剖分为网格模型，对曲梁结构采用以直代曲的方式生成有限元模型，这个过程中引入了几何误差，用等几何分析方法保证边界曲面模型精度，对其进行仿真提高计算精度。在下述算例中，变量均采用国际单位制，长度的单位采用m，力的单位采用N，应力和弹性模量的单位采用Pa。梁模型尺寸如图3所示，梁左侧边完全固定，右侧施加力44,480 N。模型的材料属性杨氏模量和泊松比分别设为2.1018E10 Pa和0.3。

图3 梁模型边界示意图Fig.3 Boundary diagram of beam model

等几何分析结果图如图4(a)所示，K细化两次后的结果如图4(c)所示，在相同材料属性和约束载荷的工况下有限元分析结果如图4(d)所示。基于等几何分析的三维悬臂梁算例与有限元结果对比，从位移对和应力对比的结果看出分布基本一致。

图4 梁分析结果Fig.4 Beam analysis results

分析位移收敛曲线如图5(a)所示，应力收敛曲线如图5(b)所示，等几何分析使用更少的单元数细化两次结果收敛。

图5 梁细化收敛曲线Fig.5 The refined convergence curve of the beam

根据表1的数据显示，IGA算法分析结果和FEM分析软件对比结果误差很小，其中位移误差为3.8%，应力误差为0.25%，仅使用单片几何模型结果已接近真值，使用多片模型分析时其收敛速度更快。这验证了三维等几何分析算法的正确性，使用等几何分析算法保证了模型的高阶连续，提高了边界的精度，且不需要划分很多网格。根据表2中IGA和FEM的节点数与单元数大小可以看出，有限元分析的刚度矩阵明显比等几何分析的刚度矩阵规模大，计算量更加烦琐，因此等几何分析的计算效率高。

表1 数值精度误差Tab.1 Numerical accuracy error

表2 结果细化收敛对比Tab.2 Comparison of results refinement convergence

5 结论(Conclusion)

等几何分析采用描述几何形状的NURBS函数作为基函数，具有几何精确特性，且离散的几何形状不随单元的稀疏而改变，这意味着即使是比较稀疏的网格划分，也能精确描述研究对象的几何形状，具有很高的数值精度。NURBS本身就具有网格，一个NURBS实体包含若干个NURBS单元，分析时这些单元成为精确描述几何形状的实体单元。NURBS单元也可以细分，基函数的次数也可提高，保证了模型的光滑，计算结果更加精确。等几何分析在设计分析优化中也具有很大的优势，使用NURBS曲线来拟合模型外轮廓以保证边界精度，将控制点坐标作为设计变量进行迭代优化，完成形状优化，有效减少优化中设计变量的数目来提高优化效率。在优化后采用等几何分析算法分析体参数化后的三维模型，避免模型数据转化，提高力学分析的计算效率和精度。

等几何分析同时存在一定的局限性，一方面，对模型单元的剖分以及模型质量要求高，一些不规则图形如点阵材料等还无法通过简单的扫描、拉伸、旋转、放样等体参数化映射方法得到，一些细微的特征如倒角、圆角和小孔等也无法完全保留;另一方面，使用直接和迭代线性代数方程解算器直接解算，对于剖面结构及壳结构通常比迭代的对角预处理共轭梯度解算更有效。但是对于NURBS实体来说，为保证收敛，细化模型多次后虽然网格精细，却无法直接用求解器计算，因为超出内存，元素的长宽比非常大，求解的方程数量也大。研究等几何分析后续需要对数值计算算法进行更深入的研究，保证在提高精度的同时计算速度快。