基于面光滑有限元的复杂三维结构拓扑优化

何智成　陈少伟　李光耀　张桂勇

1.湖南大学汽车车身先进设计制造国家重点实验室，长沙，4100822.大连理工大学工业装备结构分析国家重点实验室，大连，116023

基于面光滑有限元的复杂三维结构拓扑优化

何智成1陈少伟1李光耀1张桂勇2

1.湖南大学汽车车身先进设计制造国家重点实验室，长沙，4100822.大连理工大学工业装备结构分析国家重点实验室，大连，116023

为了增强拓扑优化计算对任意复杂模型的适应性，改进基于线性四面体有限元的拓扑优化结果，引入了一种新型高精度的基于面光滑有限元模型(FS-FEM)来进行拓扑优化，通过每次迭代时提供很好的梯度解及位移解，从而达到改善拓扑优化结果的目的。在基于面光滑有限元模型的拓扑优化中，以柔度最小作为目标函数，建立了基于固体各向同性材料惩罚插值(SIMP)的拓扑优化数学模型，该数学模型通过最优准则进行求解。多个不同载荷的拓扑优化数值算例说明，采用基于面光滑有限元进行拓扑优化，结果都能够单调收敛，且采用该方法建立的拓扑优化模型能抑制棋盘格现象。与商业软件OptiStruct的计算比较表明，该方法相比有限元方法能得到更合理的拓扑结构。

拓扑优化；光滑伽辽金；面光滑有限元；四面体单元

0　引言

在产品的开发中，保证产品性能的同时减小质量已经成为当前航空航天、机械、汽车等领域中一个重要的研究方向，因而优化技术已经广泛应用到各个领域。随着工程结构或产品越来越复杂，现有的优化技术一般都在基于离散单元的基础上进行，如拓扑优化主要基于现有的有限元技术，通过在给定的设计空间施加载荷及边界条件来寻找材料的最佳分布，使得最终设计的零部件或产品具有很好的性能以及较小的质量。在过去的几十年中，拓扑优化取得了非常大的发展。目前常用的拓扑优化方法有基于均匀化的方法[1]、基于密度的方法[2-3]、水平集法[4]、独立连续体映射模型方法[5]等。其中，变密度法具有较高的计算效率，在工程中应用最为广泛。其他结构拓扑优化方法还包括基于进化策略的进化结构法[6]等，Rozvany[7]对其进行了系统的综述。近年来，高阶单元法[8]以及无网格法[9-10]也被用于连续体结构的拓扑优化研究，研究结果表明，新型高精度数值方法能够消除棋盘格现象，使边界变得较为光滑，能很好地改进拓扑优化的结果。

在基于有限元的拓扑优化中，四面体能够对任意复杂的工程结构进行很好的网格剖分，适用于任意复杂工程结构的计算，然而由于线性四面体单元的刚度过大，在计算位移场和梯度场时存在较大的数值误差，故在拓扑优化中，不能得到很好的拓扑优化结构。另外，由于拓扑优化得到的零部件的材料分布或形状不规则、边界不光滑[11]以及存在网格依赖性(mesh-dependency)和棋盘格式(checkerboards)等数值计算缺陷[11]，难以满足工程制造的要求，故拓扑优化仍然存在着很大的挑战。现有基于离散单元的拓扑优化技术基本都是基于三维六面体单元的研究，然而六面体对复杂零件的几何形状表征能力差，不能对复杂的三维工程结构进行离散，因而，基于六面体单元的拓扑优化不能对复杂三维模型进行优化。

为了增强拓扑优化计算对任意复杂模型的适应性，改进基于线性四面体有限元的拓扑优化结果，本文通过引入一种新型高精度的基于面的光滑有限元模型(face-based smoothed finite element method, FS-FEM)[12-13]来进行拓扑优化，通过每次迭代时提供很好的梯度解及位移解来达到改善拓扑优化结果的目的。本方法通过对单元内的应变分片进行光滑操作，并形成刚度矩阵，从而使得整个模型的刚度趋于柔软，恰好部分地抵消传统线性单元过于刚硬带来的误差，从而提高了单元计算的精度，因而该方法对于拓扑优化具有很好的应用前景。本文在基于面的光滑有限元拓扑结构优化中，以柔度最小为目标函数，以基于面的光滑域作为设计的变量，并采用固体各向同性材料惩罚插值(solid isotropic material with penalization,SIMP)方法来建立拓扑优化的数学模型，通过优化准则来进行求解。最后通过三个不同载荷工况的数值算例来验证本方法的有效性。

1　三维弹性力学方程

对于连续材料的模型,其弹性介质的小变形控制方程可以使用下式进行表示:

LTσ+b=0

(1)

其中，σ为应力张量,b为体力向量, L为矩阵的微分算子。对于各向同向材料，应力和应变具有如下关系:

σ=D ε

(2)

其中， D为材料矩阵;ε 为应变张量，且与位移u具有如下几何关系：

ε=L u

(3)

一般的固体力学问题具有两类边界条件，一类是Dirichlet边界条件，另一类是Neumann边界条件。

(1)自然边界条件可以表示为

(4)

(2)位移边界条件可以表示为

(5)

2　基于面的光滑应变技术

在FS-FEM的构造中, 问题域中的应变不再采用有限元中的协调应变形式ε=L u, 而是采用光滑应变的形式。该光滑应变是通过在光滑域内对不连续的应变进行梯度光滑操作而形成的。通过对所有光滑域的组装，即

图1　基于面的三维光滑域

(6)

(7)

(8)

(9)

(10)

(11)

(12)

3　基于面的光滑有限元模型的建立

在弹性力学中，假定整个的问题域为Ω，其边界为Γ(Γ=Γu+Γt)，其广义光滑Galerkin的弱形式可以写成如下形式[13]：

∫ΓtδuTtΓdΓ=0

(13)

其中，δ为变分因子。将式(6)、式(8)代入式(13),最后的离散线性方程组可以写成如下形式:

(14)

(15)

(16)

(17)

(18)

式(18)是基于线性四面体的常应变单元来进行推导的(应变矩阵为常数)。由于应变矩阵是基于有限元的形函数进行的推导，故可以看出，三维FS-FEM的刚度矩阵为对称稀疏矩阵，可以使用通用的求解器进行线性方程组的高效率求解。

4　基于面光滑有限元的拓扑优化模型

4.1基于SIMP材料插值的拓扑优化模型

近年来，基于刚度-密度插值的模型由于实现简单，因而在拓扑优化中应用得最为广泛。其中，基于正交各向同性材料密度幂指数形式的变密度法材料密度插值理论方法是最常见的一类刚度-密度插值模型，本文也基于该类模型来实现拓扑优化。基于SIMP材料插值模型假设材料的弹性模量Ep可以表示成如下形式[11]：

E(x)=ρ(x)pE0

(19)

其中，ρ(x)为设计变量，p为幂指数惩罚函数，通过引入惩罚因子，对中间密度进行惩罚，尽量使优化过程中的密度趋近“0”(表示无材料)或“1” (表示有材料)，E0为给定各向同性材料的弹性模量，E为插值的材料弹性模量。

在基于面光滑的有限元中，由于以四面体的面光滑域作为组装整体刚度矩阵的基本单位，故可以选取基于面光滑域的密度作为设计变量。因此，在基于面光滑有限元模型的拓扑优化中，以柔度最小(即结构应变能最小，刚度最大)作为目标函数，以结构的体积比为约束条件，则基于SIMP模型的拓扑优化数学模型可以表示为

(20)

其中，目标函数c为结构的总柔度，由各个光滑域相加得到。V0、V分别为优化前和优化后设计域内的材料体积；f为体积分数；ρ为设计变量，表示面光滑域的密度，取值范围为0～1，但是为了保证数值分析的稳定性，避免出现奇异矩阵，通常选取一个非常小的值作为设计变量的下限ρmin，本文取0.01。

建立拓扑优化的数学模型后，可以采用多种优化方法对其求解，如优化准则法 (optimalitycriteria,OC)[11]、序列线性规划法(sequentiallinearprogramming,SLP)[14]以及移动近似算法(MMA)[15]等。OC算法是一种间接优化算法，它并不是直接对目标函数进行优化求解，而是利用Kuhn-Tucker条件作为最优解的满足准则，其算法收敛速度快，迭代次数少，这对于设计变量多、连续体拓扑结构的复杂优化求解是一种高效率高精度的方法[3]。基于SIMP材料插值模型的优化准则法的迭代格式如下：

(21)

其中，η(0<η<1)为阻尼因子，m为移动极限，用来控制迭代的步长。通过在迭代过程中引入这两个参数来保证设计变量迭代过程稳定，参数B通过OC准则得到，其表达式为

(22)

其中，λ为拉格朗日乘子，通过二分法得到。其迭代中止条件为

‖(ρnew-ρ)/ρnew‖≤ξ

(23)

4.2灵敏度分析

由式(22)可知，采用优化准则求解优化问题时，需要用到目标函数关于设计变量的导数，即目标函数对设计变量的灵敏度。通过对柔度求导，得到

(24)

(25)

将式(25)代入式(24)得

(26)

这样，目标函数的灵敏度就转化为求解刚度矩阵关于设计变量的灵敏度。

4.3敏度过滤及密度插值

在基于有限元的拓扑优化中，普遍存在网格依赖和棋盘格等数值不稳定的现象，尤其在低阶的单元中[11]。为了保证拓扑优化得到较好的数值解，许多学者对其进行了研究。许多研究表明，基于敏度过滤的方法能够很好地解决网格依赖和棋盘格问题[11]。敏度过滤的主要思想是在求解某个指定单元灵敏度时，取其周围相邻单元的灵敏度值加权平均。由于敏度过滤的方法没有增加额外的约束，故在计算量上不会造成大的负担，简单易用。敏度过滤是通过在迭代中改变单元的灵敏度来实现过滤的。而在基于面光滑有限元的拓扑优化中，在求解某个光滑域的灵敏度时，取相邻光滑域的灵敏度进行平均，公式为

(27)

其中，M是定义敏度过滤的半径内的所有光滑域的个数；Hf是权系数，也是卷积运算符表达式，即

(28)

其中，rmin是敏度过滤半径，|xi-xf|表示光滑域i的中心点与光滑域f中心点之间的距离(若其距离超出了过滤半径，则权重为0)。在敏度过滤技术中，灵敏度的权重是根据两个单元之间的距离来控制的，过滤半径大小的选择在一定程度上也会影响过滤效果。在过滤技术中，对优化结果的控制主要是通过调节过滤半径rmin来实现的，通过选择合适的过滤半径rmin即可得到稳定收敛的优化结果。

在FS-FEM的拓扑优化中，由于计算得到的是面光滑域的密度，笔者又引入一层密度过滤技术，将面光滑域的密度转化为单元的节点密度，则节点密度与面光滑域密度的关系为

(29)

其中，ρi为节点i的密度,n为包含节点i的面光滑域的个数，ρk为以面k为中心的光滑域的密度。

5　数值算例

本节采用3个三维拓扑优化算例来对本文提出的方法和模型进行验证。算例一为单点载荷下的三维结构拓扑优化，是拓扑优化的经典算例。由于在实际的优化中，一般结构都具有多个载荷，故算例二针对结构的多载荷工况进行了拓扑优化；另外由于在实际的应用中，结构和载荷一般都比较复杂，故算例三也针对复杂结构复杂载荷下的汽车悬架下摆臂进行了研究。在本文前两个算例的模型中，材料的弹性模量均为1Pa,泊松比为0.3，模型采用通用的线性四面体网格进行离散，算例的计算在MATLAB软件中编程实现。

5.1基于单点载荷三维悬臂梁模型

图2所示为三维悬臂梁，其尺寸为8m×4m×8m，其左端面约束，其右端面下边界中间受到单个集中载荷的作用，其力的大小为1N。该设计域被离散成1920个节点、8071个四面体单元，取体积分数的约束为0.4。采用FS-FEM对其进行拓扑优化研究。首先研究惩罚因子p对其收敛性的影响。图3为惩罚因子分别为2、3、4的收敛曲线图，可以看到当p为3和4时收敛较快，迭代43次后收敛，p值取2时，迭代次数要比惩罚因子为3和4时的迭代次数大很多。这说明惩罚因子取p≥3比较合适。

图2　三维悬臂梁模型示意图

图3　基于单点载荷拓扑优化的目标函数收敛示意图

图4a为收敛后的拓扑结构示意图，为了便于比较，采用商业计算软件OptiStruct进行计算，其结果见图4b。两者结构类似，但是根据文献[1]的二维类似问题比较，发现基于面光滑有限元的结果更接近参考解的拓扑结构。

(a)FS-FEM的优化计算结果(b)OptiStruct的优化计算结果图4　拓扑优化计算结果比较

5.2基于多载荷的三维拓扑优化设计

在实际应用中，大多数工程结构都工作在多载荷工况下，因此，对多载荷工况下连续结构拓扑优化问题进行研究具有现实的工程应用价值。图5所示为三维立方体模型，其尺寸为100m×100m×100m，其下端4个顶点约束，在顶面内受到垂直4个集中载荷的作用，其力的大小为1N。采用FS-FEM对其多载荷的拓扑结果进行优化。首先将该设计域离散成7284个节点、35 648个四面体单元，取体积分数的约束为0.4。图6为采用FS-FEM计算得到的目标函数的收敛曲线图，可以看出该模型是单调收敛的。

图5　多载荷模型示意图

图6　多载荷拓扑优化的目标函数收敛示意图

图7为收敛后的的拓扑结构示意图，为了便于比较，图8、图9所示分别为采用商业计算软件OptiStruct计算的结果以及文献[16]提供的结果。可以看出，采用FS-FEM计算的结果比OptiStruct计算的结果更加接近文献[16]的结果。

图7　FS-FEM的优化计算结果

图8　OptiStruct的优化计算结果

图9　OOlhoff等[16]计算的结果

5.3汽车悬架下摆臂的拓扑优化

近年来，轻量化是汽车的结构设计非常关注的研究方向。汽车零部件的轻量化是指在不影响零部件刚度等性能的基础上，通过设计质量小的产品达到降低汽车制造成本的目的。由于拓扑优化能够在给定的边界条件下得到最优的材料分布，因而在汽车中的应用越来越广泛。本节针对具有复杂载荷的汽车悬架下摆臂进行研究，来验证FS-FEM在复杂汽车零部件设计中的应用。图10为汽车下摆臂的设计空间示意图。由于下摆臂的左端的两个铰链与副车架相连，边界条件为约束其中一个铰链的xyz自由度，以及约束另外一个铰链的yz自由度，另外，由于下摆臂的右端与车轮相连，需要承受来自x向、y向、z向的力，因此在铰链中心的x向、y向、z向施加1kN的作用力。该下摆臂的材料为压铸铝，其材料参数如下：弹性模量为79GPa，泊松比为0.3。

图10　汽车悬架下摆臂设计空间及边界条件示意图

首先采用FS-FEM方法来对其进行拓扑优化研究，该模型离散成2372个节点、9819个四面体单元，其体积分数约束为0.4。图11为采用FS-FEM计算得到的目标函数的收敛曲线，可以看出该模型单调收敛。图12为收敛后的拓扑结构示意图，为了便于比较，采用商业计算软件OptiStruct分别对离散的四面体模型和六面体模型进行计算，其结果分别在图13、图14中给出。可以看出，采用四面体有限元计算的拓扑结构虽然为稳定的三角形结构，但是边界不光滑，优化后的拓扑结构不对称；而采用FS-FEM计算的拓扑结果与OptiStruct采用六面体的计算结果非常类似，为稳定的三角形结构，并且FS-FEM计算的拓扑边界比六面体网格的计算结果具有较好的连续性，这说明FS-FEM能很好地对复杂零部件进行拓扑优化。

图11　汽车悬架下摆臂拓扑优化的目标函数收敛示意图

图12　FS-FEM的优化计算结果

图13　基于四面体离散的OptiStruct优化计算结果

图14　基于六面体离散的OptiStruct优化计算结果

6　结论

(1)本文方法采用对问题域适应性很好的四面体网格来进行离散，对任何复杂的几何模型具有很好的适应性，适用于任意复杂三维结构的拓扑优化。

(2)基于面光滑的有限元能够提供很好的梯度解，因而进行拓扑优化计算的结构都能够单调收敛。

(3)采用本方法建立的拓扑优化模型能抑制了棋盘格现象，通过与商业软件OptiStruct的计算结果进行比较，发现相比有限元法，本文方法能得到更合理的拓扑结构，边界具有更好的连续性，具有很好的应用前景。

[1]BendsoeMP,KikuchiN.GeneratingOptimalTopologiesinStructuralDesignUsingaHomogenizationMethod[J].Comput.MethodsAppl.Mech.Eng.，1988，71(2):197-224.

[2]BendsoeMP,SigmundO.MaterialInterpolationSchemesinTopologyOptimization[J].Arch.Appl.Mech., 1999,69(9/10):635-654.

[3]罗震.基于变密度法的连续体结构拓扑优化设计技术研究[D].武汉:华中科技大学, 2005.

[4]ChallisJV.ADiscreteLevel-setTopologyOptimizationCodeWritteninMATLAB[J].Struct.Multidisc.Optim.,2010,41(3):453-464.

[5]隋允康, 杨德庆.多工况应力和位移约束下连续体结构拓扑优化[J].力学学报, 2000, 32(2): 171-179.

SuiYunkang,YangDeqing,WangBei.TopologicalOptimizationofContinuumStructurewithStressandDisplacementConstraintsunderMultipleLoadingCases[J].ActaMech.Sin.,2000,32(2): 171-179.

[6]XieYM,StevenGP.ASimpleEvolutionaryProcedureforStructuralOptimization[J].Comput.&Struct., 1993,49(5):885-896.

[7]RozvanyGIN.ACriticalReviewofEstablishedMethodsofStructuralTopologyOptimization[J].Struct.Multidisc.Optim., 2010,37(3):217-237.

[8]BruggiM.OntheSolutionoftheCheckerboardProbleminMixed-FEMTopologyOptimization[J].Comput.&Struct., 2008,86(19/20):1819-1829.

[9]ZhengJ,LongSY,XiongYB,etal.ATopologyOptimizationDesignfortheContinuumStructureBasedontheMeshlessNumericalTechnique[J].CMES:Comput.Model.Eng.&Sci., 2008,34(2):137-154.

[10]ZhouJX,ZouW.MeshlessApproximationCombinedwithImplicitTopologyDescriptionforOptimizationofContinua[J].Struct.Multidisc.Optim., 2008,36(4):347-353.

[11]Bends∅eMP,SigmundO.TopologyOptimizationTheory,MethodsandApplications[M].NewYork:Springer, 2003.

[12]Nguyen-ThoiT,LiuGR,LamKY,etal.AFace-basedSmoothedFiniteElementMethod(FS-FEM)for3DLinearandGeometricallyNonlinearSolidMechanicsProblemsUsing4-nodeTetrahedralElements[J].Int.J.Numer.Mech.Engrg.,2009,78(3):324-353.

[13]LiuGR.AGeneralizedGradientSmoothingTechniqueandtheSmoothedBilinearFormforGalerkinFormulationofaWideClassofComputationalMethods[J].Int.J.ComputationalMethods, 2008,5(2): 199-236.

[14]FujiiD,KikuchiN.ImprovementofNumericalInstabilitiesinTopologyOptimizationUsingSLPMethod[J].Struct.Multidisc.Optim., 2000,19(2):113-121.

[15]SvanbergK.TheMethodofMovingAsymptotes-aNewMethodforStructuralOptimization[J].Int.J.Numer.Meth.Engrg., 1987,24(2):359-373.

[16]OlhoffN,RonholtE,ScheelJ.TopologyOptimizationofThreeDimensionalStructuresUsingOptimumMicrostructures[J].Struct.Optim., 1998,16(1):1-18.

(编辑陈勇)

Topology Optimization Using FS-FEM for Complex Three-dimensional Models

He Zhicheng1Chen Shaowei1Li Guangyao1Zhang Guiyong2

1.State Key Laboratory of Advanced Design and Manufacturing for Vehicle Body,Hunan University,Changsha,410082 2.State Key Laboratory of Structural Analysis for Industrial Equipment,Dalian University of Technology,Dalian,Liaoning,116023

The FS-FEM was proposed recently based on tetrahedral mesh, and showed great property in solid mechanics, such as providing very good gradient solutions. The topology optimization design of the continuum structures under static load was formulated in the scheme of FS-FEM. As the face-based smoothing domains were the sub-units of assembling stiffness matrix in the discretized system of smoothed Galerkin weak form, thus the relative density of face-based smoothing domains were served as design variables. In this formulation, the minimization compliance was considered as an objective function, and the mathematic of the topology optimization model was developed using the solid isotropic material with penalization (SIMP) interpolation scheme. The topology optimization problem was then solved by the optimality criteria method. Finally, the feasibility and efficiency of the proposed method were illustrated with several 3D examples that were widely used in the topology optimization design.

topology optimization; smoothed Galerkin weak form; face-based smoothed finite element method (FS-FEM);tetrahedral mesh

2014-03-05

国家自然科学基金资助项目(11202074)；湖南大学汽车车身先进设计制造国家重点实验室自主研究课题(51375001)；工业装备结构分析国家重点实验室开放基金资助项目(GZ1403)

TH122DOI：10.3969/j.issn.1004-132X.2015.07.003

何智成，男，1983年生。湖南大学汽车车身先进设计制造国家重点实验室助理研究员。主要研究方向为汽车振动噪声、汽车轻量化。陈少伟，男，1981年生。湖南大学汽车车身先进设计制造国家重点实验室博士研究生。李光耀，男，1963年生。湖南大学汽车车身先进设计制造国家重点实验室主任、教授、博士研究生导师。张桂勇，男，1979年生。大连理工大学船舶学院教授、博士研究生导师。