一种等几何形状优化的网格更新方法

□ 李 杨 □ 张卫红 □ 蔡守宇

西北工业大学 现代设计与集成制造技术教育部重点实验室 西安 710072

结构形状优化旨在通过修改结构内外边界几何形状，改善结构特性（如降低应力集中、提高刚度、减轻结构质量），主要包括三步：参数化几何模型的建立；结构响应分析和灵敏度计算；优化算法的实现。在优化过程中，用于描述结构几何形状的参数作为设计变量，结构质量与响应（如应力、柔顺度、频率等）作为目标函数和约束。

传统的形状优化方法采用计算机辅助几何设计（CAD）进行几何建模，即利用Bezier、B样条或非均匀有理B样条（NURBS）参数化描述。由于CAD模型无法直接用于有限元分析（FEA），需要划分有限元网格后 进 行 分 析 和 优 化［1］， 几 何 模 型 （CAD） 和 分 析 模 型（FEA）相互孤立，如图 1（a）所示。

IGA方法将NURBS基函数不仅用于参数化几何建模，而且用于力学场计算［2］，实现了几何模型（CAD）、分析模型（IGA）和优化模型的无缝连接，如图1（b）所示。由于CAD模型直接用于等几何分析，几何形状的精确表示使分析精度提高。目前，基于IGA的形状优化已被成功应用于平面线弹性问题［3,4,5］、三维线弹性问题［6］、梁结构［7］和振动腔问题［8］。 基于 T 样条和Coons曲面的几何参数化方法也被应用于等几何形状优化中［9,10］。

▲图1 优化过程的比较示意图

然而，优化迭代过程中NURBS控制点的大幅变化，时常使相邻控制点距离过近或过远，甚至出现网格重叠畸形，需要在迭代过程中进行网格更新。文献［3］采用内部控制点的位置与边界控制点的位置成固定比例的方法控制网格，这种方法虽然简单但缺乏通用性。文献［9］采用虚拟位移载荷法，把边界控制点的变化量作为一个弹性辅助系统的位移边界载荷处理，通过求解相应的线性方程组得到内部控制点的变化量。方法较为通用，但是每次迭代之后相当于要额外进行一次有限元分析，效率较低。

本文借鉴文献 ［11］空气动力学中的网格调整策略，采用基于NURBS控制点间距的网格更新方法，进行IGA形状优化。首先简要介绍IGA方法、相应的形状优化模型，然后建立网格更新方法，最后给出平面线弹性形状优化算例。

1 IGA等几何分析方法

对于平面线弹性问题的IGA离散形式为：

式中：Ku=F为IGA离散的平衡方程；Ω为弹性体几何域；Γ为面力边界；K为总刚度矩阵；F为总载荷向量；u为节点位移向量；fb为体力载荷向量；fs为面力载荷向量；D 为弹性矩阵；为基函数矩阵，

上述基函数R可由双变量NURBS基函数构造得到：

式中：ξ、η 为 NURBS 节点参数；p、q 为基函数阶数；n、m 为基函数个数；ωi，j为权值。

类似于有限元分析，总刚度矩阵和总体载荷向量可以由单元刚度矩阵和单元载荷向量组装而成。由此可见，等几何分析得到的离散方程，形式上与有限元分析相同：NURBS控制点对应有限元分析中的节点，NURBS物理空间网格对应有限元分析中的网格，NURBS控制点上的位移对应有限元分析中的节点位移等。

2 基于等几何分析的形状优化

2.1 形状优化模型

形状优化就是通过改变连续体结构内外边界形状以改善结构特性（如减轻结构质量、降低应力集中等）的一种优化过程。基于IGA的形状优化问题一般可以描述为：

式中：x为与边界NURBS控制点相关的设计变量；ND为设计变量个数；f和gj分别为目标函数和约束；NC为不等式约束个数。

2.2 网格更新方法

基于IGA的形状优化在迭代过程中，常会出现设计变量的大变形。如图2（a），优化迭代之前NURBS控制网格中控制点P在Y的正方向有大变形，那么迭代之后可能会出现网格重叠［如图2（b）］。这样变形之后的网格是畸形网格，无法继续进行分析。为保证网格的协调性，采用一种基于控制点间距的网格更新方法。

▲图2 NURBS控制网格大变形示意图

考虑四分之一厚壁圆筒截面的NURBS控制点网格，如图3所示。选择实心圆代表的控制点作为设计变量，仅限制其在虚线（边界法向）的方向上扰动。为方便说明，将图中所示控制网格中的虚线上“一列”提取出来说明。 P0，P1，...，Pm为 m+1 个控制点的位置，其中，P0为边界设计变量控制点（Variable Control Point），Pm为固定控制点（Fix Control Point），其余为连接控制点（Linked Control Point）。ΔP0，ΔP1，...，ΔPm为对应的控制点变形量，di=dist（Pi，Pi-1）（1≤i≤m）为相邻两控制点之间的距离。建立如下关系：

▲图3 网格更新方法示意图

故有：从而：

下面给出式（7）的3条性质：

性质 1：当 i=m 时，ΔPm=CmΔP0=0，这与 Pm为固定控制点相符。

性质 2：若 ΔPi1=Ci1ΔP0，ΔPi2=Ci2ΔP0，则有当 i2＞i1 时，Ci2＜Ci1，所以 ΔPi1＞ΔPi2，说明:距离边界可变控制点P0越近的点，变形量越大；距离P0越远的点，变形量越小。

性质 3：当 i2＞i1时，Pi2＞Pi1， 考察扰动后两点Pi1+ΔPi1，Pi2+ΔPi2的位置关系，有：

又 ΔP0＜dsum，故（Pi2+ΔPi2）-（Pi1+ΔPi1）＞0 即 Pi2+ΔPi2＞Pi1+ΔPi1。也就是任意两控制点在网格更新之后，都保持变形前的拓扑关系，不会出现网格重叠现象。

3 数值算例

基于上述的网格更新方法，下面给出等几何形状优化的数值算例。优化算法采用移动渐近线算法（MMA）［12］。 收敛准则为：

▲图4 四分之一带孔平板初始NURBS模型

▲图5 优化前后von-Mises应力云图对比

式中：fi为第i次迭代之后的目标值；f0为目标初始值。

本文中控制点权值在优化过程中保持不变。考虑平板应力集中问题，由于对称性，仅考虑四分之一带孔平板。问题初始模型定义如图4所示，板两边受均布拉力，杨氏模量为210 GPa，泊松比为0.3，平板厚度为1 mm。在体分比99%的约束下，优化目标为最大von-Mises应力最小。初始形状下von-Mises应力云图如图5（a）所示，优化后形状和von-Mises应力云图如图5（b）所示。初始最大von-Mises应力为27.245 2 MPa，优化后为15.744 2 MPa，应力集中明显降低。

4 结论与展望

等几何分析通过NURBS将CAD、FEA和结构形状优化结合起来，省去了繁琐的模型转换，拥有诸多优点，有着极大的发展前景。此外，等几何分析对于解决一些高阶问题（如板壳问题）具有独特的优势。本文采用一种基于NURBS控制点间距的网格更新方法，成功应用于等几何形状优化。数值算例表明，此方法简单有效，下一步将结合此方法进行三维的等几何形状优化研究，并进一步尝试等几何形状优化在板壳问题中的应用。

［1］ V Braibant and C Fleury.Shape Optimal Design Using B-splines ［J］.Computer Methods in Applied Mechanics and Engineering,1984,44：247-267.

［2］ Hughes TJR,Cottrell JA,Bazilevs Y.Isogeometric Analysis：CAD,Finite Elements,NURBS,Exact Geometry and Mesh Refinement ［J］.Computer Methods in Applied Mechanics and Engineering,2005,194（39-41）： 4135-4195.

［3］ Wall WA,Frenzel MA,Cyron C.Isogeometric Structural Shape Optimization ［J］.Computer Methods in Applied Mechanics and Engineering,2008,197（33-40）： 2976-2988.

［4］ Cho S,Ha S-H.Isogeometric Shape Design Optimization：Exact Geometry and Enhanced Sensitivity［J］.Structural and Multidisciplinary Optimization,2009,38（1）： 53-70.

［5］ Xiaoping Q.Full Analytical Sensitivities in NURBS Based Isogeometric Shape Optimization ［J］.Computer Methods in Applied Mechanics and Engineering,2010,199 （29-32）：2059-2071.

［6］ Hassani B,Tavakkoli SM,Moghadam NZ.Application of Isogeometric Analysis in Structural Shape Optimization ［J］.Scientia Iranica,2011,18（4）： 846-852.

［7］ Nagy AP,Abdalla MM,Gürdal Z.Isogeometric Sizing and Shape Optimisation of Beam Structures ［J］.Computer Methods in Applied Mechanics and Engineering,2010，199（17-20）： 1216-1230.

［8］ Manh ND,Evgrafov A,Gersborg AR,et al.Isogeometric Shape Optimization of Vibrating Membranes ［J］.Computer Methods in Applied Mechanics and Engineering,2011,200（13-16）： 1343-1353.

［9］ Ha S-H,Choi K,Cho S.Numerical Method for Sshape Optimization Using T-spline Based Isogeometric Method［J］.Structural and Multidisciplinary Optimization,2010,42 （3）：417-428.

［10］ Qian X,Sigmund O.Isogeometric Shape Optimization of Photonic Crystals via Coons Patches ［J］.Computer Methods in Applied Mechanics and Engineering， 2011，200（25-28）：2237-2255.

［11］ Pandya MJ,Baysal O.Gradient-based Aerodynamic Shape Optimization Using Alternating Direction Implicit Method［J］.Journal of Aircraft,1997,34（3）： 346-352.

［12］ K.Svanberg,The Method of Moving Asymptotes： A New Method for Structural Optimization［J］.International Journal of Numerical Methods in Engineering， 1987，24：359-373.