[發明專利]一種基于拓撲優化技術的松質骨多孔結構建模方法

有效

• 申請號: 201710410036.0
• 申請日: 2017-06-03
• 公開（公告）號: CN107391784B
• 公開（公告）日: 2020-08-14
• 發明（設計）人: 劉蓉;張亞坤;郭新路;王永軒
• 申請（專利權）人: 大連理工大學
• 主分類號: G06F30/17
• 分類號: G06F30/17;G06F30/23;G06F119/14
• 代理公司: 大連星海專利事務所有限公司 21208
• 代理人: 王樹本

1.一種基于拓撲優化技術的松質骨多孔結構建模方法，其特征在于，包括以下步驟：

步驟1、建立松質骨模型，為了保證松質骨模型的真實性，基于醫學圖像建立松質骨模型，具體包括以下子步驟：

(a)根據設計要求選取股骨頭并對其進行物理處理，利用Micro-CT儀器掃描得到松質骨結構的醫學圖像；

(b)利用Mimics軟件對Micro-CT儀器掃描得到的醫學圖像進行處理，包括窗寬和窗位調整、閾值分割、區域增長，得到具有清晰骨小梁結構的松質骨圖像，并建立松質骨結構模型；

步驟2、求解松質骨模型的體積分數，將在Mimics中建立得到的松質骨模型導出為STL格式，并將其導入到Magics軟件中得到松質骨模型的體積信息，并基于此體積信息對松質骨模型的體積分數BV/TV進行求解，其中BV是松質骨實心部分的體積，TV是松質骨實心部分和空心部分的總體積；

步驟3、求解松質骨模型的力學性能，選用代表體元法對模型施加約束條件和外載荷，并基于ANSYS軟件對松質骨模型的9個力學性能參數進行求解，具體包括以下子步驟：

(a)對松質骨模型施加約束，可通過公式(1)進行描述，

式中，1,3,5行分別代表x＝0，y＝0和z＝0的面分別在x方向，y方向和z方向的位移為0；u，v，w分別是x，y，z三個方向的位移；松質骨模型代表體元在x軸方向,y軸方向和z軸方向的長度分別為a,b,c；2,4,6行分別代表x＝a，y＝b，z＝c的面的位移分別為δ_x，δ_y，δ_z，它們的值由ANSYS求解得到；

(b)求解松質骨模型的等效彈性模量，包括三個等效彈性模量，分別是x方向的彈性模量E_x，y方向的彈性模量E_y，z方向的彈性模量E_z，可通過公式(2)、(3)、(4)進行描述，

式中，p為施加的外載荷，它的大小不影響結果；a,b,c分別為松質骨模型代表體元在x軸方向,y軸方向和z軸方向的邊長；

(c)求解松質骨模型的等效泊松比，包括三個等效泊松比，分別對應于xy面的泊松比v_xy，yz面的泊松比v_yz和zx面的泊松比v_zx，可通過公式(5)、(6)、(7)進行描述，

式中，δ_xy為在x＝a的面施加外載荷p后y＝b的面的位移，δ_yz為在y＝b的面施加外載荷p后z＝c的面的位移，δ_zx為在z＝c的面施加外載荷p后x＝a的面的位移；

(d)求解松質骨模型的等效剪切模量，包括三個等效剪切模量，分別是與xy面平行的剪切模量G_xy，與yz面平行的剪切模量G_yz，與zx面平行的剪切模量G_zx，可通過公式(8)、(9)、(10)進行描述，

步驟4、提出本方法的拓撲優化模型，其中包括目標函數、約束條件以及控制方程，具體包括以下子步驟：

(a)引入拓撲優化中的材料模型，選用Solid Isotropic Material WithPenalization(SIMP)方法，可通過公式(11)進行描述，

F(x)＝ρ(x)^qF⁰ (11)

式中，F(x)為材料點x處的單元材料的彈性模量，F⁰為材料本身的彈性模量，ρ(x)為單元材料的密度，q是懲罰因子；

(b)為了能夠實時監測迭代過程中的材料模型的等效性能，采用均勻化方法通過公式(12)、(13)對迭代過程中的材料模型性能進行求解，

式中，Y為單胞域，是符合上式的廣義位移函數；v_i表示單胞中定義的連續函數；E_ijkl表示單胞中材料的彈性常數；i，j,k，l表示坐標編號，二維取1，2，三維取1，2，3，E_ijmn表示y_n處的四階彈性張量，m，n同樣是坐標編號，二維取1，2，三維取1，2，3，y_n，y_j表示不同維度的坐標，V_Ω表示設計域，通過公式(12)求解得到廣義位移函數然后通過公式(13)求解得到宏觀等效的彈性模量

(c)基于SIMP材料插值模型和均勻化方法，以松質骨結構的九個性能參數和體積分數為設計目標提出了結構建模優化列式，可通過公式(14)進行描述，

式中，ρ_e表示單元密度，具體可以參考SIMP方法；e表示單元編號，N表示單元的總數目；f(ρ_e)表示目標函數，E_x、E_y、E_z、G_xy、G_yz、G_zx、v_xy、v_yz、v_zx為多孔結構的九個設計參數，同時也是要設計要達到的性能，表示設計迭代過程中九個等效性能實時的值，V^upp表示設計的體積上限，V表示單元材料分布的總體積，F(x)＝ρ(x)^qF⁰，F(x)為材料點x處的單元材料的彈性模量，表示引入的密度插值函數，Homogenization表示均勻化方法的控制方程；

步驟5、求解多孔結構，首先確定設計域，給定一個r×r×r的正方體網格設計域，建議給定敏度過濾半徑r_min為1.5，建議最大迭代步數為100步，九個目標性能參數，體積約束，給定實心材料無孔隙條件下的彈性模量和泊松比，控制某一特定閾值以上的單元在ANSYS中顯示，所述特定閾值由具體情況進行選取，范圍在0-1之間；

步驟6、還原CAD模型，將ANSYS中顯示的多孔結構利用CAD軟件建立得到符合實際的CAD模型。