(西北工業(yè)大學(xué) 航宇材料結(jié)構(gòu)一體化設(shè)計(jì)與增材制造裝備技術(shù)國(guó)際聯(lián)合研究中心，陜西省空天結(jié)構(gòu)技術(shù)重點(diǎn)實(shí)驗(yàn)室，西安 710072)

1 引 言

循環(huán)對(duì)稱結(jié)構(gòu)是一類在極坐標(biāo)系下沿圓周方向循環(huán)對(duì)稱的特殊周期結(jié)構(gòu)，如圖1(a)所示的夾芯點(diǎn)陣薄壁圓筒結(jié)構(gòu)，其輕量化高性能設(shè)計(jì)在航空航天與機(jī)械領(lǐng)域具有重要應(yīng)用價(jià)值。同時(shí)，這類結(jié)構(gòu)也廣泛存在于骨骼和植物等自然界生物材料中[1-3]。然而，過去幾十年建立的結(jié)構(gòu)多尺度優(yōu)化方法主要針對(duì)笛卡爾坐標(biāo)系下周期結(jié)構(gòu)，如圖1(b)所示的夾芯點(diǎn)陣平板結(jié)構(gòu)。其特點(diǎn)為所有點(diǎn)陣微結(jié)構(gòu)單胞均為矩形單胞，可由某一單胞經(jīng)幾何平移操作獲得。因此，所有單胞微結(jié)構(gòu)均可假設(shè)具有相同的等效力學(xué)性能，并通過均勻化方法計(jì)算結(jié)構(gòu)的宏觀力學(xué)響應(yīng)。由此可見，均勻化方法不僅是計(jì)算單胞微結(jié)構(gòu)等效力學(xué)性能的關(guān)鍵步驟，也是關(guān)聯(lián)微結(jié)構(gòu)與宏觀性能和實(shí)現(xiàn)多尺度優(yōu)化的橋梁。目前，基于微結(jié)構(gòu)均勻化方法，相繼發(fā)展了固體各向同性懲罰方法(SIMP)[4,5]、漸進(jìn)結(jié)構(gòu)優(yōu)化方法(ESO)[6,7]和水平集方法(LSM)[8]等微結(jié)構(gòu)拓?fù)鋬?yōu)化方法，同時(shí)也促進(jìn)了功能梯度材料以及考慮尺度效應(yīng)的微結(jié)構(gòu)多目標(biāo)優(yōu)化與拓?fù)鋬?yōu)化研究[9-12]。

圖1 夾芯點(diǎn)陣

本文研究的循環(huán)對(duì)稱結(jié)構(gòu)多尺度拓?fù)鋬?yōu)化問題具有一定的代表性。與平移周期結(jié)構(gòu)相比，可以發(fā)現(xiàn)當(dāng)循環(huán)對(duì)稱結(jié)構(gòu)的半徑趨于無窮大時(shí)，該結(jié)構(gòu)退化為平移周期結(jié)構(gòu)。因此，如何建立循環(huán)對(duì)稱結(jié)構(gòu)的多尺度優(yōu)化方法，從理論上實(shí)現(xiàn)兩類周期結(jié)構(gòu)的有機(jī)統(tǒng)一，不僅具有重要的工程應(yīng)用價(jià)值，也是本文的科學(xué)意義所在。

然而，目前循環(huán)對(duì)稱結(jié)構(gòu)的拓?fù)鋬?yōu)化仍停留在單一尺度上，即直接對(duì)周期微結(jié)構(gòu)[13-18]進(jìn)行優(yōu)化。由于缺乏微結(jié)構(gòu)力學(xué)性能的均勻化等效，計(jì)算效率受到極大挑戰(zhàn)，難以實(shí)現(xiàn)大規(guī)模循環(huán)對(duì)稱結(jié)構(gòu)的多尺度優(yōu)化。Zhang等[19]嘗試將針對(duì)平移周期微結(jié)構(gòu)等效的均勻化方法直接應(yīng)用到循環(huán)對(duì)稱微結(jié)構(gòu)，并未考慮隨半徑位置以及θ方向變化時(shí)單胞微結(jié)構(gòu)性能的不同。Chatzigorgiou等[20,21]提出的改進(jìn)式均勻化方法可用于循環(huán)對(duì)稱周期微結(jié)構(gòu)均勻化等效，并通過矩陣變換得到笛卡爾坐標(biāo)系下相應(yīng)的均勻化等效性能。然而，該方法需要對(duì)不同半徑的單胞進(jìn)行均勻化，其高昂的計(jì)算成本極大限制了應(yīng)用效果，導(dǎo)致相關(guān)工作仍停留在均勻化等效層面上，循環(huán)對(duì)稱結(jié)構(gòu)多尺度拓?fù)鋬?yōu)化研究仍處于起始階段。

2 循環(huán)對(duì)稱結(jié)構(gòu)均勻化方法

2.1 循環(huán)對(duì)稱單胞均勻化等效彈性矩陣計(jì)算

圖2 循環(huán)對(duì)稱結(jié)構(gòu)

(1)

(2)

表1比較了循環(huán)對(duì)稱微結(jié)構(gòu)單胞和平移周期微結(jié)構(gòu)單胞的均勻化方法差異。

表1 循環(huán)對(duì)稱和平移周期微結(jié)構(gòu)的均勻化方法比較Tab.1 Homogenization methods of cyclic symmetrical structure and translational periodic structure

(3)

式中M為二階張量的旋轉(zhuǎn)矩陣。

(4)

2.2 循環(huán)對(duì)稱單胞均勻化等效性能的不變性

由表1可知，對(duì)于循環(huán)對(duì)稱單胞，由于微觀算子包含宏觀半徑r，單胞微結(jié)構(gòu)均勻化等效性能與其所在位置對(duì)應(yīng)的宏觀半徑變量相關(guān)，理論上應(yīng)對(duì)不同宏觀半徑的單胞分別進(jìn)行均勻化等效計(jì)算。

圖3 兩個(gè)相似扇形單胞

(5)

(6)

(7)

比較式(6,7)可得

(8)

(9)

式中本文引入單胞無量綱特征參數(shù)λ用于描述扇形單胞中心處弧長(zhǎng)與邊長(zhǎng)之比。其物理意義為，若兩相似單胞的特征參數(shù)λ相同，即使宏觀半徑不同，其均勻化等效性能仍然相同。由此可知，圖2所示的循環(huán)對(duì)稱結(jié)構(gòu)，若單胞微觀邊長(zhǎng)與其宏觀半徑呈正比例關(guān)系，則其均勻化等效性能保持相同(推導(dǎo)忽略)。

2.3 循環(huán)對(duì)稱單胞均勻化等效性能的可映射性

受文獻(xiàn)[23，24]關(guān)于平移周期單胞存在的映射關(guān)系啟發(fā)，本文基于特征參數(shù)λ建立了循環(huán)對(duì)稱單胞等效性能計(jì)算的映射關(guān)系。其重要意義為，對(duì)于不同λ的扇形單胞，無需通過式(1)每次重新劃分網(wǎng)格計(jì)算其等效性能，只需在標(biāo)準(zhǔn)單胞(λ=1)上計(jì)算求解等效性能。其中的基體材料屬性以及最終單胞等效性能通過相應(yīng)映射計(jì)算獲得，由此極大簡(jiǎn)化了不同λ取值單胞的均勻化過程。

圖4所示為標(biāo)準(zhǔn)單胞(λ=1)和任意λ取值單胞，參數(shù)分別使用上標(biāo)^和～區(qū)分。存在如下幾何關(guān)系

圖4 單胞映射變化

(10)

通過理論推導(dǎo)，任意λ取值單胞的均勻化等效性能為

(11)

式中Ti j為映射張量的分量，映射矩陣T為

(i,j=1,2)(12)

(13)

(14)

(15)

2.4 循環(huán)對(duì)稱單胞與平移周期單胞的統(tǒng)一性

考慮圖5所示的循環(huán)對(duì)稱單胞與平移周期單胞。

圖5 不同周期單胞

平移周期單胞尺寸分別滿足

(16)

(17)

(18)

即表1兩坐標(biāo)系下微觀算子等價(jià)。因此，當(dāng)r*→ ∞時(shí)，循環(huán)對(duì)稱單胞與平移周期單胞均勻化等效性能等價(jià)。

3 單胞微結(jié)構(gòu)參數(shù)化建模與等效

性能插值計(jì)算

3.1 單胞微結(jié)構(gòu)的特征驅(qū)動(dòng)參數(shù)化建模

基于文獻(xiàn)[25,26]的前期研究工作，本文首次將特征驅(qū)動(dòng)拓?fù)鋬?yōu)化方法應(yīng)用于循環(huán)對(duì)稱單胞微結(jié)構(gòu)的多尺度拓?fù)鋬?yōu)化。圖6所示為超橢圓特征，其水平集函數(shù)為

圖6 超橢圓特征水平集函數(shù)及其相關(guān)設(shè)計(jì)變量

(19)

(20)

(21)

圖7(左)所示為具有最大體積模量和最大剪切模量的兩種平移周期單胞微結(jié)構(gòu)[27]。對(duì)其分別采用相同數(shù)目的超橢圓特征進(jìn)行參數(shù)化建模，則單胞微結(jié)構(gòu)包含的超橢圓設(shè)計(jì)變量集合為s。

圖7 兩種典型的微結(jié)構(gòu)

(22)

式中(θc i,rc i,ac i,wc i和lc i)T代表第i個(gè)超橢圓設(shè)計(jì)參數(shù)。生成的單胞微結(jié)構(gòu)整體水平集函數(shù)可通過所有超橢圓特征水平集函數(shù)的布爾運(yùn)算獲得

(23)

采用以上最大體積模量和最大剪切模量單胞微結(jié)構(gòu)分別作為初始和最終微結(jié)構(gòu)構(gòu)型，則任意中間態(tài)單胞微結(jié)構(gòu)構(gòu)型可通過設(shè)計(jì)變量插值確定。

sκ=sinitial+κ(sfinal-sinitial)

(24)

圖8 由超橢圓特征組成的單胞微結(jié)構(gòu)

(25)

極坐標(biāo)系下，單胞微結(jié)構(gòu)的體積分?jǐn)?shù)計(jì)算式為

(26)

式中H(·)為海維賽德函數(shù)。

引入超橢圓寬度比例變量ω用于改變單胞體積分?jǐn)?shù)，則有

sω=ωs,ω=diag(1,1,1,ω,1)

(27)

(28)

圖9 兩種不同演變形式的參數(shù)化微結(jié)構(gòu)

(29)

(30)

3.2 參數(shù)化微結(jié)構(gòu)等效力學(xué)性能的插值計(jì)算

單胞微結(jié)構(gòu)的等效力學(xué)性能不僅取決于微結(jié)構(gòu)構(gòu)型和材料用量，還與其弧邊比特征參數(shù)λ緊密相關(guān)。為了定量地描述參數(shù)化微結(jié)構(gòu)的等效力學(xué)性能與微結(jié)構(gòu)控制參數(shù)的依賴關(guān)系和減少重復(fù)計(jì)算成本，本文基于函數(shù)擬合插值模型[21]構(gòu)造了等效力學(xué)性能與微結(jié)構(gòu)控制參數(shù)的近似表達(dá)式，具體采用三元四次多項(xiàng)式插值方式描述三個(gè)控制參數(shù)與等效力學(xué)性能的關(guān)系。

(31)

式中n=4，βI J L為對(duì)應(yīng)系數(shù)。

圖11給出了均勻化等效性能C11的剖視圖與λ=1透明截面處的截面視圖。圓圈為原始樣本點(diǎn)均勻化等效數(shù)據(jù)，均位于擬合曲面附近，具有良好的擬合精度。

圖11 均勻化等效性能C11

4 密度與插值參數(shù)κ空間分布的B樣條參數(shù)化建模與優(yōu)化模型

4.1 基于B樣條參數(shù)化的多尺度設(shè)計(jì)框架

(32)

[sinitial+κ(sfinal-sinitial)]

(33)

如圖12所示，由于相鄰單胞在公共邊界上共享相同參數(shù)變量，B樣條可以自動(dòng)保證微結(jié)構(gòu)的幾何連續(xù)性，無需施加任何額外約束。而傳統(tǒng)密度法由于不同單胞微結(jié)構(gòu)的變量離散化定義，即使對(duì)微結(jié)構(gòu)進(jìn)行細(xì)分也無法從本質(zhì)上改變相鄰微結(jié)構(gòu)之間連續(xù)性差的問題。

圖12 微結(jié)構(gòu)連接光順性的比較

4.2 循環(huán)對(duì)稱結(jié)構(gòu)多尺度拓?fù)鋬?yōu)化

研究給定體積約束下結(jié)構(gòu)柔順度最小化問題。對(duì)應(yīng)的數(shù)學(xué)模型為

min.C=UTKU

(34)

式中λm為中心線半徑處的弧邊比特征參數(shù)。

圖13為基于B樣條參數(shù)化的循環(huán)對(duì)稱結(jié)構(gòu)多尺度優(yōu)化框架。

圖13 循環(huán)對(duì)稱結(jié)構(gòu)多尺度優(yōu)化流程

5 數(shù)值算例

5.1 曲梁多尺度拓?fù)鋬?yōu)化

圖14所示曲梁左端固支，右下端加載集中力F=1 N。軸線長(zhǎng)l=600 mm，寬度h=200 mm，r0=1200/πmm，固體材料楊氏模量E=1 MPa,泊松比為0.3。假設(shè)材料用量體分比為50%，曲梁設(shè)計(jì)域在極坐標(biāo)系下劃分為20×60個(gè)單胞，并填充如圖9(左)所示的參數(shù)化微結(jié)構(gòu)，中心線半徑處λm=1。當(dāng)曲梁半徑R=r0時(shí)，曲梁扇形角為90°。隨著R的增大，曲梁的曲率逐漸變小。當(dāng)R→ ∞時(shí)，曲梁的曲率趨近于零并退化為懸臂梁，其均勻化等效性能也與傳統(tǒng)平移周期微結(jié)構(gòu)一致。

圖14 受集中力的曲梁

針對(duì)上述三種情況，將不同R的曲梁均置于相同B樣條參數(shù)域中進(jìn)行拓?fù)鋬?yōu)化，其控制點(diǎn)位置不隨R的變化而發(fā)生改變，如圖15所示。

圖15 不同半徑曲梁置于同一個(gè)B樣條參數(shù)域

為了驗(yàn)證該方法的正確性，說明平移周期結(jié)構(gòu)是循環(huán)對(duì)稱結(jié)構(gòu)的特例，取R=1×104r0，使得曲梁接近于懸臂梁的形式。此時(shí)，曲梁的弧度約為1.57×10-4，幾乎可忽略不計(jì)。多尺度優(yōu)化結(jié)果和迭代曲線如圖16和圖17所示。前者與VCUT水平集方法[29]一致，有效驗(yàn)證了循環(huán)對(duì)稱結(jié)構(gòu)多尺度拓?fù)鋬?yōu)化方法的正確性。

圖16 R=1×104r0的曲梁多尺度優(yōu)化結(jié)果

圖17 R=1×104r0時(shí)曲梁柔順度和體積分?jǐn)?shù)的迭代曲線

表2給出了不同R下曲梁的多尺度優(yōu)化結(jié)果。

可以發(fā)現(xiàn)，材料主要沿上下表面的加載路徑分布，且局部微結(jié)構(gòu)構(gòu)型與其主應(yīng)力方向相關(guān)，曲梁的上下表面構(gòu)型參數(shù)趨向于0，以橫豎桿為主。在上下主傳力路徑之間，存在一定的過渡結(jié)構(gòu)。從 圖18 可以看出，隨著半徑的增加，柔順度逐漸增加；當(dāng)半徑趨于無窮大時(shí)，逐漸收斂到懸臂梁的柔順度值。

圖18 不同R值曲梁結(jié)構(gòu)柔順度

5.2 圓盤多尺度拓?fù)鋬?yōu)化

該問題著重研究不同循環(huán)對(duì)稱單胞排列方式、弧邊比參數(shù)以及曲率半徑對(duì)結(jié)構(gòu)性能的影響。

如圖19(a)所示，圓環(huán)內(nèi)圈固定，圓環(huán)外圈施加切向載荷，外徑與內(nèi)徑之差2(R0-r0)=50 mm。假設(shè)材料體積分?jǐn)?shù)限制為40%，基體材料楊氏模量E=10 GPa，泊松比為0.3。圖19(b，c)分別給出了等間距排列單胞和等比例排列單胞。最初將圓盤劃分為6×54個(gè)單胞，優(yōu)化過程中徑向?qū)訑?shù)保持不變，圓周方向的數(shù)目與λm相關(guān)。本算例假定r0分別為15 mm,25 mm,50 mm,75 mm，λm的初始值分別為0.54,1.01,1.69,2.36。

圖19 圓盤問題

對(duì)于等間距排列單胞，位于宏觀坐標(biāo)(r,θ)的微結(jié)構(gòu)弧邊比特征參數(shù)λ與中心線弧邊比特征參數(shù)λm存在如下關(guān)系，

λ=2r/(r0+R0)λm

(35)

等間距排列單胞在不同內(nèi)徑下的優(yōu)化設(shè)計(jì)結(jié)果如圖20所示?？梢钥闯?，(a) 微結(jié)構(gòu)插值構(gòu)型參數(shù)κ始終接近1，形成X形的微結(jié)構(gòu)。(b)λm總是收斂于1附近，表明具有相同弧長(zhǎng)和邊長(zhǎng)的扇形單胞更有利于結(jié)構(gòu)優(yōu)化。(c) 無論內(nèi)徑長(zhǎng)度如何，材料均在內(nèi)徑處堆積，微結(jié)構(gòu)體積分?jǐn)?shù)從內(nèi)徑到外徑逐漸減小。 (d) 內(nèi)徑較小時(shí)，內(nèi)徑和外徑的體積分?jǐn)?shù)變化很大。隨著內(nèi)徑的增加，內(nèi)徑與外徑的體積分?jǐn)?shù)差逐漸減小。

圖20 等間距排列單胞在不同內(nèi)徑下的優(yōu)化設(shè)計(jì)結(jié)果

圖21 內(nèi)徑趨于無窮大時(shí)的結(jié)構(gòu)演化

同樣，當(dāng)r0=15 mm,25 mm,50 mm,75 mm且單胞以等比例方式排列時(shí)，可以得到如圖22所示與等間距排列單胞相似的優(yōu)化結(jié)果。該排列方式確保參數(shù)λ在圓的所有位置均等于中心線處弧邊比參數(shù)λm。

圖22 等比例排列單胞在不同內(nèi)徑下的優(yōu)化設(shè)計(jì)結(jié)果

表3比較了相同內(nèi)徑下不同單胞排列方式對(duì)柔順度的影響。柔順度結(jié)果表明，等比例排列單胞明顯優(yōu)于等間距排列單胞。當(dāng)內(nèi)徑較小時(shí)，這兩種排列結(jié)構(gòu)的柔順度差異很明顯，并在r0=15 mm時(shí)二者柔順度降幅達(dá)到了6.37%。隨著內(nèi)徑的增大，二者柔順度的差異逐漸減小。當(dāng)r0→ ∞時(shí)，兩種排列方式之間的差異將不再明顯。因此，在環(huán)向剪切載荷作用下，等比例排列比等間距排列具有更大的優(yōu)勢(shì)，在工程中對(duì)于受切向載荷的圓筒結(jié)構(gòu)可優(yōu)先考慮使用等比例排列單胞設(shè)計(jì)。

表3 不同半徑不同排列方式柔順度差異Tab.3 Differences in compliance under different radii and different arrangements

(續(xù)表)

6 結(jié) 論

本文開展了循環(huán)對(duì)稱結(jié)構(gòu)多尺度拓?fù)鋬?yōu)化方法研究，提出了用于表征等效性能映射計(jì)算的弧邊比特征參數(shù)，闡明了循環(huán)對(duì)稱單胞均勻化等效性能的不變性、可映射性以及與平移周期單胞的統(tǒng)一性，建立了微結(jié)構(gòu)單胞等效性能的三元擬合插值模型與B樣條參數(shù)化多尺度設(shè)計(jì)框架。

算例結(jié)果表明，對(duì)于受剪切載荷作用的循環(huán)對(duì)稱結(jié)構(gòu)，等比例單胞排列設(shè)計(jì)優(yōu)于等間距單胞排列設(shè)計(jì)，且當(dāng)結(jié)構(gòu)半徑趨于無窮大時(shí)，循環(huán)對(duì)稱結(jié)構(gòu)退化為平移周期性結(jié)構(gòu)。

參考文獻(xiàn)(References):

[1] Liu Z Q,Meyers M A,Zhang Z F,et al.Functional gradients and heterogeneities in biological materials:Design principles,functions,and bioinspired applications[J].ProgressinMaterialsScience,2017,88:467-498.

[2] Osorio L,Trujillo E,van Vuure A W,et al.Morphological aspects and mechanical properties of single bamboo fibers and flexural characterization of bamboo/epoxy composites[J].JournalofReinforcedPlasticsandComposites,2011,30(5):396-408.

[3] Chen W J,Zheng X N,Liu S T.Finite -element-mesh based method for modeling and optimization of lattice struc-tures for additive manufacturing[J].Materials,2018,11(11):2073.

[4] Bends?e M P,Sigmund O.Material interpolation schemes in topology optimization [J].ArchiveofApp-liedMechanics,1999,69(9-10):635-654.

[5] Sigmund O.Tailoring materials with prescribed elastic properties[J].MechanicsofMaterials,1995,20(4):351-368.

[6] Xia L,Breitkopf P.Recent advances on topology optimization of multiscale nonlinear structures[J].ArchivesofComputationalMethodsinEngineering,2017,24(2):227-249.

[7] Radman A,Huang X,Xie Y M.Topology optimization of functionally graded cellular materials[J].JournalofMaterialsScience,2013,48(4):1503-1510.

[8] Li H,Luo Z,Gao L,et al.Topology optimization for functionally graded cellular composites with metamaterials by level sets [J].ComputerMethodsinApp-liedMechanicsandEngineering,2018,328:340-364.

[9] 邱克鵬,張衛(wèi)紅.功能梯度材料結(jié)構(gòu)拓?fù)鋬?yōu)化設(shè)計(jì)研究[J].西北工業(yè)大學(xué)學(xué)報(bào),2010,28(6):851-857.(QIU Ke -peng,ZHANG Wei-hong.Improving topo -logical design of FGM (functionally graded material) structure[J].JournalofNorthwesternPolytechnicalUniversity,2010,28(6):851-857.(in Chinese))

[10] Wang Y Q,Zhang L,Daynes S,et al.Design of graded lattice structure with optimized mesostructures for additive manufacturing [J].Materials&Design,2018,142:114-123.

[11] 孫士平,張衛(wèi)紅.多相材料微結(jié)構(gòu)多目標(biāo)拓?fù)鋬?yōu)化設(shè)計(jì)[J].力學(xué)學(xué)報(bào),2006,38(5):633-638.(SUN Shi-ping,ZHANG Wei-hong.Multiple objective topology optimal design of multiphase microstructures[J].ChineseJournalofTheoreticalandAppliedMe-chanics,2006,38(5):633-638.(in Chinese))

[12] 張衛(wèi)紅,駱金威,戴高明,等.周期性多孔材料等效剪切模量與尺寸效應(yīng)研究[J].力學(xué)學(xué)報(bào),2011,43(1):144-153.(ZHANG Wei-hong,LUO Jin-wei,DAI Gao-ming,et al.Numerical predictions of effective shear modulus and size effect for periodic cellular materials [J].ChineseJournalofTheoreticalandApp-liedMechanics,2011,43(1):144-153.(in Chinese))

[13] Tang Y L,Dong G Y,Zhao Y F.A hybrid geometric modeling method for lattice structures fabricated by additive manufacturing[J].TheInternationalJournalofAdvancedManufacturingTechnology,2019,102(9-12):4011-4030.

[14] Dedè L,Borden M J,Hughes T J R.Isogeometric analysis for topology optimization with a phase field model[J].ArchivesofComputationalMethodsinEngineering,2012,19(3):427-465.

[15] Moses E,Fuchs M B,Ryvkin M.Topological design of modular structures under arbitrary loading[J].StructuralandMultidisciplinaryOptimization,2002,24(6):407-417.

[16] Gao T,Zhang W H.Topology optimization involving thermo -elastic stress loads[J].StructuralandMultidisciplinaryOptimization,2010,42(5):725-738.

[17] Zuo Z H,Xie Y M,Huang X D.Optimal topological design of periodic structures for natural frequencies[J].JournalofStructuralEngineering,2011,137(10):1229-1240.

[18] Maharaj Y,James K A.Metamaterial topology optimization of nonpneumatic tires with stress and buckling constraints[J].InternationalJournalforNumericalMethodsinEngineering,2020,121(7):1410-1439.

[19] Zhang H K,Wang Y G,Kang Z.Topology optimization for concurrent design of layer-wise graded lattice materials and structures[J].InternationalJournalofEngineeringScience,2019,138:26-49.

[20] Chatzigeorgiou G,Efendiev Y,Lagoudas D C.Homo -genization of aligned “fuzzy fiber” composites[J].InternationalJournalofSolidsandStructures,2011,48(19):2668-2680.

[21] Wang C,Zhu J H,Zhang W H,et al.Concurrent topology optimization design of structures and non-uniform parameterized lattice microstructures[J].StructuralandMultidisciplinaryOptimization,2018,58(1):35-50.

[22] Xu Z,Zhang W H,Gao T,et al.A B-spline multi-para-meterization method for multi-material topology optimization of thermoelastic structures[J].Struc-turalandMultidisciplinaryOptimization,2020,61(3):923-942.

[23] Chatzigeorgiou G,Efendiev Y,Charalambakis N,et al.Effective thermoelastic properties of composites with periodicity in cylindrical coordinates[J].InternationalJournalofSolidsandStructures,2012,49(18):2590-2603.

[24] Li Q H,Xu R,Liu J,et al.Topology optimization de -sign of multi-scale structures with alterable microstructural length-width ratios[J].CompositeStructures,2019,230:111454.

[25] Zhou Y,Zhang W H,Zhu J H,et al.Feature -driven topology optimization method with signed distance function[J].ComputerMethodsinAppliedMecha-nicsandEngineering,2016,310:1-32.

[26] Xu Z,Zhang W H,Zhou Y,et al.Multiscale topology optimization using feature -driven method[J].ChineseJournalofAeronautics,2020,33(2):621-633.

[27] Xia L,Breitkopf P.Design of materials using topology optimization and energy-based homogenization app -roach in matlab[J].StructuralandMultidisciplinaryOptimization,2015,52(6):1229-1241.

[28] Liu S T,Cheng G D,Gu Y,et al.Mapping method for sensitivity analysis of composite material property[J].StructuralandMultidisciplinaryOptimization,2002,24(3):212-217.

[29] Zong H M,Liu H,Ma Q P,et al.VCUT level set method for topology optimization of functionally graded cellular structures[J].ComputerMethodsinAppliedMechanicsandEngineering,2019,354:487-505.