選區激光熔化Voronoi 多孔結構的設計與性能預測

曾壽金，吳啟銳，何家辰，韋鐵平，葉建華

（1. 福建工程學院機械與汽車工程學院，福建福州350118；2. 蘇州大學附屬第一醫院骨科研究所，江蘇蘇州215006）

1 引 言

對損傷骨組織進行植入體的置換，一直是醫學上常用的方式。臨床上常采用力學性能較高的金屬植入體進行替換，例如不銹鋼、鈦和鈦合金、鈷鉻合金等金屬植入體［1］，用作人工膝關節、牙種植體、人工骨等［2-3］。但是人工植入體存在著遠期無菌松動和再骨折現象，研究人員普遍認為“應力屏蔽”效應［4-5］是其主要原因?！皯ζ帘巍毙侵附饘偃斯ぶ踩塍w的彈性模量與人體骨頭的彈性模量存在著巨大差異，導致應力無法從人工植入體傳遞到骨組織上。由于人體骨細胞有著應力響應的生物學特性，人工植入體一端彈性模量較大，會承受較多的應力，骨細胞得不到足夠的應力刺激時，骨細胞會發生凋亡，從而對健康骨頭產生不良影響［6-8］。

人骨小梁的孔隙率在50%～90%之間，彈性模量為0.5～20 GPa［9］，而實心316L 不銹鋼的彈性模量為210 GPa，將實心金屬設計成多孔結構能夠有效地降低金屬材料的彈性模量。因此，華南理工大學楊永強團隊［10］研究了正八面體和正六面體單元多孔結構的比表面積、孔隙率、平均孔徑對抗壓強度和彈性模量的影響。重慶大學柏龍［11］等人基于體心立方（Body-Centered Cubic，BCC）多孔結構設計出一種兼具輕質和高強性能的體心四方（Body-Centered Tetragonal，BCT）多孔結構，發現BCT 多孔結構比BCC 在力學上有著更加明顯的優勢，為多孔結構在工程領域的應用提供了力學理論依據。Simoneau［12］等人提出了一種用于醫學應用的隨機開放的多孔結構設計、制造和測試的新方法，測定了剛度和抗剪強度均符合數值預期。Yan［13］等人在研究三周期極小曲面和金剛石單元多孔結構時發現三周期極小曲面多孔結構的孔隙率與小梁骨孔隙度相當，彈性模量在0.12～1.25 GPa 內，且可根據小梁骨模量進行調整，從而減少或避免“應力屏蔽”現象，增加植入物的壽命。

選區激光熔化（Selective Laser Melting，SLM）技術［14-16］是成型金屬多孔結構重要的增材制造技術。它通過對切片后的模型數據進行路徑掃描，利用激光熔化金屬粉末，層層累積，可以成型出具有極其復雜幾何特征的零件，可以很好地滿足醫療領域個性化植入體的需求［17-19］。

天然骨組織是由外部皮質骨和內部骨小梁組成的，它們并不是密實化組織，而是復雜的無規則孔狀結構，這些孔狀結構的形態各異、性能不同［20］。如何在植入體設計階段就能考慮成型后制件的力學性能與生物相容性要求，是利用選區激光熔化技術服務于骨科醫療需求的重要課題。為此，本文通過控制不規則多孔結構的設計參數，研究了一種自下而上的多孔結構參數化設計方法，建立了利用不規則多孔結構的結構設計參數預測彈性模量、抗壓強度和孔隙率等性能參數的多元非線性數學模型，為醫用植入體的研究提供了設計依據。

2 試驗方法及多孔結構設計

2.1 試驗材料

選用316L 不銹鋼粉末，其粒徑分布為10～45 μm，對應的歐洲牌號為1.4404。利用掃描電子顯微鏡對粉末進行觀測，316L 不銹鋼粉末有很高的球形，如圖1 所示，各元素含量如表1 所示。實驗前粉末在80 ℃真空干燥10 h 以上。

表1 316L 不銹鋼粉末的化學成分Tab.1 Chemical composition of 316L stainless steel powder（wt.%）

圖1 316L 不銹鋼粉末形貌Fig. 1 Morphology of 316L stainless steel powder

2.2 試驗設備

圖2 為SLM 工作原理。試驗采用的是德國SLM Solution 公司研發的SLM-125HL 成形設備，其中包括激光系統、掃描振鏡系統、冷卻系統、密封成形缸、鋪粉機構以及控制系統。設備的成形缸尺寸為125 mm×125 mm×125 mm，最大功率為400 W，成形精度可達0.01 mm。激光熔化過程中，采用氮氣對材料進行保護。成形工藝參數為：激光功率P=250 W、掃描速度v=800 mm/s、掃描間距h=0.08 mm、鋪粉厚度d=0.03 mm，基板預熱為100 ℃，掃描策略為單向掃描。

圖2 選區激光熔化原理Fig.2 Principle of selective laser melting

2.3 不規則多孔結構設計

圖3 概率球模型Fig.3 Probability ball model

Voronoi 是幾何計算中一種常見的方法。它基于在空間或平面上隨機生成的種子點［21］，通過特定的算法將種子相互連接，形成類似于細胞的多邊形，從而實現空間或平面的分割。對于不規則性的定義目前尚未統一，本文從種子在Voronoi 空間區域的分布出發，探討了孔隙結構的不規則性。規則概率球的球面中心與不規則點的關系如圖3（a）所示，其中R為概率球的半徑，在概率球內部及其表面會按照函數關系（表2）隨機生成不規則點。Ki為相鄰概率球之間的距離（即單元距離），Li是概率球心與隨機點之間的距離（圖3（b））。I被定義為不規則度，則有：

其中N為種子點數。本研究的不規則多孔結構的空間點數均為7×7×7，共計343。

本文利用參數化建模軟件Grasshopper 對不規則多孔結構進行設計。其過程為：（1）生成規則點陣，如圖4（a）所示；（2）通過調整概率球的直徑來改變種子的不規則度，概率球直徑越大，隨機點運動的范圍越廣，種子的分散程度就越高，如圖4（b）和4（c）所示；（3）通過Voronoi 幾何算法生成完全互聯的不規則多孔結構骨架，如圖4（d）和4（e）所示，并對骨架進行包邊，得到的模型如圖4（f）所示。

表2 概率球隨機函數Tab.2 Probability sphere random functions

圖4 不規則多孔結構設計過程Fig.4 Design process of irregular porous structure

2.4 力學性能測試

SLM 制備樣件完成后如圖5（a）所示，樣品被線切割、清洗、和干燥。利用掃描電子顯微鏡觀察顯微結構，如圖5（b）～5（c）所示，多孔結構表面會有少量未熔化的粉末和掛渣，這是SLM技術成形造成的，樣件表面沒有裂紋、缺損，成型效果較好。采用中國長春機械院研發的SDS100型電液伺服疲勞試驗機測得多孔結構的彈性模量和抗壓強度，壓縮試驗滿足ISO-13314∶2011標準，壓縮速度為1 mm/min。平臺應力（σpl）被確定為20%～40%壓縮應變的算數平均值，可表示為抗壓強度。彈性模量被稱為多孔金屬材料的表觀彈性模量，由20%～70% 平臺應力（σpl）之間連接直線的斜率確定，如圖6 所示。

圖5 SLM 制備的316L 不規則多孔結構樣件及其SEM 圖像Fig.5 As-built 316L irregular porous structure by SLM and its SEM images

圖6 多孔結構的應力-應變曲線Fig.6 Stress-strain curve to determine characteristic values from compression testing of porous structure

2.5 響應面試驗

響應面法（Response Surface Method，RSM）是一種試驗綜合設計的優化方法，最早于1951 年由數學家George E.P.Box 和K.B.Wilson 提出。響應面法通過較少的試驗次數，建立了輸入變量與輸出變量之間的函數關系，通過分析多元回歸方程尋求最優結果。RSM 包括中心復合表面設計、中心復合有界設計、中心復合序貫設計以及BOX-Behnken 設 計（BOX-Behnken Design，BBD）。

圖7 三因素BBD 試驗中試驗點分布示意圖Fig.7 Distribution diagram of test points in three-factor BBD experimental

BBD 是由2n個設計與不完全區域組合設計組合而成，相同水平時會比中心復合設計擁有更少的試驗次數，其設計不包含立方體的頂點，即各變量的極值點，三因素試驗中各點分布如圖7所示。BBD 是球形設計，具有可旋轉或近似可旋轉性，沒有序貫性，所有因素不會被同時安排為高水平的試驗組合，此設計比較適用于對安全要求較高或有特別需求的試驗。

本文選擇了RSM 的BBD，試驗參數如表3所示，因素分別為孔棱直徑（D）、不規則度（I）和單元距離（K）。表中，“-1”代表低水平，“0”代表中心點，“+1”代表高水平。

表3 結構參數的水平編碼及真實值Tab.3 Horizontal encoding and true values of structural parameters

3 實驗結果及數據處理

3.1 實驗結果

SLM 成功制備樣件后，測得相關參數如表4所示。其中，E表示彈性模量，σ表示為抗壓強度，Φ表示為孔隙率。

3.2 實驗數據處理

3.2.1 數據歸一化處理

對彈性模量、抗壓強度和孔隙率3 個響應值進行歸一化處理，能夠避免因單位的不同對結果的影響。在多孔結構制備實驗的過程中，彈性模量具有望小的特性，而抗壓強度和孔隙率具有望大的特性?？箟簭姸群涂紫堵实臍w一化計算公式為：

彈性模量的歸一化計算公式為：

式中：n為響應面的目標數；i為實驗組數；xi(n)為參考序列(n)為比較序列。

表4 BBD 試驗結果Tab.4 BBD experimental results

3.2.2 灰色關聯系數計算

對于參考序列和比較序列之間的灰色關聯系數（Grey Relation Coefficient，GRC）計算為：

式中：Δ(n)被成為偏差數列，表示參考數列與比較數列的差值，即Δ(n)=|xi(n)-(n)|；Δmin和Δmax分別表示偏差數列的最小值和最大值；ξ為分別系數，通常取值為0.5。

3.2.3 主成分分析及響應值權重

在進行多目標分析的過程中，確定各響應值的權重是需要實驗數據進行降維處理。主成分分析法能夠從實驗數據中得到各因子對響應值的影響程度，是一種十分有效的數據降維處理，從而確定響應值的目標權重。

（1）構造響應值樣本矩陣：

式中：t為響應值目標數，t=3；m為實驗組數，m=17。

（2）計算相關關系系數矩陣Rjk，Rjk=[rjk]，

式中：rjk為關系系數矩陣中的各元素；Cov[xi(j)，xi(k)]為xi(j)和xi(k)的協方差；σ[xi(j)]和σ[xi(k)]分別是xi(j)和xi(k)的標準差。

（3）特征值求解得到：

式中：E為單位矩陣；λk為關系系數矩陣Rjk的特征值，特征值根據大小排列，即λ1＞λ2＞…λk＞0，k=1，2，…，n。

（4）計算主成分貢獻率ak：

式中ak為各響應值的權重。對各響應進行主成分分析，結果如表5 所示。

表5 主成分分析結果Tab.5 Results of principal component analysis

3.2.4 灰色關聯度計算

表6 BBD 實驗數據處理結果Tab.6 Result of BBD experimental data processing

3.3 灰色關聯度預測模型

本文采用多元回歸分析函數（公式（9））建立和分析輸入參數和輸出結果的模型：

其中：y為響應值，β0為截距因子，βj，βij，βjj分別是線性項、相互作用項和二次項的系數，xi，xj表示處理參數，k為因子個數，ε為殘差。通過響應面設計，建立一個基于試驗結果的不規則多孔結構參數與GRG 之間的統計預測模型。利用Design-Expert（V8.0.6.1）軟件得到GRG 的預測模型：

其中：線性效應項D為孔棱直徑，I為不規則度，K為單元距離，D×I，D×K，I×K是相互作用項，D2，I2，K2是二次項。在進行有效性篩選時，P＞F值小于0.05 被評價為顯著項。GRG 模型中各項的F值和P＞F值如表7 所示。計算得到R2為0.975 9，PredR2為0.934 6，信噪比為17.488，AdjR2為0.944 9。

表7 各項式估計值與P＞F 值表Tab.7 Estimates of each formula and P＞F values

表中，模型決定系數R2，PredR2，AdjR2均大于0.9，表明模型有較高的擬合度，P＞F值均小于0.000 1，說明該數學模型擬合方程的預測精度較高，可以準確地找到輸入與輸出之間的關系。圖8 為模型預測值與實際值的對比，可以看出，數據點分布在參考線附近，說明預測值與試驗值的差距較小，該模型的精度較高。

圖8 模型的預測值與實際值Fig.8 Predicted and actual values of proposed model

4 結果討論與分析

4.1 結構參數對彈性模量與孔隙率的影響

圖9 展示了結構參數對彈性模量與孔隙率的影響。多孔結構主要由氣相和固相兩部分構成，分別對應于孔隙率和孔棱直徑。對于多孔結構的彈性模量和孔隙率，Gibson 和Ashby 兩位學者研究做了大量研究，并發現彈性模量和孔隙率之間存在冪函數關系。Gibson-Ashby力學模型［10］如下：

孔棱直徑對彈性模量和孔隙率的影響十分顯著，隨著孔棱直徑的增大，多孔結構的孔隙率越來減小。這是因為孔棱直徑的增加，多孔內部的固相增加，在多孔結構外輪廓不變的情況下，氣相必然會減小，如圖10 所示，在孔棱直徑不斷增大的過程中，多孔結構內部孔洞逐漸減小，并趨于密實化。隨著孔隙率的減小，彈性模量會增大。不規則度和單元距離對彈性模量和孔隙率的影響較小，不規則度和單元距離與彈性模量呈負相關性，與孔隙率呈正相關性。

圖9 結構參數對彈性模量和孔隙率的影響Fig.9 Influence of structural parameters on elastic modulus and porosity

圖10 不同孔棱直徑對孔隙率的影響Fig.10 Influence of different strut diameters on porosity

4.2 結構參數對抗壓強度的影響

圖11 展示了結構參數對抗壓強度的影響。結果顯示，孔棱直徑對多孔結構的抗壓強度的影響十分顯著，抗壓強度隨孔棱直徑的增大而增強。不規則度和單元距離的增大會影響多孔結構的穩定性，導致抗壓強度變弱。如圖12 所示，在不規則度由小增大的過程中，多孔結構單元由均勻、規則的立方體逐漸變為不規則的多面體，部分元胞會發生畸變，導致局部元胞的尺寸變化程度加劇，元胞畸變較大的部位的應力集中現象較為突出，最終導致多孔結構抵抗外界載荷的能力變弱。

圖11 結構參數對抗壓強度的影響Fig.11 Influence of structural parameters on compressive strength

圖12 不規則度對抗壓強度的影響Fig.12 Influence of irregularity on compressive strength

優化因子及預測響應如表9 所示。試驗時，選用最大期望的結構參數，孔棱直徑0.3 mm，不規則度為0.5，單元距離為2 mm，驗證組應力-應變曲線如圖13 所示。在疲勞試驗機對多孔結構施加載荷的過程中，初期階段以多孔內部彈性變形為主要形式，此階段應力-應變曲線呈現出線性增長。繼續增加載荷，316L 多孔樣件進入局部變形階段，多孔結構被逐漸壓縮成密實狀態，此時應力的增長趨于平緩，直至多孔結構被壓縮成完全密實化的樣件。通過對模型預測值和試驗值的對比，驗證GRG 預測模型。誤差為預測值與試驗值的絕對差除以預測值計算。驗證實驗中灰色關聯度用第三節所述方法計算為0.789 5，誤差為1.2%，其中彈性模量為2.987 GPa，抗壓強度為210.048 MPa，孔隙率為89.43%。

表8 優化標準和目標Tab.8 Optimization criteria and targets

表9 預測結果與試驗驗證Tab.9 Prediction results and experimental verification

圖13 驗證組的應力-應變曲線Fig.13 Stress-strain curves of test group

5 結 論

本文基于Voronoi-Tessellation 原理設計的不規則多孔結構，利用參數化軟件Grasshopper實現了多孔結構設計的個性化和特殊化需求。SLM 制造的不規則多孔結構表現出良好的性能，其彈性模量為3.084～5.386 GPa，抗壓強度為218.014～378.595 MPa，孔隙率為50.85%～88.31%。與人骨的彈性模量3.65 GPa、抗壓強度119 MPa 相比，不規則多孔結構能夠滿足骨組織的性能需求。設計了BOX-Behnken 實驗，并結合灰色關聯分析確定了預測值數學模型，方差分析的決定系數R2，PredR2，AdjR2均大于0.9，預測模型具有較好的擬合性。試驗結果表明，孔棱直徑是影響多孔結構力學性能和孔隙率的主要因素，彈性模量、抗壓強度均隨著孔棱直徑的增大而增大，孔隙率則隨孔棱直徑的增大而減??；單元距離和不規則度與彈性模量和抗壓強度呈負相關性，與孔隙率呈正相關性，但影響程度較小。通過灰色關聯分析優化得到的最優結構參數如下：孔棱直徑0.3 mm，不規則度0.5，單元距離2 mm。采用最優結構參數制備不規則多孔結構，其彈性模量為2.987 GPa，抗壓強度為210.048 GPa，孔隙率為89.43%，GRG 為0.789 5，優化結果與預測結果的符合程度較高，誤差僅為1.2%。