角膜为生物软组织材料，在眼的屈光系统中起着重要作用。角膜的生物力学是指角膜在受到外界压力时相应的变形情况[1]。当其生物力学性能发生变化时，角膜会影响到眼球整体的屈光能力。因此，深入研究角膜的力学性质，有助于诊断及治疗角膜的损伤和疾病，从而预防及保护视力。在对角膜生物力学性能的探究中，由于人的角膜不容易获取，而猪的角膜不仅在取材上相对容易，其生物力学性能与人角膜的生物力学性质也较为相似[2]。因此，常常采用猪角膜作为人角膜的替代材料进行研究[3-4]。
角膜离体实验包括角膜轴向拉伸实验、角膜膨胀实验、离体全眼球膨胀实验等；膨胀实验方法[5]既可以保证角膜组织的完整性，又能模拟眼内压的压力变化，弥补了轴向拉伸实验的不足。离体全眼球膨胀实验中角膜的生理环境及角膜的边界约束与活体中的角膜更为相似，但实验过程较为复杂。ENDERSON et al[6]建立了三维球壳角膜模型，根据角膜膨胀实验数据推导出了角膜非线性力学性能的计算公式。BAO et al[7]以整体兔眼球为研究对象，进行了全眼球膨胀实验，记录和分析了当角膜的眼内压变化时，角膜的厚度及角膜迟滞、角膜阻力因子、变形幅度等12个生物力学指标的变化情况。
随着对角膜材料性能认识的深入以及有限元软件中材料本构的不断完善，采用有限元软件能更准确地对角膜实验进行仿真计算，对角膜实验进行更为详细的分析。HOELTZEL et al[8]考虑了角膜的几何和材料非线性，对角膜膨胀实验和角膜拉伸实验进行了模拟，分析角膜的原始曲率对应力分布的影响。HOLZAPFEL et al[9]通过对角膜微观结构的深入研究，提出了将角膜视为带有两族主方向的纤维材料，建立了符合角膜微观结构的本构公式，此公式也适用于动脉壁模型的计算。
本文分别从整体眼球膨胀实验及数值模拟计算两方面探讨猪角膜的生物力学性能。首先，完善眼球整体膨胀实验的实验方法，对角膜进行整体眼球膨胀实验，获得实验中猪角膜的顶点位移-眼内压之间的关系及角膜的应力-应变曲线。其次，在ABAQUS中采用考虑角膜的微观结构的Holzapfel本构模型，将角膜视为各向异性材料对角膜进行仿真计算。分析角膜组织在不同眼压作用下的力学性能，并讨论Holzapfel计算模型中其他参数的变化对角膜的力学性能的影响。
1 猪眼球整体膨胀实验
1.1 实验取材及前期准备
在当地屠宰场选取4月龄、品种为长白的健康猪眼球。对眼球表面进行处理，如图1(a)、(b)所示，只留下完好的巩膜和角膜。沿着巩膜的轮廓涂抹医用胶，用以保护眼球结构并在其外层涂抹塑形性能较好的医用石膏对巩膜进行塑形和固定，以确保在实验过程中巩膜对角膜的位移没有影响，如图1(c)、(d)所示。
图1 眼球处理后的实物图
Fig.1 Physical maps of eyeball was treated
1.2 实验方法
此实验方法主要分为以下4个步骤：
1) 首先，将激光位移传感器固定在膨胀实验装置上，并开启激光位移传感器，预热30 min，然后将激光位移传感器与笔记本电脑连接，开启LK-Navigator2软件，调节软件中的相关选项。
2) 将U型管里注满PBS溶液[10]，PBS溶液用于模拟眼球内的房水，并确保整个U型管里无气泡后，将输液管用夹子夹闭，玻璃管连同玻璃管支架固定在Instron试验机上，将自制的黑色浮标放入玻璃管内。
3) 将准备好的眼球试样固定在膨胀实验装置上，取下输液管上的夹子，输液管针头于角巩膜缘处扎入眼球内，此时，实验系统与眼球连通。
4) 用Instron5544试验机为眼球提供连续的眼内压并且控制眼球的眼内压改变速度，使角膜发生形变。具体操作如下：启动Merlin软件，设置Instron5544试验机加载速率为500 mm/min，上升高度为0～400 mm，眼球内的眼内压变化范围为0～4 kPa.
由于角膜为黏弹性材料，眼球试样首先进行3次[11]预循环实验，预实验结束后，测量角膜尺寸并找到眼球的顶点做好标记，将激光位移传感器调至与角膜顶点的垂直距离为40～60 mm范围内，开始正式实验。实验过程中采用Merlin软件控制眼球内的眼内压变化，采用LK-Navigator2软件记录角膜顶点位移随眼内压变化数据。整个实验过程确保角膜在无褶皱状态下且在无震动的环境下完成，实验在眼球取回实验室后3 h之内结束。
2 实验结果及讨论
2.1 猪眼球整体膨胀实验中角膜顶点位移与眼内压曲线
9组猪眼球整体膨胀实验得到的角膜顶点位移-眼内压曲线，如图3(a)所示，并且计算出9组实验的位移平均值，得出角膜顶点平均位移-眼内压关系曲线及误差，如图3(b)所示。
2.2 角膜的弹性模量与应力、应变
根据板壳理论，将角膜视为球壳模型用于计算猪角膜生物力学性能，通过上述实验中的角膜顶点位移与眼内压的关系得到角膜的割线弹性模量E和应变ε，进而绘制角膜顶点应力-应变曲线图。
计算公式如下：
(1)
(2)
(3)
(4)
σ=Eε.
(5)
(6)
式中：T1为角膜的中央厚度；T2为角膜缘厚度；H为角膜前房高度；Ra为角膜横半径的平均值，采用游标卡尺可以测量得出，如表1所示；由T1，H及Ra根据勾股定理可计算得到Rb及θ，t为角膜的平均厚度；R为角膜的曲率半径；p为施加的压力；r为顶点位移；ν为泊松比；β为壳的细长度模量。根据式(1)-(6)及表1测量值可得角膜的弹性模量、应变及应力。角膜顶点的应力-应变曲线及弹性模量-眼内压的平均值如图4所示。
表1 相关参数
Table 1 Related parameters
参数T1/mmT2/mmH/mmRa/mm数值0.80±0.551.03±0.944.27±0.406.97±0.52
2.3 结果及讨论
从图3可以看出，在这9组实验数据中，初始阶段角膜的位移随眼内压的变化增长较快，曲线分布较为集中。眼内压在2.04 kPa以后，9组角膜位移变化较为离散，经过一段非线性变化后，角膜顶点位移随眼内压变化较小，曲线逐渐变得平缓，整体曲线呈对数形式关系。这是由于在低眼内压阶段，角膜基质层中的基质起主要的承担载荷作用，在这一阶段角膜胶原纤维处于松弛状态[12]。在压力增加到一定阶段之后，角膜胶原纤维被拉紧，取代角膜基质承担起抵抗压力的作用，角膜的硬度较前一阶段有明显增加。在2.04～4.02 kPa范围内9组角膜顶点位移的变化均值为0～0.223 mm.当眼内压在2.04～4.02 kPa范围内，角膜弹性模量趋于线性增长趋势，其弹性模量范围为48.39～70.50 kPa，用线性方程拟合后眼内压与弹性模量的关系为y=1.709x+44.55.
3 猪角膜结构的观测及猪眼球整体膨胀实验的模拟计算
3.1 猪角膜结构的观测
对无眼内压状态的角膜进行扫描电子显微镜表征时，对角膜固定、乙醇梯度脱水、置换、干燥后，对其进行电镜表征，如图5所示。从角膜断面的板层结构可以看出角膜在厚度方向上是由多层胶原纤维板层组成；在平行于角膜表面上，胶原纤维板层与板层之间为平行排列。进行透射电子显微镜表征时，首先通过图2(a)实验设备，制备眼内压为0 kPa及4 kPa眼压状态下的眼球，将图2(b)中的输液管针头处密封，保持眼球内眼压状态；其次，将眼球进行固定后，制备角膜试样；最后，对角膜试样进行渗透及包埋后，制备尺寸为1 mm×1 mm×1 mm角膜试样，对其进行透射电子显微镜表征，结果如图6(a)-(d)所示，角膜在不受眼内压作用时胶原纤维整体及单根纤维都呈弯曲状；在眼内压4 kPa时，如图6(e)，角膜中整体板层的胶原纤维被拉直；图6(f)，在角膜胶原纤维板层交界处，胶原纤维被拉断；图6(g)、(h)，板层内部纤维拉断情况不明显。
图2 猪眼球整体膨胀实验装置及输液管针头与眼球连接示意图
Fig.2 (a)Diagram of experimental apparatus of eye globe inflation tests; (b)Contacted infusion needle with eye
图3 猪眼球整体膨胀实验角膜顶点位移-眼内压
Fig.3 Apical rise-intraocular pressure relationship with eye globe inflation tests
图4 猪眼球整体膨胀实验角膜顶点应力-应变及弹性模量眼内压平均值
Fig.4 Stress-strain and mean of elasticity modulus and intraocular press relationship of corneas with eye globe inflation tests
图5 猪角膜断面的板层结构扫描电子显微镜表征图
Fig.5 Structure of section lamellae in porcine cornea
图6 膨胀实验过程角膜透射电子显微镜表征图
Fig.6 Microscopic morphology images of the corneal during expansion experiment
3.2 猪眼球整体膨胀实验模拟计算
3.2.1模型尺寸及网格划分
利用UG三维制图软件，建立角膜计算模型如图7所示，角膜模型的几何参数与计算角膜弹性模量的几何参数相同，角膜的三维模型建好后，导入ABAQUS软件中，对其划分单元属性为3维8节点杂交缩减积分杂交实体单元的(C3D8HR)网格。网格种子尺寸为0.1 mm.
图7 角膜计算模型
Fig.7 Computational model of cornea
3.2.2本构模型及参数
通过对猪角膜的表征观测，可以发现猪角膜基质层中以纤维板层平行铺层为主且铺层具有方向性。因此，在ABAQUS中选用非线性、各向异性的HGO超弹性本构模型，对角膜进行数值模拟计算，计算公式如下：
此公式为采用应变能密度函数来描述材料的本构关系。式中：C10为基质的剪切模量；D为0，代表角膜为不可压缩性材料；k1代表胶原纤维的刚度模量；k2为常数项；κ的值代表纤维沿主方向的分散程度。参照ELSHEIKH et al[13]测量得到k2、D的数值并根据本文中的整体眼球膨胀实验，调整C10及k1的数值，如表2所示，使有限元计算结果与实验结果最为相近。计算结果如图8所示。
图8 模拟数据与实验数据曲线对比图
Fig.8 Comparison of the simulation calculation and the tests
表2 模拟参数
Table 2 Simulation parameters
参数C10/MPaDk1/MPak2κN数值0.002 6800.0454110.2422
3.2.3载荷及边界条件
在ABAQUS有限元软件中选用ABAQUS/Standard求解模块进行计算，给角膜施加眼内压时，将2.04～4.02 kPa的眼内压分成8个压力段，并赋在8个载荷步中进行计算。同时，在角膜缘处约束角膜的位移，即U1=U2=U3=0.没有约束角膜缘处的旋转情况，这一边界条件与整体眼球实验中角膜的边界约束情况一致。
3.3 模拟结果及讨论
图9分别从角膜的内表面、外表面及厚度方向得到了角膜的应力云图。从图9可以看出，在胶原纤维主方向的角膜缘处，角膜所承受的应力最大。随着眼内压的增加，角膜内的应力逐渐变大，但沿角膜中胶原纤维分布主方向发散的应力分布的特点不变；角膜外表面的应力明显小于角膜内表面的应力，且角膜外表面角膜缘处在纤维主方向的应力反而较小；从厚度方向可以得到角膜内表面到外表面应力的变化情况，随着眼内压的升高，角膜内表面应力增加较快，而外表面应力增加较慢，沿厚度方向应力分布梯度越来越大。
图9 不同眼内压下角膜内表面、外表面及厚度方向的Mises应力云图
Fig.9 Von Mises stress maps for different values of the intraocular pressure, outer surface, and the thickness direction of cornea
4 结论
1) 通过完善膨胀实验装置，设计眼球膨胀实验方案，进行了猪整体眼球膨胀实验，得到了角膜的顶点位移-眼内压之间的关系，计算得出，眼内压在2.04～4.02 kPa内，角膜的弹性模量随着眼内压的升高成线性增长，其数值范围为48.39～70.50 kPa.应变随着眼内压升高的变化范围为0.056～0.065.
2) 利用ABAQUS软件对猪眼球整体膨胀实验进行了模拟仿真，有限元计算的角膜顶点位移-眼内压曲线与实验得到的曲线吻合较好，从而通过有限元计算可以得到在不同眼内压作用下角膜的应力分布。
[1] HUSEYNOVA T,WARING G O,ROBERTS C,et al.Corneal biomechanics as a function of intraocular pressure and pachymetry by dynamic infrared signal and scheimpflug imaging analysis in normal eyes[J].American Journal of Ophthalmology,2014,157(4):885-893.
[2] ZENG Y,YANG J,HUANG K,et al.A comparison of biomechanical properties between human and porcine cornea[J].Journal of Biomechanics,2001,34(4):533-537.
[3] 杨坚,曾衍钧,黄昆,等.人与猪角膜的生物力学特性之比较[J].中国生物医学工程学报,2001(2):166-169.
YANG J,ZENG Y J,HUANG K,et al.Biomechanical properties of human and porcine cornea[J].Chinese Journal of Biomedical Engineering,2001(2):166-169.
[4] MALOLEY L A,RAZEGHINEJAD M R,HAVENS S J,et al.Pneumotonometer accuracy using manometric measurements after radial keratotomy,clear corneal incisions and lamellar dissection in porcine eyes[J].Current Eye Research,2019,45(1):1-6.
[5] WANG L,TIAN L,HUANG Y,et al.Assessment of corneal biomechanical properties with inflation test using optical coherence tomography[J].Annals of Biomedical Engineering,2018,46(2):247-256.
[6] ENDERSON K A,LSHEIKH A E,NEWSON T.Application of structural analysis to the mechanical behaviour of the cornear[J].Journal of the Royal Society Interface,2004,1(1):3-15.
[7] BAO F J,DENG M L,WANG Q M,et al.Evaluation of the relationship of corneal biomechanical metrics with physical intraocular pressure and central corneal thickness in exvivo rabbit eye globes[J].Experimental Eye Research,2015,137:11-17.
[8] HOELTZEL D A,ALTMAN P,BUZARD K,et al.Strip extensiometry for comparison of the mechanical response of bovine,rabbit,and human corneas[J].Experimental Eye Research,1992,114(2):202-215.
[9] HOLZAPFEL G A,GASSER T C,OGDEN R W.A new constitutive framework for arterial wall mechanics and a comparative study of material models[J].Journal of Elasticity,2000,61(1):1-48.
[10] 于杰.眼球体积随内压变化的实验研究及有限元数值模拟[D].太原：太原理工大学,2007.
[11] DINCEL V E,SENGELEN,SEPICI V,et al.The association of proximal femur geometry with hip fracture risk[J].Clinical Anatomy,2008,21(6):575-580.
[12] 姜黎,包芳军,张东升,等.膨胀试验测定猪眼角膜生物力学参数的研究[J].医用生物力学,2009,24(2):123-126.
JIANG L,BAO F J,ZHANG D S,et al.Determining the ex vivo biomechanical properties of porcine cornea with inflation test[J].Journal of Medical Biomechanics,2009,24(2):123-126.
[13] ELSHEIKH A,KASSEM W,JONES S W.Strain-rate sensitivity of porcine and ovine corneas[J].Acta of Bioengineering & Biomechanics,2011,13(2):25.