瞬態載荷下L4/L5 椎間盤內流固耦合效應1)

李華忠 王 博 王丹丹 劉思金 陳曉東,2) 師 彬,

* (北京理工大學宇航學院,北京 100081)

? (山東第一醫科大學(山東省醫學科學院),山東第一醫科大學附屬頸肩腰腿痛醫院,山東省骨生物力學工程實驗室,濟南 250062)

** (中國科學院生態環境研究中心,環境化學與生態毒理學國家重點實驗室,北京 100085)

引言

椎間盤(intervertebral disk,IVD)是連接相鄰椎骨的無血管的膠原結構,包括髓核(nucleus pulposus,NP)和纖維環(annulus fibrosus,AF)[1].在脊柱上傳遞和分配載荷,使脊柱具有柔韌性[2],椎間盤退變是最常見的人類運動系統慢性損傷性疾病之一,嚴重影響人類健康[3-4].生物力學變化在椎間盤退變的發生發展中起到重要作用[5],L4/L5 (第4 腰椎、第5 腰椎)椎間盤處于人體脊柱載荷傳輸的中樞位置,是椎間盤退變的高發部位[1,6-7].外部載荷通過產生椎間盤內固體基質應力?應變影響流體流動,髓核內流體的流動在維持椎間盤的機械和生物學性能中起重要作用[8].髓核流動減少會對椎間盤的營養、細胞代謝、基質金屬蛋白酶活性及多種細胞因子等因素組成的生物化學環境造成影響[5],已經被認為是引起椎間盤退變和椎間盤細胞營養物質代謝改變的主要原因[9-13].因此深入地研究椎體應力?應變[14]與椎間流動的相互作用關系,對于增進對椎間盤退變病理生理過程[13]的理解,解釋其退變機制,促進有效治療方法[15]的發展具有重要意義.

對椎間盤流動的實驗研究闡明了進出椎間盤的流體流動在保持椎間盤的性能中的重要作用[16-18].例如,van der Veen 等[19]對豬脊柱標本的研究顯示了體外模型在加載過程中發生了流體的流出;McMillan等[20]通過對人尸體腰椎間盤的加載,觀察到持續加載會減小椎間盤的流體含量;Bowden 等[21]采用磁共振擴散加權成像技術(diffusion-weighted magnetic resonance imaging,DW-MRI)研究椎間盤的流體流動,發現運動狀態的椎間盤表觀擴散系數高于久坐狀態,表明運動有助于保持椎間盤健康.這些實驗研究尚不能再現椎間盤的孔隙彈性行為,也無法在完整的動態過程中實現對流動場的測量.

有限元分析是評估椎間盤應力?應變以及流動的有效工具.建立考慮流固耦合效應的數值模型[22-23]可以反映固體變形與流體流動的相互作用關系,獲得實驗難以測定的腰椎間盤內部的固體應力?應變和流體壓力流速的分布,是研究腰椎生物力學的一條重要途徑.例如,Iatridis 等[24]發現髓核根據載荷速率可表現出不同程度的流體或固體特征,具有典型的流固耦合特征;Hassan 等[9]使用孔隙彈性模型開展有限元分析,發現終板滲透率和孔隙率降低可導致營養液的缺乏和髓核細胞凋亡;Gu 等[25]基于生物力學、電化學連續介質混合理論,建立椎間盤多孔混合物有限元模型,發現組織變性會通過降低孔隙率減少椎間盤細胞的營養供應,導致椎間盤中的細胞數量減少,并引起細胞介導的椎間盤變性.

中醫正脊治療是臨床非手術療法中治療椎間盤退變的首選手段[26],是一種公認的治療方法[27-29],目前對其生物力學效應機制還知之甚少.Deng 等[26]通過有限元模擬發現正脊手法顯著改變頸椎周圍的應力區域,減小了應力集中;田強等[30]的模擬發現纖維環后外側的位移變化可能是提拉旋轉斜扳法治療腰椎間盤突出癥有效的機制之一.目前認為椎間盤退變是力學和生物學相互作用的惡性循環[11],是機械過載、髓核失水和細胞分解等流體和固體組織相互作用的復雜生物學過程[3],對這種復雜病理生理過程的干預也是正脊治療有效的生物力學原理.因此,建立流固耦合模型以分析正脊治療中的生物力學效應,對于解釋腰椎間盤損傷退變機制和評估治療效果至關重要.

本工作測量臨床治療過程中正脊醫師對患者施加的載荷,提取該手法的典型加載曲線,進而確定用于模擬正脊手法的加載方式;根據掃描圖像和解剖學知識構建腰椎組織的幾何模型;建立基于線彈性和多孔彈性本構關系的有限元模型,開展考慮流固耦合效應的有限元分析.研究為中醫正脊治療提供了機理性認識和科學化依據,也為臨床治療中的個性化定制提供了方法和思路[31].

1 方法

1.1 載荷條件

正脊治療是指臨床醫生在脊柱上沿受控方向向目標部位施加快速、低振幅的特定大小力,導致脊柱和周圍軟組織產生瞬間變形后恢復,來治療早期椎間盤退變性疾病的方法[32-33].在臨床治療中,醫生對病人所施加的力因人而異,主要靠臨床醫生的經驗來確定.為了得到合理的載荷參數,并不失一般性,首先對正脊手法[34-35]所產生的力進行測量,進而結合文獻調研確定數值模擬采用的力的加載曲線.

使用壓力采集系統[35](上海瑞諾測控設備有限公司),以1000 Hz 的采樣頻率,在體測量志愿者在接受L4/L5 節段正脊治療中的載荷參數.測量中將薄膜傳感器貼在右側肩峰前方受力點,該部位是治療期間的主要負荷點,接受拉伸和旋轉的復合力作用[36].測量結果顯示,手法力加載可分為預加載(preload)、推動(thrust)和卸載(resolution) 3 個階段(如圖1 所示),與文獻[26,28,37-38]報道的脊柱扳動力作用規律類似.預加載階段為0 到0.6 s,推動階段為0.6～0.7 s,0.7 s 之后為卸載階段,至1.7 s 外加載荷消失.其中推動階段的時長與Herzog 等[28]的測量數據吻合(后者為0.081～ 0.15 s).

圖1 正脊手法中3 個階段的定義Fig.1 Definitions of three phases during spinal manipulation.

腰椎旋轉手法實際作用于患者肩部,對于腰椎而言,可以分解為拉伸(distraction)和旋轉(rotation)兩個作用.由于實際作用在椎體上的力難以確定,本文將拉伸作用通過在第3 腰椎(L3)表面加載軸向拉伸力實現,將旋轉作用通過在L3 表面左右端加載方向相反的力產生旋轉力矩來實現,以模擬傳遞到腰椎節段的拉伸和旋轉.

拉伸力的預加載值取為50 N,接近于蘇少亭等[39]在腰椎定點旋轉手法在體運動力學量化研究中測得的預加載力(約為43 N).拉伸力的峰值取為300 N,接近于文獻[38,40-41]中使用三維力測量系統直接測量手法操作時醫生肢體與患者之間峰值力(228 ±96) N.旋轉力矩的預加載值取為5 N·m,旋轉力矩的峰值取為22 N·m,接近于劉強等[42]實驗測量腰椎扳法對椎間盤內壓造成影響采用的25 N·m 轉矩.該峰值轉矩引起L4 棘突產生2 mm 的位移,與寧尚龍等[43]對腰椎在體運動的觀察的位移大小(1.2 ± 0.8) mm接近;引起L4 棘突旋轉的角度為3°,接近于張軍等[44]基于退變腰椎間盤模型的旋轉手法測得的椎體旋轉角度位移 2.60° ± 2.19°.

1.2 幾何模型

本文考慮L4/L5 椎間盤在不同加載形式下的流固耦合效應,為保證載荷更符合實際,幾何模型考慮L3 到L5 腰椎節段.腰椎幾何數據采集自健康中國女性志愿者(年齡42 歲,身高164 cm,體重60 kg),使用西門子16 排螺旋CT 機掃描受試者的腰椎,最終CT 圖像分辨率為0.54 mm×0.54 mm,層間距為0.625 mm,以DICOM 格式保存.在軟件Mimics 21.0(Materialise Inc.,Leuven,Belgium)內,導入DICOM圖像進行處理;在分割模塊調整灰度值范圍,利用閾值選取分離出骨性結構;運用計算3D (calculate part)模塊對圖片進行3D 計算建模,再對表面進行去三角等光滑處理、細化網格優化結構,生成.stl 文件導入3-Matic 工具.

在3-Matic 中利用設計工具,基于腰椎的上下表面的形狀和位置,結合解剖學參數構建椎間盤和終板.每個椎間盤被模擬成一個中央髓核,周圍環繞著一個環形基質纖維環,上下分別覆蓋軟骨終板和骨性終板.其中,髓核的體積大約是整個椎間盤體積的44%.根據文獻[45-46]確定韌帶起止點及橫截面積,進而創建韌帶的三維實體幾何.模型包括6 組主要韌帶:前縱韌帶、后縱韌帶、黃韌帶、棘間韌帶、棘上韌帶和關節囊韌帶.模型中將關節突關節軟骨視為一個簡化的關節間隙填充板[47],其厚度等同于從CT 數據重建的骨性關節突關節間隙.

將在3-Matic 中創建的多個幾何體分別導出為.stl 格式文件,分別導入Geomagic Wrap 2017 軟件(3D Systems Corporation,USA),在該軟件對模型表面剩余棱角進行打磨、光滑后,運用構建曲面片,進行曲面擬合,保存最后生成的NURBS 曲面光滑實體模型.將實體幾何分別導出為.igs 文件,再次輸入SolidWorks 2020 軟件(SolidWorks Corp,Dassault Systèmes,USA)檢查曲面無誤后,以.igs 格式保存.文件導入DesignModeler 2020 軟件(ANSYS Inc,USA)進行布爾運算和組裝,完成后以.stp 格式導入COMSOL 5.6 軟件 (COMSOL Inc.,USA),形成聯合體.

1.3 數學模型

本研究將皮質骨和韌帶考慮為線性彈性介質,將松質骨、纖維環和髓核考慮為多孔彈性介質,各部分的材料屬性見表1 和表2.

表1 L3-L5 節段中組織的材料屬性Table 1 Material properties of tissues in the L3-L5 segment

表2 韌帶的材料屬性Table 2 Material properties of ligaments

線性彈性介質滿足胡克定律,即[49]

其中,σ為應力,C為4 階彈性張量,εel為彈性應變.

采用Biot 多孔彈性理論描述多孔介質中流體流動和變形之間的相互作用.應力、應變和孔隙壓力的本構關系為[49]

其中,σ為Cauchy 應力張量,ε為應變張量、αB為Biot-Willis 系數、pf為流體孔隙壓力、彈性矩陣C通過恒定孔隙壓力下應力變化引起的應變來測量.

在Biot 的理論[9]中,另一個本構關系是流體含量ζ的增量與體積應變和增量孔隙壓力有關.流體孔隙壓力與孔隙基質的膨脹和流體含量的變化成正相關[49]

其中,變量M可稱為Biot-Willis 模量,是達西定律中存儲系數S的倒數.S可以用孔隙度εp,αB,流體Kf的體積模量和多孔材料基質的模量Kd計算[49]

對于基質軟的多孔介質(髓核),αB可取為1，對于基質硬的多孔介質(松質骨、纖維環)

多孔介質中的流場通過達西定律來描述,流體的質量守恒為[49]

其中,?εvol/?t 是多孔基質體積應變的變化率,ρf是流體密度.考慮重力時的達西速度為[49]

1.4 邊界條件

圖2 顯示了本文進行4 個模擬的邊界條件示意圖.模擬中將L5 椎體的下端在所有方向均固定,在L3 椎體上表面加載不同的載荷.圖2(a)為在L3 上表面加載500 N 的壓縮載荷Fv,用于驗證模型在人體直立情況下的準確性.圖2(b)為同時加載1.1 節所述的拉伸力Fd和扭轉力矩Mr(順時針),考慮旋轉拉伸同時加載(rotation and distraction,R&D).圖2(c)和圖2(d)分別考慮旋轉(rotation,R)和拉伸(distraction,D)加載.流體流動方面,纖維環的外邊界的孔隙壓力設為零[9],其他外邊界設為無流動.

圖2 腰椎載荷邊界條件示意圖Fig.2 Schematic diagram of the load and boundary conditions of the lumbar spine

2 結果和討論

2.1 模型驗證

圖3(a)顯示了壓縮載荷下椎間盤內最大壓力隨時間的變化,最大值為0.68 MPa,數值與Hassan 等[9]的模擬結果(0.53 MPa)類似.圖3(b)顯示了椎間盤的高度損失平衡至0.11 mm,與Heuer 等[50]和Hassan 等[9]報道的在相同載荷下的模擬結果類似,分別為0.16 mm 和0.15 mm.椎間盤內流速最大流速1.15 μm/s,平衡到1.00 μm/s,如圖3(c)所示,與Hassan 等[9]報道流速范圍0.50～ 5.00 μm/s 接近,與Ferguson 等[51]報道的椎間盤內最大流速數量級相同.模型所表現出來的蠕變響應證明本模型可以描述腰椎及腰椎間盤的黏彈性行為.

圖3 壓縮載荷引起椎間盤內壓力、高度損失和最大流速的變化Fig.3 Time variations of pressure,height loss and maximum flow velocity in the intervertebral disc caused by a compress load

2.2 椎間盤應力?應變

椎間盤的應力?應變可以反映外力加載下椎間盤的響應,解釋手法作用原理并預警操作風險.最大變形出現在載荷的峰值時刻,即t=0.7 s.圖4(a)顯示3 種加載方式下髓核的位移,對于旋轉拉伸、單一旋轉和單一拉伸加載,髓核的最大位移分別為536 μm,527 μm 和56 μm.前兩種加載下,髓核水平方向的位移大于垂直方向的位移.3 種加載下纖維環的最大位移分別為495 μm,486 μm 和69 μm,如圖4(b)所示.圖5 顯示了3 種加載下椎間盤的von Mises 應力分布情況.最大應力均位于纖維環前方,最大值分別為9.93 MPa,9.98 MPa 和1.53 MPa,髓核應力小于纖維環.

圖4 不同加載方式下椎間盤不同部位的位移Fig.4 Displacement of different parts of intervertebral disc under different loading modes

圖5 在0.7 s 時不同加載方式下椎間盤應力分布云圖Fig.5 Contours of intervertebral disc stress distribution at 0.7 s under different loading modes

結果表明,旋轉加載所造成的纖維環水平位移和最大應力大于單一拉伸加載.模擬結果與龐胤等[52]研究發現的腰椎椎間盤在軸向、前屈、后伸及側彎4 種工況下,纖維環形變較大且出現明顯的應力集中的結果相符.因此治療中施加的載荷需警惕損傷纖維環風險.

2.3 髓核壓力流速

流體交換是椎間盤的生物學行為特征,是髓核營養物質交換的主要途徑.髓核流速與載荷和滲透壓有關,瞬態流速取決于外部載荷,外部加載力通過引起髓核形狀和內部壓力變化引起流體流動[53].圖6為不同加載方式下髓核內最大流速的絕對值.旋轉拉伸、單一旋轉加載下髓核流速較大,這是由于旋轉加載引起較大的髓核變形.圖6 中的兩個峰值是由于髓核內的流體先流出后流入.3 種加載的髓核最大流速均出現在載荷的峰值時刻,即t=0.7 s.旋轉拉伸和單一旋轉加載產生的最大流速基本一致,遠大于單一拉伸加載.為研究載荷對髓核流動的效應,下文對t=0.7 s 時刻的髓核內壓力與流速分布進行分析.

圖6 不同加載方式在加載周期內髓核最大流速變化曲線Fig.6 Time variation of the maximum flow velocity of the nucleus pulposus in different loading modes during loading cycle

圖7(a)顯示了在旋轉拉伸、單一旋轉和單一拉伸加載下髓核表面壓力分布.壓力最大值分別是1.85 MPa,1.94 MPa 和0.03 MPa.前兩種加載下,髓核左、右側邊緣壓力大于髓核表面其他區域.髓核表面流速的絕對值和方向分布如圖7(b)所示,3 種加載的流速最大值分別為8.08 μm/s,7.94 μm/s 和1.47μm/s.

圖7 0.7 s 時不同加載方式下的髓核壓力與流速Fig.7 Nucleus pulposus pressure and flow velocity at 0.7 s under different loading modes

為研究流入髓核的流體體積,定義髓核表面的法向速度的數值為

其中,un為法向流速,u為表面流速,n為表面法向向量.流入髓核內的法向流速定義為負值,流出髓核外的法向流速為正值.

3 種加載下髓核法向速度的數值的分布如圖7(c)所示.旋轉拉伸和單一旋轉加載時,髓核右側半流動方向為流出、左側半流動方向為流入,二者流入速度最大值分別為3.18 μm/s 和3.34 μm/s,流出速度最大值分別為3.81 μm/s 和3.61 μm/s.單一拉伸力加載的平均流速在髓核后邊緣為流出,余部位為流入,流入速度和流出速度最大值分別為0.47 μm/s 和0.79 μm/s.

結果顯示,旋轉拉伸和單一旋轉加載下的髓核內流速和壓力變化都大于單一拉伸加載情況,即旋轉作用產生髓核內外壓力梯度變化并加快髓核流動.在旋轉載荷的作用下,髓核左右產生相反的法向流速,即同一時刻既有液體流入,又有液體流出.

2.4 髓核含水量變化

健康的髓核是凝膠狀、高度水合的組織[3],因而,充足的含水量是髓核維持內部靜水壓力、拉伸纖維環、支撐終板以保持軸向壓縮中椎間盤高度和剛度的主要決定因素[54],外部載荷和內部滲透壓之間的動態平衡驅動著流體流動.為研究加載力與髓核內含水量的關系,本節討論髓核表面平均法向流速隨時間的變化,以及通過時間積分,計算加載周期0～ 120 s 內髓核含水量變化.由于旋轉力加載引起髓核左右側半流動行為差異,因此將髓核從中心劃分為左右側進行研究.

圖8(a)和圖8(b)分別顯示了不同加載下髓核左側和右側的平均流速變化規律.單一拉伸下,髓核左右兩側流動方向為先流出再流入,再少量流出后達到平衡.單一旋轉下,髓核左右兩側流動方向相反,左側為先流入再流出、右側為先流出再流入.旋轉拉伸復合加載下,髓核表面平均流速特征表現為兩種單一加載力的復合效應.這樣,左側平均流速在復合加載和單一拉伸加載下的變換趨勢相反.

圖8 不同加載方式下髓核平均流速和含水量變化Fig.8 Average flow velocity and fluid volume of nucleus pulposus under different loading modes

圖8(c)和8(d)分別顯示了不同加載下髓核左側和右側的含水量變化規律.單一拉伸下,左側最大流出量(最高峰值)、流入量(最低峰值)為0.020 mm3和0.024 mm3,右側最大流出、流入量分別為0.016 mm3和0.020 mm3.單一旋轉下,左側最大流出、流入量為0.003 mm3和0.010 mm3,右側最大流出、流入量分別為0.004 mm3和0.016 mm3.旋轉拉伸復合加載下,左側最大流出、流入量為0.013 mm3和0.021 mm3,右側最大流出、流入量為0.017 mm3和0.031 mm3.最大流入和流出量的和表示加載過程中髓核與周圍組織物質交換的程度.從上述數值可知,與其他加載相比,旋轉拉伸復合加載下髓核左側平均流速和流量的變換趨勢相反,髓核右側產生了比左側更大的物質交換.

結果表明,3 種加載方式均引起髓核內外的流體交換.拉伸引起流體先流出髓核、再流入髓核.旋轉使髓核左右側流動方向相反,本例右旋加載時,髓核右側半先流出再流入,髓核左側半先流入再流出.這樣,在旋轉拉伸復合加載中,旋轉可增加同側髓核含水量的變化,降低對側髓核含水量的變化.在臨床運用中,可根據髓核病變部位指導旋轉加載方向.

3 結論

脊柱退變疾患最終的解決對策是抑制椎間盤的變性和促使椎間盤組織的再生[55-57],恢復椎間盤基質內的流體流動是啟動髓核細胞再生程序的一個重要目標,中醫正脊治療是恢復椎間流動的有效方法.本工作測量了中醫正脊治療過程中對人體施加的作用力,將該瞬態作用力加載到基于真實人體腰椎CT 掃描數據所建立的有限元模型,研究了作用力引起椎間盤應力?應變和流體運動的流固耦合效應關系.

研究通過對模型加載生理載荷驗證了椎間盤流固耦合有限元模型的有效性.對瞬態力的分解和復合模擬加載表明:旋轉引起纖維環形變較大且出現明顯的應力集中;拉伸引起流體先流出髓核、再流入髓核;旋轉可加快髓核流動,并在髓核左右產生相反的法向流速,即同一時刻既有流體流入,又有流體流出;在旋轉拉伸復合加載中,旋轉可增加同側髓核含水量的變化,降低對側髓核含水量的變化.

本文所創建的研究方法基于流固耦合物理場模型,模擬了椎間盤內應力?應變與流動的相互作用關系,解釋了外部加載引起的椎間盤固體位移形變對流場影響的物理過程,可在臨床運用中指導旋轉加載的方向.后續研究將采用更接近人體組織結構和功能狀態[47,58]的本構關系模型,進一步揭示中醫正脊治療椎間盤退變的生物力學原理,發揮本研究方法在臨床應用中指導個性化治療的優勢.