拓?fù)浣Y(jié)構(gòu)中邊緣棋盤格式的光滑處理方法

劉云龍，買買提明?艾尼

（新疆大學(xué)機(jī)械工程學(xué)院，新疆 烏魯木齊 830047）

1 引言

根據(jù)不同的優(yōu)化手段，可將結(jié)構(gòu)優(yōu)化設(shè)計(jì)分為尺寸優(yōu)化、形狀優(yōu)化和拓?fù)鋬?yōu)化[1]。而拓?fù)鋬?yōu)化則是結(jié)構(gòu)優(yōu)化設(shè)計(jì)的最高層次，其主要思想是將尋求結(jié)構(gòu)的最優(yōu)拓?fù)鋯栴}轉(zhuǎn)化為在給定的設(shè)計(jì)區(qū)域內(nèi)尋求最優(yōu)材料的分布問題，以達(dá)到結(jié)構(gòu)的最大剛度或最大基頻[2]。對于拓?fù)鋬?yōu)化設(shè)計(jì)，眾多的學(xué)者都做了許多工作。其中，文獻(xiàn)[3]在1904年提出的桁架理論可以認(rèn)為是拓?fù)鋬?yōu)化設(shè)計(jì)的起源；隨后，文獻(xiàn)[4、5]引入了材料“微結(jié)構(gòu)”的概念，指明了連續(xù)體結(jié)構(gòu)拓?fù)鋬?yōu)化的發(fā)展方向；1988年，文獻(xiàn)[6]在微結(jié)構(gòu)概念的基礎(chǔ)上，提出了均勻化拓?fù)鋬?yōu)化方法，這標(biāo)志著連續(xù)體拓?fù)鋬?yōu)化的發(fā)展進(jìn)入黃金時期，而當(dāng)前ANSYS等商業(yè)工程軟件都有基于均勻化方法的拓?fù)鋬?yōu)化模塊。

ANSYS軟件具有優(yōu)良的建模功能和強(qiáng)大的分析能力，是進(jìn)行拓?fù)鋬?yōu)化設(shè)計(jì)的有力工具之一。但是基于均勻化方法的拓?fù)鋬?yōu)化結(jié)果邊緣往往具有棋盤格式的不規(guī)則和非線性特點(diǎn)，這種棋盤格式就給拓?fù)鋬?yōu)化結(jié)果的CAD模型重構(gòu)和CAE分析帶來一定的困難[7、8]。棋盤格式的出現(xiàn)與所采用的材質(zhì)設(shè)計(jì)變量無關(guān)，而是與有限元單元有密切關(guān)系，也就是說棋盤格式是連續(xù)體結(jié)構(gòu)拓?fù)鋬?yōu)化中依賴單元類型和數(shù)量的一種固有的現(xiàn)象。因此，在拓?fù)鋬?yōu)化中，不論是采用均勻化設(shè)計(jì)方法還是采用密度法，只要用有限元法進(jìn)行拓?fù)浣Y(jié)構(gòu)優(yōu)化均會不可避免的出現(xiàn)棋盤格式。

棋盤格式問題的解決一直是很多學(xué)者關(guān)心的問題之一，并提出了一些切實(shí)有效的避免棋盤格式的求解策略。比如，采用后處理技術(shù)將棋盤格式過濾的方法，采用較為穩(wěn)定的有限元模型或改變優(yōu)化目標(biāo)函數(shù)的泛函從而使優(yōu)化過程趨于穩(wěn)定等等。這些方法都能不同程度的緩解或減弱棋盤格式現(xiàn)象，但是不能完全消除，這依然為拓?fù)鋬?yōu)化結(jié)構(gòu)的實(shí)際工程應(yīng)用帶來了一定的困難。

這里則應(yīng)用最小二乘法的線性擬合原理對邊緣棋盤格式進(jìn)行線性光滑處理，并根據(jù)邊緣光滑處理的結(jié)果重構(gòu)新的光滑拓?fù)鋬?yōu)化結(jié)構(gòu)，從而完全消除了棋盤格式。通過對新的光滑拓?fù)鋬?yōu)化結(jié)構(gòu)模型再進(jìn)行網(wǎng)格劃分和有限元分析驗(yàn)證了此方法的有效性。這種使用最小二乘法對棋盤格式進(jìn)行線性光滑處理方法，可完全消除棋盤格式并建立基本保持拓?fù)浣Y(jié)構(gòu)的新的有效拓?fù)浣Y(jié)構(gòu)模型。本方法理論簡單、容易編程實(shí)現(xiàn)、建模速度快、易于操作等優(yōu)點(diǎn)，為拓?fù)鋬?yōu)化設(shè)計(jì)提供了一種新的手段。

2 多邊緣棋盤格式的最小二乘法描述

2.1 基于均勻化方法的拓?fù)鋬?yōu)化方法

均勻化方法是一種經(jīng)典的拓?fù)鋬?yōu)化方法，有著嚴(yán)密的數(shù)學(xué)和力學(xué)理論基礎(chǔ)，由于在基于均勻化方法的拓?fù)鋬?yōu)化方面已有很多的文獻(xiàn)和詳細(xì)的報(bào)道，因此在此不再重復(fù)[6、9]。這里則是直接利用ANSYS中的基于均勻化方法的拓?fù)鋬?yōu)化方法進(jìn)行了建模和優(yōu)化分析，得到了存在棋盤格式的小梁結(jié)構(gòu)邊緣。并針對這種具有棋盤格式的邊緣，使用最小二乘法進(jìn)行了光滑處理。

2.2 基于最小二乘法的多邊緣棋盤格式光滑處理算法

最小二乘法是一個重要的數(shù)學(xué)與統(tǒng)計(jì)模型，在回歸分析、數(shù)據(jù)擬合、圖像處理等方面有著廣泛的應(yīng)用[10]。最小二乘法通過最小化誤差的平方和找到一組數(shù)據(jù)的最佳函數(shù)匹配，是用最簡單的方法求得一些絕對不可知的真值，而令誤差平方之和最小[11]。利用最小二乘法，可以對平面內(nèi)一系列非線性的散點(diǎn)進(jìn)行線性或非線性處理。本研究考慮到拓?fù)鋬?yōu)化結(jié)構(gòu)的實(shí)際工程應(yīng)用價值，盡量把小梁結(jié)構(gòu)部分簡化為線性的梁，因此本研究用線性的最小二乘法來近似小梁結(jié)構(gòu)邊緣棋盤格式的非線性散點(diǎn)（節(jié)點(diǎn)），并用一條直線方程y=ax+b來描述這一系列非線性的散點(diǎn)所組成的非線性曲線。對于復(fù)雜拓?fù)浣Y(jié)構(gòu)存在多個獨(dú)立的小梁結(jié)構(gòu)邊緣棋盤格式域，并組成多個獨(dú)立的散點(diǎn)域，其基于最小二乘法的多邊緣棋盤格式光滑處理算法介紹如下。

設(shè)在平面直角坐標(biāo)系內(nèi)有l(wèi)個散點(diǎn)域，每個散點(diǎn)域有nl個散點(diǎn)組成的散點(diǎn)集，i=1，2，3…n。在第l個散點(diǎn)域內(nèi)用一個一次函數(shù)方程式來表達(dá)這n l個散點(diǎn)，則有：

式中：y l—xl的預(yù)估值。令S為這n個散點(diǎn)的實(shí)際值與預(yù)估值的差的平方和，即：

此時，令S(al，bl)=0，并在S(al，b l)中分別對al和b l求偏導(dǎo)，則有：

聯(lián)立式（4）、式（5），且把a(bǔ)和b看作是未知數(shù)，則式（4）、式（5）的聯(lián)立式是一個兩元一次方程組，求解這個兩元一次方程組，得：

令al為斜率，bl為截距，此時可以得到一條平面直角坐標(biāo)系第l個散點(diǎn)域內(nèi)的直線。通過求解所有散點(diǎn)域的光滑直線并連接相鄰直線后即可得到完全消除棋盤格式的新的拓?fù)鋬?yōu)化結(jié)構(gòu)。

3 多邊緣棋盤格式結(jié)構(gòu)的光滑處理

3.1 拓?fù)鋬?yōu)化及其邊界條件

在平面上有一寬為3000mm，高為1000mm的板，材料選用45號鋼，即彈性模量為2.1×1011，泊松比為0.269，屈服強(qiáng)度355MPa[12]，如圖1所示。平面板的左側(cè)受固定約束；上端距左側(cè)1000mm處受一方向垂直向下，大小為1000N的力；下端最右側(cè)處受一方向垂直向上，大小為1000N的力。網(wǎng)格尺寸為50mm，采用四邊形網(wǎng)格主導(dǎo)的四邊形∕三邊形網(wǎng)格。拓?fù)鋬?yōu)化的約束條件設(shè)置為減少50%的面積。在上述條件下先在ANSYS Workbench中對模型進(jìn)行求解。

圖1 平面板的邊界條件Fig.1 Boundary Conditions of Planar Plates

為了便于拓?fù)鋬?yōu)化結(jié)果的提取，建立的項(xiàng)目圖表，如圖2所示。以便將ANSYS Workbench中的數(shù)據(jù)轉(zhuǎn)換進(jìn)入ANSYS APDL。之后右擊C2項(xiàng)目圖表，選擇Edit in Mechanical APDL以進(jìn)入ANSYS APDL環(huán)境。

圖2 轉(zhuǎn)換數(shù)據(jù)所用的項(xiàng)目圖表Fig.2 Item Charts Used to Convert Data

進(jìn)入ANSYS APDL后，利用APDL命令流對模型再進(jìn)行一次拓?fù)鋬?yōu)化計(jì)算，所需要的APDL命令流為：

即可得到約束條件為面積減少50%的拓?fù)鋬?yōu)化結(jié)果，如圖3所示。其中，紅色部分為保留材料的部分，藍(lán)色部分為去除材料的部分。可以看出，拓?fù)鋬?yōu)化結(jié)果具有明顯的多個帶棋盤格式的梁結(jié)構(gòu)。

圖3 平面板面積減少50%的拓?fù)鋬?yōu)化結(jié)果Fig.3 Topological Optimization Results for 50%Reduction of Planar Plate Area

3.2 拓?fù)鋬?yōu)化結(jié)構(gòu)的光滑處理

由第二節(jié)描述的最小二乘法可知，只要得到某條非線性曲線附近若干散點(diǎn)坐標(biāo)，就可利用最小二乘法求得一條直線方程，并用這條直線方程來近似代替非線性曲線。這里就用以上方法對復(fù)雜拓?fù)鋬?yōu)化多小梁邊緣棋盤格式結(jié)構(gòu)進(jìn)行光滑處理。其具體的光滑處理過程如下。

3.2.1 拓?fù)鋬?yōu)化結(jié)果的提取

通過ANSYS分析得到的拓?fù)鋬?yōu)化結(jié)果，如圖4所示。圖中可以看出拓?fù)鋬?yōu)化結(jié)果包含多個小梁結(jié)構(gòu)及多邊緣棋盤格式的結(jié)構(gòu)。通過詳細(xì)分析多邊緣棋盤格式將各個邊緣設(shè)置了總共22個編號，每一個編號代表一個邊緣棋盤格式散點(diǎn)域，整個拓?fù)鋬?yōu)化結(jié)構(gòu)內(nèi)總共有22個散點(diǎn)域。為了將這些散點(diǎn)域組成的棋盤格式設(shè)為一條線進(jìn)行光滑處理，首先對這些散點(diǎn)域內(nèi)的棋盤格式點(diǎn)（節(jié)點(diǎn)）進(jìn)行坐標(biāo)的提取。以圖4中的10號棋盤格式散點(diǎn)域?yàn)槔肁NSYS APDL中的Picked Entities+命令對10號棋盤格式散點(diǎn)域內(nèi)的節(jié)點(diǎn)進(jìn)行拾取，如圖5所示。

圖4 復(fù)雜拓?fù)鋬?yōu)化多小梁結(jié)構(gòu)的邊緣棋盤格式編號及散點(diǎn)域Fig.4 Edge Chessboard Format Numbering and Spot Domain of Multitrabecular Structures with Complex Topology Optimization

圖5 10號棋盤格式散點(diǎn)域內(nèi)節(jié)點(diǎn)的拾取Fig.5 Picking up Nodes in Dispersion Domain of Chessboard Format No.10

在選擇全部節(jié)點(diǎn)之后，單擊拾取器中的OK即可得到這些節(jié)點(diǎn)的坐標(biāo)。同理，利用相同的方法可以提取其余21條棋盤格式散點(diǎn)域內(nèi)的節(jié)點(diǎn)，在此不一一贅述。

在節(jié)點(diǎn)拾取的過程中通過觀察坐標(biāo)可以發(fā)現(xiàn)，1、2、5、7號線已經(jīng)是平面直角坐標(biāo)系內(nèi)的直線了，它們分別為：x=-1500、y=1000、x=1500、x=-1449.4。

3.2.2 利用最小二乘法進(jìn)行光滑處理

由圖5（a）拾取器中的信息可知，10號棋盤格式散點(diǎn)域內(nèi)共有19個節(jié)點(diǎn)，它們的坐標(biāo)分別為：（-160.57，499.72）（-110.25，499.68）（-65.786，449.70）（-116.00，449.88）（-256.12，549.52）（-206.07，549.57）（-155.88，549.57）（-19.418，400.02）（-70.070，399.95）（-302.41，599.44）（-252.38，599.44）（-22.903，350.30）（-350.02，649.36）（-300.03，649.36）（-400.01，699.36）（-350.02，699.36）（-450.01，749.36）（-400.01，749.36）（-450.01，799.36）。

將以上19個點(diǎn)帶入式（6）和式（7）得：

a10≈-0.892024562

b10≈364.9228648

以上系數(shù)a10和b10可得10號棋盤格式散點(diǎn)域光滑處理后的一條直線方程y10=-0.892024562x+364.9228648。可用這條直線方程來近似代替圖4中的10號棋盤格式散點(diǎn)域內(nèi)散點(diǎn)所組成的曲線。同理，除1、2、5、7號之外，利用相同的方法均可得到代替剩余17個棋盤格式散點(diǎn)域內(nèi)曲線的近似直線方程。

3.2.3 利用WPS表格簡化計(jì)算過程

眾所周知，WPS的表格軟件擁有強(qiáng)大的函數(shù)計(jì)算功能。在求解圖4中每條線的過程中，可以利用WPS表格來計(jì)算每條線的斜率al和截距bl。

仍以圖4中10號棋盤格式散點(diǎn)域?yàn)槔鐖D6所示。將在ANSYS APDL中提取的19個點(diǎn)的X坐標(biāo)粘貼進(jìn)A列，Y坐標(biāo)粘貼進(jìn)B列；C2單元格輸入=A2*A2，并下拉單元格至C20；D2單元格輸入=A2*B2，并下拉單元格至D20；G1單元格輸入=SUM（C2：C20）；G2單元格輸入=SUM（A2：A20）；G3單元格輸入=SUM（D2：D20）；G4單元格輸入=SUM（B2：B20）；G5單元格輸入=（19*G3-G2*G4）∕（19*G1-G2*G2）；G6單元格輸入=（G1*G4-G2*G3）∕（19*G1-G2*G2）。此時，WPS表格軟件即可自動計(jì)算出10號棋盤格式散點(diǎn)域內(nèi)所對應(yīng)的直線斜率al和截距bl，它們分別為G5和G6中的數(shù)值，如圖6所示。

圖6 利用WPS表格軟件計(jì)算直線斜率和截距Fig.6 Using WPS to Calculate Slope and Interception of Lines

使用相同方法，加之之前已得到的1、2、5、7號直線，可得22條近似直線。其中，4號直線的斜率a4為1.39167×10-33，截距b4為-2.23408×10-30，由于數(shù)值非常小，因此可以認(rèn)為4號直線為：y=0。

綜上，可得圖4中1、2、4、5、7號線的近似直線方程，如表1所示。而對于圖4中剩余17條線的近似直線方程，可以求得其截距al和斜率bl，如表2所示。

表1 1，2，4，5，7號線的近似直線的方程Tab.1 Approximate Linear Equations for Line 1，2，4，5 and 7

表2 剩余17條近似直線的斜率a l和截距b lTab.2 Slope A and Intercept B of the Remaining 17 Approximate lines

4 光滑結(jié)果的重構(gòu)與新拓?fù)浣Y(jié)構(gòu)分析

根據(jù)拓?fù)鋬?yōu)化得到的非線性非光滑結(jié)果，將3.2.3中所得的（1～5）號直線方程、（6～10）號直線方程、（11～13）號直線方程、（14～16）號直線方程、（17～19）號直線方程、（20～22）號直線方程分別聯(lián)立，求得它們的交點(diǎn)，即可得重構(gòu)線性光滑模型所需的22個關(guān)鍵點(diǎn)，如表3所示。

表3 重構(gòu)線性光滑模型所需的22個關(guān)鍵點(diǎn)坐標(biāo)Tab.3 22 Key Points Coordinates for Reconstructing Linear Smooth Model

根據(jù)表2中的22個關(guān)鍵點(diǎn)坐標(biāo)，即可在ANSYS Workbench中對拓?fù)鋬?yōu)化且進(jìn)行了光滑處理的平面板模型進(jìn)行重構(gòu)。對重構(gòu)的拓?fù)鋬?yōu)化模型進(jìn)行網(wǎng)格劃分，如圖7所示。網(wǎng)格尺寸為10mm，采用四邊形網(wǎng)格主導(dǎo)的四邊形∕三邊形網(wǎng)格。

圖7 重構(gòu)的模型Fig.7 Reconstructed Model

將3.1的邊界條件再次施加到重構(gòu)的光滑模型上，模型的左側(cè)受固定約束；上端距左側(cè)1000mm處受一方向垂直向下，大小為1000N的力；下端最右側(cè)處受一方向垂直向上，大小為1000N的力。并設(shè)置模型的彈性模量為2.1×1011，泊松比為0.269。通過有限元分析，分別對其進(jìn)行位移、應(yīng)力和應(yīng)變分析，可得的計(jì)算結(jié)果，如圖8～圖10所示。

圖8 重構(gòu)模型的位移分析Fig.8 Displacement Analysis of Reconstructed Model

圖10 重構(gòu)模型的應(yīng)變分析Fig.10 Strain Analysis of Reconstructed Model

由圖8可知，重構(gòu)模型的最大變形量為0.69344mm，最大變形分布在重構(gòu)拓?fù)鋬?yōu)化模型的最右端。但是，從圖9和圖10中可以看出，應(yīng)變和應(yīng)力都分布的比較均勻，基本消除了應(yīng)力集中。最大應(yīng)變?yōu)?.00020234，非常小，最大應(yīng)變發(fā)生在重構(gòu)模型中17、18號直線的交點(diǎn)處，如圖9所示。同樣，最大應(yīng)力為38.078MPa，最大應(yīng)力也發(fā)生在重構(gòu)模型中17、18號直線的交點(diǎn)處，這是由于交接點(diǎn)出現(xiàn)銳角而未進(jìn)行圓滑處理的緣故，雖然這樣最大應(yīng)力仍比45號鋼的屈服強(qiáng)度355MPa小得多，如圖10所示。綜上，對利用最小二乘法光滑處理后所得到的重構(gòu)拓?fù)鋬?yōu)化模型進(jìn)行有限元分析結(jié)果表明，應(yīng)力應(yīng)變分布均勻，表明最小二乘法光滑處理后重構(gòu)所得到的拓?fù)鋬?yōu)化模型基本保持了原拓?fù)浣Y(jié)構(gòu)，同時驗(yàn)證了本方法的有效性。

圖9 重構(gòu)模型的應(yīng)力分析Fig.9 Stress Analysis of Reconstructed Model

5 結(jié)論

（1）建立了多小梁邊緣棋盤格式復(fù)雜拓?fù)浣Y(jié)構(gòu)的最小二乘法光滑處理算法，并利用WPS對拓?fù)鋬?yōu)化的結(jié)果進(jìn)行計(jì)算，實(shí)現(xiàn)了光滑處理，為實(shí)際工程應(yīng)用提供了一種簡便且易于操作的多小梁邊緣棋盤格式的復(fù)雜拓?fù)浣Y(jié)構(gòu)光滑處理方法；（2）通過AN?SYS Workbench和ANSYS APDL軟件實(shí)現(xiàn)了拓?fù)鋬?yōu)化結(jié)果的提取，并將光滑處理后的拓?fù)鋬?yōu)化結(jié)果在ANSYS Workbench中進(jìn)行了重構(gòu)，實(shí)現(xiàn)了完全消除棋盤格式結(jié)構(gòu)的光滑拓?fù)鋬?yōu)化結(jié)構(gòu)，并對其進(jìn)行了相應(yīng)的有限元分析，驗(yàn)證了本方法的有效性；（3）通過新光滑拓?fù)鋬?yōu)化結(jié)構(gòu)有限元分析可知，在給定的邊界條件下和面積減少50%的拓?fù)鋬?yōu)化約束條件下，利用最小二乘法進(jìn)行邊界光滑處理得到的重構(gòu)拓?fù)鋬?yōu)化模型的應(yīng)力應(yīng)變分布均勻，消除了應(yīng)力集中，光滑處理后的重構(gòu)拓?fù)鋬?yōu)化模型基本保持了原拓?fù)浣Y(jié)構(gòu)形狀。