等强度梁理论下的多截面独立并行优化研究

徐峰祥,王得伟,邹 震,梁 锐

(1.武汉理工大学 现代汽车零部件技术湖北省重点实验室, 武汉 430070;2.武汉理工大学 汽车零部件技术湖北省协同创新中心, 武汉 430070;3.桂林航天工业学院 汽车工程学院, 广西 桂林 541004)

0 引言

在轻量化设计过程中,很多梁结构为实现减重而设计成变截面梁结构,如丁华等[1]采用连续变截面管对仪表板横梁进行轻量化,凌智勇等[2]对船用起重机臂架减重设计也采用变截面梁等,而等强度梁便是变截面梁结构的一种,其可以使所有截面上的最大应力同时达到材料的许用应力,从而提高材料的利用率[3]。此外,对于等强度梁这类异形结构的加工,一般通过拉拔或挤压成型的方法实现,可以相对容易地获得任何形状。等强度梁已在许多论文中进行了研究[4-7],但这些论文通常只考虑梁的一些最简单形状,这些形状只满足单个集中力垂向弯曲的等应力条件,无法满足其他复杂力系。

为获得等强度梁结构,主要采用2种设计方法,即理论计算方法和有限元仿真。第一种是推演等强度梁的数学模型,如孟增等[3]以矩形截面梁为例,推演了超静定梁的轴向尺寸等强度设计公式,张定等[8]给出了矩形截面和T形截面梁高的等强度设计公式等,这种方法设计出来的梁截面形状单一,考虑因素较少,难以应用于复杂的工程问题中。第二种是采用有限元软件对结构进行优化,程海鹰等[9]借助anasy软件,对仿等强度梁的关键结构尺寸进行优化设计,Mohammedali等[10]同样采用软件优化了轴向变深度的工字形梁结构,这种方法无非是软件操作,优化速度低下,优化变量数也有限。除了单独采用以上这2种设计方法外,还有采用理论计算在前、仿真优化在后的方法,如鲍明森等[11]以连续变厚度轮廓曲线对转盘进行等强度设计,而后采用有限元软件对局部应力集中进行手动调参优化,使结构整体表现为等强度,夏云[12]先建立了等截面悬臂梁数学模型,再采用软件对梁模型优化成等强度悬臂梁。

与以上传统方法获取等强度梁结构有所不同,本文中将理论计算方法融入进有限元仿真模型中,借助Abaqus软件二次开发平台设计出等强度梁结构,这种设计方法在其他结构设计上也有应用,如De等[13-14]基于Python对Abaqus进行二次开发,实现对结构的优化设计。其中,理论计算方法并非是通过公式直接推导出满应力状态的截面尺寸公式,而是构建截面优化模型,通过编程优化出面积最小且达到满应力状态的截面,有关截面优化模型的研究,主要有Ragnedda等[15]提出了一种弯曲梁截面均匀强度设计问题的优化设计方法,Croccolo等[16]通过移除离中性轴最远区域材料的优化设计方法,Liu等[17]提出了梁截面拓扑优化公式,Qin等[18]提出了一种缩短概念设计周期的两级截面优化方法,Li等[19]通过软件优化模块对梁的承载力和抗振能力进行了优化设计,Zuo等[20]提出了一种快速的拼焊板截面形状设计与优化方法。

本文中以设计任意截面形状的等强度梁为目标,在理论上建立了截面优化模型和截面应力模型,通过Abaqus-Python软件二次开发平台将理论数学模型与有限元仿真模型进行结合,构建出一套应用于工程问题优化研究的多截面独立并行优化模型。

1 一般梁模型

目前针对梁截面形状的选择多采用简单且为实心的截面,在进行等强度设计时考虑的受力也都较为简单,为研究更为一般的梁模型,本文将考虑梁结构所受的6个力系,且采用任意截面形状进行设计。将任意截面形状和受任意力系的工程梁简化为以下一般梁模型(图1),包括待设计梁截面位置、截面形状以及端面的受力状态3部分。

图1 一般梁模型

用数学模型表征以上所建三维梁模型:待设计截面位置为X=[x1,x2,…,xq],其中,q为待设计的截面数量。某一截面形心的主失与主矩记为FI(x)=[Fx,Fy,Fz],MI(x)=[Mx,My,Mz]。截面形状由数学方程描述G(y,z)=0。

当采用独立并行优化方法对梁的多个截面进行优化时,并行可提高计算速度,独立可以通过切断截面间变量的耦合提高截面参数的优化数量。切断变量间的关联产生的一些实际因素产生的误差,随后通过有限元仿真进行修正。因此,确保设计的梁具备等强度属性的关键是要在理论上建立高效的截面优化模型和应力计算模型。

2 截面优化模型

基于以往梁结构的等强度设计方法,是将梁每个截面设计成满应力状态[21-22],由于截面的面积与对应的最大应力间关系为负相关,因此,将截面面积最小为优化目标,截面最大应力为约束,建立如下优化模型。

1) 目标函数

以截面面积最小为目标。

minS

(1)

2) 设计变量

截面形状(轮廓)由一系列设计参数决定,定义截面的设计变量如下:

U={u1,u2,…,up}

(2)

式中,p为截面形状设计参数总数。

3) 约束条件

截面形状约束:

s.t.ue_min≤ue≤ue_maxe=1,2,…,p

(3)

式中,ue_min、ue_max为设计空间的边界。

强度约束:

φσmax≤[σ]

(4)

式中:[σ]为截面许用应力值,σmax为截面的最大等效应力值,φ为应力修正系数,默认为1。

截面最小厚度约束:

hmin<[h]

(5)

式中:[h]为截面最小厚度约束值。采用如下离散数学模型计算截面的最小厚度,以单连通截面为例,内外轮廓曲线离散成点的坐标如式(6)(7)所示。

(6)

(7)

式中:Yf和Zf为内轮廓各离散点的坐标,Yc和Zc为外轮廓各离散点的坐标,t为内/外轮廓离散的数量。

各点间组成的距离方阵表达如下:

h=[(Yc×1-1T×Yf)2+

(8)

式中:向量1为[1,1,…,1]1×n,n为轮廓上点的总数。则截面最薄区域厚度表达如下:

hmin=min(h)

(9)

3 截面应力模型

目前,对于任意形状截面的等效应力计算尚且缺少统一的力学模型,只有截面各应力分量的力学模型,这些模型以连续数学模型和离散数学模型为主。为了统一,本文中采用离散数学模型,将截面离散化成网格(如图2),计算各离网格节点的各应力分量及对应的等效应力值,进而求得截面最大等效应力。

图2 梁截面网格划分示意图

截面节点的应力分量记为q=[σx,τy,τz],则各点的等效应力表达如下:

(10)

式中:σ为等效应力值,σx为各点的x方向正应力,τy为各点的y方向切应力,τz为各点的z方向切应力。

其中,各点的x方向正应力主要由轴力和弯矩产生,截面所受轴力和弯矩如图1所示,根据材料力学[23]可得这些应力的数学表达:

σx=σx1+σx2+σx3

(11)

(12)

式中:σx1、σx2和σx3分别为x方向轴力、y和z方向弯矩产生的正应力,S为截面面积,Iz和Iy为z方向和y方向的惯性矩,(y,z)截面各离散点的坐标。

各点y方向和z方向的切应力主要由扭矩和剪力产生,数学表达如下:

τy=τy1+τy2

(13)

τz=τz1+τz2

(14)

式中,由扭矩产生的切应力表达如下:

(15)

(16)

式中:Jk为扭转系数,φ为截面翘曲函数,参考文献[24-25]。

由剪力产生的切应力表达如下:

(17)

(18)

式中:τiy2为第i点的y向切应力,t(yi)为第i点的z向截面厚度,S*(yi)表示截面在直线y=yi一侧面积对坐标轴的静矩,τiz2、t(yi)和S*(zi)的含义同理。

所有应力分量的求解,唯独由剪力产生的切应力无法采用矩阵计算,需要通过单点循环计算(如图3),且每次循环都得积分1次,这种求解方法非常耗时,因为矩阵计算在程序设计上通过并行计算的方法提高其计算速度[26-27],因此,可以通过建立一套剪切应力的矩阵计算模型,提高剪力切应力的计算速度。

图3 单点循环计算求解截面剪切应力

截面轮廓控制方程较为复杂,求解静矩采用积分法进行计算很不方便,且各点沿坐标轴方向的厚度计算繁琐。而将数据矩阵化则可显著提高计算截面节点切应力的效率。因此建立了如下基于截面各点静矩以及厚度的数据矩阵化的截面微方格模型。其处理流程如图4。

图4 截面微方格模型处理流程

通过截面离散化成方格,对方格进行数据化处理,计算出各行各列对应的静矩和厚度矩阵,得出截面各点对y轴和z轴的静矩和厚度矩阵,进而计算出各点的切应力值,具体模型如下。

1) 截面的离散化处理

先初始化方格矩阵,如图4左上角的图,记方格离散化矩阵为G。初始值G为全1的矩阵,矩阵中,1代表截面包络了该方格,0则相反;维度n和m的计算公式如下:

n=D/c,m=H/c

(19)

式中:D、H分别为外轮廓外接矩形宽度和高度,c为尺寸数据精度,也为单个方格尺寸,单位为mm,如尺寸宽度为10.01 mm,则c为0.01 mm;n、m表示方格集合的大小。

最后需统计各方格的位置信息,由方格中心点的坐标表示。表达式如下:

(20)

式中:Ys和Zs分别为离散点中心的横坐标和纵坐标;Ln和Lm是两组值都为1的向量,向量长度分别为n和m,(yo,zo)为截面形心。

2) 方格的数据化处理

首先,计算单个方格静矩。每个方格对y轴和z轴的静矩表达式如下:

Cy=Zs·G×c2

(21)

Cz=Ys·G×c2

(22)

然后,计算单排方格静矩。横/纵排方格对y/z轴的静矩表达式如下:

(23)

接着,计算面积静矩。横/纵排方格中心远离形心两侧所围成面积对y/z轴形成的静矩表达式分别如下:

(24)

(25)

式中:T1、T3分别为维度为 (m-z0/c)和y0/c的单位上三角矩阵,T2、T4分别为维度z0/c和 (n-y0/c)的单位下三角矩阵。

最后,计算沿y/z方向截面厚度。纵/横排网格中心横坐标沿y/z轴方向的截面厚度表达式如下:

(26)

3) 各点的切应力数据

通过以上数据,可计算每排方格中心处的切应力。每排网格中心处沿y/z轴方向的切应力表达式如下:

(27)

以上计算的是离散化方格中心的切应力值,截面上任意一点的切应力约等于该点附近方格中心的切应力值,对于截面任意一点(yi,zi),其切应力可通过对向量Uy和Uz索引获得,表达式如下:

τiy2=Uy(yi/c),τiz2=Uz(zi/c)

(28)

即剪力作用下的梁截面所有网格节点切应力值为

τy2=Uy(y/c),τz2=Uz(z/c)

(29)

4 多截面独立并行优化算法

多截面独立并行优化算法框架分为2层,外层为等强度梁设计模块,内层为截面优化模块(如图5)。内层采用理论计算公式进行理论截面优化,以实现计算时间的减少,有利于提高优化速度,但存在计算误差,外层进行结构建模和分析,采用有限元软件求解,可模拟结构的真实应力状态[28-29],如应力集中等,计算精度高,可对内层的力学模型进行修正。

图5 多截面独立并行优化算法内外层模块

内层的截面优化模块具体流程如图6所示,可计算任意形状截面的几何特征和在受力下的截面最大应力。

图6 内层截面优化模块

外层的等强度梁设计模块具体流程如图7,将三维梁模型多个待设计截面独立出来,调入内层进行优化,每个截面优化时占用1个CPU核资源,将理论模型优化出来的最优截面反馈到外层进行有限元分析,由于采用应力的理论计算存在一定的误差,需通过有限元应力解计算应力修正系数进行修正,经历多轮迭代收敛后,每个截面应力都趋于并小于约束值,进而表现出等强度特点。其中,应力修正系数等于仿真应力值与理论应力值的比值。

图7 外层等强度梁设计模块

5 工程应用

为验证本文中提出的多截面独立并行优化算法在汽车轻量化领域应用的可行性,以汽车从动桥为例,从动桥作为汽车的支撑梁,其梁上各部位所受尺寸约束和受力都较为复杂[30-31],在车桥受极限载荷工况下,采用本文中所建模型对车桥关键截面进行优化,以分析本文所建模型优化效果。其中,桥壳的极限工况具体为:受到冲击载荷的同时施加最大制动力,此时从动桥桥壳受到的转矩、垂向和横向弯矩最大。

5.1 待优化桥壳模型

图8表示从动桥壳的1/4模型,因为从动桥壳是左右对称结构,且受力也是左右对称,所以只需设计一半,从右往左依次标记待优化的截面,一共11个,这些截面主要是转角处的过渡截面、约束和载荷加载处截面以及细化结构所需的截面。

图8 从动桥壳模型

用数学模型表征以上所建桥壳模型:待设计截面位置为X=[180,200,240,280,320,360,400,440,500,550,600],单位mm。截面形状的数学方程描述为,形状为圆环型。以额定载荷为5 t的从动桥为例,对左端面和板簧连接处[550,600]mm施加约束,作用于桥壳轴承处的主失和制动器底座的主矩分别为FI(100)=[0,49 000,61 250]N,MI(320)=[20 580 000,0,0]N·mm。

5.2 待优化桥壳模型

桥壳优化问题描述为

① 优化目标:各桥壳截面面积最小;

② 优化变量:各个截面的内外半径R,r;

③ 约束条件:由于桥壳需要与轮毂、制动器和板簧底座连接,故需要对桥壳各待优化截面的外观尺寸约束,编号从小到大依次为[35,40,50,60,60,60,70,80,80,80,80]mm,最薄区域厚度不小于4 mm,最大理论等效应力为350 MPa。

优化前的桥壳截面外圆半径尺寸为约束条件的边界值,截面厚度都为9 mm,桥壳总质量为68.51 kg,桥壳的应力云图如图9,最大应力为376.45 MPa。

图9 极限工况下的等厚桥壳应力云图

5.3 待优化桥壳模型

基于所建立的截面优化模型和应力计算方法,采用Python语言对Abaqus软件进行二次开发,根据多截面独立并行优化算法编写的程序共有2层循环,在程序外层循环的应力修正迭代下,桥壳总质量和各截面最大应力随外层迭代次数的变化如图10、图11所示,最终桥壳总质量为61.38 kg。对于程序内层循环的优化迭代,对桥壳11个截面采用圆环形截面进行独立并行优化,优化模型采用差分进化算法求解,该算法搜索范围广,鲁棒性好,收敛速度快,适合于以实数编码的优化模型[32-33],以外层最后一次计算的各截面的优化迭代曲线为例,如图12。

图10 桥壳总质量随外层迭代图

图11 各截面最大应力随外层迭代图

图12 各截面面积随内层优化迭代图

外层迭代第一次和最后一次的桥壳应力云图如图13、图14所示,迭代第一次的最大应力为433.42 MPa,大于理论约束应力值,通过引入应力修正系数进行再次迭代优化后,各截面的应力都逐渐趋于理论约束应力值的下边缘,应力分布也更加均匀,最终桥壳表现出等强度或满应力的特点,如图14所示。

图13 外层迭代第一次的桥壳应力云图

图14 外层迭代最后一次的桥壳应力云图

圆环型桥壳的11个截面优化结果如图15,曲线表示不同截面位置处的设计变量最优值。

图15 圆环型桥壳各截面的最优尺寸值

当采用圆环型截面对桥壳进行优化时,11个待优化的截面共22个设计变量,采用多截面独立并行优化算法,每个核只需要优化2个变量,优于优化软件对22个变量同时进行优化的计算成本。

5.4 复杂梁截面的优化研究

为进一步研究本文中所提方法对复杂多变量梁截面的优化效果,设计了由双超椭圆方程描述的任意型截面,并采用与优化圆环型截面一样的方法对其进行优化,分析其优化效果。

截面形状的双超椭圆方程描述为

(30)

式中:A为外轮廓半宽,B为半高,a为内轮廓半宽,b为半高,这4个参数控制截面的尺寸;N1、N2、n1、n2为形状参数,这4个参数控制截面的形状。

不同参数N下的超椭圆曲线簇如图16,涵盖了从矩形到菱形、椭圆、矩形等一系列截面形状。

图16 超椭圆曲线簇示意图

相较于圆环型截面,该任意型截面的优化难度是优化变量数多,应用于前述桥壳优化中,总设计变量高达88个。应用该截面对桥壳进行优化设计,分析本文所建多截面独立并行优化算法的鲁棒性和稳定性。

桥壳优化问题的描述与圆环型桥壳优化一致,增加一条对形状参数的约束,为保证桥壳能与汽车其他零部件的配合,前6个截面的形状参数约束为1,后5个截面的形状参数约束范围为1～2。在Abaqus二次开发平台的程序优化迭代下,桥壳总质量和各截面最大应力随外层迭代次数的变化如图17、图18所示,迭代收敛的桥壳总质量为63.04 kg,优化出的桥壳结构及对应应力云图如图19所示。

图17 桥壳总质量随外层迭代图

图19 外层迭代最后一次的桥壳应力云图

各截面的尺寸变量和形状变量随截面位置的变化曲线如图20、图21。

图20 各截面纵向尺寸优化结果

图21 各截面形状调和参数优化结果

优化出的各截面形状如表1所示。可以看出,双超椭圆方程优化出的截面具有角厚边薄的特点,即结构材料远离惯性轴,经计算,截面11两个方向的单位面积抗弯截面系数分别为44.23和45.78,而优化出的圆环型截面仅有36.49,采用双超椭圆方程截面的桥壳与圆环型截面桥壳的质量相差2.63%,但抗弯曲性能提升了25%左右。

表1 不同截面位置处的最优截面形状

5.5 结果分析

将优化前的等厚度圆环型桥壳与优化后的2种截面形状桥壳进行对比,得出各类型桥壳的降重与最大等效应力降低,如表2。

表2 2种桥壳的质量和应力优化值 %

显然,优化后的2种桥壳质量和最大等效应力都有所下降,质量和应力都降低10%左右,实现了桥壳轻量化的目的。

综上,本文中所提的基于等强度理论的多截面独立并行优化算法,在工程应用中,对于桥壳这种非单一长方体组成的非均匀设计空间,能充分考虑结构的应力集中问题,进而设计出等强度的桥壳结构,此外,对于截面形状复杂的工程梁问题,该算法的优化效果同样显著。

6 结论

以设计等强度梁为目标,开发了适用于任意力系和截面形状下的多截面独立并行优化算法。随后通过建立微分格模型,统一了截面应力矩阵计算方法,并将其应用于截面优化模型。基于二次开发对理论力学模型进行修正从而得到等强度梁,并采用多截面独立并行优化算法对汽车从动桥进行了优化研究,得出如下结论。

1) 多截面独立并行优化算法解决了非均匀设计空间的应力集中、多截面变量难以收敛的问题,通过二次开发方法设计出比较理想的等强度工程梁。

2) 基于多截面独立并行优化算法对桥壳进行了轻量化设计,结果显示:桥壳整体呈现为等强度结构。圆环型桥壳在降重10.41%的基础上最大等效应力值降低10.64%。

3) 本文中所提模型设计的等强度梁,可应用于一些由复杂数学公式定义的任意型截面梁优化问题,为梁的截面形状的优化研究提供了新的解决方案。