1. 引言
青少年特发性脊柱侧凸(Adolescent Idiopathic Scoliosis, AIS)是一种涉及脊柱三维结构异常的多因素疾病，其生物力学机制的探究对于改进非手术治疗方案具有重要价值[1] [2]。在临床实践中，支具等非手术干预手段的有效性取决于对脊柱力学特性的准确评估，而有限元分析(Finite Element Analysis, FEA)凭借其无创性和良好的可重复性，已成为研究脊柱力学特征及评估矫正效果的关键技术[3]。该技术能够建立个体化的脊柱三维数字模型，用于预测侧凸脊柱在各类力学载荷作用下的变形特征、应力分布及动态力学响应。通过CT影像数据构建的四面体网格模型，可精确区分皮质骨、松质骨及韧带等结构，进而评估各节段活动范围及应力分布特征。Somtua等人通过有限元方法研究了30~70度Cobb角范围内螺钉固定系统的调整效果，并量化了三维模型中最大等效应变分布情况[4]。Ma等人采用有限元技术系统分析了脊柱侧凸矫正过程中脊髓及神经根的生物力学变化，研究显示在矫正初期，上段脊髓在各种干预条件下均存在损伤风险，而下段脊髓在推力作用下更易受损，同时牵引操作可能对中下段脊髓双侧神经组织造成潜在损伤[5]。
在有限元建模中，网格划分方式的选择直接影响模型的仿真精度与运算效能，其中六面体单元(Hexahedral)与四面体单元(Tetrahedral)的差异最为显著，这种差异不仅反映在应力场分布的仿真精度上，更显著影响数值求解过程的收敛特性[6]。网格类型的选择不仅影响仿真结果的可信度，更直接决定了生物力学研究的临床指导价值。四面体网格虽然在贴合复杂解剖形态(如椎弓根和关节突)方面具有优势，但在应力梯度较高的区域(如韧带附着点)易出现数值不稳定的问题。相反，六面体网格常需对复杂结构进行几何简化，虽能提升计算效率，却可能因过度理想化而忽略关键的局部力学信息。在有限元分析中，研究者通常会从应力及自由度两个角度去判断模型的有效性。例如Nie等人在研究中基于自由度判断[7]。Guan等人则在此基础上增加了应力收敛的测试，即网格数量与应力及位移的测试[8]。
本研究针对上述关键问题，系统比较六种实体建模方法对椎体节段生物力学性能的影响，旨在为AIS有限元建模提供兼顾仿真精度与计算效率的优化策略。
2. 材料与方法
本研究纳入一名28岁健康男性受试者(身高180 cm，体重75 kg，BMI 23.1 kg/m2)，采用联影医疗uCT 960 + 全脊柱CT扫描系统进行数据采集。扫描参数设置为：管电压120 kV，管电流250 mA，扫描层厚1 mm，共获取805层连续断层图像。本研究方案已通过机构伦理委员会审查批准(伦理批准编号：NO.XHEC-D-2023-034)。数据处理平台为HP OMEN 8台式计算机，操作系统为Windows 11专业版，配备第12代Intel(R) Core(TM) i9-12900K处理器及64 GB内存。该计算平台为后续图像重建和有限元分析提供了必要的运算性能保障。受试者在扫描前已签署知情同意书，所有操作均遵循赫尔辛基宣言的伦理准则。
本研究采用基于CT值的阈值分割技术进行骨组织提取，具体方法如下：首先设定226-3071 HU作为骨组织的灰度阈值范围，通过该阈值区间生成初始骨组织掩膜。随后采用区域生长算法实现各椎体的自动分割处理。所得骨组织的.stl文件在Geomagic软件中进行网格修复，主要包括钉状结构、自交边的识别与修补，并对椎体间接触区域进行校验。通过“精确曲面”模块中的轮廓线编辑、曲面片构建、格栅构建及曲面拟合功能，构建非均匀有理B样条(NURBS)曲面模型，并导出为.iges格式，得到初始椎骨集合模型。整体建模流程如图1所示。
Figure 1. Overview of the modeling workflow. (a) CT images of the subjects; (b) Reconstruction of STL model; (c) Entity model generation; (d1) Extraction of cortical bone; (e1) Modeling of annulus fibrosus, nucleus pulposus and articular cartilage. (c2) Finite element model of the vertebral body; (d2) Finite element model of the intervertebral disc; (e2) Finite element model of the vertebral segment
图1. 建模工作流程概览。(a) 受试者CT图像；(b) STL模型重建；(c1) 实体模型生成；(d1) 皮质骨提取；(e1) 纤维环、髓核及关节软骨建模。(c2) 椎体有限元模型；(d2) 椎间盘有限元模型；(e2) 椎体节段有限元模型
本研究设计了六种不同的有限元建模方案，各方案的关键特征及参数设置如下：
方案1：皮质骨采用实体单元建模，椎骨不区分前后部。Geomagic逆向建模完成后导入ANSYS Workbench，自动划分网格；椎间盘采用四面体单元建模。
方案2：皮质骨仍采用实体单元，椎骨不区分前后部。Geomagic建模后导入ANSYS Workbench，设置不同面网格参数，椎骨皮质骨面网格大小为1 mm，椎间盘及软骨为0.5 mm，椎间盘仍采用四面体单元。
方案3：皮质骨采用壳单元建模，厚度设为1 mm，椎骨不区分前后部。模型导入HyperMesh，皮质骨最大网格尺寸为2 mm，松质骨通过皮质骨面网格生成四面体单元，椎间盘和软骨通过solidmap工具生成六面体单元。
方案4：皮质骨采用壳单元建模，且区分椎骨前后部。Geomagic建模后导入HyperMesh，皮质骨最大网格尺寸设为1 mm，松质骨通过皮质骨生成四面体单元。椎间盘与髓核的建模方式与方案3相同。
方案5：与方案4基本一致，但软骨采用圆柱梁单元、半径2.5 mm。
方案6：与方案5基本一致，但不建立软骨，椎间盘统一一种材料参数。
研究采用弹簧单元模拟韧带组织，其材料参数详见表1 [8]-[10]。模型边界共条件设置四组，如图2所示。软骨–皮质骨界面设为绑定接触，软骨–软骨界面设为不分离接触；约束条件为固定T4椎体下
Table 1. Material properties and element types used in the FE models
表1. 有限元模型采用的材料属性与单元类型
|
弹性模量Mpa |
泊松比 |
刚度N/mm |
单元 |
皮质骨(方案1，2) |
12,000 |
0.3 |
|
Solid185 (四面体) |
皮质骨(方案3，4) |
12,000 |
0.3 |
|
Shell 28 |
椎骨后部 |
3500 |
0.3 |
|
Solid185 (四面体) |
松质骨 |
100 |
0.3 |
|
Solid185 (四面体) |
纤维环 |
4.2 |
0.4 |
|
Solid185 (六面体) |
髓核 |
1 |
0.49 |
|
Solid185 (六面体) |
关节囊软骨 |
50 |
0.3 |
|
Solid185 (六面体) |
前纵韧带 |
|
|
8.74 |
Spring |
后纵韧带 |
|
|
5.83 |
Spring |
棘间韧带 |
|
|
15.38 |
Spring |
黄韧带 |
|
|
10.85 |
Spring |
横突间韧带 |
|
|
0.19 |
Spring |
软骨关节(方案5) |
50 |
0.3 |
|
Beam |
椎间盘(方案6) |
100 |
0.3 |
|
Solid185 (六面体) |
Figure 2. Boundary conditions of the FE model. (a) Fixed constraints on the lower surface of T4; (b) A vertical compressive load is applied to the upper surface of T1; (c) Flexion/extension moment; (d) Lateral bending moment; (e) Axial rotational torque
图2. 有限元模型的边界条件。(a) T4下表面固定约束；(b) T1上表面施加垂直压缩载荷；(c) 屈曲/伸展力矩；(d) 侧向弯曲力矩；(e) 轴向旋转力矩
表面，同时在T1椎体上表面施加1000 N垂直压缩载荷，用以分析不同方案下节段的位移峰值。
在自由度(Range of Motion, ROM)验证中，对T4椎体底部进行固定(Fix)，在T1椎体顶部依次施加横向、纵向及轴向两个方向的4 N·m力矩，以模拟T1~T4节段的前屈/后伸、侧屈及轴向旋转运动。最终通过获取T1椎体顶部的可动范围来评价各建模方案的运动学表现。
3. 结果
不同方法绘制的模型如图3所示，单元与节点数量如表2所示。
Figure 3. FE models constructed using different strategies
图3. 采用不同建模策略构建的有限元模
Table 2. Number of units and nodes for different modeling methods
表2. 不同建模方法的单元与节点数量
|
单元数量 |
节点数量 |
方案1 |
47,266 |
93,280 |
方案2 |
1,130,098 |
1,745,076 |
方案3 |
76,138 |
37,205 |
方案4 |
260,603 |
86,638 |
方案5 |
255,581 |
79,342 |
方案6 |
255,581 |
79,342 |
如表3所示，网格质量评估结果显示，方案1在平均评分和最低网格质量评分方面均为最差，方案3的综合表现最佳。尽管方案2的平均评分较高，但其最低评分低于方案4。进一步分析表明，方案1与方案2的低评分主要集中在皮质骨区域，而方案3与方案4的低评分则主要分布于部分椎间盘区域，其原因在于该处六面体单元存在较低的纵横比(Aspect Ratio)，即网格扁平化问题。
在1000 N垂直载荷作用下，六种有限元建模方案呈现出不同的位移响应特征，如图4所示。各方案的最大位移量依次为：方案1 (2.96 mm)、方案2 (3.00 mm)、方案3 (2.81 mm)、方案4 (2.76 mm)、方案5 (2.91 mm)和方案6 (0.97 mm)。位移云图分析表明，所有模型的整体位移分布模式具有一致性，具体表现为：
Table 3. Mesh quality scores evaluated in ANSYS Workbench
表3. ANSYS Workbench评估的网格质量评分
最低评分 |
最高评分 |
平均分 |
标准差 |
最低评分 |
方案1 |
1 |
0.54 |
0.19 |
0.03 |
方案2 |
1 |
0.81 |
0.11 |
0.06 |
方案3 |
1 |
0.82 |
0.13 |
0.23 |
方案4 |
1 |
0.78 |
0.13 |
0.13 |
方案5 |
1 |
0.77 |
0.12 |
0.13 |
方案6 |
1 |
0.77 |
0.12 |
0.13 |
Figure 4. Displacement distribution under vertical compressive load. (a) Model 1; (b) Model 2; (c) Model 3; (d) Model 4; (e) Model 5; (f) Model 6; (g) Comparison of peak displacements among the six models
图4. 垂直压缩载荷作用下的位移分布。(a) 模型1；(b) 模型2；(c) 模型3；(d) 模型4；(e) 模型5；(f) 模型6；(g) 六种模型的峰值位移对比
(1) T1椎体呈现整体下沉位移特征，各部位位移矢量均指向下方；
(2) T2和T3椎体除主要呈现轴向位移外，还表现出前屈运动趋势，导致椎体前缘位移量较后缘显著增加(p < 0.05)，形成从前向后递减的位移梯度分布。
这种位移分布差异可能源于以下生物力学机制：
(1) T1椎体由于皮质骨的高刚度特性(弹性模量约12 GPa)，倾向于保持整体协同运动；
(2) T2~T3椎体间存在椎间盘结构(弹性模量约3.5 MPa)，其粘弹性特性为该节段提供了更大的运动自由度，从而产生相对滑移和旋转的复合位移模式。
在4 N·m力矩作用下模拟T1~T4节段前屈运动时，各建模方案呈现出不同的生物力学响应特征。T1椎体位移峰值依次为：方案1 (3.19 mm)、方案2 (3.23 mm)、方案3 (3.39 mm)、方案4 (1.21 mm)、方案5 (1.75 mm)和方案6 (1.47 mm)；对应的前屈活动度(Range of Motion, ROM)分别为5.2˚、5.3˚、5.1˚、5.8˚、6.2˚和5.5˚，如图5所示。虽然各方案间的位移峰值差异显著(最大差异达2.18 mm)，但ROM差异相对较小(最大差异仅1.1˚)，表明位移幅值与ROM之间不存在直接相关性。
Figure 5. Displacement and range of motion distribution in flexion/extension movements. (a) Displacement distribution of Model 1; (b) Displacement Distribution of Model 4; (c) Comparison of flexion/extension range of motion among six models
图5. 屈曲/伸展运动中的位移及活动度分布。(a) 模型1位移分布；(b) 模型4位移分布；(c) 六种模型屈曲/伸展活动度对比
进一步分析发现，不同方案在局部位移特征上存在明显差异。例如，方案1中椎骨棘突远端位移达到2.4 mm，而方案4仅为1.07 mm。这种差异可能源于以下因素：
(1) 材料属性分布方式的差异；
(2) 网格划分策略的不同；
(3) 边界条件的处理方式。
在本研究中方案4虽然整体位移较小，但其位移倾斜率较高，这一特征可能解释了为何该方案的ROM值(5.8˚)反而高于部分位移更大的方案(如方案3的5.1˚)。这一现象表明，在评价有限元模型性能时，需要综合考虑位移分布特征和活动度等多个生物力学指标。
Figure 6. Displacement and range of motion analysis in lateral bending motion. (a) Displacement distribution cloud map of Model 3; (b) Displacement distribution cloud map of Model 4; (c) Comparison of the lateral bending range of motion of the six groups of models
图6. 侧向弯曲运动中的位移与活动度分析。(a) 模型3位移分布云图；(b) 模型4位移分布云图；(c) 六组模型侧向弯曲活动度对比
在4 N·m力矩作用下模拟T1~T4节段左侧屈曲运动时，各有限元建模方案表现出差异比较明显的生物力学特性(图6)。具体而言，T1椎体位移峰值测量结果为：方案1 (5.79 mm)、方案2 (5.85 mm)、方案3 (4.98 mm)、方案4 (1.21 mm)、方案5 (6.47 mm)和方案6 (4.75 mm)；相应的侧屈活动度(ROM)分别为6.4˚、6.5˚、6.4˚、6.6˚、6.2˚和7.0˚。
位移分布特征分析揭示：
(1) 方案1呈现双侧横突对称性位移特征，表明其材料属性分布较为均衡；
(2) 方案4则显示出右侧横突位移(1.38 mm)显著大于左侧(1.05 mm)的非对称特性(p < 0.05)。
这种差异可能归因于：
(1) 材料参数的空间分布差异；
(2) 网格划分的非均匀性；
(3) 边界条件的非对称处理。
方案4的位移峰值最小(1.21 mm)，但其ROM值(6.6˚)却高于多数其他方案，这可能与其独特的位移分布模式(较高的位移梯度)相关。这一发现表明，在评估脊柱侧屈生物力学特性时，不仅需要考虑位移幅值，更应关注位移分布特征与整体活动度的耦合关系。
Figure 7. Displacement and range of motion characteristics in axial rotational motion. (a) Displacement distribution cloud map of Model 5; (b) Displacement distribution cloud map of Model 6; (c) Comparative analysis of the axial rotational range of motion of six groups of models
图7. 轴向旋转运动中的位移与活动度特征。(a) 模型5位移分布云图；(b) 模型6位移分布云图；(c) 六组模型轴向旋转活动度对比分析
在4 N·m扭矩作用下模拟T1~T4节段轴向左旋运动时，各有限元模型呈现出不同的旋转特性，如图7所示。具体测试数据显示：T1椎体位移峰值依次为方案1 (4.50 mm)、方案2 (4.52 mm)、方案3 (3.10 mm)、方案4 (3.28 mm)、方案5 (3.52 mm)和方案6 (3.59 mm)；相应的轴向旋转活动度(ROM)分别为7.6˚、7.8˚、7.8˚、7.3˚、7.5˚和7.3˚。
位移特征分析表明：
(1) 方案6表现出良好的双侧对称性，其左右横突位移差小于5%；
(2) 方案5则显示出明显的非对称特征，右侧横突位移较左侧增加约12% (p < 0.05)。
这种差异主要源于：
(1) 软骨组织的力学特性差异；
(2) 椎间关节接触条件的非对称设置；
(3) 旋转中心的偏移效应。
在此项研究中，方案3~6的位移峰值相近(3.10~3.59 mm)，但其ROM值仍存在0.5˚的差异，这表明除位移幅值外，旋转中心的相对位置变化也是影响轴向旋转活动度的重要因素。这一发现为理解脊柱旋转运动的生物力学机制提供了新的见解。
在计算效率方面，方案6运算速度最快，其次为方案1和方案4，方案2的计算时间最长。
4. 讨论
本研究基于六种不同的建模方法，对椎骨建模的有效性进行了验证，对比以往尸体标本的研究可以发现，6种建模方式在自由度方面均是有效的[11]-[13]。方案1与方案2主要旨在探讨网格质量与大小对整体位移的影响。从结果来看，网格尺寸对位移与自由度的影响较小，但可观察到网格尺寸越小，位移值略有增大。其可能原因在于较小的网格尺寸降低了接触面穿透的风险，能够更精确地捕捉椎间盘与皮质骨之间的接触界面，从而避免接触不稳定的现象。然而，网格尺寸过小也可能导致单元数量激增，进而显著延长运算时间。
方案1与方案2均涉及逆向重建，因此建模准备时间较长[14]。为提高建模效率，我们参考了Tomasz的研究，采用壳单元替代皮质骨，以加快三维重建过程[15]。该方法仅需划分二维表面网格，厚度通过属性定义，无需沿厚度方向进行网格划分。考虑到皮质骨厚度远小于椎体的长度与宽度，属于典型的薄壁结构，此建模策略在一定程度上具有合理性。然而，其主要局限在于壳单元忽略了厚度方向的网格划分，可能导致接触计算不稳定或出现接触穿透现象，亦无法有效捕捉厚度方向的复杂应力梯度，特别是在接触区域或应力集中区域易丢失关键力学信息。
方案4则基于Guan等人的研究，对椎骨不同分区分别赋予了不同的材料属性[8]。该方法的特点在于，尽管其位移峰值与其他方案存在一定差异，但在自由度方面表现接近，显示出较好的建模一致性。其潜在优势在于对软骨区域与椎骨后部复杂结构之间接触关系的更精细模拟。在材料参数设定方面，皮质骨为12,000 MPa，后部结构为3500 MPa，软骨为50 MPa。由于后部区域结构较为复杂，相较于壳单元建模或方案1与方案2中通过向内偏移获得的简化形态，该方法能够更真实地还原后部几何特征，同时提高接触计算的准确性。因此，该方法在处理椎骨后部与软骨接触关系方面可能具有一定优势。但其在矫形过程中的全局有效性仍有待考究。
方案5基于Bavil等人的研究，采用Beam单元简化了软骨[14]。该建模方法也较为简单，省略了软骨的建立，但是由于梁单元采用了绑定的算法，虽然可承受弯矩、剪力和轴向力，但无法模拟滑动，因此在ROM上与其余算法有些区别，可能无法抵抗非轴向移动，导致模型“软化”。整体来看仍需要调整。
方案6基于盛文倩的研究，对椎间盘采用了统一的材料参数[16]。该方案建模简易，通过椎间盘材料参数的增加，简化了软骨和椎间盘的总和作用，该方法还存在于简易的圆柱体的简化椎体。除了在压缩载荷的自由度上有着较大的差异，自由度上有较好的相似性，是一种很好的简易初实验方式。
除上述方法外，仍有多种建模策略被提出。例如，Karimi研究基于经验公式赋予皮质骨材料参数[17]。虽然该方法被很多骨科生物力学中使用，但其准确性涉及经验公式的材料数量及网格大小，且容易出现应力集中。部分学者对椎间盘材料的柔韧性进行建模并指出其对自由度有一定影响，虽然椎间盘柔韧度已经被确定会影响侧弯患者的矫正效果，但在验证阶段，尚无研究深入讨论验证其是否研究有效性[18]。在建模手段方面，还包括采用杆单元模拟椎骨[19]、忽略椎弓根的圆柱体椎骨[20]、采用全六面体进行简化建模[21]、在涉及颈椎的研究中，还会为椎间盘中引入纤维环结构并分层赋予不同材料参数等[22]。这些方法在部分文献中已获得良好的自由度验证，尽管仍需进一步探讨其生物力学合理性。部分研究通过网格无关性的应力结果验证模型有效性，但对于复杂结构而言，使用全六面体网格会显著增加建模与运算工作量。因此，大多数文献更倾向于通过自由度对比验证模型有效性，若差异在可接受范围内，即认为该模型具备一定的可行性。
5. 结论
六种建模方式均具备良好的位移预测可靠性。然而，方案6的压缩载荷与其余方案存在一定差异，需进一步探讨其潜在原因。采用六面体单元对椎间盘及软骨建模，可有效降低整体网格数量，同时保持较高的平均网格质量，从而显著提升计算效率。相比之下，将皮质骨建模为壳单元较实体单元建模更为简便，尽管可能对局部数值结果产生一定影响，但整体自由度仍保持在合理范围内，具备可接受的生物力学一致性。
基金项目
本项目由国家重点研发计划(项目编号：2023YFC2507702)资助。
NOTES
*通讯作者。