本文基于橫觀各向同性理論建立關節軟骨固液耦合雙相三維缺損及修復的有限元模型。本文通過研究鄰近修復界面的宿主軟骨的應力狀態判別其變形類型，探討致軟骨修復界面開裂的原因。研究表明，表層缺損修復時，鄰近修復界面的宿主軟骨表面節點發生壓縮變形；中間層、深層或全層缺損修復時，節點發生拉伸變形，此時軟骨徑向尺寸增加，修復界面易開裂。若采用全層缺損修復，組織工程化軟骨（TEC）的彈性模量較低（0.1 MPa、0.3 MPa）時，鄰近修復界面的宿主軟骨表層和中間層主要發生拉伸變形；TEC 的彈性模量較高（0.6 MPa、0.9 MPa）時，宿主軟骨各層均發生壓縮變形。因此，全層缺損修復時，可適當增大 TEC 的彈性模量。本文為評估軟骨組織工程修復效果提供了新的思路，或對臨床有一定的指導意義。
引用本文： 劉海英, 趙永政, 趙森, 馮晶晶, 張春秋. 關節軟骨缺損修復區應力狀態分析. 生物醫學工程學雜志, 2018, 35(5): 705-712, 719. doi: 10.7507/1001-5515.201712083 復制
引言
關節軟骨是人體關節的重要承力組織，主要功能是減小摩擦和磨損及將關節軟骨表面的載荷均勻地傳遞到軟骨下骨，對活動關節起著非常重要的支撐和保護作用。隨著國民生活水平的提高，高能、高強度運動增多，且人口老齡化日趨嚴重及肥胖人口比例上升，關節軟骨損傷成為骨科臨床常見疾病之一，若不加以有效治療最終會發展成骨性關節炎[1-2]。由于軟骨組織內沒有血供和淋巴回流為其提供營養，損傷后很難自行修復，軟骨損傷后的有效治療成為臨床上亟待解決的難題[3]。近年來，隨著組織工程技術逐漸成熟，使其有希望成為治療軟骨損傷的最有效方法[4-5]。盡管該方法在臨床上取得了一定進展，但大多短期療效較好，長期修復效果具有不確定性。主因之一是目前還無法構建出可媲美天然軟骨生物力學性能的組織工程化軟骨（tissue engineering cartilage，TEC），植入的 TEC 與宿主軟骨間力學性能的差異顯著地改變了修復界面及鄰近宿主軟骨和細胞的力學環境，常出現移植物退化變性、修復界面開裂及宿主軟骨降解等現象[6-7]。因此，鄰近修復界面宿主軟骨的力學狀態是影響修復效果的關鍵因素之一。
關節軟骨是由固相（包括：蛋白多糖、膠原纖維）和液相（水）構成的多孔結構，孔隙中充滿組織液，液相約占軟骨重量的 65% 以上，有的部位可達 80%，具有顯著的非線性粘彈特性。加載初期，液相能承擔高達 90% 的載荷從而保護膠原纖維網免于高應力的破壞，且組織液在壓力梯度作用下流動利于軟骨中營養物質的傳遞和交換。因此，學者們普遍認為固液耦合雙相模型是研究軟骨承載特性的基礎，包括纖維增強雙相模型和橫觀各向同性雙相模型。由于軟骨中纖維的直徑、含量及走向依深度變化，在軟骨的二維有限元模型中可以使用彈簧單元或者桿單元來表征膠原纖維，而纖維的空間分布規律卻沒有相應描述[8]。實際上，表層膠原纖維的走向與軟骨表面平行，如同編織的毯子，從而在軟骨表面形成了一層不透水的纖維網絡。Clark[9]據此提出了軟骨的橫觀各向同性模型，經研究驗證，此模型能夠成功地預測應力集中的位置[10]。隨著關節軟骨橫觀各向同性理論逐漸成熟，Elhamian 等[11]應用實驗和理論相結合的方法提出了關節軟骨隨深度變化的橫觀各向同性微觀力學模型，該模型進一步促進了相關研究的發展。
體外實驗研究表明，應力狀態下軟骨細胞的生物力學環境的改變致軟骨基質分解代謝的變化能顯著地影響關節軟骨的生理功能[12]。孟迪等[8]通過有限元模擬研究了缺損修復形狀對修復區的米塞斯（Mises）應力分布影響的規律。Mises 應力是基于第四強度理論（剪切應變能）的等效應力，幾乎所有的有限元分析結果中默認的都是 Mises 應力。通過 Mises 應力云圖可確定軟骨缺損修復的危險區域，是分析應力分布變化規律的重要指標，但 Mises 應力是標量，沒有力的方向性指示。因此，無法展現修復區復雜的力學環境，鄰近修復界面宿主軟骨的應力狀態可能與修復界面的整合情況密切相關，而相關研究卻鮮見報導，常用空間某一點的三個主應力（
、
、
）表示該點的應力狀態。何祝斌等[13-16]給出了在主坐標系中，空間一點正應力、剪應力的描述方程，并繪制了一點的剪應力三維圖形和偏應力空間中的正應力圖形，從圖形中能判斷不同應力狀態下材料發生變形的類型及尺寸變化的規律，為結構設計或應力、應變分析提供了一種直觀有效的方法。
本文建立了基于橫觀各向同性的關節軟骨固液耦合雙相三維缺損及修復的有限元模型。通過數值模擬和理論分析，對比研究表層、中間層、深層和全層軟骨缺損修復的情況下，鄰近修復界面的宿主軟骨表面層單元上節點的應力狀態、變形類型及材料尺寸的變化規律；并對比分析使用不同彈性模量（0.1、0.3、0.6、0.9 MPa）的 TEC 修復全層軟骨缺損后，鄰近修復界面各層宿主軟骨的應力偏量值（
、
、
）隨深度變化的規律。在三向應力狀態下，預測軟骨缺損修復界面開裂的原因具有重要的理論和臨床實際意義。
1 有限元模型構建及材料參數
?
關節軟骨位于股骨遠端和脛骨近端的關節面上，有研究表明，當軟骨缺損直徑小于 3 mm 時依靠自身的修復能力可部分或全部修復，而直徑大于 3 mm 則不能自行修復[17]。正常人體關節軟骨的厚度為 1～5 mm，由于其形狀不規則，通常采用簡化模型進行研究。本文應用有限元軟件 ABAQUS 6.10（SIMULIA，法國）建立厚 2 mm、半徑 15 mm 的圓柱形脛骨軟骨幾何模型，而圓柱形 TEC 的直徑為 4 mm，用半徑 20 mm 的剛性球模擬股骨軟骨[18]。如圖1所示，為便于施加位移載荷，在球體中心設置一個施加載荷的點；在鄰近修復界面左側的宿主軟骨內取一條路徑（如圖1中路徑一所示）、宿主軟骨表面層取一節點（如圖1中位置一所示）用于后續的應力狀態分析。
圖1
關節軟骨全層缺損修復的非圍限壓縮有限元模型截面 示意圖
Figure1.
The sectional sketch of unconfined compression finite element model for articular cartilage with full-thickness defects repair
膠原纖維是構成關節軟骨三維網狀結構的重要組成成分，能使關節適應各種機械負荷。根據膠原纖維的排列情況，關節軟骨在形態學上可分為表層、中間層和深層。表層占軟骨厚度的 10%～20%，膠原纖維較細密且走向與軟骨表面平行；中間層占軟骨厚度的 40%～60%，膠原纖維較粗且走向與軟骨表面交叉、排列不整齊；深層占軟骨厚度的 30%，纖維粗大且走向與軟骨表面垂直[19]。由于各層膠原纖維和蛋白多糖的含量、體積分數及纖維走向不同，導致關節軟骨的力學性能依深度而變化。
Elhamian 等[11]通過理論和實驗分析得到了基于橫觀各向同性理論的關節軟骨材料參數隨深度變化的曲線圖，其中軟骨模型的空間坐標系為軸向z、徑向r和圓周方向θ。各材料參數均隨深度方向的坐標值z變化：如彈性模量Ez（z）、Er（z）、Eθ（z），泊松比vrθ（z）、vrz（z）、vθz（z），切變模量Grθ（z）、Grz（z）、Gθz（z）。本文應用數學軟件 Mathematica 7.0（Wolfram Research，美國）對從曲線上提取的多組數據進行擬合，得到材料參數隨歸一化深度
值（
，h = 2 mm 為軟骨厚度）變化的非線性關系表達式，如式（1）所示。由于篇幅所限，僅以彈性模量
為例：
'/> |
軟骨有限元模型沿厚度共分為 10 層，其中 1～2 層、3～7 層、8～10 層分別為軟骨的表層、中間層、深層。將平分每層厚度的
值代入各非線性關系式可得材料參數，如表 1 所示。
表1
關節軟骨材料參數
Table1.
Material parameters of the articular cartilage
滲透率k是表征流體通過多孔材料能力的參數，其大小取決于材料的孔隙率和液體滲透方向上孔隙的幾何形狀等因素，k與應變相關，如式（2）所示[20]：
 |
式中，k0 = 1.743 × 10–15 m4/（N·s），為初始滲透率；M = 7.081 為材料滲透常數；e和e0 分別為當前孔隙率和初始孔隙率。e0 隨深度變化的關系式，如式（3）所示[21]：
 |
式中，
為軟骨表層初始孔隙比；αe = 0.413 為材料常數。TEC 的初始孔隙比取為 9[20]。軟骨底部固定且不可滲透，設軟骨側面孔壓為 0，可實現液相流動及滲出[22]。
在球體中心處施加位移載荷。首先，剛性壓頭向下移動與軟骨接觸；之后，施加z向位移載荷，確保壓縮量達到軟骨厚度的 30%；最后，令壓頭保持某一固定值，給定足夠松弛時間分析其力學行為。
2 結果
由于根據 Mises 應力云圖無法得知危險區域受拉應力還是壓應力，也無法確定該區域是主應力還是剪應力占主導，因而無法判別修復界面開裂的原因。因此研究者們更應關注應力在不同區域、不同方向上的形式。通過繪制主應力矢量圖可清晰地呈現物體內部的受力情況，如拉、壓應力的大小及方向。
2.1 關節軟骨缺損修復區的應力分布
關節軟骨具有流固耦合特性，液相在壓力梯度作用下流動持續地改變著修復區的力學環境。本文以彈性模量 0.3 MPa 的 TEC 修復全層缺損為例，數值模擬應力松弛情況下不同加載時刻軟骨缺損修復區的應力分布，如圖 2 所示。加載初期（前 2 s），在鄰近修復界面的宿主軟骨表層出現 Mises 應力及第一、二主應力（
、
）集中的現象，且
、
均為正值，在鄰近修復界面的宿主軟骨底部出現第三主應力（
）集中的現象且為負值，表明此時
、
起主導作用；當加載到 1 000 s 時，Mises 應力在修復區邊界底部的宿主軟骨上出現應力集中，此時
起主導作用。
圖2
不同加載時刻全層缺損修復區 Mises 應力和主應力（σ1、σ2、σ3）分布云圖
Figure2.
Contours of Mises stress and principal stress (σ1、σ2、σ3) in full-thickness defect repairing area at different loading time
主應力是矢量，可用箭頭表示各點主應力的方向，如圖 3 所示。
沿圓周方向（θ向）并平行于軟骨表面、
從壓頭正下方沿徑向向四周呈放射狀分布（r向）、
近似沿軸向（z向）。
圖3
軟骨缺損修復區主應力（σ1、σ2、σ3）矢量圖
Figure3.
Vector diagrams of principal stress (σ1、σ2、σ3) in defect repairing area of cartilage
2.2 修復深度不同時，鄰近修復界面的宿主軟骨表面節點的應力狀態分析
臨床研究發現，即使使用醫用生物膠粘合 TEC 與宿主軟骨，仍無法保證修復組織與鄰近宿主軟骨愈合良好，且兩者之間普遍存在裂隙。生理載荷作用下，表層軟骨變形最顯著。因此，在鄰近修復界面的宿主軟骨表層選取單元 36 693 上的節點 2 493（圖 1 中位置一）進行應力狀態分析，如圖 4 所示。研究三維應力狀態下材料的變形類型和尺寸變化的規律，探索修復界面開裂的原因。首先，從數值模擬結果中提取節點 2 493 的三個主應力（
、
、
）值。平均應力（靜水壓力）值為
= （
+
+
）/3，則各應力偏量值為
（i = 1，2，3）。節點 2 493 的各主應力和應力偏量值（
、
、
、
、
、
），如表 2 所示。
圖4
在鄰近缺損修復界面的表層宿主軟骨內選取單元和節點
Figure4.
The element and node selection of the surface layer host cartilage near the defect repairing interface
表2
不同修復深度時節點 2 493 的主應力和應力偏量
Table2.
The values of principal stress and stress deviation of node 2 493 for different repairing depth
已知在三維空間中一點的主應力為
、
、
，過該點斜面法線的方向余弦分別為l、m、n。該斜面上的全應力S、正應力 σN 及剪應力 τ，各應力值分別如式（4）、（5）、（6）所示[12]：
 |
 |
 |
設正應力在 3 個主方向上的分量分別為
、
、
，規定當正應力沿截面的外法線方向為正、沿內法線方向為負，則可得到正應力為正、負的描述方程，如式（7）所示：
 |
設剪應力在 3 個主方向上的分量分別為
、
、
，可得剪應力的描述方程，如式（8）所示：
 |
列維-米塞斯（Levy-Mises）應力應變關系，如式（9）所示[15]：
 |
式中，
為與材料性能相關的瞬時非負標量，分母為應力偏量分量，分子為相應的應變增量。公式（9）表明，一點的瞬時尺寸變化由應力偏量決定。若
，則始終有
、
，即
、
。
的正（負）由
的正（負），即
接近
和
的程度決定。當
時，
，材料為壓縮類變形；當
時，
，材料為拉伸類變形。故
的正（負）能直接判斷物體的變形類型。此外，剪應力圖形和應力偏量中的正應力圖形與物體的變形類型有內在聯系，且正應力圖形能夠判斷材料各方向上尺寸的變化趨勢[14]。
將表 2 中的數據代入式（7），可繪制節點 2 493 在應力偏量中的正應力圖形，如圖 5 所示。分別用黃（紅）色表示正應力的正（負）區域，各方向上黃色代表材料尺寸增加、紅色代表尺寸減小。可知軟骨表層缺損修復和無缺損時，節點 2 493 在三維應力狀態下發生壓縮類變形；中間層、深層、全層缺損修復時，節點 2 493 在三維應力狀態下發生拉伸類變形。
圖5
不同修復深度時，節點 2 493 在應力偏量中的正應力三維圖形
Figure5.
Three dimension normal stress graphics of the node 2 493 in stress deviation for different repairing depths
將表 2 中的數據代入式（8），可繪制節點 2 493 的剪應力圖，如圖 6 所示。軟骨表層缺損修復和無缺損時，節點 2 493 的剪應力圖形開口向上，節點發生壓縮類變形；軟骨中間層、深層和全層缺損修復時，節點的剪應力圖形開口向兩側，節點發生拉伸類變形。
上述研究結果表明，表層缺損修復時節點發生壓縮類變形，軟骨徑向尺寸減小，修復界面不易開裂；非表層缺損修復時節點發生拉伸類變形，軟骨徑向尺寸增加，修復界面有開裂的趨勢。
圖6
不同深度修復時，節點 2 493 的剪應力三維圖形
Figure6.
Three dimension shear stress graphics of node 2 493 for different repairing depths
2.3 TEC 彈性模量不同時鄰近修復界面各層宿主軟骨的應力狀態分析
圖7
TEC 彈性模量不同時，路徑一上各點的應力偏量隨歸一化深度變化的曲線
Figure7.
Curves of the stress deviator of each point on the path 1 with normalized depth for different TEC elastic?modulus
臨床上，全層缺損可使血流達到損傷區域利于修復[23]。本文通過研究不同彈性模量的 TEC 修復軟骨全層缺損后，鄰近修復界面各層宿主軟骨（圖 1 中路徑一）的應力偏量、變形類型隨深度變化的規律，以探尋合適的 TEC 材料參數。路徑一上各點的應力偏量值（
、
、
）隨深度變化的曲線圖，如圖 7 所示。TEC 彈性模量不同時，
均為正值，
均為負值，
有正、負之分，而
的正（負）反應材料的變形類型。
TEC 的彈性模量為 0.1 MPa 時，路徑一上從軟骨表層起約占軟骨厚度 65%（表層和大部分中間層軟骨）的范圍內點的
為負，即此應力狀態下軟骨發生拉伸類變形，而余下的 35%（深層軟骨）發生壓縮類變形；TEC 的彈性模量為 0.3 MPa 時，路徑一上從軟骨表層起約占軟骨厚度 47%（表層和部分中間層軟骨）的范圍內點的
為負，而余下的 53%（部分中間層和深層軟骨）發生壓縮類變形；TEC 的彈性模量較高（如 0.6 MPa、0.9 MPa）時，路徑一所有點的
均為正值，此時修復區各層宿主軟骨均發生壓縮類變形。由前文可知，當缺損修復區鄰近修復界面的宿主軟骨發生拉伸類變形時，軟骨徑向尺寸增加，該趨勢易造成修復界面開裂，因此可適當增大 TEC 的彈性模量。
3 結論
本文首次通過分析軟骨缺損修復區鄰近修復界面宿主軟骨的應力狀態來判別軟骨的變形類型，探索修復界面開裂的原因。研究發現，軟骨表層缺損修復時，鄰近修復界面的宿主軟骨表面節點發生壓縮類變形，與無缺損軟骨相同；軟骨中間層、深層和全層缺損修復時，宿主軟骨表面節點發生拉伸類變形，表明此時軟骨徑向尺寸增加，易致修復界面開裂。對于軟骨全層缺損修復，TEC 的彈性模量較低（如 0.1 MPa、0.3 MPa）時，修復界面附近的表層和大部分中間層宿主軟骨發生拉伸類變形，發生拉伸類變形的范圍隨 TEC 的彈性模量增大而減小；TEC 的彈性模量較高（如 0.6 MPa、0.9 MPa）時，鄰近修復界面各層宿主軟骨均發生壓縮類變形，利于修復。
綜上所述，臨床上若軟骨為表層缺損，即采用表層缺損修復；若缺損為表層以下（如中間層、深層）可應用手術刀具制作成全層缺損后再進行修復并適當增大 TEC 的彈性模量，從而利于缺損修復。基于此，期待本文研究成果對臨床修復具有一定的參考意義。但由于本文研究建立的模型過于理想化，無法模擬出關節軟骨缺損修復區的真實生物力學環境，因此具有一定的局限性。
引言
關節軟骨是人體關節的重要承力組織，主要功能是減小摩擦和磨損及將關節軟骨表面的載荷均勻地傳遞到軟骨下骨，對活動關節起著非常重要的支撐和保護作用。隨著國民生活水平的提高，高能、高強度運動增多，且人口老齡化日趨嚴重及肥胖人口比例上升，關節軟骨損傷成為骨科臨床常見疾病之一，若不加以有效治療最終會發展成骨性關節炎[1-2]。由于軟骨組織內沒有血供和淋巴回流為其提供營養，損傷后很難自行修復，軟骨損傷后的有效治療成為臨床上亟待解決的難題[3]。近年來，隨著組織工程技術逐漸成熟，使其有希望成為治療軟骨損傷的最有效方法[4-5]。盡管該方法在臨床上取得了一定進展，但大多短期療效較好，長期修復效果具有不確定性。主因之一是目前還無法構建出可媲美天然軟骨生物力學性能的組織工程化軟骨（tissue engineering cartilage，TEC），植入的 TEC 與宿主軟骨間力學性能的差異顯著地改變了修復界面及鄰近宿主軟骨和細胞的力學環境，常出現移植物退化變性、修復界面開裂及宿主軟骨降解等現象[6-7]。因此，鄰近修復界面宿主軟骨的力學狀態是影響修復效果的關鍵因素之一。
關節軟骨是由固相（包括：蛋白多糖、膠原纖維）和液相（水）構成的多孔結構，孔隙中充滿組織液，液相約占軟骨重量的 65% 以上，有的部位可達 80%，具有顯著的非線性粘彈特性。加載初期，液相能承擔高達 90% 的載荷從而保護膠原纖維網免于高應力的破壞，且組織液在壓力梯度作用下流動利于軟骨中營養物質的傳遞和交換。因此，學者們普遍認為固液耦合雙相模型是研究軟骨承載特性的基礎，包括纖維增強雙相模型和橫觀各向同性雙相模型。由于軟骨中纖維的直徑、含量及走向依深度變化，在軟骨的二維有限元模型中可以使用彈簧單元或者桿單元來表征膠原纖維，而纖維的空間分布規律卻沒有相應描述[8]。實際上，表層膠原纖維的走向與軟骨表面平行，如同編織的毯子，從而在軟骨表面形成了一層不透水的纖維網絡。Clark[9]據此提出了軟骨的橫觀各向同性模型，經研究驗證，此模型能夠成功地預測應力集中的位置[10]。隨著關節軟骨橫觀各向同性理論逐漸成熟，Elhamian 等[11]應用實驗和理論相結合的方法提出了關節軟骨隨深度變化的橫觀各向同性微觀力學模型，該模型進一步促進了相關研究的發展。
體外實驗研究表明，應力狀態下軟骨細胞的生物力學環境的改變致軟骨基質分解代謝的變化能顯著地影響關節軟骨的生理功能[12]。孟迪等[8]通過有限元模擬研究了缺損修復形狀對修復區的米塞斯（Mises）應力分布影響的規律。Mises 應力是基于第四強度理論（剪切應變能）的等效應力，幾乎所有的有限元分析結果中默認的都是 Mises 應力。通過 Mises 應力云圖可確定軟骨缺損修復的危險區域，是分析應力分布變化規律的重要指標，但 Mises 應力是標量，沒有力的方向性指示。因此，無法展現修復區復雜的力學環境，鄰近修復界面宿主軟骨的應力狀態可能與修復界面的整合情況密切相關，而相關研究卻鮮見報導，常用空間某一點的三個主應力（
、
、
）表示該點的應力狀態。何祝斌等[13-16]給出了在主坐標系中，空間一點正應力、剪應力的描述方程，并繪制了一點的剪應力三維圖形和偏應力空間中的正應力圖形，從圖形中能判斷不同應力狀態下材料發生變形的類型及尺寸變化的規律，為結構設計或應力、應變分析提供了一種直觀有效的方法。
本文建立了基于橫觀各向同性的關節軟骨固液耦合雙相三維缺損及修復的有限元模型。通過數值模擬和理論分析，對比研究表層、中間層、深層和全層軟骨缺損修復的情況下，鄰近修復界面的宿主軟骨表面層單元上節點的應力狀態、變形類型及材料尺寸的變化規律；并對比分析使用不同彈性模量（0.1、0.3、0.6、0.9 MPa）的 TEC 修復全層軟骨缺損后，鄰近修復界面各層宿主軟骨的應力偏量值（
、
、
）隨深度變化的規律。在三向應力狀態下，預測軟骨缺損修復界面開裂的原因具有重要的理論和臨床實際意義。
1 有限元模型構建及材料參數
?
關節軟骨位于股骨遠端和脛骨近端的關節面上，有研究表明，當軟骨缺損直徑小于 3 mm 時依靠自身的修復能力可部分或全部修復，而直徑大于 3 mm 則不能自行修復[17]。正常人體關節軟骨的厚度為 1～5 mm，由于其形狀不規則，通常采用簡化模型進行研究。本文應用有限元軟件 ABAQUS 6.10（SIMULIA，法國）建立厚 2 mm、半徑 15 mm 的圓柱形脛骨軟骨幾何模型，而圓柱形 TEC 的直徑為 4 mm，用半徑 20 mm 的剛性球模擬股骨軟骨[18]。如圖1所示，為便于施加位移載荷，在球體中心設置一個施加載荷的點；在鄰近修復界面左側的宿主軟骨內取一條路徑（如圖1中路徑一所示）、宿主軟骨表面層取一節點（如圖1中位置一所示）用于后續的應力狀態分析。
圖1
關節軟骨全層缺損修復的非圍限壓縮有限元模型截面 示意圖
Figure1.
The sectional sketch of unconfined compression finite element model for articular cartilage with full-thickness defects repair
膠原纖維是構成關節軟骨三維網狀結構的重要組成成分，能使關節適應各種機械負荷。根據膠原纖維的排列情況，關節軟骨在形態學上可分為表層、中間層和深層。表層占軟骨厚度的 10%～20%，膠原纖維較細密且走向與軟骨表面平行；中間層占軟骨厚度的 40%～60%，膠原纖維較粗且走向與軟骨表面交叉、排列不整齊；深層占軟骨厚度的 30%，纖維粗大且走向與軟骨表面垂直[19]。由于各層膠原纖維和蛋白多糖的含量、體積分數及纖維走向不同，導致關節軟骨的力學性能依深度而變化。
Elhamian 等[11]通過理論和實驗分析得到了基于橫觀各向同性理論的關節軟骨材料參數隨深度變化的曲線圖，其中軟骨模型的空間坐標系為軸向z、徑向r和圓周方向θ。各材料參數均隨深度方向的坐標值z變化：如彈性模量Ez（z）、Er（z）、Eθ（z），泊松比vrθ（z）、vrz（z）、vθz（z），切變模量Grθ（z）、Grz（z）、Gθz（z）。本文應用數學軟件 Mathematica 7.0（Wolfram Research，美國）對從曲線上提取的多組數據進行擬合，得到材料參數隨歸一化深度
值（
，h = 2 mm 為軟骨厚度）變化的非線性關系表達式，如式（1）所示。由于篇幅所限，僅以彈性模量
為例：
|
'/> |
軟骨有限元模型沿厚度共分為 10 層，其中 1～2 層、3～7 層、8～10 層分別為軟骨的表層、中間層、深層。將平分每層厚度的
值代入各非線性關系式可得材料參數，如表 1 所示。
表1
關節軟骨材料參數
Table1.
Material parameters of the articular cartilage
滲透率k是表征流體通過多孔材料能力的參數，其大小取決于材料的孔隙率和液體滲透方向上孔隙的幾何形狀等因素，k與應變相關，如式（2）所示[20]：
|
式中，k0 = 1.743 × 10–15 m4/（N·s），為初始滲透率；M = 7.081 為材料滲透常數；e和e0 分別為當前孔隙率和初始孔隙率。e0 隨深度變化的關系式，如式（3）所示[21]：
|
式中，
為軟骨表層初始孔隙比；αe = 0.413 為材料常數。TEC 的初始孔隙比取為 9[20]。軟骨底部固定且不可滲透，設軟骨側面孔壓為 0，可實現液相流動及滲出[22]。
在球體中心處施加位移載荷。首先，剛性壓頭向下移動與軟骨接觸；之后，施加z向位移載荷，確保壓縮量達到軟骨厚度的 30%；最后，令壓頭保持某一固定值，給定足夠松弛時間分析其力學行為。
2 結果
由于根據 Mises 應力云圖無法得知危險區域受拉應力還是壓應力，也無法確定該區域是主應力還是剪應力占主導，因而無法判別修復界面開裂的原因。因此研究者們更應關注應力在不同區域、不同方向上的形式。通過繪制主應力矢量圖可清晰地呈現物體內部的受力情況，如拉、壓應力的大小及方向。
2.1 關節軟骨缺損修復區的應力分布
關節軟骨具有流固耦合特性，液相在壓力梯度作用下流動持續地改變著修復區的力學環境。本文以彈性模量 0.3 MPa 的 TEC 修復全層缺損為例，數值模擬應力松弛情況下不同加載時刻軟骨缺損修復區的應力分布，如圖 2 所示。加載初期（前 2 s），在鄰近修復界面的宿主軟骨表層出現 Mises 應力及第一、二主應力（
、
）集中的現象，且
、
均為正值，在鄰近修復界面的宿主軟骨底部出現第三主應力（
）集中的現象且為負值，表明此時
、
起主導作用；當加載到 1 000 s 時，Mises 應力在修復區邊界底部的宿主軟骨上出現應力集中，此時
起主導作用。
圖2
不同加載時刻全層缺損修復區 Mises 應力和主應力（σ1、σ2、σ3）分布云圖
Figure2.
Contours of Mises stress and principal stress (σ1、σ2、σ3) in full-thickness defect repairing area at different loading time
主應力是矢量，可用箭頭表示各點主應力的方向，如圖 3 所示。
沿圓周方向（θ向）并平行于軟骨表面、
從壓頭正下方沿徑向向四周呈放射狀分布（r向）、
近似沿軸向（z向）。
圖3
軟骨缺損修復區主應力（σ1、σ2、σ3）矢量圖
Figure3.
Vector diagrams of principal stress (σ1、σ2、σ3) in defect repairing area of cartilage
2.2 修復深度不同時，鄰近修復界面的宿主軟骨表面節點的應力狀態分析
臨床研究發現，即使使用醫用生物膠粘合 TEC 與宿主軟骨，仍無法保證修復組織與鄰近宿主軟骨愈合良好，且兩者之間普遍存在裂隙。生理載荷作用下，表層軟骨變形最顯著。因此，在鄰近修復界面的宿主軟骨表層選取單元 36 693 上的節點 2 493（圖 1 中位置一）進行應力狀態分析，如圖 4 所示。研究三維應力狀態下材料的變形類型和尺寸變化的規律，探索修復界面開裂的原因。首先，從數值模擬結果中提取節點 2 493 的三個主應力（
、
、
）值。平均應力（靜水壓力）值為
= （
+
+
）/3，則各應力偏量值為
（i = 1，2，3）。節點 2 493 的各主應力和應力偏量值（
、
、
、
、
、
），如表 2 所示。
圖4
在鄰近缺損修復界面的表層宿主軟骨內選取單元和節點
Figure4.
The element and node selection of the surface layer host cartilage near the defect repairing interface
表2
不同修復深度時節點 2 493 的主應力和應力偏量
Table2.
The values of principal stress and stress deviation of node 2 493 for different repairing depth
已知在三維空間中一點的主應力為
、
、
，過該點斜面法線的方向余弦分別為l、m、n。該斜面上的全應力S、正應力 σN 及剪應力 τ，各應力值分別如式（4）、（5）、（6）所示[12]：
|
|
|
設正應力在 3 個主方向上的分量分別為
、
、
，規定當正應力沿截面的外法線方向為正、沿內法線方向為負，則可得到正應力為正、負的描述方程，如式（7）所示：
|
設剪應力在 3 個主方向上的分量分別為
、
、
，可得剪應力的描述方程，如式（8）所示：
|
列維-米塞斯（Levy-Mises）應力應變關系，如式（9）所示[15]：
|
式中，
為與材料性能相關的瞬時非負標量，分母為應力偏量分量，分子為相應的應變增量。公式（9）表明，一點的瞬時尺寸變化由應力偏量決定。若
，則始終有
、
，即
、
。
的正（負）由
的正（負），即
接近
和
的程度決定。當
時，
，材料為壓縮類變形；當
時，
，材料為拉伸類變形。故
的正（負）能直接判斷物體的變形類型。此外，剪應力圖形和應力偏量中的正應力圖形與物體的變形類型有內在聯系，且正應力圖形能夠判斷材料各方向上尺寸的變化趨勢[14]。
將表 2 中的數據代入式（7），可繪制節點 2 493 在應力偏量中的正應力圖形，如圖 5 所示。分別用黃（紅）色表示正應力的正（負）區域，各方向上黃色代表材料尺寸增加、紅色代表尺寸減小。可知軟骨表層缺損修復和無缺損時，節點 2 493 在三維應力狀態下發生壓縮類變形；中間層、深層、全層缺損修復時，節點 2 493 在三維應力狀態下發生拉伸類變形。
圖5
不同修復深度時，節點 2 493 在應力偏量中的正應力三維圖形
Figure5.
Three dimension normal stress graphics of the node 2 493 in stress deviation for different repairing depths
將表 2 中的數據代入式（8），可繪制節點 2 493 的剪應力圖，如圖 6 所示。軟骨表層缺損修復和無缺損時，節點 2 493 的剪應力圖形開口向上，節點發生壓縮類變形；軟骨中間層、深層和全層缺損修復時，節點的剪應力圖形開口向兩側，節點發生拉伸類變形。
上述研究結果表明，表層缺損修復時節點發生壓縮類變形，軟骨徑向尺寸減小，修復界面不易開裂；非表層缺損修復時節點發生拉伸類變形，軟骨徑向尺寸增加，修復界面有開裂的趨勢。
圖6
不同深度修復時，節點 2 493 的剪應力三維圖形
Figure6.
Three dimension shear stress graphics of node 2 493 for different repairing depths
2.3 TEC 彈性模量不同時鄰近修復界面各層宿主軟骨的應力狀態分析
圖7
TEC 彈性模量不同時，路徑一上各點的應力偏量隨歸一化深度變化的曲線
Figure7.
Curves of the stress deviator of each point on the path 1 with normalized depth for different TEC elastic?modulus
臨床上，全層缺損可使血流達到損傷區域利于修復[23]。本文通過研究不同彈性模量的 TEC 修復軟骨全層缺損后，鄰近修復界面各層宿主軟骨（圖 1 中路徑一）的應力偏量、變形類型隨深度變化的規律，以探尋合適的 TEC 材料參數。路徑一上各點的應力偏量值（
、
、
）隨深度變化的曲線圖，如圖 7 所示。TEC 彈性模量不同時，
均為正值，
均為負值，
有正、負之分，而
的正（負）反應材料的變形類型。
TEC 的彈性模量為 0.1 MPa 時，路徑一上從軟骨表層起約占軟骨厚度 65%（表層和大部分中間層軟骨）的范圍內點的
為負，即此應力狀態下軟骨發生拉伸類變形，而余下的 35%（深層軟骨）發生壓縮類變形；TEC 的彈性模量為 0.3 MPa 時，路徑一上從軟骨表層起約占軟骨厚度 47%（表層和部分中間層軟骨）的范圍內點的
為負，而余下的 53%（部分中間層和深層軟骨）發生壓縮類變形；TEC 的彈性模量較高（如 0.6 MPa、0.9 MPa）時，路徑一所有點的
均為正值，此時修復區各層宿主軟骨均發生壓縮類變形。由前文可知，當缺損修復區鄰近修復界面的宿主軟骨發生拉伸類變形時，軟骨徑向尺寸增加，該趨勢易造成修復界面開裂，因此可適當增大 TEC 的彈性模量。
3 結論
本文首次通過分析軟骨缺損修復區鄰近修復界面宿主軟骨的應力狀態來判別軟骨的變形類型，探索修復界面開裂的原因。研究發現，軟骨表層缺損修復時，鄰近修復界面的宿主軟骨表面節點發生壓縮類變形，與無缺損軟骨相同；軟骨中間層、深層和全層缺損修復時，宿主軟骨表面節點發生拉伸類變形，表明此時軟骨徑向尺寸增加，易致修復界面開裂。對于軟骨全層缺損修復，TEC 的彈性模量較低（如 0.1 MPa、0.3 MPa）時，修復界面附近的表層和大部分中間層宿主軟骨發生拉伸類變形，發生拉伸類變形的范圍隨 TEC 的彈性模量增大而減小；TEC 的彈性模量較高（如 0.6 MPa、0.9 MPa）時，鄰近修復界面各層宿主軟骨均發生壓縮類變形，利于修復。
綜上所述，臨床上若軟骨為表層缺損，即采用表層缺損修復；若缺損為表層以下（如中間層、深層）可應用手術刀具制作成全層缺損后再進行修復并適當增大 TEC 的彈性模量，從而利于缺損修復。基于此，期待本文研究成果對臨床修復具有一定的參考意義。但由于本文研究建立的模型過于理想化，無法模擬出關節軟骨缺損修復區的真實生物力學環境，因此具有一定的局限性。
