不同彎矩作用下腰椎正常與壓縮骨折模型的有限元分析

黃桂雄，羅世興，紀德朋

有限元分析方法是用于骨科受力分析的一種常用的方法[1-3]，因其方便及可重復性，越來越受到國內外學者的青睞。腰椎壓縮骨折常導致脊柱失穩，甚至引起脊髓損傷，常見于交通事故或高處墜落[4]。按照壓縮的椎體前后高度比，分為三度：Ⅰ度，前后比<1/3；Ⅱ度，前后比約為1/2；Ⅲ度，前后比>2/3。因為腰椎壓縮骨折后會造成脊柱不穩，有學者提出，當腰椎壓縮骨折類型為II度或III度時應行手術治療，而I度的壓縮骨折則提倡非手術治療。許多學者對腰椎壓縮骨折手術方式進行研究，但對于腰椎I度壓縮骨折是否必要行手術治療仍存在爭議[5-6]。本研究通過運用Mimics[7]、 Hypermesh[8]和Abaqus[7]軟件構建腰椎正常模型和L2壓縮骨折(I度)模型進行對比受力分析，目的是觀察在不同彎矩應力下兩種模型的椎間盤及關節突的力學分布差異。

資料與方法

1 數據來源

取一位健康男性志愿者，排除腰椎疾病，既往無腰椎骨折及結核等疾病；排除一些有可能影響到腰椎骨質的基礎性疾病。采用CT掃描志愿者腰椎,將掃描活動的CT圖像以DICOM的格式保存，以便導入相關軟件進行建模。

2 相關軟件處理

將掃描后的腰椎CT數據，選擇L1～3區域的圖片,導入Mimics 16.0軟件。進行三維重建，將在Mimics軟件中生成的腰椎正常模型，并將正常L2椎體上端切割，建立壓縮骨折模型；再分別以STL格式導入Hypermesh中進行面網格劃分,然后再以STL格式導入Abaqus軟件中劃分體網格，并建立終板、松質骨及皮質骨[9]集合；椎體、纖維環和髓核材料屬性設為各向同性；裝配兩種模型，并以Spring[10-11]單元模擬建立前縱韌帶、后縱韌帶、脊間韌帶、脊上韌帶、黃韌帶、橫突間韌帶及關節囊韌帶，并按照相關文獻[12-13]，賦予韌帶、松質骨、皮質骨及終板等物理屬性。將腰椎和纖維環、腰椎和髓核之間的關系設為綁定(Tie)，上下關節突關節面之間設為面面接觸關系，摩擦系數為0.15。接觸設置為硬接觸。

模擬：工況一，模型在一大小為500N垂直向下作用下，模擬正常人體直立上半身時所受的重力；工況二：工況一與一個彎矩大小為5Nom做前屈運動；工況三：工況一與一個彎矩大小為10Nom做前屈運動；工況四：工況一與一個彎矩大小為15Nom做前屈運動；工況五：工況一與一個彎矩大小為5Nom做后伸運動；工況六：工況一與一個彎矩大小為10Nom做后伸運動；工況七：工況一與一個彎矩大小為15Nom做后伸運動。

3 邊界條件

L3椎體和L3椎體雙側下關節突完全固定，并計算在不同工況下不同模型的L1～2椎間盤、L2～3椎間盤、L1關節突、L2關節突和L3關節突受力分布情況，結果用米塞斯應力(Von Mises)表示。

結 果

建立了正常腰椎模型和L2椎體壓縮骨折模型，包括皮質骨、松質骨、終板及相關韌帶等結構，見圖1、2，骨折數據見圖3。

圖1 正常模型圖2 骨折模型圖3 L2椎體壓縮骨折模型參數

在本實驗中，骨折模型的L1～2椎間盤最大應力值可超出正常模型的132.93%，出現在以10Nom前屈運動時；最小可超出正常模型的15.48%，出現在以15Nom做后伸運動時。骨折模型的L2～3椎間盤最大應力值可超過正常模型的66.5%，出現在以15Nom做后伸運動時；最小可超出正常模型的11.59%，出現在直立體位時(圖4～6)。骨折模型椎間盤的應力在七種工況下均大于正常模型，并且骨折模型的應力分布部位與正常模型不完全一致(表1、2)。

圖4 不同工況下L1～2椎間盤最大應力值

圖5 不同工況下L2～3椎間盤最大應力值

表1 L1～2椎間盤最大應力分布部位

表2 L2～3椎間盤最大應力分布部位

在本研究中，骨折模型的關節突在不同彎矩作用下，前屈體位時，L1關節突和L2關節突最大應力值均小于正常模型，而在后伸體位，隨著彎矩逐漸增大和關節突位置的降低，出現最大應力值逐漸增大的過程(圖6～8)。關節突最大應力的分布部位與正常模型有所差別，分布無明顯規律(表3～5)。

圖6 不同工況下L1關節突最大應力值

圖7 不同工況下L2關節突最大應力值

圖8 不同工況下L3關節突最大應力值

表3 不同工況下L1關節突最大應力分布部位

表4 不同工況下L2關節突最大應力分布部位

表5 不同工況下L3關節突最大應力分布部位

本研究結果和Chien等[14]的研究結果大致相同，以L1～2椎間盤，在10Nom前屈和后伸的計算結果云圖為例(圖9、10)。

圖9 正常模型在10Nom的彎矩下做前屈運動的云圖

圖10 骨折模型在10Nom的彎矩下做前屈運動的云圖

討 論

有限元分析[15-16]是醫學類特別是骨科受力學分析常用的研究方式。Travert等[17]已運用有限元的方法來分析材料屬性和外部載荷對脊柱強度的影響。本研究結合Mimics、Abauqs、Hypermesh軟件，建立腰椎的相關模型并進行受力分析,建立L1～3正常腰椎模型和合并L2椎體壓縮骨折的L1～3腰椎模型并進行有限元分析，觀察在不同彎矩應力下兩種模型的椎間盤及關節突的力學分布差異。

本研究結果表明，L2壓縮骨折模型在不同的前屈后伸的彎矩作用下其L1～2椎間盤和L2～3椎間盤受到應力為最高。因為應力的增加，最主要的并發癥是鄰近節段退變[18]，從而引起椎間盤突出。Saleem等[19]研究表明，腰椎間盤突出和腰椎管狹窄大多數原因是由于長期腰椎退變和抵抗應力的綜合作用的結果。綜合學者研究結果，筆者認為腰椎Ⅰ度壓縮骨折的患者脊柱的退變可能會加速,甚至將來有可能導致續發骨折[20]。

有研究表明[21]男性比女性更容易因為機械應力的增加而損害腰椎間盤。本研究結果表明，隨著前屈和后伸的彎矩增大，大部分模型的椎間盤受到的應力均增大，前屈的動作較后伸的動作增大更明顯，這就驗證了前屈的動作會增加腰椎間盤突出的概率。從事重體力活動且伴有頻繁彎腰動作容易出現腰椎間盤突出，并且胸腰段脊柱是整個脊柱的應力集中部位，容易發生骨折，尤其在骨質疏松的老年人多見[22]。

在不同彎矩的作用下，前屈運動時大部分應力位于椎間盤的前右側，這是因為人體在受向前運動彎矩的作用下，重力也向前移動，而因為慣用右手的人，脊柱的上胸椎偏向右，而腰椎則代償性凸向左側[23]，所以重心偏向右側，導致右側受應力大于左側；當在后伸的彎矩作用下，應力的分布偏向后右側的原理也如此,說明了臨床上有些患者腰椎間盤突出偏向右的原因。

本研究中，在直立狀態下腰椎壓縮骨折模型椎間盤的應力明顯大于正常模型的應力，而大部分骨折模型的關節突在前屈時應力較正常模型小；在不同的體位和不同的彎矩大小作用下，骨折模型在不同節段的關節突應力成分考慮為關節突和椎弓根受前屈到后伸的體位改變的影響，分別為關節突壓應力減小，椎弓根到拉應力增大，再轉變為關節突壓應力增大的過程；說明腰椎壓縮骨折后，關節突應力的大小和后伸的幅度和關節突的節段均有關，并且在骨折模型中，應力增大的關節突，因為力學機制的變化，有可能對關節突本身產生不利的影響。筆者認為，在直立體位時骨折模型的關節突應力最大值較正常模型的大，可能是腰椎壓縮骨折后，脊柱后凸畸形，導致脊柱重心向前，即重力的向前水平分力增加，后部受水平分力減少，相應的應力也減小。所以可推測，椎間盤和關節突的應力大小與椎體的形狀可能相關。再者，腰椎壓縮骨折時，由于壓縮的椎體橫截面積發生改變[24]，導致應力集中分布于椎體壓縮骨折皮質骨的前上半部。L1～2椎間盤的應力分布>L2～3椎間盤的原因，有可能是應力分布以L2壓縮骨折上終板前端和皮質骨交界處為中心，離該中心越遠，相同的材料屬性，應力相對越小。在不同前屈和后伸的彎矩作用下，兩種模型的椎間盤在前屈時應力增加幅度較大，說明前屈的動作對椎間盤的損害較大；后伸時椎間盤應力變化不明顯，考慮為椎體后部的楊氏模量值較椎間盤的值大，應力主要由椎體后部結構承擔，這就是應力遮擋的原理[25]。所以前屈時，椎體重心前移，脊柱骨皮質的楊氏模量值較大，承擔了大部分應力，所以腰椎后部的關節突應力變化不明顯。L1、L2和L3關節突最大應力分布情況，大部分分布在右下關節突關節面的上方，部分分布在椎弓根，這就可能導致關節突關節因為應力集中容易退變[26-27]。

本研究假設在腰椎所能承受的應力范圍內進行分析，不考慮各個材料的塑性區和斷裂區。結果表明，壓縮骨折模型腰椎間盤和部分腰椎關節突的應力較正常模型大。在生理活動下，由于高周疲勞(high-cycle fatigue，是指材料在低于其屈服強度的循環應力作用下，經10 000～100 000以上循環次數而產生的疲勞。高周疲勞的特點是作用于零件或構件的應力水平較低)[24]，有壓縮骨折病史的患者再發生骨折的可能性大。

綜上所述，有限元力學分析顯示腰椎壓縮骨折后，傷椎的椎間盤和關節突的壓應力較正常椎體增大，提示骨折后可能脊柱退變加速。相對關節突，椎間盤的應力分布較規律，關節突應力分布的不規律可能提示骨折后脊柱后柱退變加速具體部位具有不確定性。

[1] Divya VA，Edward KW,Ronald AL，et al.Bilateral pedicle screw fixation provides superior biomechanical stability in transformational lumbar interbody fusion: a finite element study [J].Spine J，2015，15(8)：1812-1822.

[2] Bisschop A,Holewijn RM,Kingma I，et al.The effects of single-level instrumented lumbar laminectomy on adjacent spinal biomechanics[J].Global Spine J，2015，5(1)：39-48.

[3] Agarwal A,Zakeri A,Agarwal AK,et al.Distraction magnitude and frequency affects the outcome in juvenile idiopathic patients with growth rods: finite element study using a representative scoliotic spine model[J].Spine J，2015，15(8)：1848-1855.

[4] Ivancic PC.Biomechanics of thoracolumbar burst and chance-type fractures during fall from height[J].Global Spine J，2014，4(3)：161-168.

[5] Dohm M,Black CM,Dacre A,et al.A randomized trial comparing balloon kyphoplasty and vertebroplasty for vertebral compression fractures due to osteoporosis [J].AJNR Am J Neuroradiol，2014，35(12)：2227-2236.

[6] Saracen A,Kotwica Z.Treatment of multiple osteoporotic vertebral compression fractures by percutaneous cement augmentation[J].Int Orthop，2014，38(11)：2309-2312.

[7] Faraipour H,Jamshidi N.Effects of different angles of the traction table on lumbar spine ligaments: a finite element study[J].Clin Orthop Surg,2017,9(4):480-488.

[8] Filardi V,Simona P,Cacciola G,et al.Finite element analysis of sagittal balance in different morphotype: forces and resulting strain in pelvis and spine[J].J Orthop,2017,14(2):268-275.

[9] Garo A,Arnoux PJ,Wagnac E, et al.Calibration of the mechanical properties in a finite element model of a lumbar vertebra under dynamic compression up to failure[J].Med Biol Eng Comput，2011，49(12)：1371-1379.

[10] Jung TG，Woo SH，Park KM,et al.Biomechanical behavior of two different cervical total disc replacement designs in relation of concavity of Articular Surfaces: ProDisc-C?VS Prestige-LP?[J].Int J Prec Eng Manufac，2013，14(5)：819-824.

[11] Clin J,Aubin CE,Lalonde N,et al.A new method to include the gravitational forces in a finite element model of the scoliotic spine[J].Med Biol Eng Comput，2011，49(8)：967-977.

[12] Goel VK，Monrone BT，Gilbertson LG，et al.Interlaminar shear stresses and laminae separation in a disc: finite element analysis of the L3-4motion segment subjected to axial compressive loads[J].Spine，1995，20(6)：689-698.

[13] Jones AC,Wilcox RK.Finite element analysis of the spine: towards a framework of verification,validation and sensitivity analysis[J].Med Eng Phys，2008，30(10)：1287-1304.

[14] Chien CY,Kuo YJ,Lin SC,et al.Kinematic and mechanical comparisons of lumbar hybrid fixation using dynesys and cosmic systems[J].Spine(Phila Pa 1976)，2014,39(15):E878-884.

[15] Imai K.Analysis of vertebral bone strength,fracture pattern,and fracture location: a validation study using a computed tomography-based nonlinear finite element analysis[J].Aging Dis，2014，6(3)：180-187.

[16] 延東，毛景松，石先明,等.不同載荷下腰1椎體內應力分布規律的有限元分析[J].中國脊柱脊髓雜志，2014，24(9)：882-887.

[17] Travert C,Jolivet E,Sapin-de Brosses E,et al.Sensitivity of patient specific vertebral finite element model from low dose imaging to material properties and loading conditions[J].Med Biol Eng Comput，2011，49(12)：1355-1361.

[18] Palepu V， Kodigudla M，Goel VK.Biomechanics of disc degeneration[J].Adv Orthop，2012，2012：1-17.

[19] Saleem S,Aslam HM,Rehmani MA,et al.Lumbar disc degeneration disease: disc degeneration symptoms and magnetic resonance image findings[J].Asian Spine J,2013，7(4)：322-334.

[20] Ren HL,Jiang JM,Chen JT,et al.Risk factors of new symptomatic vertebral compression fractures in osteoporotic patients undergone percutaneous vertebroplasty[J].Eur Spine J，2015,24(4):750-758.

[21] Wang YX，Griffith JF.Effect of menopause on lumbar disk degeneration: potential etiology[J].Radiology，2010，257(2)：318-320.

[22] Bouza C,López-Cuadrado T,Almendro N，et al.Safety of balloon kyphoplasty in the treatment of osteoporotic vertebral compression fractures in Europe: a meta-analysis of randomized controlled trials[J].Eur Spine J，2015，24(4)：715-723.

[23] 柏樹令.系統解剖學[M].2版.北京：人民衛生出版社，2013:37-44.

[24] 姚衛星.結構疲勞壽命分析[M].北京：國防工業出版社，2004：1-27.

[25] Zhao J,Chen LJ,Yu K,et al.Effects of chitosan coating on biocompatibility of Mg-6%Zn-10% Ca3(PO4)2implant[J].Trans Nonferrous Met Soc China ,2015,25(3):824-831.

[26] Rohlmann A,Lauterborn S,Dreischarf M,et al.Parameters influencing the outcome after total disc replacement at the lumbosacral junction.Part 1: misalignment of the vertebrae adjacent to a total disc replacement affects the facet joint and facet capsule forces in a probabilistic finite element analysis[J].Eur Spine J,2013,22(10):2271-2278.

[27] Park WM,Kim K,Kim YH.Changes in range of motion,intradiscal pressure,and facet joint force afterintervertebral disc and facet joint degeneration in the cervical spine[J].J Mech Sci Technol,2015,29 (7): 3031-3038.