髖關節表面接觸應力數學模型的建立及驗證

王沫楠，安賢俊

?

髖關節表面接觸應力數學模型的建立及驗證

王沫楠，安賢俊

(哈爾濱理工大學 機器人研究所，黑龍江哈爾濱，150080)

基于體外實驗的不可行性和有限元方法實現的復雜性，提出利用數學模型揭示髖關節表面接觸應力與幾何參數之間的關系進而模擬三維靜態髖關節表面生物力學特性，并以此來研究人體髖關節表面的承載區域、應力分布以及接觸應力峰值等屬性。同時，利用有限元方法建立單足靜立髖關節模型并進行模型分析，通過仿真結果與數學模型計算結果的對比，從髖關節表面接觸應力分布和接觸應力峰值2個方面對所建立的數學模型進行評價，驗證數學模型的精確性和普適性。最后，利用髖關節表面接觸應力數學模型對影響人體髖關節應力分布關鍵因素的量化分析，給出導致股骨頭產生病變的病理因素。研究結果表明：在人體正常活動范圍內，仿真結果與數學模型計算結果趨勢與數值吻合度滿足要求。

髖關節；接觸應力；數學模型；模型驗證；有限元法

髖關節是人體最大的負重關節，屬于球窩結構，起著傳導重力的作用。在人體眾多活動關節中，髖關節是最穩定的關節，這使得人類能夠完成較大范圍的動作。我們也發現，髖關節是臨床上最常見的多發病變關節，以髖關節發育不良而引起的關節炎為例，如果病理改變不能及時糾正，就會使髖臼脫離原來位置甚至造成脫位，不能以球窩狀覆蓋股骨頭，關節表面負重面積減少，應力增大，且集中在髖臼外緣，從而導致髖關節骨性關節炎發生及加重[1?2]。因此，建立髖關節表面接觸應力模擬模型，有助于及時發現患病關節病變狀況，為治療方案的確定提供重要的參考數據。髖關節生物力學特性研究方法主要有生物力學實驗法[3?5]及各種數值分析方法[6]。長期以來，人體骨骼生物力學特性都是通過對離體骨骼的實驗研究獲取的。如Taylor等[7]指出在髖關節或臀肌肌肉的作用下，股骨內側壓應力與外側張應力有著顯著差別；阿良等[8]通過MTS試驗機測試健康成人新鮮股骨標本，發現在縱向壓縮、三點彎曲、軸向扭轉狀態下，股骨干在同等載荷時各個截面所承受的力學分布不同。通過力學實驗獲取骨骼生物力學特性方法在應用方面的限制在于無法獲取個體活體組織的接觸應力分布模型。為了模擬活體組織的生物力學行為，多數研究者采用數值分析的仿真方法，有限元法是應用最廣泛的數值分析方法。有限元法起初應用在心血管系統的力學問題研究中，20世紀70年代開始應用于骨科生物力學研究[9]，80年代應用范圍普及到人體各部分骨骼[10]。國外有很多先進的研究成果，如William[11]應用有限元法建立了腿部仿真模型；美國華盛頓大學人機交互技術實驗室建立了快速有限元模型[12]。20世紀90年代末，國內出現一些利用有限元法分析組織性能的研究，有限元法通過仿真方法模擬人體組織生物力學特性，具有模型復雜，分析結果準確性高等優點，但是在計算速度方面不具優勢，同時分析過程復雜、對操作者要求較高。綜上所述，本文作者擬解決2個方面的問題：一是找到簡便易行的專門患者髖關節表面接觸應力模擬方法；二是驗證利用該模擬方法表征髖關節表面接觸應力的精確性。基于以上研究目標，本文作者提出通過建立精確三維靜態髖關節接觸應力數學模型，計算人體髖關節股骨頭表面的承載區域、應力分布以及接觸應力峰值等特性，從而揭示髖關節表面的力學特性與其幾何參數之間關系的創新思路；提出利用有限元方法建立虛擬髖關節模型，對不同載荷下關節表面接觸應力分布和接觸應力峰值進行分析，將分析結果與數學模型計算結果進行對照，驗證所構建接觸應力數學模型的精確性和普適性的研究構想；同時，利用驗證后的數學模型分析導致股骨頭產生病變的病理因素，以便為臨床治療和髖關節假體制作等提供參考依據。

1 方法

1.1 髖關節表面接觸應力數學模型的建立

1.1.1 建立坐標系

以右髖作為研究對象建立髖關節球面坐標系，將關節球面的球心置于坐標原點，和軸位于人體冠狀面上，軸位于人體矢狀面上。定義軸方向由身體右側指向左側，軸指向身體后側，軸方向垂直向上。圖1所示為所建立的髖關節模型坐標系，髖關節CE角如圖2所示。

1.1.2 建立接觸應力方程

因為骨與軟骨之間存在組織液起到潤滑作用，模型忽略摩擦力的影響，僅考慮壓應力。此時，作用于關節球面上的切向應力為零，頭臼間作用的合力就是徑向壓力。合力的作用點位于髖關節接觸面上，方向指向球心。設關節面上任意一點的球面坐標為(，，)，合力作用點的球面坐標為(，，)，其中和為方位角，而和為極角，為球面半徑。接觸面上的合力與關節球面上任意一點的接觸壓應力之間滿足以下積分關系式：

式中：d=2(sincos，sinsin，cos)sindd；=(sincos，sinsin，cos)。

1.1.3 壓力與承載區域的確定

由于關節軟骨的彈性模量與股骨和髖臼相比要小得多，將股骨頭視為剛性球體，將髖臼視為剛性的半球殼，與股骨頭同心，關節軟骨位于二者中間。在髖關節上施加載荷后，關節軟骨被擠壓變薄，但是，由于在生理載荷下關節軟骨的徑向應變與股骨頭半徑相比要小得多，變形屬于線彈性范圍，將關節軟骨看成是線彈性體。假設通過坐標原點且與載荷方向平行的方向線同髖臼半球面相交，根據理論力學原理可知，該交點與合力作用點重合，此時，關節球面上的最大應力max就是合力作用點上的接觸應力p。由此關節面上任意1點的接觸壓應力與合力作用點上的接觸壓應力p滿足以下線性關系：

式中：為接觸面上任一點與合力作用點之間的圓心角。

在討論股骨頭表面承載區域之前要引入髖臼中心邊緣角(central-edge angle，CE angle)的概念，它是股骨頭中心點的垂線與髖臼外側邊緣的夾角，如圖2所示。CE角是衡量髖關節接觸面積的重要指標，接觸面積隨CE角的減小而減小。假設髖臼為規則的半球形凹面，則接觸面積與CE角呈一次線性關系。

1.1.4 單足靜立狀態髖關節數學模型

在單足靜立狀態下，髖關節接觸力在垂直于冠狀面的分量僅為冠狀面內合力的3%，可以忽略不 計[13]。所以髖關節接觸力只有2個分量，即= (sin，0，cos)。規定角度，，，和CE從坐標系原點沿坐標軸正向看，繞坐標軸逆時針旋轉的角度為正。將坐標軸進行旋轉，使合力作用點位于關節球面的頂點上。則數學模型的積分關系在新建的坐標系下為：

，

完成上述積分，解得關節球面上合力作用點的接觸壓應力為

由上面的分析可知：當通過坐標原點且與載荷方向平行的方向線同髖臼半球面相交時，交點與關節球面的合力作用點重合，即關節面上的最大接觸壓應力就作用在該交點上，此時max=p；當通過坐標原點且與載荷方向平行的方向線同關節球面的交點落在接觸球面之外時，合力不通過球心，關節球面上接觸壓應力峰值max取承載區域邊緣上距該交點距離最近的點的應力，此時最大接觸應力max滿足方程(2)，即。

1.2 髖關節有限元模型及表面接觸應力分析

1.2.1 數據導入

以全髖關節為研究對象，對正常成年男子進行螺旋CT掃描，掃描長度為346 mm，CT機掃描層間距為1.5 mm，獲得的CT圖片以DICOM格式輸出并儲存在計算機中。將DICOM格式的CT數據導入醫學圖像處理軟件Mimics中，得到斷層影像。

1.2.2 圖像分割

CT圖像導入Mimics以后，利用閾值設定來確定圖像的編輯范圍，再利用骨組織與軟骨的灰度值差異將骨組織的數據提取出來，就可以實現最基本的分割。圖像分割后，對圖像中因掃描時出現的干擾數據進行修補和擦除，然后將其邊緣修剪整齊，生成邊界輪廓線，再將輪廓線以內的區域進行填充，從而得到較高精度的CT數據。

1.2.3 三維重建

對已分割出的2個部分進行三維重建，得到股骨和髖臼的三維模型。為了減少模型在導入Goemagic軟件后產生的自相交和釘狀物數量，在 FEA模塊中對三維模型進行處理，從而達到光順模型表面的目的。利用Mimics的角度測量功能可以對髖臼CE角進行測量。CE角為50°和10°的髖臼模型如圖3所示。

θCE/(°)：(a) 50；(b) 10

1.2.4 關節軟骨幾何建模

由于軟骨的灰度較小，在CT圖像上無法清晰顯示，不能利用CT圖像對關節軟骨進行三維重建。因此，假定關節軟骨充滿股骨頭與髖臼之間的空隙，將所有冠狀面上的髖關節CT影像中的頭臼空隙填充，利用CMF/Simulation模塊中的Boolean操作對模型進行布爾運算，股骨及髖臼模型被裁減出去，獲得關節軟骨的近似模型。CE角為40°的關節軟骨模型如圖4所示。

圖4 關節軟骨模型

1.2.5 有限元模型生成

首先，將生成的三維模型導入逆向工程三維檢測軟件Geomagic Studio 2012中，對模型進行修正優化，該軟件可以根據導入模型的幾何特征自動創建三角面片網格，并檢測模型的不合格網格，用戶可以通過手動操作對模型進行修改。將模型由多邊形格式轉換為點云格式，檢查是否存在體外孤點，在確定點云合格后將其封裝，完成對幾何模型三角面片的優化。然后，通過用若干個不同大小、不同曲率的曲面片對模型進行擬合，最終的三維幾何模型，以IGS格式保存，如圖5所示。最后將股骨、關節軟骨以及髖臼的IGS模型導入Unigraphics NX中進行裝配，得到髖關節的裝配模型，以IGS格式導出。裝配后的髖關節模型如圖6所示。

(a) 股骨；(b) 關節軟骨；(c) 髖臼

圖6 髖關節裝配模型

將已裝配的IGS格式髖關節幾何模型導入ANSYS 14.0中自動生成實體模型，在網格劃分前，要先根據骨骼的材料特性設定合適的單元屬性，這里定義皮質骨及關節軟骨均為各向同性且連續均勻的彈性材料。定義皮質骨的彈性模量為17 300 MPa，泊松比為0.29；軟骨的彈性模量為10 MPa，泊松比為0.4[14]。單元類型定義為10節點四面體實體單元(10solid187)，在設定好材料屬性和單元類型后，對模型進行自動網格劃分，為提高股骨頭處的分析精度，對整個股骨頭球體通過Remesh命令進行網格細化，得到最終的三維有限元模型。該模型單元總數為50 899，節點數為56 555，如圖7所示。

圖7 髖關節有限元模型

1.2.6 邊界條件與載荷設定

三維有限元模型模擬單足站立狀態，股骨遠端被完全固定，髖臼與關節軟骨被部分固定，只允許在冠狀面上垂直方向移動。因單腿站立時髖關節接觸力約為人體重力的3倍，設人體質量為70 kg，故加載于髖臼上的載荷大小為2 100 N，載荷方向位于平面內通過股骨頭球心，角度根據髖關節模型的CE角 設定。

1.3 模型驗證

為了驗證髖關節表面接觸應力數學模型，本研究建立了全髖關節的有限元模型，利用有限元分析法對髖關節進行表面接觸應力仿真。分別獲得有限元仿真模型表面應力分布和表面應力峰值，將此結果與通過數學模型計算獲得的表面應力分布和表面應力峰值進行對比，給出對數學模型精確性的評價。

2 結果

2.1 改變作用力方向髖關節表面接觸應力峰值分析結果

以CE=50°的髖關節作為研究對象，定義關節作用力方向從?40°~50°，仿真分析中每增加10°提取1個仿真結果，共10種工況，利用提取的仿真結果討論改變關節作用力方向對股骨頭表面應力分布的影響。通過有限元分析得到的不同關節作用力方向下的髖關節表面最大接觸應力如表1所示。從表1可以看出：髖關節表面應力峰值隨的增大而增大，且應力梯度逐漸增大，這與數學模型的計算結果一致。在=50°時，表面應力達到最大值6.626 7 MPa。

表1 不同關節作用力方向下的髖關節表面最大接觸應力

2.2 不同CE角時髖關節表面接觸應力峰值分析結果

設定載荷方向為=10°，加載于髖臼頂部，設定CE角從50°~10°遞減，每次減少10°。以分析過程中提取仿真結果，依據仿真結果討論CE角減小對股骨頭表面應力的影響。通過有限元分析得到不同CE角的髖關節表面最大接觸應力如表2所示。從表2可以看出：髖關節表面應力峰值變化趨勢為隨CE角的減小而逐漸增大，且CE角越小，應力峰值增大速度越快，在CE=10°時，達到最大值6.898 9 MPa。

表2 不同CE角的髖關節表面接觸應力峰值

2.3 髖關節表面接觸應力峰值綜合分析結果

表3所示為通過有限元分析得到的不同CE角的髖關節生物力學模型在不同載荷方向下的髖關節表面最大接觸應力。

表3 不同CE角度的髖關節生物力學模型在不同載荷方向下髖關節表面接觸應力峰值

圖8所示為不同CE角的髖關節股骨頭表面接觸應力峰值隨載荷方向的變化。其中通過有限元分析得到的數據用直線順次連接，平滑的曲線代表數學模型的理論曲線。

θCE/(°)：(a) 10；(b) 20；(c) 30；(d) 40；(e) 50

通過上述分析可以看出：對于不同CE角的髖關節，當CE?的差相等時，髖關節表面接觸應力峰值分析結果相近，這個特征與力學模型的計算結果相吻合。同時，從數學模型計算結果與仿真結果對比曲線圖可以看出，數學模型計算的理論曲線與有限元分析結果曲線無論從趨勢還是數值上都比較相近，驗證了本研究建立的數學模型的精確性。

2.4 髖關節表面接觸應力分布

將不同方向的載荷分別加載到不同CE角的髖關節有限元模型上，得到髖關節表面的Von-Mises應力圖，如圖9所示。

(a)CE=50°，=0°時股骨頭表面的應力分布；(b)CE=40°，=?30°時的應力分布；(c)CE=30°，=?40°時的應力分布；(d)CE=30°，=?50°時的應力分布；(e)CE=20°，=10°時的應力分布；(f)CE=10°，=?60°時的應力分布

圖9 股骨頭表面Von-Mises應力圖

Fig. 9 Von-Mises stress maps of Femoral head surface

從圖9(a)，(b)，(c)和(e)的分析結果可以看出：對不同CE角的髖關節有限元模型加載不同方向的載荷，股骨頭表面的最大接觸應力均發生在股骨頭表面與載荷方向交點的附近，并以放射狀向周邊分布，從中心到周邊逐漸減弱；當CE角為40°，=?30°時，應力峰值與CE角為30°，=?40°時的應力峰值相近，這與數學模型計算結果特性相互吻合。但通過圖9(c)和圖(d)的比較可以看出：CE角為30°，=?40°時的應力分布與=?50°時的應力分布差異不大；當=?50° 時，股骨頭表面的應力峰值產生的位置也未有明顯變化。當CE角為10°時增大載荷的偏轉角度至?60°，如圖10 (f)所示，股骨頭表面應力集中分布在髖臼與股骨頭的接觸邊緣。

3 討論

3.1 髖關節表面接觸應力數學模型表征生物力學特性

通過對式(4)單足靜立髖關節接觸應力數學模型的計算求解可知，影響人體髖關節頭臼接觸應力峰值的主要因素有髖關節球面的半徑、接觸合力、接觸合力的方位角以及人體髖關節CE角CE。只要對這幾個變量賦值，即可獲得單足靜立態時髖關節處的峰值應力、承載區域、應力分布等生物力學特性。通過數學模型計算得到髖關節表面接觸應力分布特點表現為以下幾點：

1) 當接觸力不變時，關節面的承載區域面積和最大接觸應力max隨CE角的變化而變化。CE角越大，承載區域面積越大，max則越小；CE角越小，承載區域面積越小，max則越大。

2) 當CE角不變時，關節面上的承載面積和接觸應力峰值max隨關節作用力的傾斜角度的變化而變化。當合力作用點位于冠狀面第一象限時，越大，承載區域面積越大，max則越小；越小，承載區域面積越小，max則越大。當合力作用點位于冠狀面第二象限時，越大，承載區域面積越小，max則越大；越小，承載區域面積越大，max則越小。

3) 當合力的方向與承載區域外邊緣所在平面垂直，即時，承載區域面積最大，覆蓋關節球面的一半，此時最大接觸壓應力最小，。

4) 當合力作用點位于承載區域邊界時，承載區域面積最小，近似看作一半圓弧，此時最大壓力max最大，。

3.2 髖關節表面接觸應力數學模型精確性評價

3.2.1 數學模型表征髖關節表面接觸應力峰值精確性評價

1) 當CE為定值時，改變作用力方向即的大小，有限元分析結果表明：當從?40°到50°逐漸增大，此時，髖關節表面接觸應力峰值隨變化逐漸增大。數學模型計算結果也表明：當合力作用點位于冠狀面第一象限時，越大，承載區域面積越大，max則越小；越小，承載區域面積越小，max則越大。當合力作用點位于冠狀面第二象限時，越大，承載區域面積越小，max則越大；越小，承載區域面積越大，max則越小。

2) 當作用力方向即的大小固定，改變CE的大小。有限元分析結果表明髖關節表面應力峰值隨CE增大而增大。這與數學模型計算結果一致。

3) 通過對髖關節表面接觸應力峰值的綜合分析，得到CE為10°，20°，30°，40°和50°共5組對比曲線，從對比曲線可以看出：髖關節表面接觸應力數學模型計算的理論曲線與有限元分析結果曲線無論在趨勢上還是數值上都比較接近，驗證了髖關節表面接觸應力數學模型的精確性以及在應用方面的普適性。

3.2.2 表征髖關節表面接觸應力分布評價

從圖9(a)，(b)，(c)和(e)的分析結果可以看出：對不同CE角的髖關節有限元模型加載不同方向的載荷，股骨頭表面的最大接觸應力均發生在股骨頭表面與載荷方向交點的附近，并呈現出以放射狀向周邊逐漸減弱分布的特點。當CE角為40°，=?30°時，應力峰值與CE角為30°，=?40°時的應力峰值相近，這與數學模型計算結果特性相互吻合。但通過圖9(c)和圖9(d)的比較可以看出：CE角為30°，當=?40°時的應力分布與=?50°時的應力分布差異不大，=?50° 時，股骨頭表面的應力峰值產生的位置也未有明顯變化。當CE角為10°時增大載荷的偏轉角度至?60°，如圖9(f)所示，髖關節表面應力集中分布在髖臼與股骨頭的接觸邊緣，這與本研究建立的數學模型計算結果差異較大。產生這種差異的原因是本文所構建的髖關節幾何模型頂部厚，周邊薄，使得股骨頭頂部與關節軟骨接觸充分，而周邊部分產生空隙；當逐漸加大載荷偏轉角度到一定值時，載荷作用位置的關節囊開始產生空隙，使得股骨頭表面與載荷方向交點不產生接觸壓應力，模型因關節囊的接觸問題而失效，股骨頭表面的應力集中在髖臼邊緣與股骨頭接觸的位置。

3.3 髖關節表面接觸應力數學模型的應用

圖10所示為利用數學模型得到的不同CE角的髖關節股骨頭表面接觸應力峰值隨載荷方向的變化情況。從圖10可以看出：在人體正常的活動范圍內(?20°＜＜10°)時，CE角越小，受力面積越小，股骨頭表面接觸應力峰值越大，且隨載荷方向的變化而增大的速度越快。因髖臼發育不良引起的生物力學屬性的變化主要表現為髖關節應力大幅度增加，負重面積減少，應力分布不均勻并集中于髖臼外緣。若髖關節長期在這種異常應力狀態下活動，則將導致髖臼的負重面變淺。當髖臼傾斜增大到一定程度時，髖關節變得不再穩定，頭臼匹配程度進一步降低，甚至導致髖關節脫位，使髖臼和股骨頭由正常狀態下的面面接觸變為點面接觸甚至點點接觸，使髖關節的關節軟骨局部反復承受過高的壓應力載荷，超過髖關節軟骨生理載荷限度，進而引起關節軟骨變形、軟骨下骨質硬化以及囊性變等退行性變[15?16]。

θCE/(°)：1—50；2—40；3—30；4—20；5—10

通過以上分析可知：對于患有髖關節發育不良(CE角小于20°)的患者，應盡量避免下肢在冠狀面上的運動，因為髖關節在冠狀面上的微小轉動都可能導致接觸表面應力的急劇變化。在手術方面，可以通過髖臼假體置換或添加髖關節外展支架的手段來提高髖臼的覆蓋率，減小股骨頭表面承受的應力。理論上CE角越大對髖關節的健康越有利，但具體手術方案要根據患者的具體情況綜合考慮而定。

4 結論

1) 數學模型能正確反應髖關節表面應力分布情況。當載荷偏轉角度α大于?40°時，數學模型計算結果和有限元模型分析結果相同條件下表現出相同的應力分布特征。當載荷偏轉角度α小于?40°時，由于有限元模型的接觸失效問題，導致兩者間有較大偏差。通常人體正常的活動范圍為?20°＜＜10°，所以，在此范圍內數學模型能正確反應髖關節表面應力分布情況。

2) 數學模型能正確計算不同作用力方向、不同CE角度髖關節表面接觸應力峰值。通過對CE為10°，20°，30°，40°和50°時不同作用力方向髖關節表面接觸應力峰值的仿真結果和數學模型模擬結果的對比，得出兩者在趨勢和數值兩方面都很接近。

3) 數學模型簡便易行且滿足精確性要求。利用髖關節表面接觸應力數學模型模擬髖關節表面應力避免了有限元從建模到分析的復雜過程，而且在數值和趨勢兩方面都能滿足臨床精確性要求。實現了在不進行體外生物力學實驗的前提下，利用簡便易行的方法獲取髖關節表面接觸應力信息的目標。

4) 可以利用該數學模型進行髖關節表面接觸應力分析，為臨床應用提供參考。通過數學模型模擬的髖關節幾何參數與表面接觸應力關系曲線，為髖關節疾病治療方案的制定和假體制造提供有價值的參考。

由于有限元模型的限制，沒有對數學模型推導出的“當合力的方向與承載區域外邊緣所在平面垂直時，股骨頭表面接觸應力峰值最小”結論進行驗證。同時，本文在臨床應用方面僅給出1個實例，數學模型在臨床上的進一步應用還有待探討。

[1] 肖海濤, 宋世鋒, 鄭南生, 等. 全髖關節置換術治療晚期非創傷性股骨頭缺血性壞死[J]. 現代生物醫學進展, 2012, 27(12): 5291?5295. XIAO Haitao, SONG Shifeng, ZHENG Nansheng, et al. Total hip arthroplasty in the treatment of advanced non-traumatic osteonecrosis of femoral head[J]. Progress in Modern Biomedicine, 2012, 27(12): 5291?5295.

[2] 李蕊, 陳濤, 樸繁林, 等. 髖關節置換與運動性股骨頭損傷[J]. 中國組織工程研究與臨床康復, 2011, 15(48): 9078?9081. LI Rui, CHEN Tao, PIAO Fanlin, et al. Hip replacement and sports femoral head injury[J]. Journal of Clinical Rehabilitative Tissue Engineering Research, 2011, 15(48): 9078?9081.

[3] Rushfeldt P D, Mann R W, Harris W H. Improved techniques for measuring in vitro the geometry and pressure distribution in the human acetabulum II: Instrumented endoprosthesis measurements of articular surface pressure distribution[J]. Journal of Biomechanics, 1981(14): 315?323.

[4] Hodge W A, Fijan R S, Carlson K L, et al. Contact pressures in the human hip joint measured in vivo[J]. Proceedings of the National Academy of Sciences of the USA, 1986, 8(3): 2879?2883.

[5] Polliack A A, Sieh R C, Craig D D, et al. Scientific validation of two commercial pressure sensor systems for prosthetic socket fit[J]. Prosthet Orthot Int, 2000, 24(1): 63?73.

[6] Maxian T A, Brown T D, Pedersen D R, et al. 32 dimensional sliding/contact computational simulation of total hip wear[J]. Clin Orthop, 1996, 33(3): 41?50.

[7] Taylor W R, Roland E, Ploeg H, et al. Determination of orthotropic bone elastic constants using FEA and modal analysis[J]. Journal of Biomechanics, 2012, 35(3): 767?773.

[8] 阿良, 吉士俊, 范慈方. 有限元分析先天性髖關節半脫位及髖臼發育不良的應力變化[J]. 中華小兒外科雜志, 2000, 21(6): 327?330. A Liang, JI Shijun, FAN Cifang. Stress change of congenital subluxation of hip and dysplasia of acetabular by finite element analysis[J]. Chinese Journal of Pediatric Surgery, 2000, 21(6): 327?330.

[9] Felson D T. The epidemiology of osteoarthritis: Prevalence and risk factors[C]//Kuettner K E, Goldberg V M. Osteoarthritic disor2 ders. Rosemont, IL: American Academy of Orthopaedic Surgeons, 1995: 13?24.

[10] Cowin S C. The mechanics properties of cancellous bone[J]. In Bone Mechanics, 1989, 8(1): 30?57.

[11] William J L, Mary C P, Glimcher M J. Electron microscopic observations of bone tissue prepared anhydrously in organic solvents[J]. Journal of Ultrastructure Research, 1977, 59(2): 185?206.

[12] McGibbon C A, Krebs D E, Trahan C A, et al. Cartilage degeneration in relation to repetitive pressure: case study of a uniateral hip hemiarthroplasty patient[J]. J Arthroplasty, 2010, 14(1): 52?58.

[13] 蘆俊鵬, 張建國. 人體顱腦三維有限元模型構建[J]. 微計算機信息, 2012, 22(8): 211?213. LU Junpeng, ZHANG Jianguo. A development of 3D FEM for cranium brain based on CT[J]. Microcomputer Information, 2012, 22(8): 211?213.

[14] 王天英, 姜福川, 王西十. 人體髖關節接觸力計算分析一例[J]. 力學與實踐, 2012, 27(2): 67?69. WANG Tianying, JIANG Fuchuan, WANG Xishi. Human hip contact force analysis[J]. Mechanics in Engineering, 2012, 27(2): 67?69.

[15] Silva M J, Keaveny T M, Hayes W C. Load sharing between the shell and centrum in the lumber vertebral body[J]. Spine, 2009, 22(2): 140?150.

[16] Hashimoto S , Nishi yama T, Hayashi S, et al. Role of p53 in human chondrocyte apoptosis in response to shear strain[J]. Arthritis Rheum, 2010, 60(8): 2340?2349.

(編輯 陳愛華)

Mathematics modeling and verification of hip joint surface contact stress

WANG Monan, AN Xianjun

(Robotics Institute, Harbin University of Science and Technology, Harbin 150080, China)

For patients, vitro experiment is not feasible and finite element method is more complex. 3D static hip joint biomechanical behaviors were expressed based on mathematics model between biomechanics and geometry parameter. Hip joint surface bearing area, stress distribution and peak contact stress were calculated from mathematical model. Then the same single static hip joint finite element model was developed and analysed. The results including femur surface stress distribution and peak contact stress from finite element simulation were compared to those from mathematical model in order to verify accuracy and universality mathematical model. Finally, femur pathologic mechanism was established based on the quantitative analyses results of mathematical model. The results show that the curve obtained from simulation is similar to the curve from mathematical model in trend and value.

hip joint; contact stress; mathematics modeling; model verification; finite element method

10.11817/j.issn.1672-7207.2015.06.015

TB12；R318.01

A

1672?7207(2015)06?2081?09

2014?06?13；

2014?08?20

教育部新世紀優秀人才支持計劃(NCET-13-0756)；國家自然科學基金(61272387)；黑龍江省杰出青年基金(JC201302)(Prooject (NCET-13-0756) supported by the New Century Excellent Talents in University; Project (61272387) supported by the National Natural Science Foundation of China; Project (JC201302) supported by the Province in Heilongjiang Outstanding Youth Fund)

王沫楠，博士，教授，從事計算機輔助醫療研究；E-mail：qqwmnan@163.com