1. 引言
近年来，近视发生年龄趋于低龄化且全球近视患病率显著增加[1]，至2050年，近视总人数将达到47.58亿，占全球人口的49.8% [2]。因此，近视防控已成为共同关注的社会问题。角膜塑形镜(简称OK镜)是目前国内外广泛应用的一种的物理近视防控方式，据相关临床数据统计，在一年的随访中，OK镜治疗对近视儿童眼轴(AL)延长的抑制作用从30%到60%不等[3] [4]。OK镜属于一种可晚间配戴的物理矫形治疗镜，其采用特殊的反几何设计，通过机械压迫、眼睑按摩以及液压作用，将角膜中央压平，中周围角膜变陡，以此纠正角膜的屈光度达到矫正近视的目的[5]。OK镜对角膜的详细作用机制至今尚未明确，角膜的生理结构和OK镜的设计参数对角膜塑形效果的影响一直是医患和研发人员关注的重点问题。
在角膜生物力学的有限元研究中，角膜模型经历了从相对简化到逐渐精细的演变过程，越来越符合真实人眼。2020年，Lihua等[6]使用椭圆体形状来定义角膜前后表面，定量分析LASIK手术后角膜生物力学性质变化；2020年，Doll等[7]根据角膜地形图数据拟合人群角膜Q值和曲率建立角膜模型，分析角膜形状对屈光力变化的影响；2023年，Wang等[8]通过患者的高度和厚度等测量数据构建了角膜的三维模型，模拟SMILE手术中的角膜变形。
近几年，有研究学者结合OK镜的理论设计参数，定量分析OK镜下角膜的空间位移和应力分布，探究其对角膜塑形的作用机制。2018年，Chang等[9]建立了角膜与RGP镜的三维模型，改变RGP镜与角膜接触处的曲率，提取相应的角膜外周的接触应力，探讨镜片的结构设计与老年人患者佩戴舒适度的相关性。2019年，Abass等[10]针对不同曲率的角膜和软性接触镜进行有限元建模，他们的模型考虑了角膜与泪液的表面张力与镜片的摩擦力，通过计算评估软性隐形眼镜适配度对光焦度的影响；2020年，利物浦大学Doll等[7]对佩戴隐形眼镜时眼的屈光变化进行了有限元建模研究，发现角膜的形态比镜片的周边定位设计对镜片的有效屈光度影响更大。2021年，Wu等[11]建立了球面设计的OK镜–角膜模型，探讨不同结构参数下的OK镜分别对不同曲率、厚度和近视度数的角膜的生物力学反应，并表明角膜的生物力学反应随着角膜参数和镜片的设计改变而发生显著变化。2023年，Zhai等[12]在Wu等[11]的基础上建立了基于非球面设计的OK镜–角膜模型，通过改变近视度数和角膜的屈光度，研究相应的OK镜下角膜前表面整体的位移和应力分布的特点。同年，胡郡琦等[13]建立球面设计的OK镜–角膜模型，通过改变角膜的形态以及角膜和塑形镜的材料参数，分析戴镜后角膜的生物力学响应。然而以往研究多将角膜简化成简单的扁椭球结构，忽略了角膜非球面非对称结构。
本文基于常见的角膜曲率、偏心率和散光度数，建立散光0.5 D下不同角膜平K (40.0, 41.0, 42.0, 43.0, 44.0, 45.0和46.0 D)的OK镜–角膜非对称结构耦合有限元模型，探究OK镜对非对称角膜生物力学的影响。
2. 材料与方法
2.1. 非对称角膜的几何模型建立
角膜的前表面轮廓近似一个扁椭圆形，其在中心有最大的曲率(或最小的半径)，而向边缘有一个更平坦的表面。角膜轮廓的非球面性通常用形状因子(p)来描述，对于扁椭圆形来说，p小于1。根据公司OK镜的理论设计参数和人眼真实散光情况，角膜被假定为形状因子p为0.75，散光度为0.5 D的扁椭圆形。在MATLAB 2020软件中根据非球面公式，建立角膜平K为43.0 D的角膜前表面，如图1所示，并保存为STL格式。
Figure 1. Schematic diagram of the anterior surface profile of an asymmetric cornea
图1. 非对称角膜前表面轮廓示意图
将STL格式的角膜前表面模型导入到Geomagic12.0 (美国Geomagic公司)中进行修整，最后将所得三维模型拟合成曲面片以IGS格式储存，利用ABAQUS 2020生成角膜前表面的网格，并由此创建相应的网格部件，使用“offset”功能生成3层角膜的整个结构化网格。然后对模型各部分赋予材料属性和接触设置，建立包括角膜上皮层、角膜基质层和角膜下皮层结构的三维有限元模型。
2.2. 非球OK镜几何模型的建立
本研究采用的OK镜后表面均采用四段弧VST非球面设计，非球面系数统一设定为0.5。其中基弧(base curve, BC)、反转弧(reverse curve, RC)、定位弧(alignment curve, AC)和周围弧(peripheral curve, PC)的弧长半径分别为3 mm、3.6 mm、4.8 mm和5.3 mm。材料为高透氧硬性材料(BOSTON XO2 nominal Dk of 141 × 10−11 cm3 O2 (cm)/((sec)(cm2) (mmHg))。根据Jesson的研究，OK镜设计中需引入“压缩因子”(亦称Jesson因子或JF因子)以减缓夜间佩戴后日间视力反弹现象。本研究中所有OK镜的压缩因子均设置为+0.75 D。基于角膜屈光度范围(40~46 D)，通过公司验配理论计算出相应的7套OK镜参数，具体数据如表1所示。
Table 1. Geometric parameters of the OK lens
表1. OK镜几何参数
角膜曲率/mm |
RBC/mm |
RRC/mm |
RAC/mm |
RPC/mm |
eBC |
BC的矢高/mm |
水平方向(平K) |
垂直方向(陡K) |
40.0 |
40.5 D |
9.22 |
7.17 |
8.26 |
10.71 |
1.387 |
0.446 |
41.0 |
41.5 D |
8.98 |
7.05 |
8.06 |
10.35 |
1.325 |
0.459 |
42.0 |
42.5 D |
8.74 |
6.93 |
7.88 |
10.00 |
1.278 |
0.473 |
43.0 |
43.5 D |
8.52 |
6.81 |
7.70 |
9.69 |
1.228 |
0.487 |
44.0 |
44.5 D |
8.31 |
6.70 |
7.53 |
9.39 |
1.181 |
0.500 |
45.0 |
45.5 D |
8.11 |
6.59 |
7.36 |
9.11 |
1.138 |
0.514 |
46.0 |
46.5 D |
7.92 |
6.49 |
7.21 |
8.85 |
1.097 |
0.527 |
通过SolidWorks软件的“方程式驱动的曲线”功能，输入各段弧的数学表达式及弧长，生成OK镜二维平面轮廓，然后通过“旋转特征”功能将该二维轮廓绕中心轴旋转，形成完整的三维实体模型，并保存为STEP格式。
2.3. 角膜—OK镜有限元模型的建立
1) 接触属性设置
本研究构建的完整有限元模型由OK镜和角膜两个主要部分组成。为准确模拟泪液层的润滑作用，模型中将OK镜与角膜相互接触的界面设置为无摩擦滑动接触，其摩擦系数设定为0。在接触对的定义中，将OK镜后表面指定为从面(slave surface)，角膜前表面定义为主面(master surface)。为了避免复杂接触分析的不收敛问题，预先增加了两个分析步(在分析步1中加载0.001 mm的位移，在分析步2中卸载前者的位移)以使有限元模型在分析时顺利建立起接触。
2) 材料属性定义
基于生物力学特性分析，本研究对角膜材料模型进行了如下假设：由于非线性材料模型在13 mmHg眼压下表现出与弹性材料模型相似的行为特征，且非线性弹性材料在微小变形范围内可呈现线性变形特性，因此将角膜假设为非均质弹性材料。在该模型中，角膜最外层(上皮层)和最内层(内皮层)的弹性模量设定为基质层的10%。OK镜片的材料参数则依据波士顿材料产品指南进行设定，具体材料参数详见表2所示。
Table 2. Material properties of the finite element model
表2. 有限元模型的材料属性
模型结构 |
泊松比 |
杨氏模量/MPa |
OK镜 |
0.30 |
1160 |
角膜最外层 |
0.49 |
0.127 |
角膜基质层 |
0.49 |
1.270 |
角膜最内层 |
0.49 |
0.127 |
3) 边界条件与载荷
眼内压和眼睑压分别设定为13 mmHg和9 mmHg，均匀地施加在角膜的内表面和OK镜的外表面。为简化模型并模拟巩膜的生物力学作用，在角膜边缘设置了与水平面呈23˚夹角的滑动支撑边界条件。简化模型如图2所示。
Figure 2. Finite element model of the OK lens and cornea
图2. OK镜与角膜的有限元模型
3. 结果
3.1. 塑形后角膜前表面的应力变化
在固定角膜偏心率(e值)为0.5的条件下，通过改变角膜曲率(40.0~46.0 D，步长1.0 D)对角膜前表面应力分布进行了系统分析。图3所示的应力云图结果表明，角膜前表面从中央区到周围区的应力变化呈现显著的非线性特征。随着角膜曲率的增加，角膜前表面中央区的应力分布趋于均匀化。
在不考虑泪液润滑作用的情况下，OK镜佩戴后的角膜前表面应力分布呈现以下特征：最大应力集中区域主要分布于角膜周围区，且在水平和垂直方向上表现出明显的非对称性分布趋势。具体而言，水平方向的最大应力值分别为2938.96 Pa (40.0 D)、3203.27 Pa (41.0 D)、3419.34 Pa (42.0 D)、3036.08 Pa (43.0 D)、3749.45 Pa (44.0 D)、3766.77 Pa (45.0 D)和3685.22 Pa (46.0 D)。相比之下，垂直方向的周围区应力分布较为均匀，未出现明显的应力峰值。
Figure 3. Stress contour plot of the anterior corneal surface for flat K values ranging from 40 to 46 D (unit: MPa)
图3. 平K为40~46 D的角膜前表面应力云图(单位：Mpa)
尽管7种不同结构角膜的前表面应力分布在整体趋势上具有相似性，但在特定位置上的应力值存在显著差异。为定量分析应力分布特征，本研究提取了角膜前表面沿路径A (角膜X轴正方向的水平子午线)的应力分布曲线，并选取曲线上4个特征位置(S1~S4)进行重点分析。图4展示了角膜e值为0.50、曲率范围40.0~46.0 D时，配戴OK镜后角膜沿路径A的应力分布曲线。
Figure 4. Stress distribution curve along path A on the anterior corneal surface for flat K values ranging from 40 to 46 D; S1 to S4 represent the four characteristic points marked by red boxes in the figure
图4. 平K为40~46 D的角膜前表面沿路径A的应力分布曲线；S1~S4为图中红色方框所标出的4处特征点位
分析结果表明，在角膜前表面顶点(X = 0 mm)至极小值S1 (X ≈ 1.44 mm)区间内，各曲率角膜的应力分布曲线基本重合。该区间内，von Mises应力(VMS)平均值在44.0 D曲率角膜中达到最大值，为507.03 Pa。这一发现揭示了角膜曲率对前表面应力分布的显著影响，特别是在角膜中央区域的应力分布特征。
在角膜中周区应力值随径向距离(Di)的增加呈现递增趋势，其分布曲线表现出典型的抛物线形态特征。应力峰值出现在Di ≈ 2.55 mm处(S2点)，各曲率角膜的应力曲线在此区域呈现轻微波动，平均应力值范围为1430.12~1469.11 Pa。
在角膜周围区应力分布表现出显著的双峰特征，且波动幅度较大。两个应力峰值分别位于Di ≈ 3.60 mm (S3点)和Di ≈ 4.8 mm (S4点)处。值得注意的是，角膜前表面von Mises应力(VMS)的最大值均出现在S4点，且具体数值分别为2649.33、3183.53、3419.34、3036.08、3646.19、3568.48和3385.94Pa。这一发现进一步证实了角膜周围区(S4点)是角膜上皮层应力集中最为显著的区域，且最大应力值与角膜曲率之间存在明显的相关性。
3.2. 塑形后角膜前表面的位移变化
在固定角膜e值为0.5的条件下，本研究通过改变角膜曲率(40.0~46.0 D，步长1.0 D)对角膜位移分布进行了系统分析。图5所示的位移云图结果表明，不同曲率角膜在配戴OK镜后的位移分布特征在中周区存在显著位移环(红色区域)，随着角膜曲率的增加，角膜位移分布整体呈现出明显的梯度变化，特别是在角膜中央区和周围区表现出不同的位移分布。这一发现为理解角膜形变机制提供了重要的可视化依据。
Figure 5. Displacement contour plot of the anterior corneal surface for flat K values ranging from 40 to 46 D (unit: µm)
图5. 平K为40~46 D的角膜前表面位移云图(单位：µm)
7种不同结构的角膜前表面表现出相似的位移分布趋势，但其位移值存在显著差异。为深入分析这种差异，本研究提取了不同曲率角膜沿路径A的位移分布曲线，如图6所示。分析结果表明，所有曲率的角膜前表面位移在靠近角膜中心区域较小，随着距角膜中心距离的增加，位移逐渐增大并在2.5 mm至3.5 mm范围内达到峰值，随后开始减小。此外在Di ≈ 2.3 mm至3.5 mm范围内，位移变化显著，不同曲率间的差异逐渐增大，然而大多在Di ≈ 2.47 mm处出现位移峰值，不同曲率的角膜位移峰值分别为28.79、28.02、27.11、26.45、26.32、26.27和25.34 μm；在距角膜中心4.8 mm以外的区域，位移变化较小且各曲率间的差异不明显。这些发现提示，角膜曲率对角膜前表面的力学响应具有重要影响，特别是在角膜中央区域，曲率小(<43 D)的角膜表现出相对更大的位移变化，这可能与OK镜与角膜曲率的适配性密切相关。
Figure 6. Stress distribution curve along path A on the anterior corneal surface for flat K values ranging from 40 to 46 D
图6. 平K为40~46D的角膜前表面沿路径A的应力分布曲线
4. 结论
本文基于有限元方法(FEM)构建了不同曲率条件下非对称角膜及相应OK镜的数值模型，旨在探究OK镜四段弧VST非球面设计参数对角膜塑形后的生物力学响应特性。研究重点关注了角膜前表面位移与应力的空间分布特征，并对特定区域的力学行为进行了对比分析。研究结果表明，角膜与OK镜相互作用后的几何形态变化趋势与临床观察结果一致[14]，主要表现为角膜中央区曲率减小(变平)和中周区曲率增大(变陡)。且中周区的最大应力和位移均出现在距角膜顶点约2.45~2.55 mm处，该区域仍位于OK镜基弧区(0~3 mm)的作用范围内。这进一步验证了OK镜通过改变角膜几何形态实现屈光矫正的生物力学原理。
研究还发现角膜前表面的应力值显著高于课题组先前的研究结果，这一差异可能与模型参数的优化及边界条件的精确设定有关。更为重要的是，研究观察到角膜前表面的应力与位移分布均呈现出明显的非中心对称特征。这种非对称性主要归因于角膜的几何结构特性：角膜并非完全对称的球面结构，其厚度和曲率在不同子午线上存在固有差异。为准确模拟这一解剖学特征，本研究在模型中设置了角膜水平子午线与垂直子午线之间的曲率差异为0.5 D。这种非对称的几何结构直接导致了角膜在OK镜作用下的应力与位移分布呈现非对称特征。这一发现不仅更真实地反映了角膜的生物力学响应特性，也为解释OK镜作用下角膜的非对称形变机制提供了重要的理论依据。
基金项目
国家自然科学基金项目(51735003)，广东省微创手术器械制造技术重点实验室( MISIMT20224)，上海市科委医学创新研究专项(21Y11911600)。
NOTES
*通讯作者。