考慮時(shí)變剛度特性的復(fù)合材料微結(jié)構(gòu)拓?fù)鋬?yōu)化設(shè)計(jì)方法1)

徐世鵬 丁曉紅 段朋云 張 橫

(上海理工大學(xué)機(jī)械工程學(xué)院,上海 200093)

引言

時(shí)變剛度是指結(jié)構(gòu)剛度隨時(shí)間發(fā)生變化,以滿足特定的使用場(chǎng)景,如在一些骨科和心血管疾病的治療時(shí),需要植入物在植入初期保持足夠大的剛度,起到有效的固定或支撐作用,隨著機(jī)體的康復(fù),需要植入物的剛度逐漸降低,避免出現(xiàn)應(yīng)力遮擋效應(yīng)、支架內(nèi)血栓等問題[1].因此,具有時(shí)變剛度的結(jié)構(gòu)是一種理想的臨時(shí)性醫(yī)療植入物.隨著材料科學(xué)的發(fā)展,生物可降解材料被應(yīng)用于臨時(shí)性植入物的設(shè)計(jì),使其具有時(shí)變剛度特性[2-3],但需要保證時(shí)變剛度特性能夠滿足骨骼或血管不同修復(fù)階段特定的生物力學(xué)需求.以骨折內(nèi)固定植入物為例,骨折愈合大致可以分為骨痂生長(zhǎng)和骨痂重塑兩個(gè)階段,不同階段具有不同的生物力學(xué)需求,同時(shí)隨著骨愈合的進(jìn)行骨痂的剛度逐漸增大,而可降解植入物的剛度逐漸減小[4-5].因此,在設(shè)計(jì)可降解植入物時(shí)需要考慮使植入物時(shí)變剛度特性與骨骼或血管不同修復(fù)階段對(duì)不同生物力學(xué)刺激的需求相匹配.

目前針對(duì)可降解植入物的設(shè)計(jì)研究主要集中在生物材料的研究,常通過調(diào)整材料配比或加工工藝等方式控制材料性能,但是單一材料通常具有相對(duì)恒定的降解速率[6],因此使用單一材料設(shè)計(jì)的植入物不能實(shí)現(xiàn)結(jié)構(gòu)剛度在不同愈合階段具有不同且特定的時(shí)變速率的設(shè)計(jì)要求.使用由多種可降解材料組成的復(fù)合結(jié)構(gòu)是一種有效的解決方案.連續(xù)體拓?fù)鋬?yōu)化方法在基于微結(jié)構(gòu)的功能梯度結(jié)構(gòu)[7]、負(fù)泊松比材料[8]、超材料[9]和復(fù)合材料的單尺度[10-11]與多尺度[12]等設(shè)計(jì)問題中已經(jīng)得到了有效應(yīng)用,但涉及可降解材料的微結(jié)構(gòu)設(shè)計(jì)的研究極少.這是因?yàn)榭紤]時(shí)變剛度特性的微結(jié)構(gòu)拓?fù)鋬?yōu)化不僅涉及到結(jié)構(gòu)的降解模型及在時(shí)間域內(nèi)的結(jié)構(gòu)性能分析,而且由于可降解結(jié)構(gòu)性能非線性強(qiáng),如何建立合理的優(yōu)化數(shù)學(xué)模型并進(jìn)行靈敏度分析,得到清晰的拓?fù)鋬?yōu)化結(jié)果也是一個(gè)難題.目前在一些拓?fù)鋬?yōu)化問題的研究中使用了時(shí)間維度的變量,如James 和Waisman[13]針對(duì)具有蠕變特性結(jié)構(gòu)的設(shè)計(jì),提出了使用時(shí)間伴隨方法考慮時(shí)間維度變量的結(jié)構(gòu)拓?fù)鋬?yōu)化設(shè)計(jì)方法,有效地降低了結(jié)構(gòu)重量.Wang 等[14]提出了同時(shí)進(jìn)行結(jié)構(gòu)拓?fù)鋬?yōu)化和增材制造順序優(yōu)化的方法,使用密度變量定義結(jié)構(gòu)布局,使用時(shí)間場(chǎng)變量定義制造順序,同時(shí)對(duì)兩組變量進(jìn)行優(yōu)化,為拓?fù)鋬?yōu)化和增材制造的結(jié)合提供了新的思路.Wu 等[15]研究了下頜骨接骨板的拓?fù)鋬?yōu)化方法,考慮了接骨板結(jié)構(gòu)拓?fù)浜凸侵厮艿鸟詈?并推導(dǎo)了時(shí)域的靈敏度,其中骨重塑過程具有時(shí)間維度的變化.但上述考慮時(shí)域的研究并未研究材料降解的問題.考慮材料降解的結(jié)構(gòu)優(yōu)化設(shè)計(jì)中,結(jié)構(gòu)的優(yōu)化和降解分別進(jìn)行,但兩者又互相影響:結(jié)構(gòu)的幾何形態(tài)和尺寸會(huì)影響整體結(jié)構(gòu)的降解速度,同時(shí)結(jié)構(gòu)的降解又會(huì)影響不同降解時(shí)間步的結(jié)構(gòu)瞬時(shí)力學(xué)性能.對(duì)于可降解結(jié)構(gòu)的拓?fù)鋬?yōu)化設(shè)計(jì),模擬結(jié)構(gòu)降解過程并量化降解對(duì)結(jié)構(gòu)性能的影響至關(guān)重要,降解過程中結(jié)構(gòu)性能隨時(shí)間變化,需要在時(shí)間域內(nèi)對(duì)降解過程進(jìn)行建模[16-18].關(guān)于降解過程的模擬,Grogan 等[19]提出了使用擴(kuò)散方程模擬材料降解的方法,Guo 等[20]提出了一種基于元胞自動(dòng)機(jī)制的模擬多孔支架降解的方法,但上述文獻(xiàn)并未研究復(fù)合結(jié)構(gòu)的降解問題,且降解過程涉及的參數(shù)都需要通過實(shí)驗(yàn)給定,不同的材料具有不同的參數(shù).

本文提出一種考慮時(shí)變剛度特性的復(fù)合材料微結(jié)構(gòu)拓?fù)鋬?yōu)化設(shè)計(jì)方法,在時(shí)間維度對(duì)微結(jié)構(gòu)降解過程進(jìn)行模擬,建立考慮材料降解的有限元模型,在不同降解時(shí)間步利用均勻化方法計(jì)算得到中間結(jié)構(gòu)的力學(xué)性能,優(yōu)化目標(biāo)是某些特定降解時(shí)間步中間結(jié)構(gòu)的剛度之和最大.

1 考慮材料降解的均勻化方法

1.1 具有特定時(shí)變剛度特性的復(fù)合材料結(jié)構(gòu)

根據(jù)Zhao 等[2]提出的理想可降解骨折內(nèi)固定植入物剛度的變化趨勢(shì),簡(jiǎn)述考慮時(shí)變剛度特性的復(fù)合材料微結(jié)構(gòu)設(shè)計(jì)構(gòu)思.骨折愈合是一個(gè)復(fù)雜的生物學(xué)過程,大致可以分為骨痂生長(zhǎng)和骨痂重塑兩個(gè)階段,如圖1 所示橫坐標(biāo)表示骨愈合的時(shí)間歷程,縱坐標(biāo)表示剛度,虛線描述了在骨愈合過程中骨痂整體剛度的變化趨勢(shì),在生長(zhǎng)階段先緩慢增大,然后以相對(duì)較快的速率增大,在重塑階段再逐漸增大.為了順應(yīng)骨愈合的規(guī)律,理想的可降解內(nèi)固定植入物的剛度應(yīng)如圖中實(shí)線所示,在骨痂生長(zhǎng)階段剛度基本保持不變,在骨痂重塑期以一定速率逐漸降低,在骨愈合后能完全降解,避免二次手術(shù).

圖1 理想骨折內(nèi)固定植入物和愈傷組織整體剛度隨時(shí)間變化趨勢(shì)Fig.1 The change trend of stiffness of ideal fracture internal fixation implant and callus with time

為了使植入物剛度達(dá)到上述變化趨勢(shì),使用兩種生物降解材料,其中材料1 剛度較大、降解速率較小,材料2 剛度較小、降解速率較大.具有特定時(shí)變剛度特性的復(fù)合材料結(jié)構(gòu)示意圖如圖2 所示,由材料1 和2 構(gòu)成,該優(yōu)化設(shè)計(jì)問題可以描述為在給定的優(yōu)化數(shù)學(xué)模型下,在微觀尺度上確定兩種材料的最優(yōu)布局,使可降解結(jié)構(gòu)達(dá)到所需的時(shí)變剛度特性.

圖2 具有特定時(shí)變剛度特性的復(fù)合材料結(jié)構(gòu)Fig.2 Periodic composite structure with specific time-changeable characteristics

為了實(shí)現(xiàn)上述設(shè)計(jì)目標(biāo),模擬微結(jié)構(gòu)降解過程并量化材料降解對(duì)微結(jié)構(gòu)性能的影響是研究的關(guān)鍵.分別使用變量xe和ye描述微結(jié)構(gòu)的材料分布和降解狀態(tài),判定依據(jù)如式(1) 和式(2) 所示,其中xe表示相對(duì)密度,表征單元e為材料1 還是材料2;ye表示單元e的殘留率,如果ye=1,則單元e未降解,如果ye=ymin則單元e完全降解,當(dāng)ye處于ymin到1 的中間值時(shí),表示單元e部分降解.

式中,xmin和ymin分別是為避免計(jì)算中數(shù)值不穩(wěn)定預(yù)設(shè)的最小值,取0.001,m為單元總數(shù)量.

1.2 結(jié)構(gòu)降解過程模擬

均勻腐蝕方法是一種常見的基于現(xiàn)象的降解模擬方法,使用平均降解率模擬結(jié)構(gòu)的降解,相關(guān)學(xué)者已通過實(shí)驗(yàn)對(duì)該方法的有效性進(jìn)行了驗(yàn)證[21-22].由于可降解材料與腐蝕環(huán)境接觸才會(huì)發(fā)生降解,將接觸邊界稱為可降解邊界,本文采用均勻腐蝕方法,為簡(jiǎn)化降解過程,做出以下假定:(1)假定降解均勻進(jìn)行且只發(fā)生在可降解邊界處,可降解邊界隨著降解過程的進(jìn)行而移動(dòng),直到所有材料完成降解;(2)假定降解過程僅受材料降解速率控制,不受其他物理參數(shù)(例如流體速度和溫度等)的影響.

降解模擬過程如圖3 所示,紅色虛線表示可降解邊界,降解過程中材料可以分為3 種類型:(1)固態(tài),表示材料未降解;(2)固液混合態(tài),表示材料部分降解;(3)液態(tài),表示材料完全降解.從圖3(a)到圖3(c),可降解邊界處的單元逐漸降解,從圖3(c)到圖3(d),可降解邊界移動(dòng)到第二層單元,從圖3(d)到圖3(f)第二層單元逐漸降解,可降解邊界移動(dòng)到下一層單元,重復(fù)上述過程,直到所有單元全部降解.

圖3 降解過程示意圖Fig.3 Schematic diagram of degradation process

在每一個(gè)降解時(shí)間步使用均勻化方法對(duì)具有中間殘留率的結(jié)構(gòu)進(jìn)行分析,獲得微結(jié)構(gòu)性能隨時(shí)間的變化.單元降解的數(shù)學(xué)表達(dá)為

式中,i為降解時(shí)間步,de為單元e的降解速率,t為單位降解時(shí)間步的時(shí)長(zhǎng).

如圖4 所示,使用k和l將單元e的位置標(biāo)記為(k,l),與單元e相鄰的4 個(gè)單元為(k-1,l),(k+1,l),(k,l-1)和(k,l+1),則相鄰單元的殘留率能夠表示為yk-1,l,yk+1,l,yk,l-1和yk,l+1.僅當(dāng)單元處于可降解邊界時(shí)該單元能夠降解,為了保證數(shù)值計(jì)算的穩(wěn)定,設(shè)定yb=0.1,假定當(dāng)單元的殘留率小于yb時(shí)可降解邊界發(fā)生移動(dòng),如當(dāng)yk-1,l＜yb時(shí),認(rèn)為單元e能夠從左側(cè)開始降解,當(dāng)yk+1,l＜yb時(shí),認(rèn)為單元e能夠從右側(cè)開始降解,上下方向同理.根據(jù)相鄰單元的降解情況,單元e的降解速率為

圖4 單元e及其相鄰單元位置標(biāo)記Fig.4 Location marking of element eand its neighbors

式中,Ω為與單元e相鄰單元的集合,α為Ω中的第α個(gè)單元,d為材料的降解速率.邏輯方程為

為了便于求導(dǎo),使用Heaviside 方程[23-24]將式(5)近似為連續(xù)方程,式(4)可以改寫為

其中H為

式中,γ為常數(shù).

由式(3)和式(6),單元e殘留率的更新方程為

降解模擬中,時(shí)間步長(zhǎng)越小模擬結(jié)果越精確,為了保證式(8) 的連續(xù)性,使用Kreisselmeier-Steinhauser (KS)方程將最大函數(shù)問題近似[25]為

式中,η為常數(shù),γ和η值越大近似誤差越小,但是同時(shí)也使問題的非線性度變大,導(dǎo)致很難使用基于梯度的優(yōu)化算法求解優(yōu)化問題[26].根據(jù)數(shù)值測(cè)試,γ=12 和η=12 能夠滿足本文的需求.

降解模擬算法可以總結(jié)為如下步驟.

第2 步:根據(jù)yk-1,l,yk+1,l,yk,l-1和yk,l+1判斷相鄰單元的降解狀態(tài);

第3 步:使用式(6)計(jì)算單元e的降解速率de;

第4 步:使用式(9)更新單元e的殘留率;

重復(fù)第1 到第4 步,直到指定時(shí)間步.

1.3 材料插值模型

為了獲得清晰的結(jié)構(gòu),使用基于SIMP 的材料插值機(jī)制[27-28],彈性模量E和降解速率d能夠表示為

式中,p為懲罰系數(shù),取3.0,下標(biāo)1 和2 分別代表可降解材料1 和2.

1.4 均勻化方法

利用均勻化方法計(jì)算單胞的等效彈性矩陣[29],通過施加周期性邊界條件對(duì)單胞進(jìn)行有限元分析,結(jié)構(gòu)平衡方程為

式中,K為單胞的整體剛度矩陣,由于材料的降解,該矩陣與時(shí)間相關(guān),U為單胞的位移場(chǎng),f為由均勻應(yīng)變場(chǎng)引起的外力.

對(duì)于可降解的單胞結(jié)構(gòu),K由單元?jiǎng)偠染仃噆e與和ye的乘積組裝得到,其中ye是關(guān)于時(shí)間t和降解速率d的函數(shù).假設(shè)殘留率y與材料剛度呈線性關(guān)系,將材料降解對(duì)力學(xué)性能的影響整合到有限元分析過程中,可降解結(jié)構(gòu)的剛度矩陣為

式中,?為單胞,?e為單胞中第e個(gè)單元的體積,B為應(yīng)變-位移矩陣,De為單元e的彈性矩陣.對(duì)于平面應(yīng)變問題,D為

式中,E為彈性模量,μ為泊松比.則f可以表示為

式中,針對(duì)二維問題,應(yīng)變?chǔ)?[ε1,ε2,ε3],ε1=[1,0,0]T,ε2=[0,1,0]T,ε3=[0,0,1]T.

當(dāng)?shù)玫轿灰浦?利用均勻化方法,均勻彈性矩陣DH可以表示為

式中,|?|為單胞的總體積,為根據(jù)式(17)通過施加單位應(yīng)變計(jì)算得到的節(jié)點(diǎn)位移場(chǎng);Ue為根據(jù)式(12)通過施加全局單位應(yīng)變計(jì)算得到的位移場(chǎng)

根據(jù)式(16),矩陣DH中的彈性張量Drs可以表示為

式中,r=1,2,3,s=1,2,3.

2 優(yōu)化數(shù)學(xué)模型及靈敏度分析

2.1 優(yōu)化數(shù)學(xué)模型

設(shè)定一個(gè)降解時(shí)間步長(zhǎng)為一天,為了使復(fù)合結(jié)構(gòu)剛度的變化趨勢(shì)接近于圖1 中實(shí)線所示的趨勢(shì),即在骨痂生長(zhǎng)期剛度下降值小,則可將在第1 和第60 降解步的復(fù)合結(jié)構(gòu)剛度之和最大作為優(yōu)化目標(biāo),以彈性張量作為剛度的等效表征量,對(duì)剛度和降解速率不同的兩種材料在微結(jié)構(gòu)中的分布進(jìn)行優(yōu)化設(shè)計(jì).復(fù)合材料可降解微結(jié)構(gòu)剛度的優(yōu)化數(shù)學(xué)模型可以表示為

其中,X為設(shè)計(jì)變量的集合,i表示降解時(shí)間步,Λ為降解時(shí)間步的集合,v1表示單胞中材料1 的最終體積,φ為材料1 的體積分?jǐn)?shù)約束值,v0表示單胞的總體積.

骨折愈合的最終目標(biāo)是使斷骨恢復(fù)到健康骨的狀態(tài),Liu 等[30-31]研究表明健康松質(zhì)骨具有明顯的軸向優(yōu)勢(shì).對(duì)比分別使用彈性張量D11,D22,D33作為優(yōu)化目標(biāo)等效表征量時(shí)的優(yōu)化結(jié)果,當(dāng)以D11或D22為等效表征量時(shí),設(shè)計(jì)結(jié)果單向剛度過大,不利于骨痂的重塑;而以D33為等效表征量時(shí),優(yōu)化結(jié)果具有較均衡的抗拉壓特性和抗彎特性.為了能夠保證骨重塑的順利進(jìn)行,以D33作為等效表征量,同時(shí)為了保證結(jié)構(gòu)在降解初期保持足夠大的剛度,將在第1 和第60 天的D33之和最大作為優(yōu)化目標(biāo),可以表示為

2.2 靈敏度分析

在每一降解時(shí)間步對(duì)中間結(jié)構(gòu)進(jìn)行有限元分析,計(jì)算目標(biāo)函數(shù)及其靈敏度.目標(biāo)函數(shù)關(guān)于設(shè)計(jì)變量的靈敏度可以表示為

目標(biāo)函數(shù)C是關(guān)于彈性張量Drs的函數(shù),根據(jù)式(13)和式(18),Drs可以表示為關(guān)于x,y和t的函數(shù)

彈性張量Drs在ti時(shí)刻關(guān)于xe的導(dǎo)數(shù)可以表示為

通過鏈?zhǔn)角髮?dǎo)法則,式(23)可以表示為

根據(jù)式(10)和式(11),E(xe)和d(xe)關(guān)于xe的靈敏度能夠表示為

根據(jù)式(9),y關(guān)于xe的靈敏度可以表示為

材料1 的體積可以表示為

材料1 體積v1關(guān)于xe的靈敏度可以表示為

為克服棋盤格效應(yīng),采用基于Heaviside 靈敏度過濾方法使密度變量趨近于0 或1[32].

2.3 優(yōu)化設(shè)計(jì)流程

考慮時(shí)變剛度特性的復(fù)合材料微結(jié)構(gòu)拓?fù)鋬?yōu)化設(shè)計(jì)流程如圖5 所示,主要步驟如下.

圖5 考慮時(shí)變剛度特性的復(fù)合材料微結(jié)構(gòu)拓?fù)鋬?yōu)化流程Fig.5 Multi material microstructure topology optimization process considering time-changeable characteristics

步驟1:構(gòu)造微結(jié)構(gòu)的初始構(gòu)型,定義數(shù)學(xué)模型參數(shù),給定材料1 體積約束和過濾半徑,設(shè)定x和y的初始值,并定義可降解邊界等條件;

步驟2:根據(jù)式(9)進(jìn)行降解過程模擬,更新殘留率;

步驟3:根據(jù)式(16),在考慮材料降解的情況下對(duì)微結(jié)構(gòu)進(jìn)行均勻化分析;

步驟4:根據(jù)式(24)進(jìn)行靈敏度分析;

重復(fù)步驟2～4,直至達(dá)到規(guī)定的降解時(shí)間步;

步驟5:使用MMA 方法更新設(shè)計(jì)變量;

重復(fù)步驟2～5,直至滿足設(shè)計(jì)約束或者目標(biāo)函數(shù)收斂為止.當(dāng)滿足下列兩個(gè)條件之一時(shí),迭代終止:(1)連續(xù)兩次迭代的目標(biāo)函數(shù)變化小于收斂判據(jù);(2)迭代步數(shù)達(dá)到指定值,取100.

3 數(shù)值算例

本節(jié)通過典型數(shù)值算例對(duì)所提出拓?fù)鋬?yōu)化方法的有效性進(jìn)行驗(yàn)證.選定材料1 為可降解鐵基合金,具有高剛度但降解速率低,材料2 為可降解鎂基合金,降解速率高但剛度低,根據(jù)文獻(xiàn)[5],材料屬性如表1 所示.

表1 材料屬性[5]Table 1 Material properties[5]

可降解邊界的位置和數(shù)量對(duì)降解過程中結(jié)構(gòu)整體性能的變化具有重要影響,需要對(duì)具有不同可降解邊界的設(shè)計(jì)問題進(jìn)行研究.首先研究具有兩個(gè)可降解邊界的情況.如圖6 所示,建立幾何尺寸為L(zhǎng)×W=4 mm × 4 mm 的微結(jié)構(gòu)單胞,使用4 節(jié)點(diǎn)正方形平面應(yīng)變單元將其離散為單元數(shù)量為60 × 60 的有限元模型,過濾半徑為2,假定結(jié)構(gòu)的上側(cè)和下側(cè)為可降解邊界,左側(cè)和右側(cè)為不可降解邊界.

圖6 單胞幾何模型和兩個(gè)可降解邊界Fig.6 Unit cell model and two degradable interface

3.1 體積分?jǐn)?shù)約束對(duì)優(yōu)化結(jié)果的影響

將單胞初始構(gòu)型設(shè)置為均勻形態(tài),即將所有單元的初始相對(duì)密度設(shè)置為φ,φ分別取0.4,0.5,0.6,0.7.優(yōu)化目標(biāo)為第1 天和第60 天的D33之和最大.圖7 是不同φ時(shí)的優(yōu)化結(jié)果和3 × 3 結(jié)構(gòu),由圖可知,當(dāng)φ等于0.4,0.6 和0.7 時(shí)材料1 都形成了類似X 型的支撐結(jié)構(gòu),并隨著φ的增大更多的材料1 向X 區(qū)域集中,這是因?yàn)椴牧? 的剛度高,形成的支撐結(jié)構(gòu)能夠使優(yōu)化結(jié)果在同一φ情況下達(dá)到剛度最大的優(yōu)化目標(biāo).φ等于0.5 時(shí)材料1 形成了類似于正八邊形的支撐結(jié)構(gòu),同樣能夠能使優(yōu)化結(jié)果在同一φ情況下具有較好的剛度.從3 × 3 的宏觀構(gòu)型上看,不同體積分?jǐn)?shù)下設(shè)計(jì)結(jié)果的基本構(gòu)型類似,材料1 都是交叉分布.不同的優(yōu)化結(jié)果除了主要的支撐結(jié)構(gòu)之外,材料1 都在上下側(cè)有較多的分布,這是因?yàn)樯舷聜?cè)為初始的可降解邊界,而材料1 的降解速率低,在上下側(cè)分布較多的材料1 能夠保證在降解初期結(jié)構(gòu)剛度以較低速率減小,這與優(yōu)化目的相符.

圖7 不同φ的設(shè)計(jì)結(jié)果:(a)初始構(gòu)型,(b)優(yōu)化結(jié)果,(c) 3 × 3 結(jié)構(gòu)Fig.7 Optimal design with different φ:(a) initial configuration,(b) topology optimization result,(c) 3 × 3 structure

不同φ情況下優(yōu)化迭代歷程如圖8 所示,優(yōu)化目標(biāo)都達(dá)到了收斂狀態(tài),且φ對(duì)單胞優(yōu)化結(jié)果的D33具有較大的影響,因?yàn)棣帐菍?duì)材料1 的體積分?jǐn)?shù)上限進(jìn)行約束,而材料1 具有高剛度的材料特性,φ越大單胞優(yōu)化結(jié)果的剛度越大,因此可以根據(jù)不同的需求,通過調(diào)整φ的大小設(shè)計(jì)具有不同初始剛度的可降解植入物.圖9 為φ等于0.5 時(shí),單胞構(gòu)型和體積分?jǐn)?shù)隨迭代歷程,隨著迭代的進(jìn)行逐漸得到清晰的材料分布,體積分?jǐn)?shù)也保持穩(wěn)定.

圖8 不同φ的優(yōu)化迭代歷程Fig.8 Iterative process of optimization with different φ

圖9 φ為0.5 時(shí)單胞構(gòu)型和材料1 體積分?jǐn)?shù)隨迭代歷程變化趨勢(shì)Fig.9 The unit cell configuration and volume fraction vary with the iterative process when φ=0.5

不同φ情況下單胞優(yōu)化結(jié)果降解過程在第1～70 天D33的變化趨勢(shì)如圖10 所示,大致都表現(xiàn)為在1～60 天以相對(duì)較慢的速率降低,在60～70 天以相對(duì)較快的速率降低.不同優(yōu)化結(jié)果在第1～60 天和第60～70 天D33的降低量(第1～60 天D33的降低量為第1 天D33與第60 天D33的差值,其他時(shí)間段降低量同理) 如表2 中所示,不同優(yōu)化結(jié)果在第60～70 天D33天的降低量都達(dá)到了第1～60 天降低量的50%以上,即不同優(yōu)化結(jié)果的剛度都能在規(guī)定的天數(shù)內(nèi)以相對(duì)較低的速率減小,說明了所提出的拓?fù)鋬?yōu)化方法對(duì)結(jié)構(gòu)的時(shí)變剛度特性起到了調(diào)控作用.

圖10 不同φ的優(yōu)化結(jié)果降解過程中D33 變化趨勢(shì)Fig.10 Variation trend of D33 during degradation with different φoptimization results

表2 不同φ優(yōu)化結(jié)果降解過程中D33(GPa)在不同時(shí)間段的降低量Table 2 The variation of D33 (GPa) in different time periodsduring the degradation of different optimization results

3.2 與無時(shí)變剛度特性調(diào)控優(yōu)化結(jié)果對(duì)比

為了進(jìn)一步驗(yàn)證所提出設(shè)計(jì)方法的有效性,以φ等于0.5 的情況為例,與無時(shí)變剛度特性調(diào)控(以第1 天D33最大為優(yōu)化目標(biāo))的優(yōu)化結(jié)果、僅使用材料1 或2 的結(jié)構(gòu)在降解過程中結(jié)構(gòu)性能的變化進(jìn)行對(duì)比.圖11(a)和圖11(b)分別為有和無時(shí)變剛度特性調(diào)控優(yōu)化結(jié)果,后者材料1 形成了X 型支撐結(jié)構(gòu),但在上下側(cè)初始可降解邊界附近材料1 的分布較少,這是因?yàn)闊o時(shí)變剛度特性調(diào)控的設(shè)計(jì)僅以第1 天D33最大為優(yōu)化目標(biāo).從上述兩種優(yōu)化結(jié)果的結(jié)構(gòu)特征對(duì)比,體現(xiàn)出了所提出優(yōu)化方法對(duì)時(shí)變剛度特性調(diào)控的特點(diǎn).四種結(jié)構(gòu)在第1～301 天的降解歷程如圖12 所示,其中(a)列和(b)列分別為僅使用材料1 和僅使用材料2 結(jié)構(gòu)的降解歷程,材料的降解量與給定的材料降解速率一致,證明了降解模擬的準(zhǔn)確性,(c)列和(d)列分別為考慮時(shí)變剛度特性和無時(shí)變剛度特性調(diào)控單胞優(yōu)化結(jié)果的降解歷程,無時(shí)變剛度特性調(diào)控單胞優(yōu)化結(jié)果中由于在上下邊緣附近材料1 的分布較少,X 上下開口內(nèi)材料2 降解較快,導(dǎo)致結(jié)構(gòu)在降解初期的降解速率較快.

圖11 有/無時(shí)變剛度特性調(diào)控的優(yōu)化結(jié)果比較:(a) 考慮時(shí)變剛度特性,(b) 無時(shí)變剛度特性調(diào)控Fig.11 Optimization result:(a) considering time-changeable characteristics,(b) non-time-changeable characteristic regulation

圖12 不同結(jié)構(gòu)降解歷程:(a)僅使用材料1,(b)僅使用材料2,(c)考慮時(shí)變剛度特性,(d)無時(shí)變剛度特性調(diào)控Fig.12 Degradation of different structures:(a) using material 1 only,(b) using material 2 only,(c) considering time-changeable characteristics,(d) non-time-changeable characteristic regulation

4 種結(jié)構(gòu)降解過程在第1～70 天中D33的變化如圖13 所示,僅使用材料1 的結(jié)構(gòu)D33的初始值較大,且下降緩慢,僅使用材料2 的結(jié)構(gòu)D33的初始值較小且下降較快,與材料屬性差異相符.無時(shí)變剛度特性調(diào)控的單胞優(yōu)化結(jié)果在第1～60 天D33從34.24 GPa 下降到26.28 GPa,降低量為7.96 GPa,在第60～70 天D33從26.28 GPa 下降到23.07 GPa,降低量為3.21 GPa,而考慮時(shí)變剛度特性的單胞優(yōu)化結(jié)果在第1～60 天D33能夠以相對(duì)緩慢的速率降低,降低量為4.27 GPa,之后變化速率加快,第60～70 天D33的降低量為2.79 GPa,與優(yōu)化目的保持一致.通過上述對(duì)比進(jìn)一步驗(yàn)證了所提出的考慮時(shí)變剛度特性拓?fù)鋬?yōu)化設(shè)計(jì)方法的有效性.

圖13 不同結(jié)構(gòu)降解過程中D33 變化趨勢(shì)Fig.13 Variation trend of D33 during degradation with different structures

3.3 初始構(gòu)型對(duì)優(yōu)化結(jié)果的影響

由上述分析可知,材料1 的分布情況對(duì)結(jié)構(gòu)的時(shí)變特性有較大的影響,因此,研究單胞初始構(gòu)型中材料1 的分布對(duì)優(yōu)化結(jié)果的影響.分別以均勻形態(tài)(初始構(gòu)型1)、上下邊附近為材料1,其余為均勻形態(tài)(初始構(gòu)型2)、水平中線附近為材料1,其余為均勻形態(tài)(初始構(gòu)型3)為初始構(gòu)型,以第1 和第60 天D33之和最大為優(yōu)化目標(biāo),φ取0.5.初始構(gòu)型1,2,3 及對(duì)應(yīng)的拓?fù)鋬?yōu)化結(jié)果和3 × 3 結(jié)構(gòu)分別如圖14(a)～14(c)所示,由初始構(gòu)型1 和3 得到的拓?fù)鋬?yōu)化結(jié)果較為相似,形成了類似八邊形的材料1 分布,并且四角都有枝狀的材料1 分布,區(qū)別在于由材料2 形成的四邊形結(jié)構(gòu)尺寸不同,且后者左右側(cè)中部材料2 形成了局部分布.由初始構(gòu)型2 得到的拓?fù)鋬?yōu)化結(jié)果與上述兩者有較大區(qū)別,材料1 形成了類似菱形的支撐結(jié)構(gòu),并且在上下邊附近有較多的材料1 分布,這有利于在降解初期使結(jié)構(gòu)性能緩慢下降.

圖14 不同初始構(gòu)型優(yōu)化設(shè)計(jì):(a)初始構(gòu)型1,(b)初始構(gòu)型2,(c)初始構(gòu)型3Fig.14 Optimal design of different initial configurations:(a) initial 1,(b) initial 2,(c) initial 3

不同初始構(gòu)型得到的單胞優(yōu)化結(jié)果在第1～70天降解過程中D33的變化趨勢(shì)如圖15 所示,都表現(xiàn)為在第1～60 天以相對(duì)較慢的速率下降,60 天之后下降速率變大,與優(yōu)化目標(biāo)保持一致.其中由初始構(gòu)型2 得到的優(yōu)化設(shè)計(jì)結(jié)果D33的初始值雖略小,但在第1～60 天降低量最小,為3.40 GPa.因此,以上下側(cè)可降解邊界附近分布材料1,其他區(qū)域?yàn)榫鶆蛐螒B(tài)的初始構(gòu)型2 更有利于達(dá)到前期保持一定結(jié)構(gòu)剛度的目的.

圖15 不同初始構(gòu)型優(yōu)化結(jié)果降解過程中D33 變化趨勢(shì)Fig.15 Variation trend of D33 during degradation of optimization results of different initial configurations

3.4 宏觀結(jié)構(gòu)分析

上述案例僅分析了單胞結(jié)構(gòu)在降解過程中結(jié)構(gòu)性能的變化趨勢(shì),接下來由以初始構(gòu)型2 為單胞初始形態(tài),φ等于0.5,以第1 和第60 天D33之和最大為優(yōu)化目標(biāo)的單胞優(yōu)化結(jié)果構(gòu)成的宏觀結(jié)構(gòu)為例,通過與由無時(shí)變剛度特性調(diào)控優(yōu)化、僅使用材料1、僅使用材料2 的單胞組成的宏觀結(jié)構(gòu)對(duì)比,驗(yàn)證所提出設(shè)計(jì)方法對(duì)由單胞組成的宏觀結(jié)構(gòu)時(shí)變剛度特性的調(diào)控能力.如圖16 所示,由1 × 6 個(gè)單胞構(gòu)成幾何尺寸為L(zhǎng)×W=24 mm × 4 mm 的宏觀結(jié)構(gòu),左右兩側(cè)固支,中心點(diǎn)受到豎直向下的載荷,大小為500 N.

圖16 宏觀結(jié)構(gòu)幾何尺寸與邊界條件Fig.16 Geometric dimensions and boundary conditions of macrostructure

假設(shè)骨愈合周期為18 個(gè)月,由考慮時(shí)變剛度特性和無時(shí)變剛度特性調(diào)控單胞優(yōu)化結(jié)果組成的宏觀結(jié)構(gòu)在第1～541 天的降解歷程分別如圖17(a) 和圖17(b)所示,對(duì)比可知,考慮時(shí)變剛度特性設(shè)計(jì)的結(jié)構(gòu)在上下邊附近由材料1 形成的交叉結(jié)構(gòu),延緩了降解初期結(jié)構(gòu)的降解速率.4 種宏觀結(jié)構(gòu)在第1～541天降解過程中加載點(diǎn)在加載方向的位移變化如圖18所示,位移越大說明結(jié)構(gòu)剛度越小,在降解過程中僅使用材料1 結(jié)構(gòu)的剛度降低幅度很小,而僅使用材料2 結(jié)構(gòu)的剛度則快速降低,無時(shí)變剛度特性調(diào)控優(yōu)化結(jié)構(gòu)的剛度在降解過程的降解速率逐漸加快,且呈現(xiàn)出加速速率逐漸增大的趨勢(shì).考慮時(shí)變剛度特性優(yōu)化結(jié)構(gòu)的剛度則先以相對(duì)緩慢且穩(wěn)定的速率下降,然后加速下降.這說明所提出的設(shè)計(jì)方法達(dá)到了通過對(duì)單胞的優(yōu)化設(shè)計(jì),從而對(duì)由其組成的宏觀結(jié)構(gòu)的時(shí)變剛度特性進(jìn)行調(diào)控的目的.

圖17 有/無時(shí)變剛度特性調(diào)控設(shè)計(jì)的宏觀結(jié)構(gòu)降解歷程:(a) 考慮時(shí)變剛度特性,(b) 無時(shí)變剛度特性調(diào)控Fig.17 Degradation of different macro-structures:(a) considering timechangeable characteristics,(b) non-time-changeable characteristic regulation

圖17 有/無時(shí)變剛度特性調(diào)控設(shè)計(jì)的宏觀結(jié)構(gòu)降解歷程:(a) 考慮時(shí)變剛度特性,(b) 無時(shí)變剛度特性調(diào)控 (續(xù))Fig.17 Degradation of different macro-structures:(a) considering timechangeable characteristics,(b) non-time-changeable characteristic regulation (continued)

圖18 不同宏觀結(jié)構(gòu)降解過程中加載點(diǎn)在加載方向的位移變化趨勢(shì)Fig.18 Variation trend of displacement of loading point in loading direction in the degradation process of different macro-structures

3.5 可降解邊界條件對(duì)優(yōu)化結(jié)果影響

可降解邊界的位置和數(shù)量對(duì)結(jié)構(gòu)降解有重要的影響,為了分析可降解邊界條件對(duì)優(yōu)化結(jié)果的影響,并驗(yàn)證所提出的考慮時(shí)變剛度特性拓?fù)鋬?yōu)化方法對(duì)具有不同降解邊界問題的處理能力,僅將上側(cè)設(shè)置為可降解邊界,下側(cè)、左側(cè)、右側(cè)為不可降解邊界與具有上下側(cè)兩個(gè)可降解邊界的優(yōu)化結(jié)果進(jìn)行對(duì)比.以均勻形態(tài)為單胞初始構(gòu)型,優(yōu)化目標(biāo)為第1 和第60 天的D33之和最大,φ取0.5,具有上下側(cè)兩個(gè)可降解邊界和只有上側(cè)一個(gè)可降解邊界的單胞優(yōu)化結(jié)果和3 × 3 結(jié)構(gòu)分別如圖19(a)和圖19(b)所示.由于僅上側(cè)為可降解邊界,單胞優(yōu)化結(jié)果在上側(cè)分布了較多的材料1,憑借材料1 較慢的降解速率,能夠使結(jié)構(gòu)性能在較長(zhǎng)時(shí)間內(nèi)以相對(duì)較小的速率降低,雖然形成了較為復(fù)雜的單胞構(gòu)型,但仍具有良好的連接性.兩種優(yōu)化結(jié)果在第1～70 天降解過程中D33變化如圖20 所示,都呈現(xiàn)出在第1～60 天以相對(duì)緩慢的速率降低,在第60～70 天加速降低的趨勢(shì).D33初始值基本相等,但是僅有一個(gè)可降解邊界優(yōu)化結(jié)果D33在第1～60 天降低速率更為平緩,D33降低量為1.81 GPa.因此可以通過調(diào)整可降解邊界的數(shù)量和位置對(duì)結(jié)構(gòu)的時(shí)變剛度特性進(jìn)行有效的調(diào)控.

圖19 不同可降解邊界條件優(yōu)化結(jié)果Fig.19 Optimization results of different degradable interface condition

圖20 不同可降解邊界條件優(yōu)化結(jié)果降解過程中D33 變化趨勢(shì)Fig.20 Variation trend of D33 during degradation of optimization results with different degradable interface condition

4 結(jié)論

本文提出了一種考慮時(shí)變剛度特性的復(fù)合材料微結(jié)構(gòu)拓?fù)鋬?yōu)化方法.使用剛度和降解速率不同的兩種可降解材料,建立考慮時(shí)間維度材料降解的有限元模型,利用均勻化方法計(jì)算不同降解時(shí)間步中間結(jié)構(gòu)的力學(xué)性能,以特定降解時(shí)間步的結(jié)構(gòu)剛度之和最大為優(yōu)化目標(biāo),使用SIMP 方法對(duì)兩種材料的布局進(jìn)行優(yōu)化設(shè)計(jì),得到了最優(yōu)且輪廓清晰的材料分布.通過設(shè)計(jì)不同材料的分布,實(shí)現(xiàn)了對(duì)結(jié)構(gòu)時(shí)變剛度特性調(diào)控的目的.以典型算例對(duì)微結(jié)構(gòu)構(gòu)型和時(shí)變剛度特性進(jìn)行了研究,結(jié)果表明:

(1)相對(duì)于僅使用一種生物可降解材料和設(shè)計(jì)時(shí)無時(shí)變剛度特性調(diào)控的單胞優(yōu)化結(jié)果,考慮時(shí)變剛度特性單胞優(yōu)化結(jié)果在降解過程中D33在第1～60 天能夠以相對(duì)緩慢的速率降低,保證結(jié)構(gòu)在降解初期維持一定的剛度.

(2)材料1 的體積分?jǐn)?shù)約束φ、初始構(gòu)型及可降解邊界的位置和數(shù)量影響單胞拓?fù)浣Y(jié)構(gòu).隨著材料1 的體積分?jǐn)?shù)約束φ增大,材料1 更趨向于集中在X 型的支撐結(jié)構(gòu)上,微結(jié)構(gòu)的剛度增大;在其他條件相同的情況下,在設(shè)計(jì)域上下側(cè)可降解邊界附近分布材料1 的初始構(gòu)型得到了降解初期結(jié)構(gòu)剛度下降最小的優(yōu)化結(jié)果;僅有1 個(gè)降解邊界的設(shè)計(jì)結(jié)果,由于降解邊界少,降解緩慢,D33在第1～60 天降低速率更為平緩.

(3)考慮時(shí)變剛度特性優(yōu)化的單胞組成的宏觀結(jié)構(gòu),其剛度能夠在一定時(shí)間內(nèi)保持相對(duì)緩慢且穩(wěn)定的降低速率,而后降低速率逐漸增大,基本復(fù)合骨愈合所需的內(nèi)固定植入物的剛度變化特性需求.

本文提出的具有時(shí)變剛度調(diào)控的復(fù)合材料微結(jié)構(gòu)拓?fù)鋬?yōu)化設(shè)計(jì)方法,通過對(duì)兩種生物降解材料的布局優(yōu)化,可實(shí)現(xiàn)特定的結(jié)構(gòu)時(shí)變剛度特性.目前增材制造技術(shù)的發(fā)展能夠完成微結(jié)構(gòu)中的精細(xì)化結(jié)構(gòu)的加工,并且能夠達(dá)到足夠的加工精度和可接受的成本.后續(xù)將開展結(jié)合骨愈合理論的多尺度復(fù)合結(jié)構(gòu)的設(shè)計(jì),以達(dá)到特定的時(shí)變剛度特性,滿足骨折愈合不同階段對(duì)生物力學(xué)的需求,并進(jìn)行面向增材制造的可降解植入物拓?fù)鋬?yōu)化設(shè)計(jì),實(shí)現(xiàn)優(yōu)化設(shè)計(jì)-制造一體化研究,進(jìn)而取得更好的臨床效果.