徐世鵬 丁曉紅 段朋云 張 橫

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

引言

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

式中,i為降解時間步,de為單元e的降解速率,t為單位降解時間步的時長.

如圖4 所示,使用k和l將單元e的位置標記為(k,l),與單元e相鄰的4 個單元為(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ù)值計算的穩(wěn)定,設(shè)定yb=0.1,假定當(dāng)單元的殘留率小于yb時可降解邊界發(fā)生移動,如當(dāng)yk-1,l＜yb時,認為單元e能夠從左側(cè)開始降解,當(dāng)yk+1,l＜yb時,認為單元e能夠從右側(cè)開始降解,上下方向同理.根據(jù)相鄰單元的降解情況,單元e的降解速率為

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

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

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

其中H為

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

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

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

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

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

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

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

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

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

1.3 材料插值模型

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

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

1.4 均勻化方法

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

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

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

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

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

式中,針對二維問題,應(yīng)變ε=[ε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)變計算得到的節(jié)點位移場;Ue為根據(jù)式(12)通過施加全局單位應(yī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è)定一個降解時間步長為一天,為了使復(fù)合結(jié)構(gòu)剛度的變化趨勢接近于圖1 中實線所示的趨勢,即在骨痂生長期剛度下降值小,則可將在第1 和第60 降解步的復(fù)合結(jié)構(gòu)剛度之和最大作為優(yōu)化目標,以彈性張量作為剛度的等效表征量,對剛度和降解速率不同的兩種材料在微結(jié)構(gòu)中的分布進行優(yōu)化設(shè)計.復(fù)合材料可降解微結(jié)構(gòu)剛度的優(yōu)化數(shù)學(xué)模型可以表示為

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

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

2.2 靈敏度分析

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

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

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

通過鏈式求導(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è)計流程

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

圖5 考慮時變剛度特性的復(fù)合材料微結(jié)構(gòu)拓撲優(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)進行降解過程模擬,更新殘留率;

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

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

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

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

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

3 數(shù)值算例

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

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

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

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

3.1 體積分數(shù)約束對優(yōu)化結(jié)果的影響

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

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

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

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

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

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

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

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

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

圖11 有/無時變剛度特性調(diào)控的優(yōu)化結(jié)果比較:(a) 考慮時變剛度特性,(b) 無時變剛度特性調(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)考慮時變剛度特性,(d)無時變剛度特性調(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的初始值較小且下降較快,與材料屬性差異相符.無時變剛度特性調(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,而考慮時變剛度特性的單胞優(yōu)化結(jié)果在第1～60 天D33能夠以相對緩慢的速率降低,降低量為4.27 GPa,之后變化速率加快,第60～70 天D33的降低量為2.79 GPa,與優(yōu)化目的保持一致.通過上述對比進一步驗證了所提出的考慮時變剛度特性拓撲優(yōu)化設(shè)計方法的有效性.

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

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

由上述分析可知,材料1 的分布情況對結(jié)構(gòu)的時變特性有較大的影響,因此,研究單胞初始構(gòu)型中材料1 的分布對優(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)化目標,φ取0.5.初始構(gòu)型1,2,3 及對應(yīng)的拓撲優(yōu)化結(jié)果和3 × 3 結(jié)構(gòu)分別如圖14(a)～14(c)所示,由初始構(gòu)型1 和3 得到的拓撲優(yōu)化結(jié)果較為相似,形成了類似八邊形的材料1 分布,并且四角都有枝狀的材料1 分布,區(qū)別在于由材料2 形成的四邊形結(jié)構(gòu)尺寸不同,且后者左右側(cè)中部材料2 形成了局部分布.由初始構(gòu)型2 得到的拓撲優(yōu)化結(jié)果與上述兩者有較大區(qū)別,材料1 形成了類似菱形的支撐結(jié)構(gòu),并且在上下邊附近有較多的材料1 分布,這有利于在降解初期使結(jié)構(gòu)性能緩慢下降.

圖14 不同初始構(gòu)型優(yōu)化設(shè)計:(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的變化趨勢如圖15 所示,都表現(xiàn)為在第1～60 天以相對較慢的速率下降,60 天之后下降速率變大,與優(yōu)化目標保持一致.其中由初始構(gòu)型2 得到的優(yōu)化設(shè)計結(jié)果D33的初始值雖略小,但在第1～60 天降低量最小,為3.40 GPa.因此,以上下側(cè)可降解邊界附近分布材料1,其他區(qū)域為均勻形態(tài)的初始構(gòu)型2 更有利于達到前期保持一定結(jié)構(gòu)剛度的目的.

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

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

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

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

圖17 有/無時變剛度特性調(diào)控設(shè)計的宏觀結(jié)構(gòu)降解歷程:(a) 考慮時變剛度特性,(b) 無時變剛度特性調(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)降解過程中加載點在加載方向的位移變化趨勢Fig.18 Variation trend of displacement of loading point in loading direction in the degradation process of different macro-structures

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

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

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

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

4 結(jié)論

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

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

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

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

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