黃桂雄，羅世興，紀(jì)德朋

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

資料與方法

1 數(shù)據(jù)來源

取一位健康男性志愿者，排除腰椎疾病，既往無腰椎骨折及結(jié)核等疾病；排除一些有可能影響到腰椎骨質(zhì)的基礎(chǔ)性疾病。采用CT掃描志愿者腰椎,將掃描活動(dòng)的CT圖像以DICOM的格式保存，以便導(dǎo)入相關(guān)軟件進(jìn)行建模。

2 相關(guān)軟件處理

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

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

3 邊界條件

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

結(jié) 果

建立了正常腰椎模型和L2椎體壓縮骨折模型，包括皮質(zhì)骨、松質(zhì)骨、終板及相關(guān)韌帶等結(jié)構(gòu)，見圖1、2，骨折數(shù)據(jù)見圖3。

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

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

圖4 不同工況下L1～2椎間盤最大應(yīng)力值

圖5 不同工況下L2～3椎間盤最大應(yīng)力值

表1 L1～2椎間盤最大應(yīng)力分布部位

表2 L2～3椎間盤最大應(yīng)力分布部位

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

圖6 不同工況下L1關(guān)節(jié)突最大應(yīng)力值

圖7 不同工況下L2關(guān)節(jié)突最大應(yīng)力值

圖8 不同工況下L3關(guān)節(jié)突最大應(yīng)力值

表3 不同工況下L1關(guān)節(jié)突最大應(yīng)力分布部位

表4 不同工況下L2關(guān)節(jié)突最大應(yīng)力分布部位

表5 不同工況下L3關(guān)節(jié)突最大應(yīng)力分布部位

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

圖9 正常模型在10Nom的彎矩下做前屈運(yùn)動(dòng)的云圖

圖10 骨折模型在10Nom的彎矩下做前屈運(yùn)動(dòng)的云圖

討 論

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

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

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

在不同彎矩的作用下，前屈運(yùn)動(dòng)時(shí)大部分應(yīng)力位于椎間盤的前右側(cè)，這是因?yàn)槿梭w在受向前運(yùn)動(dòng)彎矩的作用下，重力也向前移動(dòng)，而因?yàn)閼T用右手的人，脊柱的上胸椎偏向右，而腰椎則代償性凸向左側(cè)[23]，所以重心偏向右側(cè)，導(dǎo)致右側(cè)受應(yīng)力大于左側(cè)；當(dāng)在后伸的彎矩作用下，應(yīng)力的分布偏向后右側(cè)的原理也如此,說明了臨床上有些患者腰椎間盤突出偏向右的原因。

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

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

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

[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] 延?xùn)|，毛景松，石先明,等.不同載荷下腰1椎體內(nèi)應(yīng)力分布規(guī)律的有限元分析[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] 柏樹令.系統(tǒng)解剖學(xué)[M].2版.北京：人民衛(wèi)生出版社，2013:37-44.

[24] 姚衛(wèi)星.結(jié)構(gòu)疲勞壽命分析[M].北京：國防工業(yè)出版社，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.