基于等幾何裁剪分析的拓?fù)渑c形狀集成優(yōu)化

第一作者傅曉錦男，博士，教授，1964年11月生

基于等幾何裁剪分析的拓?fù)渑c形狀集成優(yōu)化

傅曉錦1，龍凱2，周利明1，闕春蘭1，葉航1

(1. 上海電機(jī)學(xué)院機(jī)械學(xué)院，上海200245； 2. 華北電力大學(xué)新能源電力系統(tǒng)國家重點(diǎn)實(shí)驗(yàn)室，北京102206)

摘要：提出了將設(shè)計(jì)和分析、拓?fù)渑c形狀優(yōu)化集成的思想，探索了基于等幾何裁剪分析的拓?fù)渑c形狀集成優(yōu)化設(shè)計(jì)算法，該方法統(tǒng)一了結(jié)構(gòu)優(yōu)化的計(jì)算機(jī)輔助設(shè)計(jì)、計(jì)算機(jī)輔助工程分析和優(yōu)化設(shè)計(jì)的模型，基于B樣條的等幾何裁剪分析既能準(zhǔn)確表達(dá)幾何形狀，又可以用裁剪面分析方便處理任意復(fù)雜拓?fù)鋬?yōu)化問題，由裁剪選擇標(biāo)準(zhǔn)確定合理的拓?fù)浣Y(jié)構(gòu)變動(dòng)方向，結(jié)構(gòu)變動(dòng)時(shí)無需重新劃分網(wǎng)格，設(shè)計(jì)結(jié)果突破初始設(shè)計(jì)空間的限制，還可方便優(yōu)化形狀。建立了等幾何裁剪靈敏度分析的計(jì)算方法，給出了等幾何裁剪分析拓?fù)渑c形狀集成優(yōu)化算法，通過典型實(shí)例表明所用方法的正確性和有效性。

關(guān)鍵詞：形狀優(yōu)化； 靈敏度分析； 拓?fù)鋬?yōu)化； 等幾何分析； 裁剪B樣條； 有限元分析

基金項(xiàng)目：國家自然科學(xué)基金(11202078);上海市教委科研創(chuàng)新重點(diǎn)項(xiàng)目(13ZZ145);上海市自然科學(xué)基金(11ZR1413800)

收稿日期：2014-02-24修改稿收到日期：2014-04-21

中圖分類號(hào):O241；TH122

文獻(xiàn)標(biāo)志碼：A

DOI:10.13465/j.cnki.jvs.2015.07.026

Abstract:An idea of integration of computer aided design, computer aided engineering analysis, topology and shape optimization design was introduced and a kind of optimization strategy for integation of topology and shape optimization design of continuum structure based on isogeometric trimmed surface analysis was explored. The proposed method offers a unified model of computer aided design, computer aided engineering and optimization design in structure optimization. The isogeometric trimmed surface analysis based on B splines has the capability of both expressing the geometry shape accurately and solving arbitrarily complex topology optimization problem. The trimmed criteria selected determine the reasonable direction of topology changes. It does not need to remesh during structural optimization process. The design results can be beyond the limitations of initial design space. It is also convenient in shape optimization. In the integration method, the calculation approach for sensitivity analysis of isogeometric trimmed surface, and the algorithm for integation of topology and shape optimization design based on isogeometric trimmed analysis were developed. A numerical example illustrates the correctness and effectiveness of the method.

Integration of topology and shape optimization design of continuum structure based on isogeometric trimmed surface analysis

FUXiao-jin1,LONGKa2,ZHOULi-ming1,QUEChun-lan1,YEHang1(1. School of Mechanical Engineering, Shanghai Dianji University, Shanghai 200245,China 2.State Key Laboratory for Alternate Electrical Power System with Renewable Energy Sources, North China Electric Power University, Beijing 102206,China)

Key words:shape optimization; sensitivity analysis; topology optimization; isogeometric analysis; trimmed B-spline; finite element analysis

計(jì)算機(jī)輔助結(jié)構(gòu)優(yōu)化設(shè)計(jì)是一項(xiàng)多學(xué)科相交叉的新技術(shù)。它的研究領(lǐng)域涉及到力學(xué)、機(jī)械、結(jié)構(gòu)工程、數(shù)學(xué)和信息技術(shù)，具有巨大的發(fā)展?jié)摿1-3]。我國程耿東院士和郭旭在布局優(yōu)化研究中提出求解奇異最優(yōu)解的松弛算法，是具有重要意義的貢獻(xiàn)[4]。Michael等[5]提出一種的水平集方法(Level set method)，優(yōu)化問題中結(jié)構(gòu)的邊界用嵌入到高維尺度函數(shù)中的水平集模型來表示，模型在描述復(fù)雜結(jié)構(gòu)的拓?fù)浼斑吔缱兓矫婢哂休^好的靈活性。目前，多數(shù)研究者先計(jì)算機(jī)輔助設(shè)計(jì)(Computer Aided Design，CAD)，之后，進(jìn)行計(jì)算機(jī)輔助工程(Computer Aided Engineering，CAE)分析，幾何和分析模型會(huì)有出入，CAD/CAE模型難以融合[6-7]，優(yōu)化方法的最優(yōu)解與初始產(chǎn)生空洞的數(shù)量與位置有很大關(guān)系[8-9]。其中存在二個(gè)明顯的問題：第一，設(shè)計(jì)空間與初始網(wǎng)格有關(guān)，設(shè)計(jì)空間為固定網(wǎng)格，設(shè)計(jì)結(jié)果與設(shè)計(jì)空間相關(guān)，若限制設(shè)計(jì)空間，有可能得不到解。另一方面，若網(wǎng)格劃分過細(xì)，會(huì)使計(jì)算量大大增加。第二，將計(jì)算結(jié)果返回到CAD系統(tǒng)的后處理工作量也大，這是由于常規(guī)優(yōu)化方法設(shè)計(jì)和分析的模型不是CAD系統(tǒng)中樣條數(shù)據(jù)，需要大量的后處理工作。第一個(gè)缺點(diǎn)可采用無網(wǎng)格解決。克服第二個(gè)缺點(diǎn)需要統(tǒng)一設(shè)計(jì)和分析的模型[10-11]。Hughes等[12]提出了一種等幾何分析的新型有限元方法，該方法采用Non-Uniform Rational B-Splines (NURBS)作為有限元離散的插值函數(shù)，將控制頂點(diǎn)作為網(wǎng)格節(jié)點(diǎn)，使CAD和CAE之間通信既方便容易，又有效[13-17]。在Hughes等[18-25]幾何分析的啟發(fā)下，很多學(xué)者在結(jié)構(gòu)優(yōu)化領(lǐng)域結(jié)合等幾何分析進(jìn)行了廣泛的研究，優(yōu)化技術(shù)迅速發(fā)展，創(chuàng)造了不少成果，具有較好的參考價(jià)值。但以上文獻(xiàn)的拓?fù)鋬?yōu)化方法多采用細(xì)胞或單元，然而，細(xì)胞表示的拓?fù)涫遣灰?guī)則或模糊邊界輪廓，所以，試圖研究一種能確定的光滑的材料邊界的方法，最近，水平集拓?fù)鋬?yōu)化已被廣泛的研究[26]，水平集拓?fù)鋬?yōu)化中，由雅可比方程確定進(jìn)化水平集函數(shù)表示材料邊界，設(shè)計(jì)靈敏度分析計(jì)算形狀速度，一般，水平集拓?fù)鋬?yōu)化有一個(gè)嚴(yán)重的不足，在優(yōu)化過程中，新的內(nèi)部空洞不能自動(dòng)產(chǎn)生，所以，初始的水平集會(huì)包含一些用戶定義的內(nèi)部空洞，這樣，該拓?fù)鋬?yōu)化方法的問題是會(huì)依靠初始空洞的數(shù)量和位置，在優(yōu)化過程中，合理的策略應(yīng)該是基于拓?fù)潇`敏度和應(yīng)變能量密度引入新的空洞。本文在前期集成優(yōu)化研究的基礎(chǔ)上[27-28]，探索等幾何裁剪靈敏度分析方法，利用裁剪面分析任意復(fù)雜拓?fù)浣Y(jié)構(gòu)，無需將其分成多塊，從而降低由于分塊引起的巨大計(jì)算量，有限元模型與幾何實(shí)體重合，統(tǒng)一了CAD、CAE和優(yōu)化設(shè)計(jì)的模型。

1幾何造型

非均勻有理B樣條(Non-Uniform Rational B-Splines，NURBS)為曲線曲面幾何造型主要方法[29-31]。

1.1節(jié)點(diǎn)矢量與B樣條

節(jié)點(diǎn)是參數(shù)軸上的分割點(diǎn)，一個(gè)非遞減的參數(shù)坐標(biāo)的集合定義為節(jié)點(diǎn)矢量，記為

Ξ=[ξ0,ξ1,…ξn+p]

式中，ξi∈R是第i個(gè)節(jié)點(diǎn)，p是B樣條基函數(shù)的次數(shù)，n是基函數(shù)的數(shù)量。

定義B樣條基函數(shù)

(1)

(2)

用n個(gè)p次B樣條基函數(shù)Ni,p(ξ)(i=1,2,…,n)和n個(gè)控制頂點(diǎn)Pi(i=1,2,…,n)確定p次B樣條曲線：

(3)

1.2NURBS描述

(4)

NURBS曲線可表示為

(5)

由節(jié)點(diǎn)矢量ξ=[ξ0,ξ1,…ξn+p]和節(jié)點(diǎn)矢量η=[η0,η1,…ηm+q]，n×m個(gè)控制頂點(diǎn)Pi,j(i=1,2,…,n;j=1,2,…,m)，構(gòu)成的B樣條曲面：

(6)

由節(jié)點(diǎn)矢量ξ=[ξ0,ξ1,…ξn+p]、節(jié)點(diǎn)矢量η=[η0,η1,…ηm+q]和ζ=[ζ0,ζ1,…ζm+q]，n×m×l個(gè)控制頂點(diǎn)Pi,j,k(i=1,2,…,n;j=1,2,…,m;k=1,2,…,l)構(gòu)成的B樣條體：

S(ξ,η,ζ)=

(7)

1.3基于NURBS的等幾何分析

建立準(zhǔn)確的幾何模型是等幾何分析的關(guān)鍵，CAD(計(jì)算機(jī)輔助設(shè)計(jì))和CAE(計(jì)算機(jī)輔助工程分析)是同一個(gè)模型，即等幾何，用有限元方法進(jìn)行分析。用B樣條基函數(shù)代替多項(xiàng)式基函數(shù)，控制頂點(diǎn)視同有限元節(jié)點(diǎn)，可得：

(8)

式中，u為如位移、應(yīng)力和應(yīng)變等控制變量的物理場(chǎng)，d可以為控制頂點(diǎn)的位移、應(yīng)力和應(yīng)變等。

2裁剪面分析

常規(guī)的等幾何分析時(shí)，分析模型是分解成幾個(gè)B樣條小塊，所以，若拓?fù)浣Y(jié)構(gòu)復(fù)雜，分析較為困難，若借助計(jì)算機(jī)圖形學(xué)的裁剪操作，即使復(fù)雜模型也容易實(shí)現(xiàn)。裁剪單元按照單元的頂點(diǎn)數(shù)量分為三類，如圖1所示。復(fù)雜的裁剪單元與這三類裁剪單元不匹配，那么繼續(xù)四分細(xì)化該單元，直到滿足這三類的其中一類為止。裁剪單元繼續(xù)分為三角形和曲邊三角形，曲形三角形有一邊為B樣條裁剪曲線。

圖1　三類裁剪單元 Fig.1 Three kind of trimmed elements

因參數(shù)空間比物理空間簡單，故單元分解在未裁剪的參數(shù)空間操作。單元?jiǎng)偠染仃嚍?/p>

(9)

B=MΓG

(10)

其中，Jξ為由參數(shù)空間向物理空間轉(zhuǎn)換的變換矩陣，表示為

(11)

J=JTkJTn=

(12)

如圖2所示，在有限元中，線性變換Tkt能用線性三角形形狀函數(shù)表示。

圖2　三角形單元變換 Fig.2 Transformations in triangular elements

為方便積分，采用線性變換Tkt

Tkt:{u,v}→{ξ,η}

(13)

如圖3所示，正積分三角形單元的二頂點(diǎn)(ξ1,η1)和(ξ2,η2)與裁剪矩形單元的頂點(diǎn)一致，另一點(diǎn)(ξ3,η3)位于裁剪曲線上，可求裁剪曲線交點(diǎn)及其參數(shù)。假定，點(diǎn)(ξ3,η3)中，ξ3為已知，η3為未知。

(1)由已知參數(shù)用下式可求λ0

(2)由下式求得未知參數(shù)

圖3　在三角形單元中求裁剪曲線交點(diǎn)及其參數(shù) Fig.3 Finding intersection and its corresponding parameter of trimming curve in triangular elements

正三角形單元的雅可比矩陣為

(14)

圖4為由高斯積分單元轉(zhuǎn)化為物理空間曲邊單元，物理空間單元用ΩPH表示；相應(yīng)參數(shù)空間單元用ΩPM表示。

圖4　曲邊三角形積分單元變換 Fig.4 Transformations in curved triangular integration elements

(1)由已知參數(shù)用下式可求λ1

(2)由下式求得未知參數(shù)

假定，點(diǎn)(ξ3,η3)，ξ3為已知，η3為未知。

(1)由已知參數(shù)用下式可求λ2

(2)由下式求得未知參數(shù)

圖5　在曲邊三角形單元中求裁剪曲線交點(diǎn)及其參數(shù) Fig.5 Finding intersection and its corresponding parameter of trimming curve in triangular elements

在圖4中，線性變換Ti，Tj和Tkt為

(16)

(17)

在αβ坐標(biāo)系φ(λ)=[φα(λ),φβ(λ)]T裁剪曲線能用逆矩陣R-1變換

φ(λ)=R-1·C(λ)=A-1(C(λ)-Z)

(18)

曲邊積分三角形單元的雅可比矩陣為

J=JTiJTjJTktJTn=

(19)

式中，

(20)

圖6　積分單元分解和單元的積分點(diǎn) Fig.6 Integration elements decomposition and integration points in the elements

圖6顯示積分單元分解和積分點(diǎn)。單元按圖5規(guī)則分解，在曲邊三角形積分單元中的積分點(diǎn)均布在邊界，曲邊三角形積分單元的剛度矩陣由式(19)和式(9)計(jì)算可得。

3裁剪靈敏度分析

本文取邊界控制頂點(diǎn)的坐標(biāo)為設(shè)計(jì)變量，建立優(yōu)化模型

MinΨ(X,u(X))

s.t.hj(X,u(X))=0j=1,…,p1

(21)

3.1靈敏度計(jì)算

目標(biāo)函數(shù)方程

(22)

式中，F(xiàn)為線性方程的力矢量，u為位移矢量。

目標(biāo)函數(shù)對(duì)設(shè)計(jì)變量的靈敏度為

(23)

對(duì)式(9)的單元?jiǎng)偠染仃嚽髮?dǎo)

(24)

對(duì)B的導(dǎo)數(shù)

(25)

對(duì)于矩形單元：

對(duì)于正三角形單元：

對(duì)于曲邊三角形單元：

(26)

用式(25)和式(26)代入式(24)求得單元?jiǎng)偠让舳龋纱丝傻每倓偠让舳?/p>

(27)

式中，ne為單元數(shù)量。

形狀優(yōu)化時(shí)希望結(jié)構(gòu)的最大應(yīng)力最小，取最大Von-Mises應(yīng)力為研究對(duì)象

(28)

Von-Mises應(yīng)力對(duì)控制頂點(diǎn)坐標(biāo)的導(dǎo)數(shù)為

(29)

計(jì)算極值點(diǎn)，求得最大應(yīng)力

(30)

式中，d為在u全局內(nèi)最大應(yīng)力點(diǎn)的控制頂點(diǎn)位移，A為由全局位移矢量轉(zhuǎn)為局部位移矢量的變換矩陣，A在該處值為1，其它地方為0。局部處靈敏度為

(31)

式中，v為下列伴隨方程的解。

(32)

4優(yōu)化方法

本部分介紹優(yōu)化方法和策略，確定設(shè)計(jì)變量，定義選擇標(biāo)準(zhǔn)，提出相應(yīng)算法。

4.1設(shè)計(jì)變量

在二維優(yōu)化問題，首先，未裁剪表面的邊界表面控制頂點(diǎn)坐標(biāo)為設(shè)計(jì)變量，因?yàn)閮?nèi)部表面控制頂點(diǎn)不影響邊界形狀，所以，內(nèi)部表面控制頂點(diǎn)不作為設(shè)計(jì)變量。一旦空洞生成，裁剪曲線的控制頂點(diǎn)坐標(biāo)也作為設(shè)計(jì)變量了。控制頂點(diǎn)P的法向量n定義

(33)

式中，(ξP,ηP)是參數(shù)空間最接近裁剪點(diǎn)的控制頂點(diǎn)的坐標(biāo)，t為切向量，k為單位向量。目標(biāo)函數(shù)設(shè)計(jì)靈敏度與設(shè)計(jì)表面控制頂點(diǎn)法向移動(dòng)有關(guān)，目標(biāo)函數(shù)設(shè)計(jì)靈敏度為

(34)

4.2產(chǎn)生空洞的標(biāo)準(zhǔn)

對(duì)于式(21)的優(yōu)化模型，定義選擇標(biāo)準(zhǔn)為目標(biāo)函數(shù)對(duì)約束的靈敏度

(35)

體積作為約束條件時(shí)，SC通常為正。當(dāng)空洞產(chǎn)生時(shí)，體積就會(huì)變小。如果SC為負(fù)，沒有改變約束，不可能再使目標(biāo)函數(shù)性能提高。

SC有二種，一種是設(shè)計(jì)模型內(nèi)產(chǎn)生空洞的選擇標(biāo)準(zhǔn)SC洞，另一種是設(shè)計(jì)模型內(nèi)外邊界形狀改變的選擇標(biāo)準(zhǔn)SC形。目標(biāo)函數(shù)的形狀靈敏度和設(shè)計(jì)邊界的控制頂點(diǎn)約束都要計(jì)算相應(yīng)的SC形。比較二個(gè)選擇標(biāo)準(zhǔn)SC形和SC洞能確定設(shè)計(jì)模型內(nèi)部空洞是否產(chǎn)生。

若設(shè)計(jì)模型中，約束不變，設(shè)計(jì)控制頂點(diǎn)的變動(dòng)而引起邊界攝動(dòng)，內(nèi)部空洞產(chǎn)生，則

(36)

(37)

由式(36)和式(37)可得

(38)

圖7　在邊界變動(dòng)時(shí)表面設(shè)計(jì)控制頂點(diǎn)和選擇標(biāo)準(zhǔn) Fig.7 Surface design control points and their selestion criteria on boundary movement

4.3算法

(1)定義模型

(2)裁剪面分析

(3)靈敏度分析

(5)更新設(shè)計(jì)變量

(6)內(nèi)部空洞是否合并?若合并返回(2)

(7)裁剪面分析

(8)收斂否?若不收斂，返回(3)

(9)收斂得優(yōu)化結(jié)果

5應(yīng)用實(shí)例

實(shí)例1

如圖8所示為多載荷梁的結(jié)構(gòu)優(yōu)化計(jì)算模型，結(jié)構(gòu)尺寸為200 mm×100 mm，彈性模量E=2.01×105MPa，泊松比ν=0.3，定義整個(gè)區(qū)域?yàn)樵O(shè)計(jì)域。在下端分別施加100 N、600 N、100 N三個(gè)集中力載荷，并在左下角約束了6個(gè)自由度，右下角約束了y方向的移動(dòng)自由度。

圖8　優(yōu)化模型 Fig.8 Optimization model

●表示設(shè)計(jì)面的控制頂點(diǎn)，■表示初始面的控制頂點(diǎn) 圖9　物理空間初始面的控制頂點(diǎn)和設(shè)計(jì)面的控制頂點(diǎn) Fig.9 Initial control pointsand design surface control points in physical area

如圖9所示優(yōu)化模型為二次B樣條面，分為20×10個(gè)單元。在物理空間，以體積作為約束，以小方塊表示初始表面的控制頂點(diǎn)，以小圓點(diǎn)表示設(shè)計(jì)面邊界(頂端邊界)的控制頂點(diǎn)，首先以設(shè)計(jì)面邊界的控制頂點(diǎn)為設(shè)計(jì)變量，為了避免網(wǎng)格扭曲變形，二末端控制頂點(diǎn)允許沿垂直方向移動(dòng)，其余控制頂點(diǎn)允許沿法向移動(dòng)，其余邊界無設(shè)計(jì)面控制頂點(diǎn)，但可由裁剪曲線按需要改變邊界。按照上述優(yōu)化方法進(jìn)行優(yōu)化，進(jìn)行靈敏度分析，設(shè)計(jì)面邊界的控制頂點(diǎn)(即設(shè)計(jì)變量)就會(huì)按計(jì)算的靈敏度方向移動(dòng)，其余表面控制頂點(diǎn)會(huì)相應(yīng)變動(dòng)，盡管在小的區(qū)域控制頂點(diǎn)會(huì)較密集，但自適應(yīng)使其不會(huì)相交。通過改變表面控制頂點(diǎn)坐標(biāo)進(jìn)行優(yōu)化，直到由拓?fù)鋵?dǎo)數(shù)確定內(nèi)部產(chǎn)生空洞位置為止，得如圖10所示物理空間中表面的優(yōu)化結(jié)果。

圖10　物理空間中表面的優(yōu)化結(jié)果 Fig.10 Optimal result in physical area

一旦空洞生成，裁剪曲線控制頂點(diǎn)的坐標(biāo)也成為設(shè)計(jì)變量。下一步重點(diǎn)分析裁剪曲線控制頂點(diǎn)的靈敏度。

裁剪曲線方程為

(39)

式中，λ為曲線參數(shù)，ng為裁剪曲線控制頂點(diǎn)數(shù)量，裁剪曲線第k個(gè)控制頂點(diǎn)坐標(biāo)為Pk=[pk,qk]。

圖11　交點(diǎn)靈敏度計(jì)算時(shí)裁剪曲線 控制頂點(diǎn)變動(dòng)及其相應(yīng)曲線參數(shù) Fig.11 Movement of trimming curve control point for sensitivity compute of intersections and their corresponding curve parameters

圖11說明由于裁剪曲線控制頂點(diǎn)攝動(dòng)，二個(gè)交點(diǎn)及其相應(yīng)參數(shù)變動(dòng)情況。攝動(dòng)前，二交點(diǎn)分別為(ξ2b,η2b)和(ξ3b,η3b)，其相應(yīng)曲線參數(shù)分別為λ2b和λ3b。裁剪曲線控制頂點(diǎn)垂直攝動(dòng)后，二交點(diǎn)分別變?yōu)?ξ2a,η2a)和(ξ3a,η3a)，其相應(yīng)曲線參數(shù)也分別變?yōu)棣?a和λ3a。交點(diǎn)及其參數(shù)對(duì)曲線控制頂點(diǎn)攝動(dòng)的靈敏度分別為

(40)

應(yīng)變位移矩陣B對(duì)曲線控制頂點(diǎn)的靈敏度用Γ和G的靈敏度表示為

(41)

其中，Γ的靈敏度表示為

(42)

式中，

(43)

G對(duì)曲線控制頂點(diǎn)的靈敏度表示為

(44)

第k個(gè)基函數(shù)對(duì)曲線控制頂點(diǎn)的一階導(dǎo)數(shù)用參數(shù)(ξ,η)靈敏度和基函數(shù)的二階導(dǎo)數(shù)表示

(45)

(46)

由式(14)可得雅可比靈敏度

(47)

由式(17)計(jì)算曲線三角形積分單元的(ξ,η)靈敏度

(48)

由式(40)計(jì)算交點(diǎn)的靈敏度

(49)

由式(16)計(jì)算式(48)中參數(shù)α和β的靈敏度

(50)

由式(18)可得

(51)

其中

(52)

由式(39)可得裁剪曲線靈敏度

(53)

由式(15)可得曲線參數(shù)λ靈敏度

(54)

(55)

(56)

(57)

(58)

(59)

對(duì)裁剪曲線求偏導(dǎo)

(60)

它對(duì)曲線控制頂點(diǎn)的靈敏度為

(61)

由式(56)到式(61)可計(jì)算曲邊三角形積分單元的式(55)雅可比行列式靈敏度。

圖12 在參數(shù)空間中的模型Fig.12Modelintheparametricarea圖13 在物理空間中裁剪曲線控制頂點(diǎn)Fig.13Controlpointsinthephysicalarea圖14 在物理空間中選定裁剪曲線的詳細(xì)信息Fig.14Detailedinformationofselectcurveinthephysicalarea

圖15 物理空間中未裁剪表面的優(yōu)化結(jié)果Fig.15Optimaluntrimmedresultinthephysicalarea圖16 物理空間表面控制頂點(diǎn)和逐步細(xì)化積分單元Fig.16Controlpointsofdecomposedintegrationelementsinthephysicalarea圖17 裁剪面分析應(yīng)力云圖(Von-Mises應(yīng)力場(chǎng))Fig.17Stresscontoursoftrimmingcurve

圖18 裁剪優(yōu)化結(jié)果Fig.18Optimaltrimmedresult圖19 形狀優(yōu)化結(jié)果Fig.19Shapeoptimalresult圖20 形狀優(yōu)化應(yīng)力云圖(Von-Mises應(yīng)力場(chǎng))Fig.20Stresscontoursofshapeoptimization

圖16所示為在物理空間中的模型，顯示適應(yīng)表面控制頂點(diǎn)的敏度，圖17所示為裁剪面分析應(yīng)力云圖(Von-Mises應(yīng)力場(chǎng))，圖18所示為裁剪優(yōu)化結(jié)果。

然后，再進(jìn)行形狀優(yōu)化靈敏度分析，形狀優(yōu)化時(shí)希望結(jié)構(gòu)的最大應(yīng)力最小，取最大Von-Mises應(yīng)力為研究對(duì)象,求Von-Mises應(yīng)力對(duì)控制頂點(diǎn)坐標(biāo)的導(dǎo)數(shù)，計(jì)算極值點(diǎn)，得優(yōu)化形狀，圖19 所示為形狀優(yōu)化結(jié)果，圖20形狀優(yōu)化應(yīng)力云圖。

實(shí)例2

如圖21所示為單載荷梁的結(jié)構(gòu)優(yōu)化計(jì)算模型，結(jié)構(gòu)尺寸為200 mm×100 mm，彈性模量E=2.01×105MPa，泊松比ν=0.3，定義整個(gè)區(qū)域?yàn)樵O(shè)計(jì)域。在右邊緣下端施加一個(gè)100 N的載荷，并約束了左邊緣6個(gè)自由度。

如圖22所示優(yōu)化模型為二次B樣條面，分為20×10個(gè)單元。在物理空間，以體積作為約束，以小方塊表示初始表面的控制頂點(diǎn)，以小圓點(diǎn)表示設(shè)計(jì)面邊界(頂端邊界)的控制頂點(diǎn)，首先以設(shè)計(jì)面邊界的控制頂點(diǎn)為設(shè)計(jì)變量，為了避免網(wǎng)格扭曲變形，二末端控制頂點(diǎn)允許沿垂直方向移動(dòng)，其余控制頂點(diǎn)允許沿法向移動(dòng)，其余邊界無設(shè)計(jì)面控制頂點(diǎn)，但可由裁剪曲線按需要改變邊界。按照上述優(yōu)化方法進(jìn)行優(yōu)化，進(jìn)行靈敏度分析，設(shè)計(jì)面邊界的控制頂點(diǎn)(即設(shè)計(jì)變量)就會(huì)按計(jì)算的靈敏度方向移動(dòng)，其余表面控制頂點(diǎn)會(huì)相應(yīng)變動(dòng)，盡管在小的區(qū)域控制頂點(diǎn)會(huì)較密集，但自適應(yīng)使其不會(huì)相交。通過改變表面控制頂點(diǎn)坐標(biāo)進(jìn)行優(yōu)化，直到由拓?fù)鋵?dǎo)數(shù)確定內(nèi)部產(chǎn)生空洞位置為止，空洞生成后，裁剪曲線控制頂點(diǎn)的坐標(biāo)也成為設(shè)計(jì)變量，更新設(shè)計(jì)變量，外部邊界和內(nèi)部邊界同時(shí)優(yōu)化，物理空間和參數(shù)空間的拓?fù)浣Y(jié)構(gòu)同時(shí)相應(yīng)改變。

圖21 優(yōu)化模型Fig.21Optimizationmodel圖22 物理空間初始面的控制頂點(diǎn)和設(shè)計(jì)面的控制頂點(diǎn)Fig.22Initialcontrolpointsanddesignsurfacecontrolpointsinphysicalarea圖23 在參數(shù)空間中的模型Fig.23Modelintheparametricarea

圖24　在物理空間中裁剪曲線控制頂點(diǎn) Fig.24 Control points in the physical area

圖25 在物理空間中選定裁剪曲線的詳細(xì)信息Fig.25Detailedinformationofselectcurveinthephysicalarea圖26 物理空間表面控制頂點(diǎn)和逐步細(xì)化積分單元Fig.26Controlpointsofdecomposedintegrationelementsinthephysicalarea圖27 裁剪面分析應(yīng)力云圖(Von-Mises應(yīng)力場(chǎng))Fig.27Stresscontoursoftrimmingcurve

圖26所示為在物理空間中的模型，顯示適應(yīng)表面控制頂點(diǎn)的敏度，圖27所示為裁剪面分析應(yīng)力云圖(Von-Mises 應(yīng)力場(chǎng))，圖28所示為裁剪優(yōu)化結(jié)果。

然后，再進(jìn)行形狀優(yōu)化靈敏度分析，形狀優(yōu)化時(shí)希望結(jié)構(gòu)的最大應(yīng)力最小，取最大Von-Mises應(yīng)力為研究對(duì)象,求Von-Mises應(yīng)力對(duì)控制頂點(diǎn)坐標(biāo)的導(dǎo)數(shù)，計(jì)算極值點(diǎn)，得優(yōu)化形狀，圖29 所示為形狀優(yōu)化結(jié)果，圖30形狀優(yōu)化應(yīng)力云圖。

圖28 裁剪優(yōu)化結(jié)果Fig.28Optimaltrimmedresult圖29 形狀優(yōu)化結(jié)果Fig.29Shapeoptimalresult圖30 形狀優(yōu)化應(yīng)力云圖(Von-Mises應(yīng)力場(chǎng))Fig.30Stresscontoursofshapeoptimization

6結(jié)論

(1)本文探討了將設(shè)計(jì)和分析、拓?fù)渑c形狀優(yōu)化雙集成的基于等幾何裁剪分析的連續(xù)體結(jié)構(gòu)拓?fù)渑c形狀集成優(yōu)化設(shè)計(jì)方法，該方法將在CAD系統(tǒng)表達(dá)幾何模型控制頂點(diǎn)和基函數(shù)的樣條信息用于CAE，CAD、CAE和優(yōu)化設(shè)計(jì)是同一個(gè)模型，實(shí)現(xiàn)其無縫集成。

(2)用B樣條和裁剪曲線表示設(shè)計(jì)模型的內(nèi)外邊界，用裁剪面分析能方便處理復(fù)雜樣條面，改變拓?fù)浣Y(jié)構(gòu)，設(shè)計(jì)結(jié)果不限定在初始的設(shè)計(jì)空間內(nèi)。將未裁剪表面控制頂點(diǎn)坐標(biāo)為設(shè)計(jì)變量，一旦空洞生成，裁剪曲線控制頂點(diǎn)坐標(biāo)也作為設(shè)計(jì)變量，無需額外的設(shè)計(jì)變量參數(shù)。為了避免未裁剪面的幾個(gè)網(wǎng)格發(fā)生畸變，將除了角點(diǎn)外的設(shè)計(jì)面控制頂點(diǎn)沿法線方向移動(dòng)，確定產(chǎn)生合理裁剪曲線的選擇標(biāo)準(zhǔn)使拓?fù)浣Y(jié)構(gòu)往更合理方向變動(dòng)，結(jié)構(gòu)變動(dòng)時(shí)無需重新劃分網(wǎng)格。

(3)建立了等幾何裁剪曲線和表面的設(shè)計(jì)控制頂點(diǎn)靈敏度計(jì)算方法，優(yōu)化過程中，內(nèi)部的拓?fù)浣Y(jié)構(gòu)是靈活可變的，提出了內(nèi)部空洞產(chǎn)生和合并算法，常規(guī)的優(yōu)化方法，設(shè)計(jì)空間相關(guān)是一個(gè)不足，提出的樣條優(yōu)化很好解決了這個(gè)問題，在拓?fù)鋬?yōu)化中，對(duì)常規(guī)網(wǎng)格相關(guān)載荷也較難處理，用提出的方法容易處理，在數(shù)字分析和優(yōu)化中，轉(zhuǎn)到CAD模型時(shí)，樣條信息是一樣的，消除后處理的工作量，同時(shí)優(yōu)化拓?fù)浣Y(jié)構(gòu)和形狀，典型實(shí)例表明所用方法的合理性和有效性。該方法容易擴(kuò)展設(shè)計(jì)空間，方便分析處理，高度集成優(yōu)化設(shè)計(jì)，高效可行，具有工程應(yīng)用價(jià)值。

參考文獻(xiàn)

[1]Huang X, Xie Y M. Evolutionary Topology Optimization of Continuum Structures: Methods and Application[M]. London: Springer, 2010.

[2]Bendsoe M P, Sigmund O. Topology optimization: theory, methods, and applications[M]. New York:Springer-Verlag, 2003.

[3]Qian X, Sigmund O.Isogeometric shape optimization of photonic crystals via Coons patches[J]. Computer Methods in Applied Mechanics and Engineering, 2011,200(25-28):2237-2255.

[4]程耿東，許林. 基于可靠度的結(jié)構(gòu)優(yōu)化的序列近似規(guī)劃算法[J].計(jì)算力學(xué)學(xué)報(bào)，2006,23(6):641-646.

CHENG Geng-dong, XU Lin.Sequential approximate programming approach to reliability based structural optimization[J].Chinese Journal of Computational Mechanics 2006,23(6):641-646.

[5]Wang M Y, Wang X, Guo D. A level set method for structural topology optimization[J]. Computer Methods in Applied Mechanics and Engineering, 2003, 192(1-2): 227-246.

[6]Uhm T K, Kim K S, Seo Y D, et al. A locally refinable T-spline finite element method for CAD/CAE integration[J]. Structural Engineering and Mechanics, 2008,30(2):225-245.

[7]Wall W A, Frenzel M A, Cyron C. Isogeometric structural shape optimization[J]. Computer Methods in Applied Mechanics and Engineering, 2008,197(33-40): 2976-2988.

[8]Dedè L, Borden M J, Hughes T J R.Isogeometric analysis for topology optimization with a phase field model[J].Archives of Computational Methods in Engineering,2012,19(3):427-465.

[9]Hassani B, Khanzadi M, Tavakkoli S M.An isogeometrical approach to structural topology optimization by optimality criteria[J]. Structural and Multidisciplinary Optimization, 2012,45(2):223-233.

[10]Lipton S, Evans J A, Bazilevs Y, et al. Robustness of isogeometric structural discretizations under severe mesh distortion[J]. Computer Methods in Applied Mechanics and Engineering 2010, 199(5-8): 357-373.

[11]Hassani B, Tavakkoli S M, Moghadam N Z.Application of isogeometric analysis in structural shape optimization[J].Scientia Iranica.2011,18(4):846-852.

[12]Hughes T J R, Cottrell J A, Bazilevs Y.Isogemetric analysis:CAD, finite element, NURBS,exact geometry and mesh refinement[J]. Computer Methods in Applied Mechanics and Engineering,2005,194(39-41):4135-4195.

[13]Cottrell J A, Reali A, Bazilevs Y, et al. Isogeometric analysis of structural vibrations[J]. Computer Methods in Applied Mechanics and Engineering, 2006, 195(41-43):5257-5296.

[14]Auricchio F, Da Veiga L B, Hughes T J R, et al.Isogeometric collocation methods[J]. Mathematical Models and Methods in Applied Sciences, 2010, 20(11):2075-2107.

[15]Bazilevs Y, Da Veiga L B, Cottrell J,et al. Isogeometric analysis: approximation, stability and error estimates for h-refined meshes[J]. Mathematical Models and Methods in Applied Sciences, 2006, 16(7):1031-1090.

[16]Hughes T J R, Reali A, Sangalli G. Efficient quadrature for NURBS based isogeometric analysis[J]. Computer Methods in Applied Mechanics and Engineering, 2010,199(5-8):301-313.

[17]Hughes T J R, The finite element method: linear static and dynamic finite element analysis[M]. Dover Publications, Mineola, New York:2000.

[18]Roh H Y, Cho M, Integration of geometric design and mechanical analysis using B-spline functions on surface[J]. International Journal for Numerical Methods in Engineering. 2005, 62(14):1927-1949.

[19]Li K, Qian X. Isogeometric analysis and shape optimization via boundary integral[J].Computer Aided Design,2011,43(11):1427-1437.

[20]Hesch C, Betsch P.Isogeometric analysis and domain decomposition methods[J].Computer Methods in Applied Mechanics and Engineering,2012,213-216:104-112.

[21]Qian X, Full analytical sensitivities in NURBS based isogeometric shape optimization[J]. Computer Methods in Applied Mechanics and Engineering.2010, 199(29-32):2059-2071.

[22]Seo Y D, Kim H J, Youn S K. Shape optimization and its extension to topological design based on isogeometric analysis[J]. International Journal of Solids and Structures. 2010, 47(11-12):1618-1640.

[23]Cea J, Garreau S, Guillaume P,et al. The shape and topological optimizations connection[J]. Computer Methods in Applied Mechanics and Engineering.2000,188(4):713-726.

[24]Manh N D, Evgrafov A, Gersborg A R,et al.Isogeometric shape optimization of vibrating membranes[J]. Computer Methods in Applied Mechanics and Engineering, 2011,200(13-16):1343-1353.

[25]Navarrina F, Gómez H,París J,et al.Isogeometric shape sensitivity analysis[J].WIT Transactions on the Built Environment,2012,125:119-130.

[26]Allaire G, Jouve F, Toader A M. Structural optimization using sensitivity analysis and a level set method[J]. Journal of Computational Physics. 2004,194(1):363-393.

[27]傅曉錦.連續(xù)體結(jié)構(gòu)綜合優(yōu)化設(shè)計(jì)[J].機(jī)械工程學(xué)報(bào)，2012,48(1): 128-134.

FU Xiao-jin. Integration optimization design of continuum structure[J].Chinese Journal of Mechanical Engineering, 2012,48(1): 128-134.

[28]傅曉錦,盧炎麟.拓?fù)浜托螤罴蓛?yōu)化設(shè)計(jì)[J].計(jì)算機(jī)集成制造系統(tǒng)，2008,14(6): 1041-1048.

FU Xiao-jin,LU Yan-lin. Integration of topology and shape optimization[J].Computer Integrated Manufacturing Systems, 2008,14(6): 1041-1048.

[29]施法中. 計(jì)算機(jī)輔助幾何設(shè)計(jì)與非均勻有理B樣條[M].北京：高等教育出版社，2001.8.

[30]Piegl L, Tiller W, The NURBS Book (2nd Edition)[M], Springer-Verlag Berlin Heidelberg, Germany: 1997.

[31]朱心雄. 自由曲線曲面造型[M].北京：科學(xué)出版社,2000.

[32]Cho S, Ha S H. Isogeometric shape design optimization: exact geometry and enhanced sensitivity[J]. Structural and Multidisciplinary Optimization. 2009, 38(1):53-70.

[33]Choi K K, Kim N H. Structural sensitivity analysis and optimization 2:nonlinear systems and applications[M]. New York: Springer, 2005.