连续纤维增强复合材料结构基频最大化设计1)

程长征 卞光耀 王选,2) 龙凯 李景传 吴乔国

*(合肥工业大学工程力学系,合肥 230009)

†(合肥工业大学土木工程结构与材料安徽省重点实验室,合肥 230009)

**(华北电力大学新能源电力系统国家重点实验室,北京 102206)

引言

随着现代工业的发展,很多新型材料被研制并成功应用于不同领域,其中纤维增强复合材料就是最常用的一种,已在汽车、建筑以及航空航天工业等领域得到广泛应用[1].纤维增强复合材料由基体材料和增强材料两部分组成.与单独组分材料和传统的金属材料相比,纤维增强复合材料在强度、刚度、抗断裂、抗疲劳、热稳定性等诸多方面具备更优良的性能.另一方面,当结构所受的外激励频率与自身的固有频率接近时,就会发生共振现象,共振会导致机器破坏、桥梁倒塌等严重危害.因此,本文开展连续纤维增强复合材料结构基频最大化设计研究具有重要意义.

基频最大化是结构动力学拓扑优化中被广泛关注的优化问题之一.Ma 等[2]采用均匀化方法较早研究了动力学拓扑优化问题.Pedersen[3]通过弓入新的材料刚度插值格式避免了基于SIMP (solid isotropic material with penalization) 方法的频率优化中在低密度区出现的局部模态现象.为了适应实际制造工艺要求,Niu 等[4]采用均匀化方法实现了蜂窝材料结构在宏观和微观两个尺度上同步优化以最大化结构基频,Huang 等[5]基于BESO(Bi-directional evolutionary structural optimization)方法研究了体积约束下的结构一阶频率最大化问题.Xia 等[6]基于水平集方法研究了结构在无阻尼自由振动下的频率最大化问题,得到了光滑的拓扑构型边界,在此基础上,他们进一步考虑了均匀边界腐蚀对动力学优化问题中基频的影响[7].高兴军和马海涛[8]提出了一种将移频与虚假模态识别相结合的通用方法,以消去连续体结构动力拓扑优化中局部模态的不利影响.动力学拓扑优化方面的研究还包括最大化高阶本征频率或最大化两阶本征频率的带宽[9]、最大化受迫振动下的动柔度[10-13]等.此外,还可以将上面度量结构的振动效应的指标作为约束函数来建模[14-15].文献[16]从不同优化方法角度对结构动力学拓扑优化问题做了较为全面的综述.

纤维增强复合材料结构的纤维角度优化问题也是被大家关注的热点问题之一.Pedersen[17]基于OC准则方法(optimality criteria based method) 实现了层压膜结构的厚度和纤维角度优化.Luo 和Gea[18]基于能量的方法以应变能最小化或基频最大化为目标,实现了三维板壳结构的纤维角度优化问题.Setoodeh等[19]基于SIMP 方法研究了单工况和多工况载荷作用下二维复合结构的拓扑构型和纤维角度的同步优化问题,以最大化结构的刚度.Stegmann 和Lund[20]基于多相材料插值思想提出了一种新颖的DMO(discrete material optimization method)实现了一般的层压复合壳的纤维优化.Niu 等[21]将DMO 法拓展到考虑最小声辐射的层压复合板结构的离散纤维角度优化问题中.Gao 等[22]将复合材料结构的离散纤维优化问题转化为可用连续体拓扑优化方法处理的材料选择问题,提出了一种有效的二值参数化列式.段尊义等[23]将连续化惩罚策略与Heaviside 惩罚函数弓入传统DMO 模型中,提高了复合材料单层板结构的纤维角度优化的收敛率.他们还将所提方法拓展到纤维角度和结构拓扑构型一体化优化问题中,以最大化结构的刚度[24].Xia 和Shi[25]基于Shepard 插值方法[26]构造了一个连续分布的纤维场,研究纤维增强复合材料结构的刚度最大化问题.Yan 等[27]以最大结构刚度为目标,基于双向渐进结构优化法提出了一种实现正交材料纤维方向和结构拓扑构型的同步优化.Xu 等[28]实现了考虑载荷不确定性的纤维增强复合结构的纤维方向的鲁棒性优化.Papapetrou等[29]将SIMP 方法和水平集方法结合,实现了纤维角度和结构拓扑的同步优化,以最大化结构的刚度.Luo 等[30]针对纤维增强复合材料的刚度最大化问题,提出了离散-连续的参数化列式,可有效减小优化问题陷入局部优化解的风险.

可以看出,现有文献大部分关注单一均质材料结构的动力学拓扑构型的优化、或单纯考虑纤维角度优化问题,另外,关于纤维增强复合材料结构的拓扑构型和纤维方向的同步优化方面的研究大部分关注刚度最大化优化问题,而在动力学背景下考虑结构拓扑构型和纤维方向同步优化的研究较少.因此,本文旨在提出一种针对连续纤维增强复合材料结构无阻尼自由振动下的基频最大化问题的拓扑优化方法,实现结构拓扑构型与纤维方向的同步优化,以最大化结构的第一阶本征频率.

1 连续纤维增强复合材料结构优化模型

1.1 复合材料结构自由振动有限元分析

考虑二维正交各向异性材料线弹性结构无阻尼自由振动分析的离散型有限元方程表示如下

其中,ui为第i阶特征频率ωi所对应的特征向量.K和M分别为结构的总体刚度阵和质量阵,由对应的单元刚度阵和质量阵组装形成.和分别由式(2)和式(3)确定

其中,h为结构的厚度,ρ0为材料的质量密度,N为双线性等参元形函数组成的矩阵.B为应变位移矩阵.和分别为单元e所对应的单元质量阵和刚度阵.θe为单元e所对应的表征纤维方向的角度设计变量,即建立在单元e上的纤维主方向和垂直于纤维主方向上的局部坐标X′-Y′与全局坐标X-Y之间的旋转角度,如图1 所示.

图1 坐标转换示意图Fig.1 Schematic illustration of coordinate transformation

式(3)中的C为全局坐标下的正交本构材料矩阵,其表达式为

其中,C0为局部坐标下的正交本构材料矩阵,由沿纤维轴向的弹性刚度Ex,纤维横向弹性刚度Ey,剪切模量Gxy,泊松比Vxy,Vyx五个材料常数组成,其表达式为

式中,Ex,Ey,Vxy,Vyx满足VyxEx=VxyEy.

式(4)中的T(θe)为C与C0之间的转换矩阵,其表达式由下式定义

1.2 材料插值格式

为了实现结构拓扑构型与纤维角度的同步优化,这里采用SIMP 方法来实现结构拓扑构型的优化,即给每一个单元弓入一个密度设计变量xe,然后采用SIMP 材料插值模型建立单元刚度阵和质量阵与单元密度设计变量xe之间的关系

其中,Me和Ke分别为弓入表征拓扑构型的密度变量之后的单元质量阵和刚度阵,p和q分别为刚度阵和质量阵的惩罚因子,和分别由式(2) 和式(3)确定.

众所周知,在用SIMP 求解动力学优化问题时,材料低密度区域通常会出现局部模态[3,8-9,31],这是由于在低密度区域(如xe≤0.1),基于q=1 惩罚之后的质量阵与基于p=3 惩罚之后的刚度阵比值过高造成的[3,9].避免局部模态的常用方法有:(1)修改刚度矩阵的惩罚方式[3,5]; (2) 修改质量矩阵的惩罚方式[9-10,12].在这里,为了避免低密度区域中的局部模态,采用文献[9]中的方法,修改质量阵的惩罚方式

这里设置刚度矩阵惩罚因子p=3,对于质量矩阵,当单元密度xe≥0.1 时,设置惩罚因子为1; 当单元密度xe＜0.1 时,设置惩罚因子q=6.通过式(8) 和式(9)可知,通过修改质量阵的惩罚方式可以很好地避免底密度区域质量阵与刚度阵的比值过高的情况,从而有效避免局部模态.

1.3 优化模型

结构拓扑优化旨在满足某些约束条件下寻找最优的材料分布,以期获得最佳的结构性能.与文献[32-34] 关注的优化问题不同,本文研究连续纤维增强复合材料结构在无阻尼自由振动下的基频最大化优化问题,以实现纤维方向角度与结构拓扑构型的同步优化,其对应的数学列式可表示如下

在此优化模型中,有两类设计变量,其中x为表征结构拓扑构型的密度设计变量集合,θ为表征纤维方向的角度设计变量集合,Ne为用来离散设计域的单元个数.ωi为第i阶特征频率,ui为其对应的特征向量.K和M分别为弓入两类设计变量之后结构的总体刚度阵和质量阵.Ve和V0分别为单元e和整个设计域的体积,fv为设计领域中可用实体材料的体积分数比.xmin为很小的正数,作为密度设计变量的下限以防止刚度矩阵奇异,本文取xmin=10-3.本文设置每个纤维角度变量θe的范围为[-2π,2π],宽的纤维角度范围允许角度旋转到最佳位置,可以避免较窄的角度变量范围导致的潜在的局部解问题[35].

2 灵敏度分析与灵敏度过滤

动力学优化问题可能存在重特征根问题,特别是那些设计自由度多、依赖的设计参数多的结构,在优化的过程中可能会出现两阶特征根相等的情况,重特征根会导致优化问题不连续,进而导致多重特征根关于设计变量的敏度难于计算的困难.目前也有较多优秀工作讨论了这一问题,并给出了有效的解决方案[9,36].重特征根问题不是这篇文章的关注点,因此,这里假设所关注的特征根均为单根,即特征根均互不相等,这样可以推导单特征值关于两类设计变量的灵敏度.

2.1 特征值对密度设计变量的灵敏度

式(10)中无阻尼自由振动分析的有限元方程两边同时对密度设计变量xe求导,并化简可得

对式(11)两边同时左乘,再经过简单化简可得特征值λi=关于设计变量xe的灵敏度

其中,∂K/∂xe和∂M/∂xe与材料插值格式有关,具体表达式可根据式(8)和式(9)计算得到

2.2 特征值对角度设计变量的灵敏度

类似于求特征值关于密度设计变量的灵敏度,可以推导特征值关于角度设计变量θe的灵敏度,其表达式如下

值得注意的是,这里假设质量矩阵与纤维方向无关.式(15)中刚度阵K对角度设计变量θe的导数∂K/∂θe可表示为

关于体积约束关于设计变量的敏度可直接求导获得,这里省略其细节推导.

2.3 灵敏度过滤

由于有限元网格离散的影响,基于单元密度的拓扑优化方法经常遭受棋盘格和网格依赖性等数值不稳定性问题.为了获得清晰的黑白设计,本文采用式(17)对式(12)中定义的第一阶特征值(最小特征值)λ1关于密度设计变量xe的敏度进行过滤

其中权重函数Hei定义为

式中,dist(e,i) 为单元e的中心和单元i的中心之间的欧式距离,Nei是指单元中心落在以单元e为中心、半径为rmin的圆形邻域内的单元集合.以过滤后的敏度代替式(12)参入优化迭代.

上述优化问题可以用不同的优化算法来求解,如序列线性规划法(SLP)[37],移动渐近线算法(MMA)[38]等,本文采用MMA 算法求解上述优化问题.数值测试表明角度设计变量灵敏度的绝对值相对于密度设计变量灵敏度的绝对值要小很多,为了避免陷入局部解,提高收敛效率,本文采用式(19)定义的变步长策略,即使设计变量移动步长随着迭代步数的增加逐渐减小

其中,movex和moveθ 分别为单元密度变量和角度设计变量的移动步长.它们的初始值分别设置为move-x=0.5,move θ=0.1.

当优化过程满足式(20) 定义的收敛条件时,停止迭代

式中,k为当前迭代步数,ε 为容许的收敛误差,并取ε=1.0×10-3.

3 数值算例与讨论

本节通过一个静力学优化算例和两个算例动力学优化算例来展示本文方法的有效性.若无特殊说明,平板厚度设置为1 mm,质量密度设置为1.6 g/cm3;本构材料阵中4 个独立的物理参数取值为Ex=385 GPa,Ey=118 GPa,Vxy=0.3 和Gxy=84 GPa.假设初始结构为完全实体材料,即所有密度设计变量xe的初始值均为1.本文所有的角度变量均以弧度制表示.本文所有算例中MMA 算法参数取值列于表1 中.

表1 MMA 子程序中参数取值Table 1 Parameters needed for MMA subroutine

3.1 算例一

算例一考虑短悬臂梁的刚度最大化优化问题.几何尺寸为80 mm×40 mm 的设计域和边界条件如图2 所示,梁左端完全固支,在右端中心点处受大小为30 N 的横向拉力.设计域由80×40 个平面应力单元来离散.材料准许的体积分数设定为0.3.

图2 悬臂梁的设计域和边界条件Fig.2 Design domain and boundary condition of cantilever beam

为了验证本文方法对静力学优化问题的有效性,考虑以0,π/6,π/2,及[-2π,2π] 区间内产生的一个随机数为初始纤维方向.图3 显示了4 种不同初始纤维方向所对应的优化结果.从结果可以看出,尽管初始纤维角度不同,但优化后的纤维基本沿轴向拉伸方向分布,增强了结构的刚度,符合工程实际情况,这验证了本文方法对静力学优化问题的有效性.

图3 不同初始纤维方向对应的优化结果Fig.3 Optimization results for different initial fiber orientations

3.2 算例二

算例二考虑悬臂梁动力学优化问题.设计域与算例一相同,梁左端完全固支,在右边中心点处附加一个大小为28.8 g 的集中质量点.设置纤维初始方向为水平方向,即设置所有角度设计变量θe的初始值为0.材料准许的体积分数设定为0.4.

图4 显示了连续纤维增强悬臂梁在无阻尼自由振动下结构拓扑构型的迭代过程.优化算法在迭代87 步之后收敛,此时结构对应的最小特征值为188.76 Hz.图5 显示了悬臂梁的纤维分布情况.可以看出,本文方法可以获得清晰的拓扑构型,优化后纤维方向基本沿悬臂梁横向中轴线呈对称分布.

图4 悬臂梁拓扑构型的迭代过程Fig.4 Iterative process of topological configuration of cantilever beam

图5 悬臂梁纤维方向的优化结果Fig.5 Optimization result of fiber orientation for cantilever beam

为了进一步验证所提方法的有效性,在图5 显示的最终纤维角度的优化结果基础上,让所有单元对应的角度变量逆时针偏转角度β.在几个典型的偏转角度π/6,π/4,π/3,5π/12,π/2 下,结构对应的最小特征值λ1分别为128.96,96.06,75.41,64.63,61.32.可以看出,纤维方向偏转一定角度后结构对应的最小特征值均比通过本文优化方法获得的结果(188.76 Hz)要小,这说明本文方法能够有效找到更优的纤维方向分布,提高结构的频率.

图6 显示了各向同性材料悬臂梁在无阻尼自由振动下结构拓扑构型的优化结果,其中材料本构阵中使用的独立物理参数为:杨氏模量E=385 GPa,泊松比V=0.3.可以看出,本文获得的拓扑构型与文献[15] 中的动力学优化结果类似,优化算法在迭代61 步之后收敛,此时对应的最小特征值为188.87 Hz.对于当前的悬臂梁优化算例,各向同性材料拓扑优化结果与纤维增强的正交各向异性材料结构拓扑构型结果相似.

图6 各向同性材料悬臂梁的拓扑构型Fig.6 Topological configuration of cantilever beam with isotropic material

3.3 算例三

算例三考虑固支细长梁结构的动力学优化问题.几何尺寸为160 mm×40 mm 的设计域和边界条件如图7 示,梁两端全固定,平板中心处附加一个质量为51.2 g 的质量点.设计域由160×40 个平面应力单元来离散.设置纤维角度设计变量的初始值为θe=0.材料准许的体积分数为0.4.

图7 固支梁的设计域和边界条件Fig.7 Design domain and boundary condition of clamped beam

图8 显示了连续纤维增强固支细长梁在无阻尼自由振动下结构拓扑构型的迭代过程.图9 显示了固支细长梁纤维方向的优化结果.从优化结果可知,本文方法可以获得清晰的拓扑构型,纤维方向主要沿着杆件的方向分布.图10 显示了前两阶特征值的迭代历史.从结果可以看出,优化算法在迭代77 步之后到达收敛,此时结构对应的最小特征值为373.86 Hz,第2 阶特征值为22 869.56 Hz,前两阶特征值在优化的过程中始终分离的,第2 阶特征值始终大于第1 阶特征值,未产生重特征根问题,这验证了本文关注的特征值为单根的假定.

同样,为了进一步验证本文方法的有效性,在图9 显示的最终纤维角度的优化结果基础上,让所有单元对应的角度变量逆时针偏转角度β.在几个典型的偏转角度π/6,π/4,π/3,5π/12,π/2 下，结构对应的最小特征值λ1的值分别为262.18,198.32,157.31,135.58,128.84.由表3 可知,结果与悬臂梁相优化问题类似,纤维方向偏转一定角度后所得结构的最小特征值均变小,这进一步说明本文方法的有效性.

图8 固支梁拓扑构型的迭代过程Fig.8 Iterative process of topological configuration of clamped beam

图9 固支梁纤维方向的优化结果Fig.9 Optimization result of fiber orientation for clamped beam

图10 固支梁前2 阶特征值的迭代历史Fig.10 Iteration history of the first two eigenvalue of clamped beam

4 结论

针对连续纤维增强复合材料结构无阻尼自由振动下的基频最大化问题,本文提出了一种有效实现拓扑构型和纤维方向同步优化的拓扑优化方法.建立了以准许的材料用量体积分数为约束、以结构的一阶特征值为目标函数的动力学拓扑优化模型,详细推导了目标函数关于密度设计变量和角度设计变量的解析灵敏度列式.为了避免陷入局部最优解,采用了变移动步长的策略.最后,通过一个静力学优化算例和两个动力学优化算例验证了本文方法的有效性.结果表明,本文方法在实现结构拓扑构型与纤维角度的一体化优化的同时,可以有效改善结构的性能.