考慮時變剛度特性的復合材料微結構拓撲優化設計方法1)

徐世鵬 丁曉紅 段朋云 張 橫

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

引言

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

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

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

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

1.1 具有特定時變剛度特性的復合材料結構

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

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

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

圖2 具有特定時變剛度特性的復合材料結構Fig.2 Periodic composite structure with specific time-changeable characteristics

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

式中,xmin和ymin分別是為避免計算中數值不穩定預設的最小值,取0.001,m為單元總數量.

1.2 結構降解過程模擬

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

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

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

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

式中,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.僅當單元處于可降解邊界時該單元能夠降解,為了保證數值計算的穩定,設定yb=0.1,假定當單元的殘留率小于yb時可降解邊界發生移動,如當yk-1,l＜yb時,認為單元e能夠從左側開始降解,當yk+1,l＜yb時,認為單元e能夠從右側開始降解,上下方向同理.根據相鄰單元的降解情況,單元e的降解速率為

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

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

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

其中H為

式中,γ為常數.

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

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

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

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

第2 步:根據yk-1,l,yk+1,l,yk,l-1和yk,l+1判斷相鄰單元的降解狀態;

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

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

重復第1 到第4 步,直到指定時間步.

1.3 材料插值模型

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

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

1.4 均勻化方法

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

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

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

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

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

式中,針對二維問題,應變ε=[ε1,ε2,ε3],ε1=[1,0,0]T,ε2=[0,1,0]T,ε3=[0,0,1]T.

當得到位移之后,利用均勻化方法,均勻彈性矩陣DH可以表示為

式中,|?|為單胞的總體積,為根據式(17)通過施加單位應變計算得到的節點位移場;Ue為根據式(12)通過施加全局單位應變計算得到的位移場

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

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

2 優化數學模型及靈敏度分析

2.1 優化數學模型

設定一個降解時間步長為一天,為了使復合結構剛度的變化趨勢接近于圖1 中實線所示的趨勢,即在骨痂生長期剛度下降值小,則可將在第1 和第60 降解步的復合結構剛度之和最大作為優化目標,以彈性張量作為剛度的等效表征量,對剛度和降解速率不同的兩種材料在微結構中的分布進行優化設計.復合材料可降解微結構剛度的優化數學模型可以表示為

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

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

2.2 靈敏度分析

在每一降解時間步對中間結構進行有限元分析,計算目標函數及其靈敏度.目標函數關于設計變量的靈敏度可以表示為

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

彈性張量Drs在ti時刻關于xe的導數可以表示為

通過鏈式求導法則,式(23)可以表示為

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

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

材料1 的體積可以表示為

材料1 體積v1關于xe的靈敏度可以表示為

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

2.3 優化設計流程

考慮時變剛度特性的復合材料微結構拓撲優化設計流程如圖5 所示,主要步驟如下.

圖5 考慮時變剛度特性的復合材料微結構拓撲優化流程Fig.5 Multi material microstructure topology optimization process considering time-changeable characteristics

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

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

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

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

重復步驟2～4,直至達到規定的降解時間步;

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

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

3 數值算例

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

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

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

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

3.1 體積分數約束對優化結果的影響

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

圖7 不同φ的設計結果:(a)初始構型,(b)優化結果,(c) 3 × 3 結構Fig.7 Optimal design with different φ:(a) initial configuration,(b) topology optimization result,(c) 3 × 3 structure

不同φ情況下優化迭代歷程如圖8 所示,優化目標都達到了收斂狀態,且φ對單胞優化結果的D33具有較大的影響,因為φ是對材料1 的體積分數上限進行約束,而材料1 具有高剛度的材料特性,φ越大單胞優化結果的剛度越大,因此可以根據不同的需求,通過調整φ的大小設計具有不同初始剛度的可降解植入物.圖9 為φ等于0.5 時,單胞構型和體積分數隨迭代歷程,隨著迭代的進行逐漸得到清晰的材料分布,體積分數也保持穩定.

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

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

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

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

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

3.2 與無時變剛度特性調控優化結果對比

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

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

圖12 不同結構降解歷程:(a)僅使用材料1,(b)僅使用材料2,(c)考慮時變剛度特性,(d)無時變剛度特性調控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 種結構降解過程在第1～70 天中D33的變化如圖13 所示,僅使用材料1 的結構D33的初始值較大,且下降緩慢,僅使用材料2 的結構D33的初始值較小且下降較快,與材料屬性差異相符.無時變剛度特性調控的單胞優化結果在第1～60 天D33從34.24 GPa 下降到26.28 GPa,降低量為7.96 GPa,在第60～70 天D33從26.28 GPa 下降到23.07 GPa,降低量為3.21 GPa,而考慮時變剛度特性的單胞優化結果在第1～60 天D33能夠以相對緩慢的速率降低,降低量為4.27 GPa,之后變化速率加快,第60～70 天D33的降低量為2.79 GPa,與優化目的保持一致.通過上述對比進一步驗證了所提出的考慮時變剛度特性拓撲優化設計方法的有效性.

圖13 不同結構降解過程中D33 變化趨勢Fig.13 Variation trend of D33 during degradation with different structures

3.3 初始構型對優化結果的影響

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

圖14 不同初始構型優化設計:(a)初始構型1,(b)初始構型2,(c)初始構型3Fig.14 Optimal design of different initial configurations:(a) initial 1,(b) initial 2,(c) initial 3

不同初始構型得到的單胞優化結果在第1～70天降解過程中D33的變化趨勢如圖15 所示,都表現為在第1～60 天以相對較慢的速率下降,60 天之后下降速率變大,與優化目標保持一致.其中由初始構型2 得到的優化設計結果D33的初始值雖略小,但在第1～60 天降低量最小,為3.40 GPa.因此,以上下側可降解邊界附近分布材料1,其他區域為均勻形態的初始構型2 更有利于達到前期保持一定結構剛度的目的.

圖15 不同初始構型優化結果降解過程中D33 變化趨勢Fig.15 Variation trend of D33 during degradation of optimization results of different initial configurations

3.4 宏觀結構分析

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

圖16 宏觀結構幾何尺寸與邊界條件Fig.16 Geometric dimensions and boundary conditions of macrostructure

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

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

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

圖18 不同宏觀結構降解過程中加載點在加載方向的位移變化趨勢Fig.18 Variation trend of displacement of loading point in loading direction in the degradation process of different macro-structures

3.5 可降解邊界條件對優化結果影響

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

圖19 不同可降解邊界條件優化結果Fig.19 Optimization results of different degradable interface condition

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

4 結論

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

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

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

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

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