人体脊椎腰骶段的三维有限元仿真模拟0724

1概述有限元方法是根据变分法原理求解数学上可描述的物理问题的一种数值计算方法，其基本思想是用简单问题代替复杂问题然后再求解，它将求解域看成是由许多称为有限元的小的互连子域组成，对每一单元假定一个合适的近似解，然后推导求解出这个域总的满足条件，从而得到问题的近似解。与以往的数值计算方法不同的是，有限元方法能够解决结构、材料性质和载荷情况相对较为复杂的问题。由于人体脊椎在结构形状、材料特性以及承载能力等方面都比较复杂，以往的研究手段都难以获得全域性信息，而有限元方法可用任意形状的网格来分割区域，布置节点，尤其是对脊椎这样复杂的结构有较好的适应性。本文将在Simpleware软件的辅助下建立人体脊椎腰骶段L3-L4段的三维有限元网格模型，然后给出实际情况中有限元应力应变的计算方法，最后利用ANSYS软件对得到的有限元网格模型进行受力仿真模拟，并比较验证了腰椎间盘突出情况下脊椎的应力应变分布。2有限元网格模型的建立本文采用的数据是65帧连续的DICOM格式的人体脊椎矢状位图像，其中图像序列间隔为1.25mm，像素间隔也是1.25mm，分辨率为256x256，如图1所示。基金项目：国家自然科学基金资助重点项目；首先将DICOM图像导入到专业的医学图像处理软件Simpleware中，然后对二维DICOM数据进行手动分割，得到人体脊椎腰骶L3-L4段椎体外部、椎体内部、椎间盘、间盘内部的二维轮廓，通过三维重建得到与真实椎体高度相似的三维模型，最后对模型进行网格划分，进而得到人体脊椎腰骶段L3-L4段的三维有限元体网格模型，为下一步的有限元力学分析提供仿真模型，图2和图3分别是三维有限元模型的单元和节点表示。图1人体脊椎矢状位图像人体脊椎腰骶段的三维有限元建模与仿真摘要：人体脊椎的有限元仿真计算能够精确描述其在被施加集中力载荷后的应力应变情况，进而帮助验证病变手术结果的有效性。本文首先建立了人体脊椎腰骶段L3-L4段的三维有限元网格模型，然后给出了有限单元应力应变的计算方法，最后利用有限元分析软件ANSYS模拟了正常情况和椎间盘突出情况下脊椎的受力情况，并分析了结果的准确性和有效性。关键词：脊椎，有限元，三维模型TheThree-dimensionalFiniteElementSimulationofHumanSpinalLumbarsegmentAbstract:Thefiniteelementcalculationofthehumanspinecandescribethedeformationofthespineaccurately,whichcanbeusedtotesttheeffectivenessoftherelevantoperation.Inthispaper,thethree-dimensionalfiniteelementmodelofhumanL3-L4lumbarsegmentisbuiltwiththesoftwareSimplewarefirstly,thenthecalculatingmethodofstressanddeformationofthefiniteelementisprovided.Finally,thestressanddisplacementinthenormalcaseandtheprotrusionoflumbarvertebraldisccaseissimulatedandcomparedwiththesoftwareANSYSwhichisusedforthefiniteelementanalysisspecially,andtheaccuracyandeffectivenessoftheresultareanalysedattheendofthepaper.Keywords:Spine,FiniteElement,Three-dimensionalModel图2单元表示图3节点表示在有限元网格模型输出时，需要赋予各结构材料不同的科学参数：弹性模量（杨氏模量）和泊松比，从而将椎骨各部分简化为各项同性的弹性材料。其中椎体外表面单元用皮质骨描述，内部单元用松质骨描述，椎间盘纤维环基质也可以简化为各向同性弹性材料，现将其设置为可伸缩的弹簧单元，而将髓核假定为一种不可压缩的固体材料。图3中浅蓝色节点表示的是脊椎的L3段皮质骨，红色节点表示的是脊椎的L4段皮质骨，深蓝色节点表示的是脊椎L3、L4段内的松质骨，黄色节点表示的是两段椎骨间的纤维环，紫色节点则表示的是纤维环内的髓核。表1给出了输出的有限元网格模型各结构材料所包含的节点、单元数量和材料参数值。表1各结构单元数量和材料参数3有限元的应力应变计算将得到的网格模型导入到专业有限元分析软件ANSYS中，为了简化计算本文仅采用SOLID45这个八节点六面体模型来模拟所有的四面体和六面体计算单元。3.1八节点六面体单元的应力应变计算在实际的计算中，需要将八节点六面体单元映射为正方体单元，二者之间的坐标变换关系可以用局部坐标表示的形函数][N来表示：TTzyxzyxNzyx]...][[][888111，其中1281281280000......00[]0000......000000......00NNNNNNNNNN，,,(1,2,...,8)iiixyzi表示实际单元的节点坐标。用u，v，w分别表示节点沿x轴，y轴和z轴方向的位移变化，若令[][]TUuvw，111888[][...]Tuvwuvw，则八节点六面体单元的应变位移为：[][][]UN。将节点的应力状态用六个分量表示：[][]Txxyyzzxyyzxz，其中xxyyzz是正应力，xxyyzz是剪应力；节点的应变状态也用六个分量表示：[][]Txxyyzzxyyzxz，其中xxyyzz是正应变，xxyyzz是剪应变，则节点应力可以通过由弹性模量E和泊松比构成的三维弹性矩阵[D]和节点应变的乘积得到：[][][]D，其中组成部分节点数六面体数四面体数单元总数弹性模量泊松比皮质骨L32090168475299159838120000.3皮质骨L42262577625594963711120000.3松质骨12553641822785292031000.2纤维环661626581338516043920.45髓核1835915244233571.00.499100012121210001212121000121212[]1100000210000021000002ED节点的应变和位移的关系为：xxux，yyvy，zzwz，xyuvyx，yzvwzy，xzuwzx，则有[][][]B，其中81818188118811881100...0000...0000...00[]0...00...00...0NNxxNNyyNNzzBNNNNyxyxNNNNzyzyNNNNzxzx根据虚位移原理就可以计算出三维八节点等参单元的刚度矩阵，即：[][][][]TVKBDBdxdydz，其中，V表示单元的体积。对于负荷矩阵的计算，因为所采用的是集中力载荷，所以只需要把该负荷放在某个适当的节点，并使它的负荷分量与节点的坐标方向保持某种一致，就可以获得相应的负荷矩阵。3.2面面接触和边界条件模拟结构之间的接触分析对于有限元的受力分析也是至关重要的，尤其是对于人体脊椎这样的复杂结构，某个椎管的受力都会引起其他椎管的连锁反应。ANSYS软件支持刚体柔体的面面接触单元，刚性面被当做目标面，柔性体的表面被当做接触面，本文采用TARGE170来模拟三维目标面，采用CONTA173来模拟三维接触面。其中，CONTA173表示的是三维四节点的低阶四边形单元，位于三维实体的表面。在之前得到的三维有限元网格模型中，各个组成部分的接触部分是由节点表示的，而非连接表面。在有限元理论中，边界条件与自由度约束是对应的，椎体的三维运动可以用6个自由度表示，包括3个线变量确定椎体上的某点的坐标，3个角变量确定椎体的空间方向，所以在有限元分析时需要添加自由度约束，取模型表面的关键点来施加集中力载荷，进而得到模型应力应变分布云图。4结果与讨论4.1L3-L4的应力应变分布通过施加集中力载荷来模拟正常的生理应力应变情况，在L4段定义位移自由度约束为0，在L3表面定义集中力载荷，方向为FZ，大小为500N。图4描述了L3-L4应力之后的形状变化，图中白色线框表示的是受力之前L3-L4的形状，而蓝色表示受力之后L3-L4的形状，从中可以看出椎管向左有明显的倾斜；而图5和图6则描述了L3-L4应力之后的位移变化分布，图中红色表示应变位移最大处，值为0.563867，蓝色表示的是应变位移最小处，值为0。图4向右侧弯的形状变化图5向右侧弯应变分布后面观图6向右侧弯应变分布右面观4.2椎间盘突出情况模拟验证为了验证有限元计算方法结果的有效性，本文模拟了腰椎间盘突出时的受力情况，腰椎间盘突出症是指椎间盘的纤维环和髓核组织突出，压迫和刺激神经根所引起的一系列症状，图7是实际的椎间盘突出症状图，图8是模拟的三维有限元网格模型。表2比较了正常情况和间盘突出情况下各组成部分有限单元的数量，可以发现其主要区别在于纤维环和髓核的单元数量。图7腰椎间盘突出症状图8腰椎间盘突出模拟模型表2正常情形和椎间盘突出情形下的数据比较在正常腰椎L3-L4的节点15407、15567、15761处设置边界条件约束，自由度值为0，在节点36630处施加集中力载荷，方向为FX，大小为-500，力矩为MX，值为-15；在间盘突出的模型15545、15704、15907节点处设置边界条件约束，自由度值为0，在节点3834处施加集中力载荷，方向为FX，大小为-500，力矩为MX，值为-15，两个模型上施加边界条件约束和集中力载荷的点是一一对应的，图9给出了在后伸情况下两个模型的受力仿真比较。（a）正常后伸形变L3L4松质骨纤维环髓核正常情形598386371129203160433357异常情形598946364029211158815018（b）异常后伸形变（c）正常应变分布（d）异常应变分布图9两个模型的应变情况比较表3和表4给出了两个模型在后伸情况下的最大位移数据和平均应力数据比较，分析数据可得，在间盘突出的情况下腰椎的应力和应变都要比正常情况下的严重，这样所造成的结果就是当病人弯腰时间盘向后的严重突出会导致脊神经被压迫，从而引起一系列的症状。表3最大位移数据比较表4平均应力数据比较5总结与展望本文利用Simpleware软件建立了人体脊椎腰骶段L3-L4段的有限元网格模型，然后给出了该模型基于有限元分析软件ANSYS的受力仿真模拟的结果，最后模拟了椎间盘突出情况下的受力情况，比较了两种情况下的应力应变数据，从而验证了方法的有效性。然而，在有限元模型建立的过程中，我们只是模拟了真实情况下人体脊椎的各个组成部分，与真实情况会有些偏差，继而影响有限元分析结果的准确性，在以后的研究工作中，希望能在这方面有所改进。X方向Y方向Z方向最大位移正常情况0.0890240.1529810.0477600.607375异常情况0.1095540.1945010.0527900.692577X方向Y方向Z方向VonMises正常情况-7.578-18.084-14.628324.966异常情况-56.979-20.653-63.533414.327