首页| 行业资讯| 企业名录| 周边产品| 数字城市| 增强现实| 工业仿真| 解决方案| 虚拟医疗| 行业仿真| 图形处理| 军事战场
用户登录/注册 ×
资讯首页
行业资讯 >> 专业文献
缝针刺入角膜的力学建模与虚拟仿真
时间:2019-11-21    评论:0
    来源:第三维度
    作者:杨洋、刘笑宇
    单位:北京航空航天大学机械工程及自动化学院

    摘要:利用微力材料试验机以不同的刺入速度 , 对缝针刺入角膜组织产生的微小力进行一系列的力学测量试验.在此基础上, 以超弹性和粘弹性理论为基础构建了角膜组织的材料本构模型.考虑眼外肌和眼压的作用 , 根据自下而上的建模方式建立角膜的几何模型.基于接触问题中的粘—滑理论, 使用点—面和面—面接触的方式建立缝针刺入角膜的接触模型 , 并利用单元生死技术实现手术中缝针刺入角膜过程的有限元仿真.仿真的结果可为角膜移植机器人的精确控制提供重要的依据, 并进一步保证机器人执行角膜移植手术的安全性和可靠性.

    目前已经研制出的医疗机器人可以执行不同的手术任务 , 然而对于角膜移植这样的显微外科手术而言 , 完全或部分依靠机器人操作还是具有相当的难度.对于一个对象为直径不到 10 mm , 厚度大约0.7 mm 的角膜植片, 机器人在手术过程中若有细微的误差就会将其撕破, 甚至会由于用力不当而伤害到病人的眼球, 造成不良的后果[1].为了避免这样的事故发生, 有必要在机器人实施手术之前进行缝针刺入角膜的力学建模并仿真整个过程, 从而可获得眼角膜内部所有位置的力学信息 , 为角膜移植机器人进行手术精准控制提供依据.

    现代外科手术最常见的作业形式就是缝针刺入人体器官软组织.在这个过程中, 深入组织内部的缝针与非均匀的细胞物质相互碰撞、 接触, 并且没有任何的视觉反馈[2].医生仅可以通过缝针造成的触觉感知和组织表面的变形情况来定性地判断手术的缝合质量[3].然而目前的手术机器人还远远不具备这些定性判断的能力.针对这一问题, 近年来许多学者对活体细胞的刺入控制方面进行了深入的研究[4 , 5].其中值得一提的是 , 在虚拟手术训练系统的开发中 , 研究者可以用计算机模拟出缝针刺入组织的过程, 并将计算得到组织表面的变形情况和受力信息反馈给模拟器[6].另外 , Simon 曾把与软组织相似的复合材料 PVC 作为试验试件, 进行了刺入力的研究试验 , 确定了这种材料的力学特性和刺入过程中的力的分布情况[7].

    本文首先针对新鲜的猪眼球进行缝针刺入角膜的力学测量试验 , 得到了刺破点力的大小及针轴力和刺入位移的关系.接着, 根据角膜材料测试试验中的单向拉伸和应力松弛试验结果建立角膜的材料模型, 并通过 APDL 的参数化建模方式建立其几何模型, 同时利用非线性弹簧单元模拟眼外肌的作用 , 并施加相应的眼压.最后, 采用粘—滑接触模型为通过 APDL 语言编写的过程控制流 , 综合采用点 —面接触模型, 面—面接触模型和单元生死技术对整个过程进行有限元仿真.

    1 针刺角膜的微小力测量试验

    1.1 试验材料和方法

    缝针刺入角膜的微小力测量试验系统由试验机 、 新鲜猪眼球样品 、 不同型号的手术缝针和固定台等组成.其中 , 试验机采用中国科学院力学研究所的微力材料试验机 —Instron 5848 M icro Tester ,它的载荷分辨率为 50 mN , 位置精度为 1 nm.手术缝针通过自制的针柄夹持器固定在试验机上, 夹持器带动缝针沿竖直方向移动, 行程为 100 mm , 速率为 0.3 至 25 mm/ s, 试验系统如图 1 所示.

图 1 针刺角膜的试验系统

    相关研究结果表明猪角膜的力学特性与人最为接近[8], 因此采用 4 h 之内宰杀的 20 个猪的新鲜眼球作为试验对象.为了保持生物组织的活性, 将取下的眼球存放在盛有少量庆大霉素的 4 ℃冰盒内以备用.试验前, 把眼球嵌入固定器皿内以保证眼球相对于工作台静止.值得注意的是 , 固定在试验台上的眼球由于失水影响其眼压往往小于正常范围(11 —21 mmHg(1 mmHg =133.322 Pa))[9].为了确保试验的一致性 , 需要向眼球的玻璃体腔内注射生理盐水, 使眼压达到平均值 17 mmHg.试验选取的医用缝针为眼科 1 号针, 2 号针和 5 号针 , 它的平均直径分别为 0.1 , 0.2 和 0.5 mm.临床手术研究表明 , 手术缝针刺入器官组织的速度一般为 0.4 —15mm/ s.这里我们分别采用 0.5 , 1 , 10 和15 mm/ s的刺入速度进行了 4 组试验.

    1.2 试验结果及分析

    这里就选取眼科 2 号缝针 , 使其以 1 mm/ s 的速度刺入角膜.通过试验获得针轴力 —位移关系曲线.把相同条件下 10 个眼球试样的统计关系用Matlab 描绘成图 2 的关系曲线.

    可以看出, 缝针刺入眼角膜的过程可分为 3 个阶段:第一个阶段为缝针从 0 点开始接触角膜表面直至 A 刺破点.此阶段中 , 由于未刺破角膜表面 ,针轴力完全表现为表面接触刚度力.这里 0 —A 按照线性变化处理, 其斜率为 2.64 N/mm ;A 点处的刺破力为 0.366 N , 对应的位移为 0.142 mm.刺入角膜组织后, 针尖受到的力急剧下降至 B(0.146 ,0.239), 此时角膜的形变表面也开始回复.B —C 阶段由于缝针刚刺入角膜, 摩擦阻力并不明显, 而缝针两刃产生的切割力使得该段受力表现出明显的不稳定波动.从 C 点开始, 针轴力主要表现为摩擦阻力并呈现出稳定线性上升的趋势 , 其斜率为 0.321N/mm.

图 2 针轴力—位移关系曲线

    2 缝针刺入角膜的力学建模

    2.1 角膜材料的生物力学建模

    角膜材料模型是依据角膜试片的单向拉伸和应力松弛材料特性测量试验为基础 , 综合考虑采用超弹性和粘弹性模型建立角膜材料的本构模型 , 通过曲线拟合的手段得到相应的数学描述.根据拟合的模型参数, 在有限元软件 ANSYS 中进行角膜特性试验的仿真 , 模拟出在拉力的作用下角膜的变形和应力的影响 , 以及该影响随时间的变化过程.

    2.1.1 超弹性本构方程的建立

    生物组织具有的主要特征是大变形和不可压缩性 , 并且在一定的时间内发生连续的形变回复[10].因此, 在建立材料模型的时候不能简单的利用虎克定律去处 理, 而Mo oney 和 Riv lin 等[11] 建立的超弹性理论可以很好的描述生物材料这种大变形和不可压缩的特点.这里采用以试验为基础的唯象理论建立角膜材料的超弹性 本构 模 型, 描 述角 膜 在 外 力 作用 下 的 行为[12 , 13].唯象理论是 M ooney(1940 年)提出一种描述超弹性行为的本构理论[14], 它是从试验出发建立描述橡胶力学行为的数学表达式, 其应变储能函数具有以下形式:


    其中 Ci, j 称为 M oo ney 常数 , 反应的是材料自身的特性 , 弹性储能公式为 I1 , I2 , I3 三个应力不变量的函数, 对于角膜的单向拉伸试验而言 , λ1 =λ,λ2 =λ3 =λ-1/2, 此时 :


    应变储能函数W 的系数可根据试验的图像性质判断 , 这里采用 3 个系数的储能函数来拟合试验数据, 得到 :


    根据最小二乘法利用公式(3)对其进行拟合处理, 得到 C10 =-3.622 , C01 =3.963 , C11 =5.614 ,在 MA T LAB 下得到角膜单向拉伸状态下的应力 —伸长率拟合曲线见图 3 所示.

图3 拟合超弹性模型

    在 ANSYS 环境中,选择角膜的单元类型为 SO LID183 , 建立与角膜单向拉伸试验相同的边界条件, 利用 M A TLAB 曲线拟合得到的超弹性材料系数模拟角膜试样的拉伸过程, 得到如图 4 所示的角膜试片的应力—应变曲线.

图4 ANSYS 模拟超弹性

    2.1.2 粘弹性本构方程的建立

    对物质的力学形态而言, 各种生物组织均由固体和液体两部分构成.在力学或生物流变学中, 能较好的反应生物组织内部固体和液体相互作用的理论之一是粘弹性理论[15].利用M axw ell 单元和弹簧单元的并联构建角膜粘弹性本构模型[16], 基于这种形式的粘弹性模型具有十分简单的数学表达式, 非常利于数值模拟计算.为了反映眼角膜应力松弛的应力—时间曲线关系 , 眼角膜粘弹性本构模型可以用 4 个 M axw ell 单元和一个弹簧单元表达为如下形式:


    其中 Ei 代表弹簧的弹性系数 , ηi 代表阻尼器的粘性系数.在 M A TLAB 软件中, 将松弛试验的数据点用公式(4)进行曲线拟合如图 5 , 求得弹性系数 Ei和阻尼器的粘性系数 ηi 分别为 :E1 =0.1204 , E2 =0.1743 , E3 =0.2026 , E4 =0.2023 , E5 =0.2996 ;η1 = 0.0385 , η2 = 0.6887 , η3 = 7.8043 , η4 =75.824.

图5 拟合粘弹性模型

    在 ANSYS 中 , 设置与应力松弛试验同样的边界条件 , 模拟应力松弛试验的整个过程, 得到角膜试样的应力随时间的变化过程, 如图 6 所示.

图6 ANSYS 模拟应力松弛

    利用 ANSYS 软件的二次开发功能 , 综合超弹性和粘弹性模型来模拟角膜的应力松弛试验, 得到图 7 所示的应力—时间关系图, 它反应的是角膜材料的非线性粘弹性行为[17].

图 7 ANSYS模拟角膜非线性特征

    由于拉伸的过程时间短暂, 而松弛又是一个漫长的过程 , 所以时间轴采用指数形式.从图 7 可以看到, 在最高点之前, 角膜试样被拉伸, 曲线表现的是角膜的超弹性行为的非线性特征;在曲线的最高点之后 , 角膜试样停止拉伸 , 开始转为松弛阶段.图像的特征表明其应力随时间几乎呈线性衰减.

    2.2 角膜的几何建模

    在实际的临床应用中, 角膜的几何信息是通过OCT 断层图像扫描仪采集到的, 这些信息包括 :角膜表面的曲率半径 , 角膜的厚度 、 宽度等 , 其采集界面如图 8 所示.

图 8 OCT采集界面

    在分析接触问题时, 网格的稀疏程度直接影响到接触作用的精度, 缝针与角膜直接接触的部分需要很高精度的网格划分.然而 , 过密的网格应用于整个模型必然会增加运算时间的消耗.最佳的办法就是将缝针与角膜接触部分的局部网格细化 , 而其他部分的网格则只需保留一定的精度.另一方面 ,由于模型将采用单元生死的方法来控制缝针的行进位移和角膜的有限元分析载荷, 因此在建立有限元模型后需要准确的知道节点的编号和位置.

    综合上述的考虑 , 角膜有限元模型的建立采用自下而上的建模方式, 即根据角膜的几何形状利用APDL 参数化语言建立一系列的节点, 然后将这些节点组合成单元模型.这些节点和单元都是以直接的方式生成的, 避免了自上而下建模方式中节点和单元无法控制的问题.在缝针刺入的位置进行网格细化, 而远离接触部位的网格逐渐变稀如图 9 所示.

图9 自上而下方式建立角膜到节点

    整个角膜的单元集共分为 3 个部分 :细化单元, 中间单元和粗化单元.每个单元集都要用到APDL 编写的过渡网格连接进行连接.过渡网格的连接采用三角单元形式 , 而 SOLID183 派生的三角形单元模型恰好可以满足要求 , 如图 10 所示.

图10 过渡网格单元

    2.3 边界条件

    2.3.1 眼外肌的作用 由于眼球并不是完全的固

    定在眼眶中, 而是通过眼外肌相连接.眼球靠眼外肌的收缩和松弛而产生协调的运动, 其中包括一定范围内的旋转和平移运动.在进行角膜移植等眼科手术时, 为了便于医生操作的准确性 , 需要注射药物以限制眼球的运动.然而 , 经过测试发现, 眼球在径向上仍存在着移动, 通过挤压眼球的试验 , 得到眼球在压力的作用下的纵向位移与作用力的关系[18].为了在建模过程中反映这种影响, 采用非线性弹簧单元 COM BINE39 来模拟眼外肌的非线性特征.COM BIN39 是一个具有非线性功能的弹簧单元, 可对此单元输入广义的力—变形曲线的关键点来定义非线性的程度.定义曲线上的各点(D1 ,F1)(D2 , F2)…(D15 , F15)来代表力 —平动位移关系曲线见图 11.

图 11 COMBIN39 单元的输入数据

    2.3.2 眼内压的施加 眼内压简称 “眼压” , 是指眼球内容物对眼球壁产生的压力 , 它的主要作用是维持眼球形态 、 保持正常生理功能.正常情况下 ,成年人 处于 正常 生理状 态的 眼压保 持在 11 —21 mmHg范围内.眼内压的存在导致角膜组织内部应力的产生.在实际的临床分析中, 需要将测量得到的眼压以压强的方式施加于角膜的内表面 , 施加后的眼压和眼外肌如图 12 所示.

图 12 眼压和眼外肌的施加

    3 缝针刺入角膜的有限元仿真

    3.1 接触分析

    接触问题是一种高度非线性行为, 它有两种基本类型 :刚 ─柔接触和柔 ─柔接触.在缝针刺入角膜的接触模型中 , 由于缝针的刚度比角膜的刚度大的多(钢的弹性模量约为 200 GPa , 而角膜的拟弹性模量只有 17.6 M Pa), 所以选择刚—柔接触类型.

    有限元的计算中存在 3 种接触方式:点 ─点接触, 点─面接触和面─面接触, 不同接触方式适用于不同的分析类型.其中 , 缝针的针尖和角膜的接触采用点 —面接触, 而缝针的针杆和角膜的接触定义为面—面接触.点 —面接触采用二维 Co nta175单元模拟刚性点和柔性面的接触.在面—面接触中, 缝针针杆的外表面被当做刚性的 “目标面” , 采用二维单元 Target169 模拟 ;角膜的外表面以及杀死单元后的内表面看成是柔性体的 “接触面” , 采用Conta172 单元来模拟, 具体的描述如图 13 所示.

图 13 接触单元的类型

    3.2 角膜单元的生死过程及算法设计

    所谓 “单元生死” 是指在模型中加入或删除材料 , 使模型中相应的单元 “生” 或 “ 死” , 达到“单元死” 的效果[19].这个过程并非将 “ 杀死” 的单元从模型中删除, 而是将其刚度矩阵乘以一个很小的因子.在一些情况下 , 单元的 “生死” 状态可以根据 ANSYS 的计算数值决定 , 如应力 , 应变等.

    在模拟缝针刺入角膜的过程中, 通过逐步杀死接触的角膜单元 , 使缝针不断的前行, 直至将角膜完全穿破.缝针和角膜组织之间的作用按照粘 —滑接触模型建立.这种模型认为相对接触滑动的两个物体表面之间存在两种状态 :粘着和滑动[20].如果相互运动产生的阻力小于或等于临界力, 则呈现粘着状态 , 接触表面彼此间不存在相对运动 ;反之则处于滑动状态[21].控制单元生死的判断准则是以缝针刺入角膜的微小力测试试验为依据 , 根据力 —位移曲线判断是否杀死接触的角膜单元而使得缝针继续前行 , 其求解过程的描述如图 14 所示 :

图 14 单元生死控制的描述
 (a)缝针未接触角膜;(b)缝针接触第一个节点;(c)缝针与第一个节点粘着; 
(d)杀死接触单元;(e)缝针与第二个节点粘着;(f)杀死单元, 移至下一点 

    缝针作为刚体的目标面 , 当给定一个缝针行进的步长, 其针尖接触到角膜的第一个节点时发生点—面接触 , 如图 14(b).随着缝针的继续前行,接触的角膜表面发生形变, 如图 14(c).由于接触的第一个节点是刺破点, 它的临界力为试验测得的刺破力 0.366 N.在每个载荷步中 , 缝针每前行一个位移便计算它所受到的反作用力.当计算值大于临界力时, 接触状态由粘着变为滑动, 将接触的失效单元杀死.同时缝针针尖脱离该节点, 指向沿刺入方向下一个单元上的节点上, 如图 14(d).此时 ,经过程序的渗透检验之后建立它们的面—面接触模型.并且通过给定的临界力判断接触状态是否改变 , 如图 14(e).由于此时缝针已经刺入角膜 , 临界力的判断准则为图 2 C 点之后的斜率0.321N/mm.

    当计算值大于临界值时, 接触状态由粘着变成滑动,此时将会再次杀死接触的角膜单元, 如图14(f).

    3.3 缝针刺入角膜的仿真结果分析

    在 ANSYS 环境下建立缝针刺入角膜的有限元模型并对其进行仿真 , 具体的分析结果如下:当缝针接触角膜之后 , 针尖与角膜表面单元开始发生接触作用, 此时角膜开始发生变形, 并在针尖的接触部位产生了应力集中.随着缝针位移的不断行进,集中部位的应力值不断增加 , 当缝针从接触角膜表面开始行进至 1.46 mm 时角膜被刺破, 此时针尖处的接触应力为 7.536 M Pa , 而角膜内表面的最大应力为 2.336 M Pa , 见图 15.

图 15 缝针刺破角膜时刻的角膜应力图

    缝针在刺入角膜的过程中, 角膜内表面产生的集中应力主要分布在缝针行进位移的延长线附近(如图 16 所示).从医学角度来讲, 角膜内表面造成的应力损伤不仅很难愈合 , 而且就目前而言不能通过有效的手段来加以控制.我们进行仿真的目的就是要根据模拟的结果判断内表面的应力是否达到或超过角膜组织所承受的应力极限, 从而保证术后的角膜不会造成组织损伤.根据生物力学理论 , 通过生物组织的应力 —应变特性曲线可以判断生物组织的损伤和失效情况.

图 16 角膜内表面的集中应力

    根据 角 膜 拉 伸 试 验[18], 得 到 应 力 极 限 为3.0 M Pa , 即约定如果角膜内表面的集中应力超过该值 , 仿真程序就停止运算.由图 16 可以看出,角膜被刺穿的瞬间 , 内表面的最大应力并没有超过组织损伤的应力极限.

    当接触的角膜达到判别标准之后便将被刺破的单元杀死.可以看到图 17 中, 角膜在刺破之后有一个明显的回弹趋势 , 接触应力降低至 4.83 M Pa ,角膜内表面应力也随之降低, 角膜的应力分布比起刺破状态有了相当的改善.

图 17 刺破单元之后角膜的应力情况

    随着缝针的继续行进, 由于角膜材料模型的松弛效应 , 尽管作用在针轴上的力逐渐增加 , 但是角膜的变 形 却在 不 断 的恢 复 , 当缝 针 的 位 移为2.65 mm时, 角膜的内表面的集中应力达到了判断极限 3.0 M Pa , 见图 18.如果缝针再继续前行 , 角膜内部将出现损伤.经过计算 , 在该状态下缝针刺入角膜的深度为 0.843 mm , 约为角膜总厚度的85 %.临床手术中医生使用缝针缝合角膜的一般深度不允许超过角膜厚度的 3/4 , 因此, 基本上与利用有限元方法模拟缝针刺入角膜的仿真结果一致.

图 18 角膜内表面的集中应力超出极限值

    4 结论

    (1)利用微力材料试验机进行了缝针刺入角膜的微小力测量试验 , 得到了正常眼压下医用眼科缝针以 1 mm/ s 的速度刺入新鲜猪眼球的位移—针轴力关系曲线, 并归纳出刺入角膜的 3 个力学阶段.

    (2)以超弹性和粘弹性模型为基础, 构建了角膜的材料模型.并根据得到的数学描述 , 在 ANSYS 中进行角膜特性的模拟试验, 描述出在外力作用下角膜的应力和应变的响应曲线.

    (3)采用的是自下而上的方式 , 通过 APDL 语言建立角膜的有限元几何模型, 并利用非线性弹簧模拟眼外肌的影响, 同时考虑眼压的作用.

    (4)采用粘 —滑接触模型为基础, 利用点—面接触模型 , 面—面接触模型和单元生死技术对整个过程进行有限元仿真.仿真结果表明 , 角膜的大变形和松弛特性得到了很好的体现, 仿真结果也与实际的临床手术保持一致.

    参考文献(略)

标签:眼科医疗手术有限元
上一篇:基于Unity3D的气相色谱仪虚拟仿真实验系统的构建下一篇:半实物仿真与航天工程
网友评论:缝针刺入角膜的力学建模与虚拟仿真
留名: 验证码:
最新评论
查看全部评论0
暂无评论
您可能还需要关注一下内容:
·裸眼3D胸腔镜与2D胸腔镜在微创食管癌根治术中的临床对比研究
·人民日报:虚拟现实 由虚向实
·水轮机叶片热模压成型过程仿真及残余应力分析
·我军海上医疗救护模拟仿真训练系统问世
·利用虚拟仿真技术改进医疗设备设计
·虚拟仿真提升植入式医疗设备的性能和安全性
·加拿大CAE在广东成立飞行模拟器合资公司
·虚拟手术中血管与布料切割训练模块的设计与实现
·使用Unity3D和HTC VIVE实现下颌骨虚拟手术
·汽车正面耐碰撞性有限元仿真分析
☏ 推荐产品

小宅 Z5 2018青春版
商家:小宅

杰瑞特运动平台
商家:杰瑞特智能

Dikalis眼动追踪
商家:赢富仪器

魔神 Hawk
商家:魔神运动分析

5DT Binoculars
商家:四维宇宙

Christie DS+750
商家:四维宇宙

Zalman M220W
商家:四维宇宙

全息360°
商家:四维宇宙

PD F10 AS3D
商家:四维宇宙

PHANTOM系列
商家:四维宇宙
☞ 外设导航
☏ 企业名录
【宁波】宁波维真显示科技股份有限公司
【潍坊】歌尔股份有限公司
【上海】霍尼韦尔(中国)有限公司
【北京】科视Christie-中国
【北京】北京华如科技股份有限公司
【北京】北京乐卡仕技术有限公司
【广州】广州弥德科技有限公司
【上海】刃之砺信息科技(上海)有限公司
【北京】北京度量科技有限公司
【北京】北京小鸟看看科技有限公司
关于本站联系我们融资计划免责声明网站建设广告服务咨询策划行业推广
北京第三维度科技有限公司 版权所有 京ICP备15051154号-3
2008-2020 Beijing The third dimension Inc. All Rights Reserved.
Email:d3dweb@163.com  QQ:496466882
Mob:13371637112(24小时)
关注虚拟现实
关注第三维度