組織液-纖維環(huán)流固耦合的椎間盤軟組織損傷演變預(yù)測方法技術(shù)

本發(fā)明專利技術(shù)公開了一種組織液

【技術(shù)實現(xiàn)步驟摘要】
組織液
?
纖維環(huán)流固耦合的椎間盤軟組織損傷演變預(yù)測方法


[0001]本專利技術(shù)屬于骨科生物力學(xué)領(lǐng)域，涉及考慮椎間盤組織液機械作用的基礎(chǔ)上的組織液
?
纖維環(huán)流固耦合的椎間盤軟組織損傷演變預(yù)測方法。

技術(shù)介紹

[0002]椎間盤由中央髓核、纖維環(huán)和上下軟骨終板組成，可以吸收緩沖腰椎活動時的機械荷載。其中纖維環(huán)繞髓核連接上下軟骨終板，主要由蛋白聚糖和多層的定向膠原纖維網(wǎng)格構(gòu)成，為椎間盤提供支撐并實現(xiàn)高拉伸強度，是椎間盤的主要承重部件。
[0003]隨著年齡的增長，椎間盤內(nèi)部會發(fā)生微結(jié)構(gòu)和化學(xué)成份的變化，這些變化會影響其承受日常載荷的能力。椎間盤發(fā)生損傷退變是引起下腰痛的主要原因，由于體外尸體試驗并不能完全反映椎間盤在人體內(nèi)的機械環(huán)境，有限元被用于研究人體腰椎間盤損傷的一種新方法，通過數(shù)值計算得到椎間盤應(yīng)力、應(yīng)變和孔隙壓力的分布，并與體外試驗互相驗證。然而現(xiàn)有的椎間盤損傷模型都存在一定局限性：
[0004]首先是模型忽略了組織液的流入流出，因此難以捕捉到椎間盤的蠕變行為；
[0005]其次，采用不可壓縮的單一固相模型不能反應(yīng)出液相(組織液)的機械貢獻(xiàn)。

技術(shù)實現(xiàn)思路

[0006]針對以上模型的局限性，本專利技術(shù)給出了一個新的解決方案。在椎間盤多孔彈性模型的基礎(chǔ)上，提出纖維環(huán)的多孔
?
損傷數(shù)值計算模型，即將椎間盤纖維環(huán)視為由不可壓縮的液體和可壓縮的各向異性固相組成。新的建模思路將多孔模型和超彈性各向異性損傷模型相結(jié)合，多孔模型部分采用傳統(tǒng)的滲流
?
應(yīng)力耦合分析，得到椎間盤纖維環(huán)在固相有效應(yīng)力下的真實損傷演化，從而可以更準(zhǔn)確的描述椎間盤纖維環(huán)生物力學(xué)損傷響應(yīng)。利用這種方法可以捕捉到椎間盤組織的蠕變行為，并在空間上和時間上對特定工況下的椎間盤纖維環(huán)的膠原纖維損傷進(jìn)程進(jìn)行預(yù)測。
[0007]為了實現(xiàn)上述技術(shù)目的，本專利技術(shù)采用如下技術(shù)手段：
[0008]一種組織液
?
纖維環(huán)流固耦合的椎間盤軟組織損傷演變預(yù)測方法，結(jié)合多孔模型和超彈性各向異性損傷模型對椎間盤的力學(xué)響應(yīng)進(jìn)行數(shù)值模擬，實際模擬計算過程通過ABAQUS有限元軟件實現(xiàn)，該方法包括以下步驟：
[0009]1)在每一計算步開始時，根據(jù)滲流
?
應(yīng)力耦合分析方法，計算液相壓力p以及液相流動帶來的單元體積比J；將有限元積分點處的變形梯度F解耦為體積變化部分F
vol
和等容變化部分用解耦后的變形梯度分別計算固相體積部分應(yīng)變能密度Ψ
vol
和各向同性等容應(yīng)變能密度Ψ
iso
；
[0010]2)用子程序UEXTERNALDB導(dǎo)入外部纖維數(shù)據(jù)，儲存在纖維方向數(shù)組中，在單元計算時，與纖維方向數(shù)組中的單元編號進(jìn)行匹配，讀取當(dāng)前單元的兩族膠原纖維方向并計算未損傷的各向異性應(yīng)變能密度函數(shù)
[0011]3)根據(jù)當(dāng)前未損傷的各向異性應(yīng)變能密度函數(shù)對損傷準(zhǔn)則進(jìn)行判斷，計算纖維環(huán)中的各向異性損傷系數(shù)D
i
，并基于考慮損傷的固相應(yīng)變能密度函數(shù)，在空間參考下計算考慮損傷的固相柯西應(yīng)力σ
e
以及剛度張量
[0012]4)根據(jù)應(yīng)力和剛度張量公式編寫UMAT子程序來表征纖維環(huán)的材料損傷行為，對特定工況下某一時間點的椎間盤損傷位置和損傷值進(jìn)行預(yù)測。
[0013]所述步驟1)中滲流
?
應(yīng)力耦合分析的計算利用ABAQUS有限元軟件中的Soil分析步，將總柯西應(yīng)力σ看成是固相有效應(yīng)力σ
e
和液相壓力pI之和,其中 I為二階單位張量。并根據(jù)Darcy定理由孔隙壓力梯度計算組織液的等效凈流速場 v，從而得到每點處液體流出導(dǎo)致的單元體積變化：
[0014]σ＝σ
e
+pI
[0015]其中，μ
f
為動力粘度系數(shù)；v是等效凈流速場，表示液體相對固相的流動速度，k為應(yīng)變依賴的滲透率系數(shù)，由初始孔隙比e0和初始滲透率k0確定：
[0016][0017]其中M為修正參數(shù)，根據(jù)試驗數(shù)據(jù)擬合確定，e是現(xiàn)時孔隙比，與單元體積比的關(guān)系為：J＝(1+e)/(1+e0)。
[0018]所述步驟1)中，高斯積分點處的變形梯度解耦公式為其中F
vol
＝J
1/3
I，應(yīng)變能密度的各向同性部分采用Neo
?
Hooken本構(gòu)模型：
[0019][0020]其中D1和C
10
是材料參數(shù)，是等容部分右柯西
?
格林變形張量的第一不變量，J是體積比，表示現(xiàn)時體元大小與初始體元大小之比，同時也是變形梯度F的行列式值。
[0021]所述步驟2)中，未損傷纖維環(huán)的各向異性應(yīng)變能密度函數(shù)采用 Gasser
?
Holzapfel
?
Ogend模型進(jìn)行描述，兩族膠原纖維的初始方向分別由兩個纖維方向向量m0和n0表示，未損傷的各項異性部分應(yīng)變能密度函數(shù)由結(jié)構(gòu)張量表示，未損傷的各項異性部分應(yīng)變能密度函數(shù)由結(jié)構(gòu)張量和確定：
[0022][0023]其中，分別是等容部分的右柯西
?
格林變形張量的第四、六不變量；
[0024]k1表示纖維模量，
[0025]k2是指數(shù)系數(shù)，為無量綱量；
[0026]損傷折減后的固相各向異性應(yīng)變能密度函數(shù)為：
[0027][0028]所述步驟2)中，確定兩族膠原纖維方向的具體流程為：
[0029]21)首先將劃分的10層環(huán)狀基質(zhì)網(wǎng)格的外表面提取出來，得到共10層二維網(wǎng)格面，輸出二維網(wǎng)格面以及對應(yīng)層基質(zhì)單元的inp文件，為后續(xù)膠原纖維方向的確定作準(zhǔn)備；
[0030]22)處理inp文件得到每個單元在對應(yīng)面上的兩個對角線方向，將該單元的對角線方向作為纖維環(huán)在該單元處的膠原纖維方向，作為外部數(shù)據(jù)儲存在新的inp文件中，第一列代表單元編號，后續(xù)3列是該單元處膠原纖維方向在空間上的x,y,z分量；
[0031]23)在UMAT子程序中通過UEXTERNAL子程序?qū)w維方向進(jìn)行讀取，作為全局變量儲存在數(shù)組中，在計算流程中，對單元編號進(jìn)行識別，提取數(shù)組中的纖維方向并進(jìn)行歸一化。
[0032]步驟3)中所述纖維環(huán)中的各向異性損傷系數(shù)D
i
的計算具體流程為：
[0033]31)根據(jù)等容變形梯度計算當(dāng)前步下各向異性部分的
[0034]32)計算損傷判斷準(zhǔn)則中的指標(biāo)并與歷史最大指標(biāo)值進(jìn)行比較，若損傷面發(fā)生移動，根據(jù)下式計算損傷折減系數(shù)以及纖維環(huán)中的各向異性損傷系數(shù)D
i
[0035][0036][0037]其中，β
f
是膠原纖維的損傷演化規(guī)律控制參數(shù)；和分別表示損傷初始時刻和完全損傷時刻的對應(yīng)損傷指標(biāo)。
[0038]所述步驟3)中固相有效應(yīng)力是基于應(yīng)變能密度函數(shù)計算的，根據(jù)有限變形理論，在材料參考系下，纖維環(huán)的PK2固相應(yīng)力S...

【技術(shù)保護(hù)點】

【技術(shù)特征摘要】
1.一種組織液
?
纖維環(huán)流固耦合的椎間盤軟組織損傷演變預(yù)測方法，其特征在于，結(jié)合多孔模型和超彈性各向異性損傷模型對椎間盤的力學(xué)響應(yīng)進(jìn)行數(shù)值模擬，實際模擬計算過程通過ABAQUS有限元軟件實現(xiàn)，該方法包括以下步驟：1)在每一計算步開始時，根據(jù)滲流
?
應(yīng)力耦合分析方法，計算液相壓力p以及液相流動帶來的單元體積比J；將有限元積分點處的變形梯度F解耦為體積變化部分F
vol
和等容變化部分用解耦后的變形梯度分別計算固相體積部分應(yīng)變能密度Ψ
vol
和各向同性等容應(yīng)變能密度Ψ
iso
；2)用子程序UEXTERNALDB導(dǎo)入外部纖維數(shù)據(jù)，儲存在纖維方向數(shù)組中，在單元計算時，與纖維方向數(shù)組中的單元編號進(jìn)行匹配，讀取當(dāng)前單元的兩族膠原纖維方向并計算未損傷的各向異性應(yīng)變能密度函數(shù)3)根據(jù)當(dāng)前未損傷的各向異性應(yīng)變能密度函數(shù)對損傷準(zhǔn)則進(jìn)行判斷，計算纖維環(huán)中的各向異性損傷系數(shù)D
i
，并基于考慮損傷的固相應(yīng)變能密度函數(shù)，在空間參考下計算考慮損傷的固相柯西應(yīng)力σ
e
以及剛度張量4)根據(jù)應(yīng)力和剛度張量公式編寫UMAT子程序來表征纖維環(huán)的材料損傷行為，對特定工況下某一時間點的椎間盤損傷位置和損傷值進(jìn)行預(yù)測。2.根據(jù)權(quán)利要求1所述組織液
?
纖維環(huán)流固耦合的椎間盤軟組織損傷演變預(yù)測方法，其特征在于，所述步驟1)中滲流
?
應(yīng)力耦合分析的計算利用ABAQUS有限元軟件中的Soil分析步，將總柯西應(yīng)力σ看成是固相有效應(yīng)力σ
e
和液相壓力pI之和，其中I為二階單位張量，并根據(jù)Darcy定理由孔隙壓力梯度計算組織液的等效凈流速場v，從而得到每點處液體流出導(dǎo)致的單元體積變化：σ＝σ
e
+pI其中，μ
f
為動力粘度系數(shù)；v是等效凈流速場，表示液體相對固相的流動速度，k為應(yīng)變依賴的滲透率系數(shù)，由初始孔隙比e0和初始滲透率k0確定：其中M為修正參數(shù)，根據(jù)試驗數(shù)據(jù)擬合確定，e是現(xiàn)時孔隙比，與體積比的關(guān)系為：J＝(1+e)/(1+e0)。3.根據(jù)權(quán)利要求1所述組織液
?
纖維環(huán)流固耦合的椎間盤軟組織損傷演變預(yù)測方法，其特征在于，所述步驟1)中，高斯積分點處的變形梯度解耦公式為其中F
vol
＝J
1/3
I，應(yīng)變能密度的各向同性部分采用Neo
?
Hooken本構(gòu)模型：其中，D1和C
10
是材料參數(shù)；是等容部分右柯西
?
格林變形張量的第一不變量；J是體積比，表示現(xiàn)時體元大小與初始體元大小之比，同時也是變形梯度F的行列式值。4.根據(jù)權(quán)利要求1所述組織液
?
纖維環(huán)流固耦合的椎間盤軟組織損傷演變預(yù)測方法，其特征在于，所述步驟2)中，未損傷纖維環(huán)的各向異性應(yīng)變能密度函數(shù)采用Gasser
?
Holzapfel
?
Ogend模型進(jìn)行描述，兩族膠原纖維的初始方向分別由兩個纖維方向向量m0和n0表...

【專利技術(shù)屬性】
技術(shù)研發(fā)人員：郭小明，許萬玟，
申請(專利權(quán))人：東南大學(xué)，
類型：發(fā)明
國別省市：