本發(fā)明涉及關(guān)節(jié)軟骨兩相模型的建立方法。
背景技術(shù)：
：關(guān)節(jié)軟骨主要由水和復(fù)合有機(jī)質(zhì)組成。復(fù)合有機(jī)質(zhì)主要包括膠原纖維和蛋白多糖，間隙流體中主要含有水，這些流體的流動(dòng)不僅對軟骨的力學(xué)性能有著重要影響，還與無血管組織的營養(yǎng)物質(zhì)傳輸有著密切關(guān)系。同時(shí)，關(guān)節(jié)軟骨使關(guān)節(jié)具有極低的摩擦系數(shù)，起到良好的潤滑作用。關(guān)節(jié)軟骨的潤滑作用對動(dòng)關(guān)節(jié)運(yùn)動(dòng)和承受載荷的能力具有極其重要的意義。因此，研究關(guān)節(jié)軟骨的力學(xué)特性，認(rèn)識(shí)其在運(yùn)動(dòng)過程中的應(yīng)力、應(yīng)變狀況及其在運(yùn)動(dòng)中各力學(xué)量的變化規(guī)律，對醫(yī)學(xué)研究、臨床診斷和新型生物材料的開發(fā)，具有非常重要的理論意義和實(shí)用價(jià)值。在關(guān)節(jié)軟骨力學(xué)建模的研究初期，人們將關(guān)節(jié)軟骨看作是各相同性線彈性單相材料。在一段時(shí)間里，研究人員接受了粘彈性的單相關(guān)節(jié)軟骨模型，然而，通過這種本構(gòu)關(guān)系得到的仿真結(jié)果總是與軟骨的實(shí)際情況和實(shí)驗(yàn)數(shù)據(jù)有很大的差異。目前，關(guān)節(jié)軟骨被視為兩相模型。但是，這些研究中大都將軟骨的固相看作是線彈性材料，線彈性模型只對小變形有效，針對關(guān)節(jié)軟骨在力的施加過程中會(huì)產(chǎn)生較大變形的問題，將軟骨的固相看作是線彈性材料的模型就會(huì)與實(shí)際情況有很大的差異，甚至失效。技術(shù)實(shí)現(xiàn)要素：本發(fā)明為了解決將軟骨的固相視作線彈性材料的模型針對關(guān)節(jié)軟骨產(chǎn)生較大變形的時(shí)候與實(shí)際情況有很大的差異、甚至失效的問題?；诔瑥椥怨腆w相特性關(guān)節(jié)軟骨兩相模型的建立方法，包括以下步驟：步驟1、幾何模型和網(wǎng)格劃分：以關(guān)節(jié)為研究對象，通過獲得關(guān)節(jié)的CT數(shù)據(jù)，將得到的數(shù)據(jù)以DICOM格式輸出并儲(chǔ)存在計(jì)算機(jī)中；通過MIMICS中圖像分割功能將CT數(shù)據(jù)中關(guān)節(jié)軟骨的區(qū)域從其他組織中分離出來，通過三維重建得到需要研究的人體組織的幾何模型；從DICOM數(shù)據(jù)到生成網(wǎng)格模型并獲取模型的單元編號(hào)和節(jié)點(diǎn)坐標(biāo)作為形變計(jì)算輸入的原始數(shù)據(jù)，獲取過程如圖1所示；步驟2、建立基于混合物理論的關(guān)節(jié)軟骨兩相模型及基于v-p變量的控制方程，包括以下步驟：步驟2.1、建立關(guān)節(jié)軟骨的兩相模型：將關(guān)節(jié)軟骨視為由超彈性固體和理想流體組成的兩相介質(zhì)混合物，且兩相具有獨(dú)立的運(yùn)動(dòng)規(guī)律，用s表示固體相，f表示液體相，那么φs表示固體相體積分?jǐn)?shù)，φf表示液體相體積分?jǐn)?shù)；在初始構(gòu)型中，兩相的初始體積分?jǐn)?shù)為φs0、φf0，兩相的初始密度為且均勻分布；由飽和下的體積分?jǐn)?shù)式得到：φs+φf＝1(1)質(zhì)量平衡方程為：▿·(φsvs+φfvf)=0---(2)]]>其中，表示矢量微分算子符號(hào)；vs是固體相速度，vf是液體相速度；在擬靜態(tài)問題中，動(dòng)量平衡方程轉(zhuǎn)化為靜力平衡方程，加速度a等于零；此研究中忽略外部體積力的作用，在兩相混合物中，由間隙液體流過多孔固體所產(chǎn)生的摩擦阻力形成的動(dòng)量交換為Ps、Pf，且有Ps＝-Pf＝K′(vf-vs)(3)式中，K′為擴(kuò)散阻力系數(shù)；K′=(φf)2κ---(4)]]>式中，κ為軟骨滲透率；固體相動(dòng)量為Ps，液體相動(dòng)量Pf；令σs和σf分別為固相和液相的柯西應(yīng)力張量，關(guān)節(jié)軟骨在當(dāng)前構(gòu)形的動(dòng)量平衡方程為：▿·σs+Ps=0---(5)]]>▿·σf+Pf=0---(6)]]>由軟骨的不可壓縮性和力學(xué)特性得到關(guān)節(jié)軟骨彈性固體相和非粘性流體相的力學(xué)本構(gòu)方程為：σs＝-φspI+σe(7)σf＝-φfpI(8)其中，I是單位矩陣，p為壓力；σe為固體相的有效應(yīng)力也叫做與固相的變形量相一致的彈性應(yīng)力張量；經(jīng)過以上推導(dǎo)得出由質(zhì)量平衡方程(2)、動(dòng)量平衡方程(5)(6)和兩相混合物的本構(gòu)方程(7)(8)組成了關(guān)節(jié)軟骨的兩相模型；根據(jù)第二Piola-Kirchhoff應(yīng)力張量的定義可知第二Piola-Kirchhoff應(yīng)力張量S與彈性應(yīng)力張量σe的關(guān)系為：σe＝J-1FSFT＝FSFT(9)其中，F(xiàn)為變形梯度矩陣；J為彈性體積比，是變形梯度矩陣F的雅可比行列式；S為第二Piola-Kirchhoff應(yīng)力張量，該參數(shù)在步驟2.3中求解，通過S的聯(lián)結(jié)而得到關(guān)節(jié)軟骨的超彈性兩相介質(zhì)模型；在有限變形中，考慮到物體受力前后空間位置的改變，令t＝0時(shí)刻物體在空間所占據(jù)的區(qū)域V0為初始構(gòu)形；在當(dāng)前時(shí)刻t，物體所占據(jù)的空間區(qū)域V表示為當(dāng)前構(gòu)形；以初始構(gòu)形為參考構(gòu)型，利用物質(zhì)描述法可得出體元、面元在有限變形中的變化，即當(dāng)前構(gòu)形中的單元物質(zhì)的體積與參考構(gòu)型中的體積滿足以下關(guān)系：dV＝JdV0(10)步驟2.2、建立超彈性模型：關(guān)節(jié)軟骨的超彈性模型由Helmholtz應(yīng)變能函數(shù)來表示，即：Ψ=α0eα1(I1-3)+α2(I2-3)I3β---(14)]]>α0、α1、α2、β是材料的力學(xué)參數(shù)，來源于數(shù)學(xué)模型與實(shí)驗(yàn)數(shù)據(jù)的擬合結(jié)果；I1、I2、I3分別是C的第一、第二、第三主不變量；對于不可壓縮材料而言，I3＝det(C)＝1，簡化后的超彈性模型為：Ψ=α0eα1(I1-3)+α2(I2-3)---(15)]]>步驟2.3、求解第二Piola-Kirchhoff應(yīng)力張量S：關(guān)節(jié)軟骨的超彈性固體相，應(yīng)力應(yīng)變不再滿足線性關(guān)系，彈性矩陣與應(yīng)變張量有關(guān)，第二Piola-Kirchhoff應(yīng)力張量S和Green應(yīng)變張量E都是以初始構(gòu)形為參考構(gòu)型的應(yīng)力和應(yīng)變張量，在一個(gè)瞬時(shí)變形過程中，本構(gòu)關(guān)系為：dS＝DdE；Green應(yīng)變張量E為：E=E11E12E13E21E22E23E31E32E33---(16)]]>S=∂Ψ∂E=Ψ2α1+4α2(1+E22+E33)-8α2E12-8α2E13-8α2E122α1+4α2(1+E11+E33)-8α2E23-8α2E13-8α2E232α1+4α2(1+E11+E22)---(17)]]>彈性張量D表示為：D=D1111D1122D1133D1123D1113D1112D2211D2222D2233D2223D2213D2212D3311D3322D3333D3323D3313D3312D2311D2322D2333D2323D2313D2312D1311D1322D1333D1323D1313D1312D1211D1222D1233D1223D1213D1212---(18)]]>式中，D中的元素Si″j″表示S中的第i″行第j″列的元素，Ek″l″表示E中的第k″行第l″列的元素；令：2α1+4α2(1+E22+E33)＝G2α1+4α2(1+E11+E22)＝H2α1+4α2(1+E11+E33)＝B得到對稱矩陣D，該矩陣中的對稱部分用sym表示：D=ΨG24α2+GB4α2+GH-8α2E23G-8α2E13G-8α2E12GB24α2+HB-8α2E23B-8α2E13B-8α2E12BH2-8α2E23H-8α2E13H-8α2E12H64α22E232-8α264α22E13E2364α22E12E23sym64α22E132-8α264α22E13E1264α22E122-8α2---(19)]]>步驟2.4、建立基于v-p變量的控制方程本研究中選取速度v、壓力p為未知量建立基于v-p變量的有限元方程，進(jìn)而表達(dá)關(guān)節(jié)軟骨基于v-p變量的控制方程；為了消除流體相的速度，先將式(3)和式(8)帶到式(6)中，通過式(4)得到：φf(vf-vs)=-κ▿p---(20)]]>在飽和條件式(1)下，將式(20)帶入式(2)中，得到只包含固體相速度的關(guān)系式：▿·(vs-κ▿p)=0---(21)]]>接下來由動(dòng)量守恒方程(5)(6)和本構(gòu)方程(7)(8)聯(lián)立可得：▿·(σe-pI)=0---(22)]]>由式(21)和式(22)就構(gòu)成了基于v-p變量的控制方程；步驟3、建立關(guān)節(jié)軟骨力學(xué)平衡方程并進(jìn)行有限元計(jì)算：步驟3.1、建立邊界條件：流體通量定義為則基于v-p變量的控制方程的邊界條件為：us(x,t)=u‾s,x∈Γus---(23)]]>σ(x,t)·n*=T‾,x∈ΓT---(24)]]>p(x,t)=p‾,x∈Γp---(25)]]>-κ▿p·n*=Q‾,x∈ΓQ---(26)]]>式中，x指在某一時(shí)刻t時(shí)當(dāng)前構(gòu)型中的空間點(diǎn)；為已知的固相位移，為總的載荷力，為已知的壓力，n*是邊界Γ的外法線向量；Γp和ΓQ分別代表液體相體積域Ωf的邊界，和ΓT分別代表固體相體積域Ωs的邊界，且邊界滿足條件：步驟3.2、構(gòu)造形函數(shù)：定義有限元四面體單元的形函數(shù)為：N＝[Ni′·I3×3Nj′·I3×3Nm′·I3×3Np′·I3×3](44)其中，Ni′=16Ve(ai′+bi′x′+ci′y′+di′z′)---(45)]]>式中，Ve是四面體單元的體積，ai′、bi′、ci′、di′為與Ni′對應(yīng)的系數(shù)；Nj′、Nm′、Np′與Ni′類似，有Nj′=16Ve(aj′+bj′x′+cj′y′+dj′z′)Nm′=16Ve(am′+bm′x′+cm′y′+dm′z′)Np′=16Ve(ap′+bp′x′+cp′y′+dp′z′)---(46)]]>步驟3.3、建立關(guān)節(jié)軟骨力學(xué)平衡方程：基于有限元求解的單元類型為四面體單元，構(gòu)造插值函數(shù)；通過基于v-p變量的控制方程(21)、(22)以及邊界條件(23)-(26)，采用伽遼金加權(quán)余量法建立關(guān)節(jié)軟骨力學(xué)平衡方程；步驟3.4、求解關(guān)節(jié)軟骨力學(xué)平衡方程，實(shí)現(xiàn)關(guān)節(jié)軟骨兩相模型的仿真。本發(fā)明具有以下優(yōu)點(diǎn)：1、目前的關(guān)節(jié)軟骨被視為兩項(xiàng)模型。但是，這些研究中大都將軟骨的固相看作是線彈性材料，線彈性模型只對小變形有效。針對關(guān)節(jié)軟骨在力的施加過程中會(huì)產(chǎn)生較大變形的問題，本發(fā)明提出基于超彈性固體相特性建立關(guān)節(jié)軟骨兩相模型，實(shí)現(xiàn)了關(guān)節(jié)軟骨兩相模型的建立和仿真。2、本發(fā)明針對關(guān)節(jié)軟骨的結(jié)構(gòu)特點(diǎn)，將關(guān)節(jié)軟骨描述為由不可壓縮的非粘性液體相和具有超彈性、橫觀各向同性的固體相所組成的二相混合物，并且在變形的條件下還會(huì)引起軟骨組織滲透性的改變。3、依據(jù)Helmholtz應(yīng)變能函數(shù)定義固體相的超彈性特性，液體相定義為理想流體，該模型能夠描述關(guān)節(jié)軟骨的非線性、不可壓縮性和滲透性。通過伽遼金加權(quán)余量法獲得基于v-p變量的有限元系統(tǒng)平衡方程，通過有限差分法對平衡方程進(jìn)行數(shù)值計(jì)算。4、本發(fā)明針對關(guān)節(jié)軟骨在力的施加過程中產(chǎn)生較大變形的情況，本發(fā)明建立的關(guān)節(jié)軟骨兩相模型也能夠準(zhǔn)確的反應(yīng)實(shí)際關(guān)節(jié)軟骨的情況，與實(shí)際關(guān)節(jié)軟骨的情況的相比，誤差低于13％。附圖說明圖1為從DICOM數(shù)據(jù)到生成網(wǎng)格模型并獲取模型的單元編號(hào)和節(jié)點(diǎn)坐標(biāo)的流程示意圖。圖2為單軸壓縮的應(yīng)力-速度實(shí)驗(yàn)曲線和有限元仿真結(jié)果對比圖。具體實(shí)施方式具體實(shí)施方式一：基于超彈性固體相特性關(guān)節(jié)軟骨兩相模型的建立方法，包括以下步驟：步驟1、幾何模型和網(wǎng)格劃分：以關(guān)節(jié)為研究對象，通過獲得關(guān)節(jié)的CT數(shù)據(jù)，將得到的數(shù)據(jù)以DICOM格式輸出并儲(chǔ)存在計(jì)算機(jī)中；通過MIMICS中圖像分割功能將CT數(shù)據(jù)中關(guān)節(jié)軟骨的區(qū)域從其他組織中分離出來，通過三維重建得到需要研究的人體組織的幾何模型；從DICOM數(shù)據(jù)到生成網(wǎng)格模型并獲取模型的單元編號(hào)和節(jié)點(diǎn)坐標(biāo)作為形變計(jì)算輸入的原始數(shù)據(jù)，獲取過程如圖1所示；步驟2、建立基于混合物理論的關(guān)節(jié)軟骨兩相模型及基于v-p變量的控制方程，包括以下步驟：步驟2.1、建立關(guān)節(jié)軟骨的兩相模型：將關(guān)節(jié)軟骨視為由超彈性固體和理想流體組成的兩相介質(zhì)混合物，且兩相具有獨(dú)立的運(yùn)動(dòng)規(guī)律，用s表示固體相，f表示液體相，那么φs表示固體相體積分?jǐn)?shù)，φf表示液體相體積分?jǐn)?shù)；在初始構(gòu)型中，兩相的初始體積分?jǐn)?shù)為φs0、φf0，兩相的初始密度為且均勻分布；由飽和下的體積分?jǐn)?shù)式得到：φs+φf＝1(1)質(zhì)量平衡方程為：▿·(φsvs+φfvf)=0---(2)]]>其中，表示矢量微分算子符號(hào)；vs是固體相速度，vf是液體相速度；在擬靜態(tài)問題中，動(dòng)量平衡方程轉(zhuǎn)化為靜力平衡方程，加速度a等于零；此研究中忽略外部體積力的作用，在兩相混合物中，由間隙液體流過多孔固體所產(chǎn)生的摩擦阻力形成的動(dòng)量交換為Ps、Pf，且有Ps＝-Pf＝K′(vf-vs)(3)式中，K′為擴(kuò)散阻力系數(shù)；K′=(φf)2κ---(4)]]>式中，κ為軟骨滲透率；固體相動(dòng)量為Ps，液體相動(dòng)量Pf；令σs和σf分別為固相和液相的柯西應(yīng)力張量，關(guān)節(jié)軟骨在當(dāng)前構(gòu)形的動(dòng)量平衡方程為：▿·σs+Ps=0---(5)]]>▿·σf+Pf=0---(6)]]>由軟骨的不可壓縮性和力學(xué)特性得到關(guān)節(jié)軟骨彈性固體相和非粘性流體相的力學(xué)本構(gòu)方程為：σs＝-φspI+σe(7)σf＝-φfpI(8)其中，I是單位矩陣，p為壓力；σe為固體相的有效應(yīng)力也叫做與固相的變形量相一致的彈性應(yīng)力張量；經(jīng)過以上推導(dǎo)得出由質(zhì)量平衡方程(2)、動(dòng)量平衡方程(5)(6)和兩相混合物的本構(gòu)方程(7)(8)組成了關(guān)節(jié)軟骨的兩相模型；根據(jù)第二Piola-Kirchhoff應(yīng)力張量的定義可知第二Piola-Kirchhoff應(yīng)力張量S與彈性應(yīng)力張量σe的關(guān)系為：σe＝J-1FSFT＝FSFT(9)其中，F(xiàn)為變形梯度矩陣；J為彈性體積比，是變形梯度矩陣F的雅可比行列式；S為第二Piola-Kirchhoff應(yīng)力張量，該參數(shù)在2.3中求解，通過S的聯(lián)結(jié)而得到關(guān)節(jié)軟骨的超彈性兩相介質(zhì)模型；在有限變形中，考慮到物體受力前后空間位置的改變，令t＝0時(shí)刻物體在空間所占據(jù)的區(qū)域V0為初始構(gòu)形；在當(dāng)前時(shí)刻t，物體所占據(jù)的空間區(qū)域V表示為當(dāng)前構(gòu)形；以初始構(gòu)形為參考構(gòu)型，利用物質(zhì)描述法可得出體元、面元在有限變形中的變化，即當(dāng)前構(gòu)形中的單元物質(zhì)的體積與參考構(gòu)型中的體積滿足以下關(guān)系：dV＝JdV0(10)步驟2.2、建立超彈性模型：關(guān)節(jié)軟骨的超彈性模型由Helmholtz應(yīng)變能函數(shù)來表示，即：Ψ=α0eα1(I1-3)+α2(I2-3)I3β---(14)]]>α0、α1、α2、β是材料的力學(xué)參數(shù)，來源于數(shù)學(xué)模型與實(shí)驗(yàn)數(shù)據(jù)的擬合結(jié)果；I1、I2、I3分別是C的第一、第二、第三主不變量；對于不可壓縮材料而言，I3＝det(C)＝1，簡化后的超彈性模型為：Ψ=α0eα1(I1-3)+α2(I2-3)---(15)]]>步驟2.3、求解第二Piola-Kirchhoff應(yīng)力張量S：關(guān)節(jié)軟骨的超彈性固體相，應(yīng)力應(yīng)變不再滿足線性關(guān)系，彈性矩陣與應(yīng)變張量有關(guān)，第二Piola-Kirchhoff應(yīng)力張量S和Green應(yīng)變張量E都是以初始構(gòu)形為參考構(gòu)型的應(yīng)力和應(yīng)變張量，在一個(gè)瞬時(shí)變形過程中，本構(gòu)關(guān)系為：dS＝DdE；Green應(yīng)變張量E為：E=E11E12E13E21E22E23E31E32E33---(16)]]>S=∂Ψ∂E=Ψ2α1+4α2(1+E22+E33)-8α2E12-8α2E13-8α2E122α1+4α2(1+E11+E33)-8α2E23-8α2E13-8α2E232α1+4α2(1+E11+E22)---(17)]]>彈性張量D表示為：D=D1111D1122D1133D1123D1113D1112D2211D2222D2233D2223D2213D2212D3311D3322D3333D3323D3313D3312D2311D2322D2333D2323D2313D2312D1311D1322D1333D1323D1313D1312D1211D1222D1233D1223D1213D1212---(18)]]>式中，D中的元素Si″j″表示S中的第i″行第j″列的元素，Ek″l″表示E中的第k″行第l″列的元素；令：2α1+4α2(1+E22+E33)＝G2α1+4α2(1+E11+E22)＝H2α1+4α2(1+E11+E33)＝B得到對稱矩陣D，該矩陣中的對稱部分用sym表示：D=ΨG24α2+GB4α2+GH-8α2E23G-8α2E13G-8α2E12GB24α2+HB-8α2E23B-8α2E13B-8α2E12BH2-8α2E23H-8α2E13H-8α2E12H64α22E232-8α264α22E13E2364α22E12E23sym64α22E132-8α264α22E13E1264α22E122-8α2---(19)]]>步驟2.4、建立基于v-p變量的控制方程本研究中選取速度v、壓力p為未知量建立基于v-p變量的有限元方程，進(jìn)而表達(dá)關(guān)節(jié)軟骨基于v-p變量的控制方程；為了消除流體相的速度，先將式(3)和式(8)帶到式(6)中，通過式(4)得到：φf(vf-vs)=-κ▿p---(20)]]>在飽和條件式(1)下，將式(20)帶入式(2)中，得到只包含固體相速度的關(guān)系式：▿·(vs-κ▿p)=0---(21)]]>接下來由動(dòng)量守恒方程(5)(6)和本構(gòu)方程(7)(8)聯(lián)立可得：▿·(σe-pI)=0---(22)]]>由式(21)和式(22)就構(gòu)成了基于v-p變量的控制方程；步驟3、建立關(guān)節(jié)軟骨力學(xué)平衡方程并進(jìn)行有限元計(jì)算：步驟3.1、建立邊界條件：流體通量定義為則基于v-p變量的控制方程的邊界條件為：us(x,t)=u‾s,x∈Γus---(23)]]>σ(x,t)·n*=T‾,x∈ΓT---(24)]]>p(x,t)=p‾,x∈Γp---(25)]]>-κ▿p·n*=Q‾,x∈ΓQ---(26)]]>式中，x指在某一時(shí)刻t時(shí)當(dāng)前構(gòu)型中的空間點(diǎn)；為已知的固相位移，為總的載荷力，為已知的壓力，n*是邊界Γ的外法線向量；Γp和ΓQ分別代表液體相體積域Ωf的邊界，和ΓT分別代表固體相體積域Ωs的邊界，且邊界滿足條件：步驟3.2、構(gòu)造形函數(shù)：定義有限元四面體單元的形函數(shù)為：N＝[Ni′·I3×3Nj′·I3×3Nm′·I3×3Np′·I3×3](44)其中，Ni′=16Ve(ai′+bi′x′+ci′y′+di′z′)---(45)]]>式中，Ve是四面體單元的體積，ai′、bi′、ci′、di′為與Ni′對應(yīng)的系數(shù)；Nj′、Nm′、Np′與Ni′類似，有Nj′=16Ve(aj′+bj′x′+cj′y′+dj′z′)Nm′=16Ve(am′+bm′x′+cm′y′+dm′z′)Np′=16Ve(ap′+bp′x′+cp′y′+dp′z′)---(46)]]>步驟3.3、建立關(guān)節(jié)軟骨力學(xué)平衡方程：基于有限元求解的單元類型為四面體單元，構(gòu)造插值函數(shù)；通過基于v-p變量的控制方程(21)、(22)以及邊界條件(23)-(26)，采用伽遼金加權(quán)余量法建立關(guān)節(jié)軟骨力學(xué)平衡方程；步驟3.4、求解關(guān)節(jié)軟骨力學(xué)平衡方程，實(shí)現(xiàn)關(guān)節(jié)軟骨兩相模型的仿真。具體實(shí)施方式二：本實(shí)施方式所述步驟3.3建立關(guān)節(jié)軟骨力學(xué)平衡方程的具體過程如下：有限元求解的單元類型為四面體單元，選擇形函數(shù)作為微分方程和邊界條件的權(quán)函數(shù)，令加權(quán)余量積分表達(dá)式中的權(quán)函數(shù)分別為NJ和-NJ，得到相應(yīng)的加權(quán)余量表達(dá)式，通過高斯散度定理和總牽引力的定義得到簡化后的加權(quán)余量積分表達(dá)式為：∫Ω{▿NJT:σe-NJT▿·NJT-κ(▿p·▿NJT)}dΩ=∫ΓTNJT·T‾dΓ+∫ΓQNJT·Q‾dΓ---(29)]]>將問題域離散化成單元形式，對速度和壓力進(jìn)行插值：v＝Nve(30)p＝Npe(31)其中，νe代表四面體單元速度，pe代表四面體單元壓力；則在一個(gè)有限元四面體單元網(wǎng)格上的加權(quán)余量表達(dá)式的矩陣形式為：enen{0GeTGeHevepe+Me0}=enenTQ---(32)]]>式中，en為一個(gè)單位向量，由于en是非零常向量，因此可以直接消掉；Ge=∫Ω0-NJT▿·NdΩ0,GeT=∫Ω0-(▿·N)NJTdΩ0,He=∫Ω0-▿·NJTκ▿·NdΩ0,]]>Me=∫Ω0tr[(▿0NJ)TFS]dΩ0,T=∫ΓT0NJT·T‾JΓdΓ0,Q=∫ΓQ0NJT·Q‾JΓdΓ0;]]>▿NJT:FSFT=tr[(▿0NJ)TFS],▿0=▿·F;]]>令：Y=0GeTGeHe,v=vepe,M=Me0,f=TQ---(33)]]>得到最終的關(guān)節(jié)軟骨力學(xué)平衡方程：Yv+M＝f(34)式中，Y是一個(gè)容量矩陣，M是在固體相中與彈性應(yīng)力張量相關(guān)的節(jié)點(diǎn)力的非線性彈性向量，f是外加載荷力向量。其它步驟及參數(shù)與具體實(shí)施方式一相同。具體實(shí)施方式三：本實(shí)施方式所述步驟3.4所述的求解關(guān)節(jié)軟骨力學(xué)平衡方程的具體過程如下：通過有限差分法對關(guān)節(jié)軟骨力學(xué)平衡方程(34)進(jìn)行數(shù)值計(jì)算：首先在時(shí)域上采用增量法，將時(shí)間t離散化為若干個(gè)時(shí)間點(diǎn)，即：0＝t0<t1<…<tn<tn+1<…<tN＝t(35)時(shí)間增量表示為：△tn+1＝tn+1-tn(36)采用完全拉格朗日法，令初始時(shí)刻t0＝0的構(gòu)形作為參考構(gòu)形，并且該參考構(gòu)形不隨時(shí)間的變化而發(fā)生變化；在tn+1時(shí)刻式(34)寫成：Yn+1vn+1+Mn+1＝fn+1(37)通過Newton-Raphson方法，M的線性化形式由以下遞推形式獲得：Mn+1i=Mn+1i0,Mn+1i=Mn+1i-1+Kn+1i-1·un+1i---(38)]]>式中，Kn+1是固相的切線剛度矩陣，Kn+1是的變分；i代表迭代次數(shù)；由梯形法確定迭代位移的遞推關(guān)系為：Δun+1i=un+1i-un+1i-1=ωΔtΔvn+1i---(39)]]>式中，u為節(jié)點(diǎn)位移，為tn+1時(shí)刻第i次迭代的節(jié)點(diǎn)位移；ω是指定的時(shí)間積分參數(shù)，0≤ω≤1，且當(dāng)ω≥1/2時(shí)，為隱式積分，方程無條件穩(wěn)定；速度的增量表示為：vn+1i=vn+1i-1+Δvn+1i---(40)]]>在速度的第i次迭代中，方程(34)可寫成：Yn+1i-1vn+1i+Mn+1i=fn+1i-1---(41)]]>由(38)(39)(40)(41)能夠得到：(Yn+1i-1+ωΔtKn+1i-1)Δvn+1i=fn+1i-1-Yn+1i-1vn+1i-1-Mn+1i-1---(42)]]>在△tn+1時(shí)間步長內(nèi)，對式(42)反復(fù)進(jìn)行迭代，直到速度增量滿足收斂準(zhǔn)則；在一個(gè)時(shí)間步長內(nèi)，迭代的收斂準(zhǔn)則為：|Δvi||vi|<TOL,TOL=10-4---(43)]]>在當(dāng)前步長內(nèi)迭代收斂發(fā)生時(shí)，將得到的解作為初始值更新到下一個(gè)時(shí)間步長內(nèi)，并且進(jìn)行重復(fù)迭代直到收斂再次發(fā)生；反復(fù)迭代并進(jìn)行數(shù)值計(jì)算，最終實(shí)現(xiàn)關(guān)節(jié)軟骨兩相模型的仿真。反復(fù)迭代并進(jìn)行數(shù)值計(jì)算的過程在VisualStudio平臺(tái)C語言環(huán)境下編程實(shí)現(xiàn)。其它步驟及參數(shù)與具體實(shí)施方式一或二相同。具體實(shí)施方式四：本實(shí)施方式步驟3.4中所述的反復(fù)迭代并進(jìn)行數(shù)值計(jì)算的迭代求解流程如下：(a)在△t＝tn+1-tn增量上定義已知初始速度vn和初始位移un，且有vn+10=vn,un+10=un+Δtvn]]>(b)設(shè)置迭代次數(shù)：i＝1(c)使用來計(jì)算方程(42)中的矩陣和其余未知量，解算這個(gè)線性代數(shù)方程得到第一個(gè)速度增量(Yn+10+ωΔtKn+10)Δvn+11=fn+10-Yn+10vn+10-Mn+10]]>(d)結(jié)束第一次迭代，更新位移和速度向量：vn+11=vn+10+Δvn+11,un+11=un+10+ωΔtΔvn+11]]>(e)更新迭代次數(shù)到i(f)使用更新矩陣和向量，解線性方程得到速度增量(Yn+1i-1+ωΔtKn+1i-1)Δvn+1i=fn+1i-1-Yn+1i-1vn+1i-1-Mn+1i-1]]>(g)結(jié)束當(dāng)前迭代，更新速度和位移向量：vn+1i=vn+1i-1+Δvn+1i,un+1i=un+1i-1+ωΔtΔvn+1i]]>回到步驟(e)繼續(xù)進(jìn)行新的迭代，直到速度的增量滿足收斂準(zhǔn)則。其它步驟及參數(shù)與具體實(shí)施方式一至三之一相同。具體實(shí)施方式五：本實(shí)施方式步驟3.2中所述的Ve、ai′、bi′、ci′、di′的具體形式如下Ve=161xi′′yi′′zi′′1xj′′yj′′zj′′1xm′′ym′′zm′′1xp′′yp′′zp′′---(47)]]>ai=xj′′yj′′zj′′xm′′ym′′zm′′xp′′yp′′zp′′,bi=-1yj′′zj′′1ym′′zm′′1yp′′zp′′ci=-xj′′1zj′′xm′′1zm′′xp′′1zp′′,di=-xj′′yj′′1xm′′ym′′1xp′′yp′′1---(48)]]>式中，(x′i′，y′i′，z′i′)、(x′j′，y′j′，z′j′)、(x′m′，y′m′，z′m′)、(x′p′，y′p′，z′p′)分別為四面體單元的四個(gè)頂點(diǎn)坐標(biāo)。其它步驟及參數(shù)與具體實(shí)施方式一至四之一相同。具體實(shí)施方式六：本實(shí)施方式步驟2.1所述的彈性體積比J＝1，由于關(guān)節(jié)軟骨被視為不可壓縮成型材料，所以J＝1。其它步驟及參數(shù)與具體實(shí)施方式一至五之一相同。具體實(shí)施方式七：本實(shí)施方式步驟2.1中所述的軟骨滲透率其中L是常數(shù)，下標(biāo)“0”代表未變形前的初始配置。關(guān)節(jié)軟骨不僅具有超彈性，而且在變形的條件下還會(huì)引起組織滲透性的改變。其它步驟及參數(shù)與具體實(shí)施方式一至六之一相同。具體實(shí)施方式八：本實(shí)施方式所述中的φs0＝0.2，κ0＝2.519×10-15m4/(Ns)，L＝0.0848。其它步驟及參數(shù)與具體實(shí)施方式一至七之一相同。具體實(shí)施方式九：本實(shí)施方式步驟2.2中所述的α0＝0.1084Mpa，α1＝0.592Mpa，α2＝0.0846Mpa。其它步驟及參數(shù)與具體實(shí)施方式一至八之一相同。本發(fā)明針對關(guān)節(jié)軟骨在力的施加過程中產(chǎn)生較大變形的情況，本發(fā)明建立的關(guān)節(jié)軟骨兩相模型也能夠準(zhǔn)確的反應(yīng)實(shí)際關(guān)節(jié)軟骨的情況。如圖2所示，通過與相同條件下的實(shí)驗(yàn)結(jié)果對比，有限元分析結(jié)果在數(shù)值和趨勢上都與實(shí)驗(yàn)數(shù)據(jù)曲線吻合，圖中“實(shí)驗(yàn)結(jié)果”為實(shí)際關(guān)節(jié)軟骨的實(shí)驗(yàn)數(shù)據(jù)，“仿真結(jié)果”為本發(fā)明的建立的基于超彈性固體相特性關(guān)節(jié)軟骨兩相模型的仿真數(shù)據(jù)。通過仿真結(jié)果與實(shí)驗(yàn)曲線的對比發(fā)現(xiàn)：基于超彈性固體相特性關(guān)節(jié)軟骨兩相模型能夠精確描述關(guān)節(jié)軟骨的力學(xué)特性，與實(shí)際關(guān)節(jié)軟骨的情況的相比，誤差低于13％，驗(yàn)證了模型的精確性和有限元程序的有效性。當(dāng)前第1頁1 2 3