等強度梁理論下的多截面獨立并行優化研究

徐峰祥,王得偉,鄒 震,梁 銳

(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) 本文中所提模型設計的等強度梁,可應用于一些由復雜數學公式定義的任意型截面梁優化問題,為梁的截面形狀的優化研究提供了新的解決方案。