人工骨支架結(jié)構(gòu)的單元設(shè)計(jì)及分析

楊 輝,齊建昌

(東南大學(xué)機(jī)械工程學(xué)院，江蘇 南京 211189)

在組織工程學(xué)領(lǐng)域，支架起著支撐及引導(dǎo)組織再生的功能，理想的支架要與真實(shí)骨的力學(xué)性能相匹配。三維周期性極小曲面(TPMS)[1-2]，由于其孔徑大小、孔的形狀、空隙率可控以及孔與孔之間相互連通等優(yōu)點(diǎn)，為設(shè)計(jì)者提供了較大的自由度，故用TPMS設(shè)計(jì)支架引起學(xué)者的廣泛關(guān)注。

文獻(xiàn)[3]用TPMS單元庫(kù)中的IWP分別于P、G融合，生成結(jié)構(gòu)上更加復(fù)雜的支架，并在CATIA CAE模塊模擬支架的受力情況，分析支架的應(yīng)力分布；文獻(xiàn)[4]使用薄板樣條徑向基函數(shù)設(shè)計(jì)非均質(zhì)化支架，計(jì)算G單元在不同空隙率下的楊氏模量。

文獻(xiàn)[5]在Rhino中設(shè)計(jì)12種多面體單元，并在ABAQUS中分析單元在彈性形變范圍內(nèi)的受力情況，比較不同單元的楊氏模量；文獻(xiàn)[6]使用模擬和實(shí)驗(yàn)的方法分析TPMS單元的塑性變形過(guò)程，重點(diǎn)研究了P、D單元在體積分?jǐn)?shù)為30%和60%時(shí)的應(yīng)力應(yīng)變曲線；文獻(xiàn)[7]在ANSYS中分析Kelvin單元的彈性變形過(guò)程，主要研究了單元的體積分?jǐn)?shù)與楊氏模量之間的關(guān)系，并用實(shí)驗(yàn)的方法驗(yàn)證了模擬結(jié)果。

文獻(xiàn)[8]在ABAQUS中分析G單元的彈性變形過(guò)程，主要研究了空隙率為55%和70%兩個(gè)模型，并在Fluent中模擬了支架內(nèi)部的流動(dòng)性；文獻(xiàn)[9]主要研究P、IWP單元的空隙率與曲面厚度以及曲面半徑之間的關(guān)系，并在ABAQUS中分析曲面厚度、曲面半徑取不同值時(shí)，單元的楊氏模量的變化規(guī)律。

本文在分析TPMS單元的基礎(chǔ)上，根據(jù)距離場(chǎng)設(shè)計(jì)保持其拓?fù)溥B通的人工單元，并在ABAQUS中分析TPMS偏置單元和人工單元的楊氏模量。

1 三維周期性極小曲面

1.1 TPMS的三角函數(shù)表示形式

TPMS中的P、D、G、IWP 4種經(jīng)典類型曲面的三角函數(shù)表示形式為[1]：

(1)

式中：X=2πx，Y=2πy，Z=2πz，其中x∈R，y∈R，z∈R。規(guī)定φ<0表示隱式曲面的內(nèi)部，φ>0表示隱式曲面的外部，則式(1)中的φG,φIWP的隱式方程需要反向。

1.2 單元的體積分?jǐn)?shù)與空隙率

TPMS曲面在三維空間中是復(fù)雜幾何及拓?fù)溥B通的結(jié)構(gòu)，其將三維空間分割為內(nèi)外兩部分，并且內(nèi)外空間的體積相等。將曲面裁剪成實(shí)體可以更好地觀察其幾何拓?fù)湫螒B(tài)，比如曲面P(φP=0)對(duì)所在空間進(jìn)行分割，其內(nèi)部空間為桿系結(jié)構(gòu)(φP≤0)，外部空間成復(fù)雜的孔洞結(jié)構(gòu)(φP≥0)。

定義1：給定TPMS曲面φ，TPMS實(shí)體單元定義為{φ(x,y,z)≤0|?x,y,z∈R,0≤x≤1,0≤y≤1,0≤z≤1}。

圖1展示了TPMS實(shí)體。設(shè)TPMS實(shí)體單元內(nèi)部所圍成的體積為V，其最小包圍盒體積為VBox，則單元空隙率定義為：

(2)

圖1 TPMS實(shí)體單元

2 保持拓?fù)溥B通的人工單元設(shè)計(jì)

2.1 TPMS單元問(wèn)題描述

式(1)所描述的TPMS曲面可進(jìn)行內(nèi)外偏置(φ=t)，常量t為負(fù)數(shù)時(shí)進(jìn)行內(nèi)偏置，常量t為正數(shù)時(shí)進(jìn)行外偏置，TPMS曲面偏置后通過(guò)對(duì)空間的分割得到減薄或加厚的TPMS實(shí)體。規(guī)定當(dāng)t≠0為偏置單元，用Po,Do表示。

當(dāng)TPMS曲面進(jìn)行內(nèi)偏置時(shí)，隨著偏置量的增大，曲面向中心線方向壓縮的過(guò)程中會(huì)出現(xiàn)自相交，對(duì)應(yīng)的實(shí)體會(huì)被截?cái)啵瑥亩鵁o(wú)法保證連通性，很難獲得高空隙率 (p>90%) 的支架。

如圖2(a)～(c)所示，P單元常量t為0時(shí)，空隙率為50%，桿系交匯處凸起較明顯，在制作時(shí)會(huì)浪費(fèi)材料；當(dāng)常量t為-0.8時(shí)，空隙率為73%，桿系半徑明顯減?。划?dāng)常量t為-1.0時(shí)，桿系出現(xiàn)自相交現(xiàn)象，無(wú)法保證連通性。圖2(d)～(f)展示了D單元隨著空隙率的增大，出現(xiàn)了桿系斷裂的現(xiàn)象。

2.2 TPMS單元的拓?fù)浞治?/h3>

TPMS原始單元(t=0)很難觀察其拓?fù)溥B通性，為了方便觀察，將原始單元進(jìn)行內(nèi)偏置得到減薄單元。如圖3所示，各個(gè)桿系用圓柱表示，記錄圓柱中心線兩端點(diǎn)坐標(biāo)；桿系匯交處用小球表示，并記錄球心的坐標(biāo)。球心坐標(biāo)、圓柱中心線端點(diǎn)坐標(biāo)以及坐標(biāo)間的連通關(guān)系用于設(shè)計(jì)人工單元。

圖2 空隙率逐漸增大出現(xiàn)的桿系自相交和斷裂

圖3 TPMS單元的拓?fù)溥B通圖

由圖3中的球和圓柱的連通規(guī)律可知，P單元為六連通結(jié)構(gòu)，D單元為四連通結(jié)構(gòu)，G單元為三連通結(jié)構(gòu)，IWP單元為八連通結(jié)構(gòu)。

2.3 人工單元設(shè)計(jì)

TPMS的拓?fù)涑Ｒ?jiàn)于真實(shí)骨的結(jié)構(gòu)中，市場(chǎng)迫切需要一種既保持TPMS單元的拓?fù)溥B通性，又能得到高空隙率的支架。2.2節(jié)具體分析了單元的拓?fù)溥B通性，并記錄了關(guān)鍵數(shù)據(jù)點(diǎn)的坐標(biāo)，借助這些信息就能夠得到TPMS桿系結(jié)構(gòu)的中心線分布，再以中心線為軸，求出光滑過(guò)渡的等距面，就能得到滿足設(shè)計(jì)需求的人工單元。

由距離場(chǎng)設(shè)計(jì)人工單元的步驟如下：

1)在空間中等距離布點(diǎn)，并記錄點(diǎn)的坐標(biāo)。

2)由TPMS的拓?fù)溥B通圖得到其中心線分布。

3)計(jì)算點(diǎn)到各中心線的距離，記為di。

4)由距離函數(shù)F(d)求得各點(diǎn)處的總場(chǎng)強(qiáng)。

(3)

5)MarchingCubes算法顯示曲面。

圖4給出了由TPMS單元設(shè)計(jì)人工單元的主要過(guò)程，用Pt表示人工單元。

圖4 由TPMS單元設(shè)計(jì)人工單元的主要過(guò)程

3 單元分析

將由P單元組成的四面體網(wǎng)格模型導(dǎo)入到ABAQUS CAE中，圖5展示了有限元分析過(guò)程，首先給模型賦予材料屬性[10]，彈性模量E=2GPa，泊松比為0.3；然后設(shè)置邊界條件，限制底面的6個(gè)自由度，頂面位移為壓縮方向的1%，提交分析直至達(dá)到平衡狀態(tài)，求得此時(shí)網(wǎng)格頂部節(jié)點(diǎn)處的反作用力F。設(shè)網(wǎng)格模型的最小包圍盒橫截面面積為A，壓縮方向的初始長(zhǎng)度為L(zhǎng)，壓縮量為ΔL，應(yīng)變用ε表示，則有效楊氏模量Ef定義為：

(4)

其中應(yīng)變?chǔ)哦x為：

(5)

圖5 有限元分析的主要過(guò)程

3.1 Ef與單元陣列的關(guān)系

選擇單元尺寸為1mm、空隙率為75%的P單元進(jìn)行分析，將P的偏置單元(Po)和人工單元(Pt)分別裝配為1，2×2×2，3×3×3，4×4×4，5×5×5的結(jié)構(gòu)。將上述四面體網(wǎng)格模型導(dǎo)入ABAQUS CAE中，底面固定，頂面壓縮，壓縮量為Y方向長(zhǎng)度的1%，記錄平衡狀態(tài)下頂部節(jié)點(diǎn)處的反作用力F，根據(jù)式(4)求得該結(jié)構(gòu)的有效楊氏模量Ef。圖6展示了空隙率為75%時(shí)，P單元的有效楊氏模量(單位為MPa)與單元陣列的關(guān)系。

圖6 Ef與單元陣列的關(guān)系

由圖6可知，當(dāng)空隙率為75%時(shí)，P的人工單元的有效楊氏模量大于偏置單元?？赡艿脑蚴牵谠摽障堵氏?，偏置單元在桿系匯交處凸起較明顯，單元拼接處的桿系半徑過(guò)小，使結(jié)構(gòu)出現(xiàn)薄弱環(huán)節(jié)，經(jīng)優(yōu)化后的人工單元，材料繞桿系中心線分布更加均勻，故抵抗彈性變形能力更強(qiáng)；隨著單元陣列數(shù)目的變化，P的偏置單元和人工單元的有效楊氏模量出現(xiàn)微小的波動(dòng)，若考慮造型過(guò)程中空隙率的微小偏差以及網(wǎng)格質(zhì)量差異對(duì)分析結(jié)果造成的影響，可近似認(rèn)為，單元的有效楊氏模量為一定值。根據(jù)上述分析，可得出如下結(jié)論：1)當(dāng)空隙率為75%時(shí)，P的人工單元抵抗彈性變形能力優(yōu)于偏置單元；2)當(dāng)空隙率為75%時(shí)，P的偏置單元和人工單元的有效楊氏模量均與單元陣列數(shù)目無(wú)關(guān)。

3.2 Ef與空隙率的關(guān)系

通過(guò)曲面偏置可得到不同空隙率單元，取單元陣列為2×2×2，單元尺寸為1mm，空隙率為30%～75%(步長(zhǎng)為5%)共10個(gè)模型進(jìn)行分析，將模型導(dǎo)入ABAQUS CAE中，底面固定，頂面壓縮，壓縮量為0.02，記錄平衡狀態(tài)下頂面節(jié)點(diǎn)處的反作用力F，根據(jù)式(4)求得結(jié)構(gòu)的有效楊氏模量Ef。圖7展示了單元陣列為2×2×2時(shí)，P單元的有效楊氏模量(單位為MPa)與空隙率的關(guān)系。

由圖7可知，隨著空隙率的增大，P的偏置單元(Po)和人工單元(Pt)的有效楊氏模量減小，且偏置單元的下降速度更快。當(dāng)空隙率小于55%時(shí)，P的偏置單元的有效楊氏模量大于人工單元；當(dāng)空隙率大于55%時(shí)，P的偏置單元的有效楊氏模量小于人工單元，出現(xiàn)此現(xiàn)象的原因在于，隨著空隙率的增大，P的偏置單元桿系半徑減小速度快于人工單元，而單元拼接處是該結(jié)構(gòu)最薄弱的位置。根據(jù)以上數(shù)據(jù)，可得出如下結(jié)論：1)隨著空隙率的增大，P的偏置單元和人工單元的有效楊氏模量均減??；2)隨著空隙率的增大，P的偏置單元的有效楊氏模量減小的速度大于人工單元。

圖7 Ef與空隙率的關(guān)系

4 結(jié)束語(yǔ)

TPMS單元作為支架設(shè)計(jì)的造型單元有其獨(dú)特的優(yōu)勢(shì)，比如空隙率可控，拓?fù)溥B通規(guī)律與骨小梁類似等，但是很難得到高空隙率的支架。而本文提出的應(yīng)用距離場(chǎng)設(shè)計(jì)保持TPMS拓?fù)溥B通的人工單元的方法，能夠得到更高空隙率的支架。由于條件限制，文中僅做了理論分析和有限元模擬，缺乏試驗(yàn)數(shù)據(jù)做支撐，因而后續(xù)工作將盡可能與相關(guān)企業(yè)及其他高校合作，將試驗(yàn)數(shù)據(jù)補(bǔ)充完整。

[1] SCHOEN A H. Infinite periodic minimal surfaces without self- intersection[M]. Cambridge: National Aeronautics and Space Administration,1970.

[2] JUNG Y, CHU K T, TORQUATO S. A variational level set approach for surface area minimization of triply-periodic surfaces[J]. Journal of Computational Physics,2007,223(2):711-730.

[3] YOO D J, KIM K H. An advanced multi-morphology porous scaffold design method using volumetric distance field and beta growth function[J]．International Journal of Precision Engineering and Manufacturing,2015,16(9):2021-2032.

[4] YOO D. Heterogeneous minimal surface porous scaffold design using the distance field and radial basis functions[J]．Medical Engineering & Physics,2012,34(5):625-639.

[5] WETTERGREEN M A,BUCKLEN B S,STARLY B,et al. Creation of a unit block library of architectures for use in assembled scaffold engineering[J]．Computer-Aided Design,2005,37(11): 1141-1149.

[6] KADKHODAPOUR J, MONTAZERIAN H, RAEISI S. Investigating internal architecture effect in plastic deformation and failure for TPMS-based scaffolds using simulation methods and experimental procedure[J]．Materials Science & Engineering C,2014,43(43):587-597.

[7] ZARGARIAN A, ESFAHANIAN M, KADKHODAPOUR J, et al. Effect of solid distribution on elastic properties of open-cell cellular solids using numerical and experimental methods[J]．Journal of the Mechanical Behavior of Biomedical Materials,2014,37(37C):264-273.

[8] OLIVARES A L, ELIA M, PLANELL J A, et al. Finite element study of scaffold architecture design and culture conditions for tissue engineering[J]．Biomaterials,2009,30(30):6142-6149.

[9] ALMEIDA H A, BARTOLO P J. Design of tissue engineering scaffolds based on hyperbolic surfaces:structural numerical evaluation[J]．Medical Engineering & Physics,2014,36(8):1033-1040.

[10] SUN W, DARLING A, STARLY B, et al. Computer-aided tissue engineering: overview, scope and challenges[J]．Biotechnology & Applied Biochemistry,2004,39(1):29-47.