0 引言
根据产品设计阶段的不同，可以将结构优化分为3个层次：拓扑优化、形状优化和尺寸优化[1]。拓扑优化是在给定的材料设计域内，通过改变拓扑结构来确定最佳的产品雏形，本质上是最佳传力路径的确定。形状优化位于拓扑优化的下一阶段，对局部形状进行优化，以达到减小局部应力、减重的目的[2]。尺寸优化是对离散尺寸变量的优化，以确定出产品最终的精确尺寸模型。在优化边界描述上，早期方法采用节点坐标作为设计变量，如CASSIS等[3]研究了大坝的形状优化问题，设计变量取有限元网格节点，优化算法选用序列线性规划法。BRAIBANT等[4]提出利用B样条来优化边界，并分析了单元数量和B样条阶次对优化结果的影响。QIAN[5]提出了采用NURBS样条描述优化边界的方法，更为精确地描述了边界形状。
基于体参数化模型的等几何分析方法最早由HUGHES等[6] 于2005 年提出。与其他方法相比，等几何分析方法应用到产品设计具有较多优势[7]。目前，等几何分析的研究主题主要有：等几何分析计算基础理论，包括计算精度、收敛性、奇异性等问题[8-9]；等几何分析方法的改进，如边界元法、配点法等[10-11]；基于各种样条基函数来实现等几何分析[12-13]；面向各种实际应用展开等几何分析的应用[14-15]。基于等几何分析方法的等几何优化方法近年来得到了较多重视[16]，这受益于样条在表达形状方面的优越性。目前，基于等几何分析的形状优化研究最多，且已取得了较大进展[17-18]。有学者在拓扑优化方面也进行了一些尝试[19-20]，而有关基于等几何分析方法的尺寸优化方面的研究报道较少[21-22]，主要原因是样条表达的形体尺寸优化不是非常直观，同时用样条来表达的形体大多是自由形体，自由形体模型本身的参数化造型也是较为困难的问题。
由于等几何优化采用NURBS样条或B样条作为基函数，这会提高计算域的连续性，因此，灵敏度分析成为了等几何优化问题中的重要问题[5]。BRAIBANT等[23]将B样条用于设计域的描述中，采用等几何分析方法进行结构响应的计算，并推导出了二维线弹性问题的灵敏度求解方程。BLANCHARD等[24]将有体积约束的最优化问题通过拉格朗日方程转化为无约束最优化问题进行求解，并对该无约束方程进行了灵敏度求解，给出了弹性优化问题的验证实例。张升刚等[25]在等几何分析框架中推导出了壳体形状优化中优化方程对控制点坐标的灵敏度方程。陈龙等[26]提出将B样条控制点作为设计变量以及将几何模型边界尺寸作为设计参数来表示出控制点，从而推出了尺寸参数相对于控制点的灵敏度，且给出了三维形状优化实例。
上述方法大多将控制点作为设计变量进行优化，并未对权重给予足够的重视，当结构具有圆锥曲线的形状表达需求时就无能为力，而权重在模型结构的表达能力上至关重要。基于上述考虑，本文在推导权重灵敏度的基础上，提出了两种形状优化方式：一种为同时考虑控制点和权重的全局形状优化，一种为基于尺寸优化结果再利用权重进行局部形状优化。本文阐述了等几何分析的基础理论，给出了优化方程并推导出了以控制点和权重同时为优化变量的灵敏度方程，最后通过实例对两种形状优化方法的可行性进行了验证。
1 NURBS表达和等几何优化方法
1.1 NURBS曲线和曲面
NURBS样条是B样条的广义形式。次数为p的NURBS曲线C(s)和两个方向上次数分别为p和q的NURBS曲面S(s,t)定义分别如下：
式中，s、t为参数；Ni,p(s)表示X方向上第i个控制点的p次B样条基函数；Mj,q(t)表示Y方向上第j个控制点的q次B样条基函数；n、m分别为X、Y方向上控制点个数；wi、wi,j为权重;Pi、Pi,j为控制点。
与B样条不同，NURBS具有合理的形状和可变权重，且每个控制点都有相应的权重值，因此，NURBS保留了B样条的所有性质且自身具有一些特有的性质。NURBS具有局部修改性，可以移动控制点Pi或改变权因子wi，仅影响区间s∈[si,si+p+1)上的那部分曲线形状。 由于其合理的形状和权重，故NURBS可以精确地表示圆形、椭圆形等圆锥曲线。
1.2 等几何分析方法
为解决现有CAD系统和CAE系统的无缝集成问题，HUGHES等[6]提出了等几何分析方法，该方法流程见图1。
图1 等几何分析方法流程图
Fig.1 Flow chart of IGA method
根据最小势能及虚功原理[27]，可得到离散化的总体平衡方程：
Ku=F
(3)
式中，u为待求解的控制点位移矢量；K为刚度矩阵；F为外力矩阵。
本文主要研究的是二维结构，已知二维刚度矩阵可表示为
式中，Ω为整体结构的参数域；B为应变矩阵；RI(I=1,2,…,H)为有限元分析中的形函数，在等几何分析中表示NURBS基函数；H为基函数个数；x、y分别为控制点的横坐标和纵坐标；D为弹性系数矩阵；Ep为弹性模量；ν为泊松比；J为雅可比矩阵，可从物理域转换至参数域。
1.3 优化问题及优化流程阐述
在一定面积约束下，以最小结构柔度为优化目标，该结构优化问题的数学表达式如下：
式中，f(·)表示柔度方程；b为优化设计变量矩阵，包括控制点和权重；u(b)表示受设计变量影响的位移矢量函数；S为模型的面积；S*为模型的最大面积值；bimin、bimax分别为设计变量bi的最小值和最大值。
本文采用灵敏度分析的结构优化方法建立样条模型并进行了优化问题的前处理和定义，对优化变量、约束和目标等进行了定义，并采用移动渐近线法(method of moving asymptotes,MMA)[28]对模型进行优化迭代从而得到了最终优化结果，优化流程见图2。
图2 基于等几何分析的形状优化方法流程图
Fig.2 Flow chart of shape optimization method based on IGA
根据NURBS曲线的局部修改性，可以移动控制点或改变权重来实现形状的局部修改。如图3所示， 其中w1=1，w2=0.707 1，w3=0.1。可以看出，若权重增大，则NURBS曲线靠近第②个控制点；若权重减小，则NURBS曲线远离第②个控制点；NURBS曲面同理。由此可知，权重对样条几何模型有很好的形状修改效果。
图3 不同权重的NURBS曲线变化图
Fig.3 Diagram of changing from different weights of NURBS curve
特别指出的是，为了达到指定的边界形状，NURBS曲线权重的初始值应设为1[29]。在优化过程中，本文将每个权重函数都限定在[0,1]范围内，即权重满足规范性。
2 基于灵敏度方程的优化方法
2.1 灵敏度分析推导过程
灵敏度分析就是分析输出响应对设计变量变化的敏感性，灵敏度的数值可反映出各设计变量对输出响应变化的贡献率。在形状优化中，根据目标函数得到其灵敏度可表示为
为了计算关于设计变量的柔度形状灵敏度，应计算单元刚度矩阵的灵敏度，其表达式如下：
式中，ig(ig=1,2，…，ng)为高斯点编号；ng为高斯点个数；Wig为第ig个转换后高斯积分权系数；e(e=1,2,…，Ac)为单元编号；Ac为单元总数目。
为计算单元刚度矩阵的灵敏度，需要计算单元应变矩阵和单元雅可比矩阵的灵敏度[5]。此时，与设计变量是控制点位置的情况不同，若以权重为设计变量，则NURBS基函数也是权重的函数。约束函数中的模型面积S对设计变量的导数可表示为
为求出
，需求出单元应变矩阵Be中形函数的偏导数，可根据复合函数的求导法，利用对参数域的偏导可求得
式(9)和式(10)可转换为
已知N可表示出单元应变矩阵Be，由式(11)可得N=J-1M，同时雅可比矩阵可由NURBS曲面公式推导得到：
针对本文提出的局部形状优化方法和全局形状优化方法，需对设计变量为控制点及设计变量为控制点和权重两种情况分别进行灵敏度求导。
2.2 以控制点为设计变量的灵敏度推导
当设计变量为控制点P(x,y)时，N对控制点求导可得
-J-1(M
)N=-NN
(13)
根据二阶行列式求导法则，可求出设计变量为控制点P(x,y)时，单元雅可比矩阵|J|e对控制点的导数为
2.3 以控制点和权重为设计变量的灵敏度推导
当设计变量为权重和控制点P(xw,yw)时，N对权重和控制点求导可得
由式(15)可知，需要求出NURBS基函数对权重的导数，设
单元雅可比矩阵|J|e对控制点和权重的导数可表示为
(17)
综上，得到单元刚度矩阵对设计变量的灵敏度后进而组装为总体刚度矩阵，即可得到总体刚度矩阵对设计变量的灵敏度矩阵，则可求解以柔度最小即刚度最大的目标函数灵敏度分析问题。
3 基于权重灵敏度分析的全局形状优化
首先考虑同时变动控制点及其权重的全局形状优化方法。以扳手为例，在扳手右上端施加FH=10 N的均布牵引力及其所受约束，如图4所示。
图4 扳手边界条件示意图
Fig.4 Diagram of boundary conditions for the wrench
构建扳手的几何模型时，两个方向的节点集合分别为
和
次数分别为p=2和q=1，初始几何参数如下：扳手长度aw=25 mm，扳手宽度bw=12 mm，槽口宽度c=4 mm。材料参数为：弹性模量Ep=150 GPa，泊松比ν=0.3。边界控制点经过孔斯插值后可得到由11×6个控制点组成的扳手模型，如图5所示。
图5 扳手控制点示意图
Fig.5 Schematic diagram of control points for the wrench
优化目标为最小柔度即最大刚度，优化的设计变量为扳手顶面上A、B、C、D、E共5个控制点的垂直位置，后面6个控制点的垂直位置与D点相同。约束条件为面积S≤S*=9 000 mm2，并设置对称约束。得到以控制点为设计变量的最终优化几何形状如图6所示。
(a)扳手控制点图
(b)扳手边界形状图
(c)扳手等参线图
图6 扳手控制点形状优化图(方法1)
Fig.6 Diagram of control points shape optimization for the wrench(method one)
将扳手顶面上的控制点及权重同时作为设计变量，得到的最终优化几何形状如图7所示。
(a)扳手控制点图
(b)扳手边界形状图
(c)扳手等参线图
图7 扳手控制点形状优化图(方法2)
Fig.7 Diagram of control points shape optimization for the wrench(method two)
将采用两种方法得到的目标柔度曲线(图8)的对比结果绘制于表1, 其中方法1是以控制点为设计变量的形状优化结果，方法2是以控制点及权重同时为设计变量的形状优化结果。由表1可知，加上权重的优化后，扳手面积基本不变，且扳手形状优化后的柔度更小。
(a)以控制点为设计变量
(b)以控制点和权重同时为设计变量
图8 扳手的柔度优化过程
Fig.8 Flexibility optimization processes for the wrench
表1 扳手优化对比
Tab.1 Optimization comparison for the wrench
优化对比柔度面积(mm2)方法10.120 78 958方法20.073 98 992
4 基于尺寸优化结果的局部形状优化
4.1 尺寸优化流程
尺寸优化是经典的优化技术，通常也称为参数优化技术，即改变模型参数值时模型发生变化，但无需重新划分网格。在等几何分析中可对几何模型的各个参数(如板件厚度、杆梁截面、尺寸、弹性元件刚度等)用样条几何的参数表征，进而优化模型，如图9所示。
图9 尺寸优化流程图
Fig.9 Flow chart of size optimization
4.2 尺寸优化设计变量的确定
与结构优化的优化模型相同，同样以柔度最小为优化目标，以几何模型的尺寸为设计变量，用NURBS样条几何模型中的控制点表征尺寸信息，即
P(xi,yi)=P(αi)
(18)
其中，αi为第i个控制点的尺寸信息。为表示出灵敏度分析中的
和
需求出
其表达式如下：
将式(19)先代入式(6)再代入式(5)，即可求得柔度最小下的最优尺寸变量结果。
4.3 尺寸优化结果
以工字梁为例，将工字梁截面的7个尺寸参数作为优化设计变量来进行尺寸优化。由工字梁四分之一截面图(图10)可知，工字梁截面的7个尺寸变量分别为：截面高度h，腿宽度b，腰宽度d，腿高度t，腿端圆弧半径r1，内圆弧半径r2，内弧倾斜度θ。
图10 工字梁四分之一截面图
Fig.10 One fourth of the cross section from the I-beam
将图10中的9个控制点用尺寸参数表示，即
构建工字梁的几何模型时，两个方向的节点集合分别为
和
次数分别为p=2和q=2，弹性模量Ep=150 GPa，泊松比ν=0.3。经过孔斯插值后可得到由25×7个控制点组成的工字梁几何模型，下底边固定，在上顶边施加均布力FN=10 N。约束条件为工字梁面积S≤S*=141 800 mm2，如图11所示。
图11 工字梁优化边界条件示意图
Fig.11 Initial optimization boundary conditions of I-beam
经过尺寸优化后的工字梁如图12所示，各参数的优化收敛过程如图13所示，工字梁尺寸优化前后的各参数对比如表2所示。
(a)控制点图 (b)边界形状图 (c)等参线图
图12 工字梁尺寸优化图(方法3)
Fig.12 Diagram of size optimization for the I-beam (method three)
(a)截面高度 (b)腿宽度
(c)腰宽度 (d)腿高度
(e)内圆弧半径 (f)腿端圆弧半径
(g)内弧倾斜度
图13 各参数的优化收敛过程
Fig.13 The optimal convergence process of each parameter
优化过程中，柔度变化如图14所示，可以看出，曲线收敛速度快，柔度数值变化大。
表2 工字梁尺寸优化前后对比
Tab.2 Comparison of I-beam sizes before and after optimization
优化对比柔度面积(mm2)h(mm)b(mm)d(mm)t(mm)r2(mm)r1(mm)θ(rad)优化前3.07141 8001 00068045.076.065330.170优化后2.72141 50098066045.478.860280.192
图14 工字梁柔度优化过程
Fig.13 Flexibility optimization process of I-beam
4.4 基于权重的形状优化方法的数值算例
由于权重具有局部支撑性，因此对于局部应力过大的部位可利用这一特性进行局部形状优化。在相同面积约束下，在优化尺寸以变动控制点达到目标柔度最小的基础上，再以权重为设计变量进行局部形状优化，从而可以得到更小的柔度。
图15 工字梁局部控制点优化示意图
Fig.15 Schematic diagram of local control points for I-beam
对尺寸优化结果之后的工字梁模型进行等几何分析可知，在工字梁内圆弧处的应力集中。将内圆弧周围的控制点权重作为设计变量(图15)，即为A、B、C、D、E、F、G、H的平行与垂直方向及I、J两点的垂直方向位置，不改变约束面积且保持控制点和权重的对称约束，得到的最终工字梁模型的形状优化结果如图16所示，优化过程如图17所示。两种优化方式的对比结果见表3，其中方法3是以控制点为设计变量的尺寸优化结果；方法4是以控制点和权重同时为设计变量的局部形状优化结果。由表3的对比结果可知，加入权重的形状优化结果可使工字梁的柔度更小，优化后的边界形状更加光顺。
(a)控制点图 (b)边界形状图 (c)等参线图
图16 工字梁尺寸优化图(方法4)
Fig.16 Diagram of size optimization for the I-beam (method four)
图17 工字梁目标柔度优化过程图
Fig.17 I-beam target flexibility optimization process diagram
表3 工字梁优化对比
Tab.3 Optimization comparison of I-beam
优化对比柔度面积(mm2)方法32.72141 500方法42.56140 700
5 结论
(1)本文以柔度最小为形状优化目标、面积为约束，进行了以权重和控制点同时为设计变量的灵敏度分析推导。
(2)提出了以权重和控制点共同为设计变量的全局形状优化方式，以及在以尺寸为设计变量的尺寸优化基础上的权重局部形状优化方式。
(3)通过实例对比分析了几何模型的形状优化结果，可以得到两种形状优化方法是切实可行的，且优化结果更佳，权重的优化使得几何模型形状更加光顺。
本文在扳手算例中的形状优化主要针对边界点进行，因此无法控制模型内部控制点的变化，导致其等参线分布不均，进而会影响到等几何分析的结果。此后可在此基础上进行内部控制点的优化。同时本文在优化过程中限制了权重的取值范围[0,1], 后续可扩大权重因子的取值范围。此外，本文算例的模型较为简单，后续可将本文的算法拓展到复杂模型上。
[1] KLARBRING A. An Introduction to Structural Optimization[M]. Berlin: Springer, 2008:15-35.
[2] SOKOLOWSKI J, ZOLESIO J P. Introduction to Shape Optimization[M].Berlin: Springer, 1992:229-263.
[3] CASSIS J H, SCHMIT L A. Optimum Structural Design with Dynamic Constraints[J]. Journal of the Structural Division,1976,102(10):2053-2071.
[4] BRAIBANT V, FLEURY C. Shape Optimal Design Using B-splines[J]. Computer Methods in Applied Mechanicsand Engineering,1984, 44(3):247-267.
[5] QIAN X . Full Analytical Sensitivities in NURBS Based Isogeometric Shape Optimization[J]. Computer Methods in Applied Mechanics and Engineering, 2010, 199(29/32):2059-2071.
[6] HUGHES T J R, COTTRELL J A, BAZILEVS Y. Isogeometric Analysis: CAD, Finite Elements, NURBS, Exact Geometry and Mesh Refinement[J]. Computer Methods in Applied Mechanics & Engineering, 194(39/41):4135-4195.
[7] 吴紫俊, 黄正东, 左兵权, 等. 等几何分析研究概述[J]. 机械工程学报,2015,51(5):114-129.
WU Zijun, HUANG Zhengdong, ZUO Bingquan, et al.Perspectives on Isogeometric Analysis[J]. Journal of Mechanical Engineering,2015,51(5):114-129.
[8] XU G, KWOK T, WANG C C L. Isogeometric Computation Reuse Method for Complex Objects with Topology-consistent Volumetric Parameterization[J].Computer-Aided Design,2017,91:1-13.
[9] CHARAWI L A. Isogeometric Overlapping Schwarz Preconditioners for the Bidomain Reaction-diffusion System[J]. Computer Methods in Applied Mechanics and Engineering,2017,319:472-490.
[10] GU Jinliang, ZHANG Jianming,LI Guangyao. Isogeometric analysis in BIE for 3-D potential problem[J]. Engineering Analysis with Boundary Elements,2012,36(5):858-865.
[11] COOX L, ATAK O, VANDEPITTE D, et al. An Isogeometric Indirect Boundary Element Method for Solving Acoustic Problems in Open-boundary Domains[J]. Computer Methods in Applied Mechanics and Engineering, 2017, 316:186-208.
[12] DORFELMR, BERT J, BERND S. Adaptive Isogeometric Analysis by Local H-refinement with T-splines[J]. Computer Methods in Applied Mechanics and Engineering, 2010,199(5):264-275.
[13] BAZILEVS Y, CALO J A, COTTRELL J A, et al. Isogeometric Analysis Using T-splines[J]. Computer Methods in Applied Mechanics and Engineering, 2010, 199(5):229-263.
[14] LU Jia. Isogeometric Contact Analysis: Geometric Basis and Formulation for Frictionless Contact[J]. Computer Methods in Applied Mechanics and Engineering,2011, 200(5):726-741.
[15] LORENZIS L D, WRIGGERS P, HUGHES T J R. Isogeometric Contact: a Review[J]. GAMM-Mitteilungen,2014,37(1):85-123.
[16] 刘宏亮, 祝雪峰, 杨迪雄. 基于等几何分析的结构优化设计研究进展[J]. 固体力学学报,2018, 39(3):248-267.
LIU Hongliang, ZHU Xuefeng, YANG Dixiong. Research on Structural Shape Optimization Based on the Isogeometric Analysis Method[J].Acta Mechanica Solida Sinica,2018,39(3):248-267.
[17] QIAN X, SIGMUND O. Isogeometric Shape Optimization of Photonic Crystals via Coons Patches[J]. Computer Methods in Applied Mechanics and Engineering,2011,200(25):2237-2255.
[18] KOSTAS K V, GINNIS A I, POLITIS C G, et al. Shape-optimization of 2D Hydrofoils Using an Isogeometric BEM Solver[J]. Computer-Aided Design,2017,82(13):79-87.
[19] QIAN X. Topology Optimization in B-spline Space[J].Computer Methods in Applied Mechanics and Engineering,2013,265:15-35.
[20] PARK J, SUTRADHAR A. A Multi-resolution Method for 3D Multi-material Topology Optimization[J].Computer Methods in Applied Mechanics and Engineering,2015,285:571-586.
[21] HERREMA A J, WIESE N M, DARLING C N, et al. A Framework for Parametric Design Optimization Using Isogeometric Analysis[J]. Computer Methods in Applied Mechanics and Engineering, 2017,316:944-965.
[22] DING C S, CUI X Y, LI G Y. Accurate Analysis and Thickness Optimization of Tailor Rolled Blanks Based on Isogeometric Analysis[J]. Structural & Multidisciplinary Optimization,2016,54(4):871-887.
[23] BRAIBANT V, FLEURY C. Shape Optimal Design Using B-splines[J]. Computer Methods in Applied Mechanics and Engineering,1984,44(3):247-267.
[24] BLANCHARD L, DUVIGNEAU R, VUONG A V, et al. Shape Gradient Computation in Isogeometric Analysis for Linear Elasticity[J]. Numerical Analysis,2014(2):361-367.
[25] 张升刚，王彦伟.等几何壳体分析与形状优化[J].计算力学学报,2014,31(1):115-119.
ZHANG Shenggang, WANG Yanwei. Isogeometric Shell Analysis & Shape Optimization[J].Chinese Journal of Computational Mechanics, 2014,31(1):115-119.
[26] 陈龙, 张高朋. 基于灵敏度分析的三维结构等几何形状优化[J]. 中国机械工程, 2017,28(14):1723-1729.
CHEN Long，ZHANG Gaopeng.Three Dimensional Structural Shape Optimization Based on Sensitivity Analysis Using Isogeometric Analysis[J].China Mechanical Engineering, 2017, 28(14): 1723-1729.
[27] SVANBERG K. The Method of Moving Asymptotes: a New Method for Structural Optimization[J]. International Journal for Numerical Methods in Engineering,1987,24(2):359-373.
[28] 潘振宇,何钢,夏婷，等.平面线弹性问题的等几何形状优化方法[J].计算机辅助工程,2014,23(2):27-30.
PAN Zhengyu, HE Gang, XIA Ting, et al. Isogeometric Shape Optimization Method of Plane Linear Elastic Problem[J].Computer Aided Engineering,2014,23(2):27-30.
[29] PIEGL L A, TILLER W. The NURBS Book[M]. Berlin: Springer,1997:47-79.
