基于等幾何分析的非均質(zhì)股骨近端模型的受力分析

產(chǎn)啟平，陳 龍

(上海理工大學(xué) 機(jī)械工程學(xué)院，上海 200093)

基于等幾何分析的非均質(zhì)股骨近端模型的受力分析

產(chǎn)啟平，陳 龍

(上海理工大學(xué) 機(jī)械工程學(xué)院，上海 200093)

人工股骨置換手術(shù)是目前治療股骨骨折比較成熟的一種方法。但由于人體骨骼的復(fù)雜性，采用單一材料均勻分布進(jìn)行力學(xué)性能分析時(shí)難以得到準(zhǔn)確結(jié)果。文中采用3種不同的材料賦值方法進(jìn)行材料非均質(zhì)建模以提高模型建模精度，采用等幾何分析方法以提高性能分析精度。從股骨幾何點(diǎn)云和材料灰度信息數(shù)據(jù)出發(fā)，用體參數(shù)化的建模理論，得到人體股骨的B樣條體參數(shù)化模型，然后采用中心點(diǎn)插值、節(jié)點(diǎn)插值、高斯積分點(diǎn)3種插值法賦予模型10種不同的材料，得到非均質(zhì)模型。根據(jù)有限元基礎(chǔ)理論推導(dǎo)出剛度矩陣，并施加邊界和約束條件，最終通過求解線性方程組得到靜力學(xué)分析結(jié)果。實(shí)驗(yàn)表明，文中建模方法精確高效，得到的結(jié)果合理且更接近實(shí)際。

等幾何分析; 非均質(zhì); 體參數(shù)化; 股骨; 靜力學(xué)分析

目前治療股骨骨折的常見方法是人工股骨置換手術(shù)，鑒于股骨的重要性，所以對置換假體的性能要求也非常高。置換假體不僅要接近真實(shí)形態(tài)，而且要承受的住日常生活所施加的載荷，使其避免出現(xiàn)斷裂、錯(cuò)位等現(xiàn)象。因此置換前，需對股骨進(jìn)行精確的模型構(gòu)建和力學(xué)性能分析。

對于股骨的分析需要建立精確的三維模型和準(zhǔn)確的材料賦值，這也是進(jìn)行股骨力學(xué)分析的基礎(chǔ)和關(guān)鍵。由于股骨結(jié)構(gòu)的復(fù)雜性，幾何形態(tài)的不規(guī)則性，研究表明直接用CAD、CAE軟件進(jìn)行建模很困難[1]。目前國內(nèi)外對于人體股骨的重建有多種方法：(1)三維坐標(biāo)儀建模法[2]；(2)通過掃描圖片利用CAD逆向三維重建出模型[3]；(3)組織切片法[4]。另一方面，目前很多關(guān)于股骨的研究是將其材料假設(shè)為均勻分布的密質(zhì)骨和松質(zhì)骨，更有將其設(shè)置為一種材料的研究，針對股骨的復(fù)雜性，這樣的賦值會對研究產(chǎn)生嚴(yán)重的誤差。本文結(jié)合真實(shí)股骨的材料分布，采用3種插值方法對股骨模型進(jìn)行10種不同材料的賦值，實(shí)現(xiàn)了股骨材料建模。分析方面目前大多數(shù)采用的是有限元方法[5-7]。由于傳統(tǒng)有限元分析中CAD設(shè)計(jì)模型和CAE分析模型往往是在不同的平臺下完成，導(dǎo)致CAE模型的分析結(jié)果不能夠準(zhǔn)確反應(yīng)原CAD模型的特征，而且對于越復(fù)雜的產(chǎn)品模型，誤差越大。

HUGHES[8-11]等在傳統(tǒng)有限元方法基礎(chǔ)上提出了等幾何分析方法(Isogeometric Analysis，IGA)。等幾何分析方法基于有限元分析方法的等參單元思想，將模型設(shè)計(jì)中用于表達(dá)幾何模型的非均勻有理B樣條(NURBS)的基函數(shù)作為形函數(shù)，實(shí)現(xiàn)了計(jì)算機(jī)輔助設(shè)計(jì)(CAD)和計(jì)算機(jī)輔助工程(CAE)的無縫結(jié)合。

本文采用等幾何分析方法對股骨進(jìn)行分析。由于幾何模型可直接用于分析，不需要經(jīng)過轉(zhuǎn)換、處理等操作，顯著提高了求解精度和計(jì)算效率。

1 股骨體參數(shù)化模型構(gòu)建

體參數(shù)化模型的構(gòu)建，是通過對已知邊界面進(jìn)行體插值，實(shí)現(xiàn)與模板模型的一一映射，體參數(shù)化模型是等幾何分析的基礎(chǔ)。體參數(shù)化模型的構(gòu)建過程如圖1所示。

圖1 體參數(shù)化模型構(gòu)建

(1)點(diǎn)云數(shù)據(jù)預(yù)處理。點(diǎn)云分片要求將點(diǎn)云模型分解為數(shù)目有限、具有相似特征且互不重疊的子區(qū)域。在參照特征邊界點(diǎn)的前提下，采用交互式分割法，將點(diǎn)云模型分割為6個(gè)點(diǎn)云面片。

每個(gè)完成分割的點(diǎn)云面片有4條邊界，采用散亂點(diǎn)云最小二乘B樣條曲面擬合算法來實(shí)現(xiàn)6個(gè)點(diǎn)云面片的曲面擬合。其步驟為：1)根據(jù)分片曲面的散亂點(diǎn)及其邊界曲線構(gòu)造基面；2)將散亂數(shù)據(jù)點(diǎn)投影到已構(gòu)造好的參數(shù)曲面上進(jìn)行投影參數(shù)化；3)采用最小二乘法 實(shí)現(xiàn)散亂點(diǎn)云的曲面初始重建；4)將初始重建曲面作為基面迭代。再將曲面拼接起來形成一個(gè)封閉的空間。

(2)B樣條體構(gòu)建。在對六曲面擬合完成以后，只是得到體參數(shù)化模型的表面控制點(diǎn)，通過對六個(gè)表面的曲面片進(jìn)行插值得到初始內(nèi)部控制點(diǎn)，為提高網(wǎng)格頂點(diǎn)的質(zhì)量，需要對這些控制點(diǎn)進(jìn)行微調(diào)，避免網(wǎng)格點(diǎn)的自交，從而得到最終的體內(nèi)控制點(diǎn)。在所有控制點(diǎn)確定之后，即可得到B樣條體的可視化模型。根據(jù)B樣條理論，三維NURBS體的表達(dá)式如下

(1)

其中，S(ξ,η,ζ)指用空間中的點(diǎn)(ξ,η，ζ)表示的三維參數(shù)體；Pi,j,k(i=0,1,…,n;j=0,1,…,m;k=0,1,…,l)表示控制點(diǎn)；Ri,j,k(ξ,η,ζ)表示3變量有理基函數(shù)，其表達(dá)式為

(2)

基于上述理論，然后按照前述方法構(gòu)建了股骨近端參數(shù)化模型如圖2所示，可以看出，基于體參數(shù)化方法構(gòu)建的人體股骨三維模型和真實(shí)的人體股骨模型結(jié)構(gòu)外觀相吻合。

圖2 股骨近端體參數(shù)化模型

2 股骨體參數(shù)化模型材料賦值

目前國內(nèi)外關(guān)于骨骼的分析材料的選取，只考慮密質(zhì)骨和松質(zhì)骨兩種材料的簡化，甚至只考慮單均一材料的影響。傳統(tǒng)這種材料賦值方法遠(yuǎn)遠(yuǎn)不能夠滿足分析的要求。根據(jù)張國棟[12-13]等人關(guān)于骨骼材料賦值的研究和綜合考慮計(jì)算效率，將股骨按照圖像的灰度值完成10種材料特性的賦值(見表1)，其中泊松比[14]為0.3，彈性模量為Mimics提供經(jīng)驗(yàn)公式

Density=-13.4+1 017×Grayvalue (3)

對于非均質(zhì)材料，希望B樣條的密度函數(shù) 能夠應(yīng)用于各種分析中，對于分析所用的密度可以用不同的方法估算得到。在本文提供了3種方法來估算有限單元中，B樣條表達(dá)的密度分布。

(5)

圖3 B樣條系數(shù)與單元中心

(2)基于節(jié)點(diǎn)的近似值。密度由單元節(jié)點(diǎn)估算而來。即單元密度為

(6)

圖4 網(wǎng)格與節(jié)點(diǎn)

(3)基于高斯積分點(diǎn)的近似值。密度由單元積分點(diǎn)處的值估算得到。則單元密度為

(7)

圖5 NURBS離散化與積分點(diǎn)

3 等幾何分析方法

等幾何分析和傳統(tǒng)有限元最本質(zhì)的區(qū)別是在于基函數(shù)的選擇，其他方面很相似， 根據(jù)有限元方法中的虛功原理及最小勢能原理推導(dǎo)出單元的剛度矩陣

(8)

(9)

其中，D0是涉及單元材料的彈性矩陣；ωj是積分點(diǎn)j處的權(quán)重；E(ρ(pj))是積分點(diǎn)j處的彈性模量。用SIMP插值算法，彈性模量可以表示為

E(ρ(x))=Emin+ρ(x)s(E0-Emin)

(10)

其中，E0是固體的剛度；s是懲罰系數(shù)；ρ(x)是單元的物理密度；Emin是最小彈性模量，避免出現(xiàn)剛度矩陣奇異。

如果ρ是常數(shù)，所有積分點(diǎn)的密度被認(rèn)為和單元中心處的密度一樣，即ρ(pj)=ρ(pc)。因此單元的剛度矩陣可以被簡化為

(11)

(12)

其中，Nk,j表示節(jié)點(diǎn)k在單元的第j個(gè)積分點(diǎn)處估算的形函數(shù)。

(13)

得到單元?jiǎng)偠染仃嚭?，借鑒傳統(tǒng)有限元中整體剛度矩陣組裝的方法進(jìn)行組裝得到整體剛度矩陣K，從而得到總體平衡方程

Kd=F

(14)

其中載荷矢量

(15)

4 股骨模型的分析

基于上述原理與流程，分析人體正常站立位時(shí)股骨近端的縱向受壓情況。本文簡化邊界條件，將股骨上端施加600 N豎直向下的壓力，股骨下端定為完全約束，即底面上所有的節(jié)點(diǎn)的自由度為0，分別采用3種材料賦值方法對模型賦值進(jìn)行分析計(jì)算，計(jì)算結(jié)果如圖所示。

圖7 應(yīng)力云圖

圖8 單一材料分析結(jié)果

圖6分別為3種賦值方法得到的各點(diǎn)的位移變化情況，圖7分別為3種方法得到的應(yīng)力分布特征。從上述分析結(jié)果可以看出，3種賦值方法得到的分析結(jié)果基本相同，應(yīng)力主要集中在股骨的上下兩端，其3種賦值方法所得最大應(yīng)力分別為13.6 MPa，14.1 Mpa和12.9 MPa，與文獻(xiàn)[17]中股骨縱向受壓仿真得到的結(jié)果非常接近。圖8為采用單一材料分析得到的結(jié)果，顯然應(yīng)力集中的地方不是很明顯，且最大應(yīng)力為9.8 MPa，與本文仿真結(jié)果相差較大。因此本文能更好的顯示出股骨在載荷作用下的應(yīng)力分布特征。

5 結(jié)束語

本文采用體參數(shù)化建模方法，對股骨模型進(jìn)行了精確、快速重建，并通過掃描得到的灰度值，采用3種不同的方法對股骨進(jìn)行材料賦值，得到非均質(zhì)股骨模型，使得所建模型更加逼真，從而提高了分析的可靠性。對股骨縱向受力分析采用等幾何分析方法來代替?zhèn)鹘y(tǒng)的有限元方法，從而使得設(shè)計(jì)模型與分析模型一致，避免了因模型轉(zhuǎn)換帶來的誤差，進(jìn)一步提高了計(jì)算精度。

[1] 胡小春,郭松青,葉銘.基于CT圖像建立人體足部骨骼三維有限元模型的研究[J].合肥工業(yè)大學(xué)學(xué)報(bào),2006,28(9): 1189-1191.

[2] 張美超,廖進(jìn)民,李敏,等.激光三維掃描系統(tǒng)重建下頜骨[J].第一軍醫(yī)大學(xué)學(xué)報(bào),2004,24(7):756-757.

[3] 胡鑫,習(xí)俊通,顧冬云,等.髖關(guān)節(jié)CAD模型的建立與幾何形態(tài)分析[J].北京生物醫(yī)學(xué)工程,2003,22(2):116-120.

[4] 李安安,劉謙,龔輝,等.“虛擬中國人男性一號”高精度骨骼系統(tǒng)的三維建模[J].中國臨床解剖學(xué)雜志,2006,24(3): 292-294.

[5] 姜海波,葛世榮.基于CT掃描人體股骨的有限元分析[J].工程力學(xué),2007,24(10)：156-159.

[6] 羅承剛,龔憲生.近端股骨的非均勻及各向異性有限元模擬[J].重慶大學(xué)學(xué)報(bào)，2004 (2):20-23.

[7] 劉海,謝力勤,徐孝根,等.給定荷載下的骨骼重建的有限元分析[J].清華大學(xué)學(xué)報(bào),2002,42(8):1079-1082.

[8] T J R Hughes,J A Cottrell,Y Bazilevs. Isogeometric analysis:CAD,finite elements,NURBS,exact geometry and mesh refinement[J].Computer Methods in Applied Mechanics Engineering,2005(194):4135-4195.

[9] Cottrell J A,Hughes T J R,Bazilevs Y.Isogeometric analysis: toward integration of CAD and FEA[M].New York:Wiley Press,2009.

[10] Kim H J,Seo Y D,Youn S K. Isogeometric analysis for trimmed CAD surfaces[J].Computer Methods in Applied Mechanics and Engineering, 2009(198):2982-2995.

[11] Noor A K.Evaluation of element stiffness matrices on CDC STAR-100 computer[J].Computers and Structures,1978(9):151-161.

[12] 張國棟,廖維靖,陶圣祥,等.股骨頸有限元分析的賦材料屬性方法探討及有效性驗(yàn)證[J].中國組織工程研究與臨床康復(fù),2009,13(52) :8436-8441.

[13] Lengsfeld M,Schmitt J,Alter P,et al. Comparison of geometry-based and CT voxel-based finite element modelling and experiment validation[J].Medicine Engineering Physics,1998,20(7):515-522.

[14] 高建偉,陳龍.基于等幾何分析的股骨模型靜力學(xué)分析[J].計(jì)算機(jī)仿真,2015, 32(5): 340-343.

[15] 陳樹廣.基于病例CT圖像的健康股骨頭建模研究[J].計(jì)算機(jī)仿真,2009,26(9): 172-175.

[16] 韓樹洋,葛世榮.人工全髖關(guān)節(jié)置換術(shù)對天然股骨生物力學(xué)行為的影響[J].醫(yī)用生物力學(xué),2010,12(6):471-474.

[17] 陳研妤,應(yīng)祖光,朱振康,等.股骨縱向受壓時(shí)的應(yīng)力分布特征分析[J].固體力學(xué)學(xué)報(bào),2016,37(5):74-78.

Stress Analysis of Heterogeneous Femoral Proximal Model by Using Isogeometric Analysis

CHAN Qiping,CHEN Long

(School of Mechanical Engineering,University of Shanghai for Science and Technology,Shanghai 200093,China)

Artificial femoral replacement surgery is a method for the treatment of femoral fractures. However it is difficult to get a currate results by using. Because of the human skeleton is very complex, uniform distribution or single material simulation. In this paper, three different interpolation methods were used to assign the 10 different materials to the model, and the femur was modeled and analyzed by geometric analysis. The geometrical information expressed in the form of point cloud and the material data expressed in grayscale information were used as the basis of geometric modeling and material modeling. The parametric geometry of the femur was obtained by the parametric construction theory of body model. Then, we used the central point interpolation, the node interpolation and the Gaussian integral point interpolation method to assign the model 10 different materials to obtain the heterogeneous model. According to the principle of virtual work and the minimum potential energy principle, the element stiffness matrix and the overall stiffness matrix were deduced, and the boundary condition was imposed. Finally, the statics analysis results were obtained by solving the linear equations. Experiments showed that the modeling method was accurate and efficient, and the heterogeneous femur model obtained by the three interpolation methods was reasonable and closer to the actual result.

isogeometric analysis; heterogeneity; volume parameterization; femur; static stress analysis

2017- 03- 29

國家自然科學(xué)基金(51475309)

產(chǎn)啟平(1991-)，男，碩士研究生。研究方向：CAD/CAE。陳龍(1978-)，男，博士，副教授。研究方向：產(chǎn)品計(jì)算設(shè)計(jì)、產(chǎn)品快速設(shè)計(jì)、圖像處理技術(shù)等。

TN04；TP391.9

A

1007-7820(2018)02-015-05