一种基于生物力学的嵌甲预测方法与流程

本发明属于嵌甲相关技术领域，更具体地，涉及一种生物力学的嵌甲预测方法。

背景技术：

通俗地讲，嵌甲就是趾甲长到肉里去，造成嵌甲的原因之一是由于不良的修剪趾甲习惯所造成，就是通常人们为了适应趾腹的形状，而将趾甲修剪成圆弧形，将趾甲两侧边缘修剪的太深太低。这时如果再穿一双尖头的鞋子，就会将甲沟组织挤向趾甲，最终导致嵌甲。嵌甲给人们带来了极大不便和痛苦。

嵌甲是一种严重影响患者生活质量的足部疾病，但由于其形成机制还明确，导致嵌甲的治疗方法还存在争议，治疗效果差强人意。趾甲的现有模型主要是板壳模型与小应变模型，没有涉及软组织的生长机制，不能够反映甲板与甲皱的耦合作用。

技术实现要素：

针对现有技术的以上缺陷或改进需求，本发明提供了一种基于生物力学的嵌甲预测方法，其建立了体现甲板与甲皱相互耦合作用的趾甲有限元模型，基于有限生长运动学和开放系统热力学理论来考虑甲皱的生长，并基于骨骼适应性生长理论研究甲板重建，通过abaqus二次开发对趾甲有限元模型进行数值实现，其中通过接触和约束的设定实现甲板与甲皱之间的耦合作用，进而基于所述趾甲有限元模型对嵌甲的形成及生长进行模拟及分析，以实现对嵌甲的预测，且为嵌甲形成和治疗提供理论依据。

为实现上述目的，按照本发明的一个方面，提供了一种基于生物力学的嵌甲预测方法，所述嵌甲预测方法包括以下步骤：

(1)基于有限生长运动学方程及开放系统热力学方程确定甲皱的应变能函数及生长驱动方程，进而构建甲皱生长模型；同时，基于骨骼自适应生长理论构建甲板重建模型，所述甲板重建模型考虑了湿度及应力影响；所述甲皱生长模型及所述甲板重建模型构成趾甲的材料模型；其中，所述甲板重建模型将应变能密度作为力学激励，且其引入了死区效应；

(2)基于趾甲的样本数据确定趾甲的材料及生长参数，并构建甲板的三维模型及甲皱的三维模型，通过接触和约束设置将甲板的三维模型及甲皱的三维模型耦合以得到趾甲三维模型，继而将所述材料模型加入所述趾甲三维模型得到趾甲有限元模型；

(3)将待预测趾甲的参数输入到所述趾甲有限元模型中，所述趾甲有限元模型基于接收到的参数数据对趾甲的生长过程进行模拟，进而实现对嵌甲的预测。

进一步地，所述应变能函数的公式为：

式中，i1^e和i3^e分别是弹性部分右柯西格林应变张量的第一不变量和第三不变量；λ为拉梅第一常数；μ为拉梅第二常数。

进一步地，所述甲板重建模型采用的内部重建控制方程为：

式中，β是重建系数；k是一个稳态值；s是死区边界阈值；u为应变能密度；ρ为甲板密度。

进一步地，应变能密度的公式为：

式中，εij是应变；σij是柯西应力。

进一步地，甲板局部的弹性模量与表观密度之间的关系为：

ei＝ξρiⁿ

式中，ξ和n是常数；ρi为甲板局部密度。

进一步地，甲板重建模型的本构方程为：

σij＝eijklεij

式中，eijkl为刚度矩阵；εij为应变。

进一步地，通过将内部重建控制方程代入本构方程，继而得到甲板重建模型。

进一步地，趾甲的弹性模量的计算公式为：

式中，ρ为甲板密度；ρ0为初始密度；ψ为湿度。

进一步地，步骤(2)中趾甲的材料及生长参数包括甲板的弹性模量、拉梅第一常数、拉梅第二常数、最大肌节长度比值、适应速度常数、组织沉积非线性常数、临界应力、初始密度、最大密度、弹性模量与密度的相关系数、重建系数、死区边界阈值及稳态值。

总体而言，通过本发明所构思的以上技术方案与现有技术相比，本发明提供的基于生物力学的嵌甲预测方法主要具有以下有益效果：

1.基于有限生长运动学方程及开放系统热力学方程确定甲皱的应变能函数及生长驱动方程，进而构建甲皱生长模型；同时，基于骨骼自适应生长理论构建甲板重建模型，所述甲板重建模型考虑了湿度及应力影响；所述甲皱生长模型及所述甲板重建模型构成趾甲的材料模型，所建立的模型能够较好地反应甲板与甲皱的相互耦合作用，继而提高了准确性。

2.将所述材料模型加入所述趾甲三维模型得到趾甲有限元模型，将待预测趾甲的参数输入到所述趾甲有限元模型中，所述趾甲有限元模型基于接收到的参数数据对趾甲的生长过程进行模拟，进而实现对嵌甲的预测，由此能够较好地研究嵌甲的形成机制，且可准确地对嵌甲进行预测，为趾甲的生长和演化特性研究提供了数据支持。

3.通过建立的甲皱生长模型推导了考虑生长的甲皱本构方程和焦曼率形式的空间切线刚度矩阵，采用umat对甲板和甲皱生长重建模型进行了算法实现，将趾甲模型在有限元软件abaqus中模拟保证了理论模型的正确性，适用性较好。

4.通过纳米压痕实验测出甲板的弹性模量，结合实验和临床数据对模型参数进行确定，选择了合适的参数，进而保证了模型预测的准确性。

附图说明

图1是本发明提供的基于生物力学的嵌甲预测方法的流程示意图；

图2是图1中的基于生物力学的嵌甲预测方法涉及的甲板重建的流程示意图；

图3是图1中的基于生物力学的嵌甲预测方法涉及的甲皱有限生长的流程示意图；

图4是本发明通过纳米压痕仪获得的甲板的弹性模量示意图；

图5是本发明涉及的趾甲弹性模量随湿度变化的拟合曲线。

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。此外，下面所描述的本发明各个实施方式中所涉及到的技术特征只要彼此之间未构成冲突就可以相互组合。

请参阅图1、图2、图3、图4及图5，本发明提供的基于生物力学的嵌甲预测方法，所述嵌甲预测方法主要包括以下步骤：

步骤一，基于有限生长运动学方程及开放系统热力学方程确定甲皱的应变能函数及生长驱动方程，进而构建甲皱生长模型；同时，基于骨骼自适应生长理论构建甲板重建模型，所述甲板重建模型考虑了湿度及应力影响；所述甲皱生长模型及所述甲板重建模型构成趾甲的材料模型；其中，所述甲板重建模型将应变能密度作为力学激励，且其引入了死区效应。

具体地，甲皱生长模型的构建过程为：

(1)基于软组织生长理论建立甲皱的有限生长运动学方程。

具体地，在有限生长运动学方程中，表征变形的关键运动参数是变形梯度张量，其为：其中，χ表示材料构型坐标与相应空间构型坐标之间的映射，表示χ相对于未变形坐标x的变形梯度。变形梯度张量f分解为弹性部分fe和生长部分fg，f＝fe·fg。

假设甲皱在每个方向上的生长特征相同，引入一个标量生长参数θ表征一系列的肌节沉积产生的纵向生长，生长变形梯度可以表示为fg＝θi，右柯西格林应变张量c和左柯西格林应变张量b可表示成：c＝f^t·f，b＝f·f^t。

本实施方式中，应变能函数的表达式为：

式中，i1^e和i3^e分别是弹性部分右柯西格林应变张量的第一不变量和第三不变量；λ为拉梅第一常数；μ为拉梅第二常数；用一对共轭的第二类p-k应力张量s和右柯西格林应变表征应力应变关系：

第二类p-k应力张量弹性部分为：

kirchhoff应力为：

τ＝fsf^t＝fe·se·fe^t

由kirchhoff应力可得柯西(cauchy)应力：

得到本构方程，kirchhoff应力表达式为：

(2)将所述有限生长运动学方程嵌入到开发系统热力学方程中，进而得到甲皱生长模型。

具体地，本步骤中主要确定质量平衡方程及生长控制方程。其中，设定中间构型的密度为ρg，对应的体积单元为dvg，空间构型的密度为ρt，对应的体积单元为dv，质量单元初始的密度为ρ0^*，生长后的质量单元可表示为：

由于从材料构型到中间构型和从中间构型到空间构型都是纯弹性映射dm＝ρgdvg，dm＝ρtdv，其中dv＝jdv，dvg＝jgdv，dv＝jedvg可得到：

ρ0＝jρt＝jgρg

则有材料构型的局部质量平衡方程为：

空间速度梯度l：

弹性速度梯度le和生长速度梯度lg分别为：

生长速度梯度lg可表示成：

没有一般性的损失，质量流r被认为是微不足道可以忽略的，则有：

推导出中间构型的质量平衡方程为：

质量变化的生长有三种类型，分别为密度不变的生长、体积保持不变的生长及密度和体积同时改变的生长。对于体积不变的生长jg＝1，对于密度不变的生长ρ0^*＝ρg为常数，质量源表示为：

在甲皱生长模型中，将甲皱组织的生长简化成纯体积生长，即密度不变的生长，当生长变形梯度fg和它的导数已知，质量源可以直接确定：

mandel应力m弹性部分：me＝ce·se。

由于me和生长速度梯度lg是共轭的，trme与kirchhoff应力的迹trτ等价，在密度不变的生长中质量源引起的伸长率进化的方程可以表示成如下形式：

基于上式引入如下应力驱动甲皱生长的进化方程：其中，pcrit是临界载荷，k(θ)是θ的尺度函数，η是生长适应速度相关的常数，θ^max是甲皱组织在各方向的肌节长度相对未生长时的最大比值，γ是甲皱组织沉积的非线性常数。

甲板重建模型的构建过程中，在甲板的力激励生长模型中引入类似骨骼生长的适应性生长理论，将应变能密度作为力学激励，应变能密度为：

式中，εij是应变；σij是柯西应力。同时，还在甲板重建模型中引入死区效应，内部重建控制方程为：

式中，β是重建系数，k是一个稳态值，s是死区边界阈值。内部重建与甲板的表观密度相关，甲板局部的弹性模量和表观密度之间的关系可由如下方程表示：

ei＝ξρiⁿ

其中，ξ和n是常数，可将甲板当做各向同性材料进行处理，本构方程如下：

σij＝eijklεij

通过将内部重建控制方程代入本构方程，继而得到甲板重建模型。

本实施方式中，考虑了湿度影响的甲板区域，受湿度影响的甲板弹性模量由应变能密度和湿度共同决定，将手指甲的弹性模量与湿度的关系式按照比例换算成脚趾甲的弹性模量与湿度的关系，趾甲弹性模量为：

其中，等效密度ρ^*为：

将考虑湿度影响的甲板区域单独定义成一个材料属性，在abaqus/cae中进行设置。先由未考虑湿度影响的模型计算出当前载荷下甲皱q点每个增量步的生长参数θ，在usdfld中将每步计算得到的状态变量按照顺序定义成数组a，将a(kinc)赋值给状态变量供湿度影响的甲板umat子程序调用。

步骤二，基于趾甲的样本数据确定趾甲的材料及生长参数，并构建甲板的三维模型及甲皱的三维模型，通过接触和约束设置将甲板的三维模型及甲皱的三维模型耦合以得到趾甲三维模型，继而将所述材料模型及所述趾甲三维模型进行耦合以得到趾甲有限元模型。

具体地，结合临床数据及甲板纳米压痕实验确定趾甲的材料及生长参数：甲板的弹性模量为7.9gpa，拉梅第一常数λ＝2.857mpa，第二常数μ＝0.714mpa，最大肌节长度比值设为θmax＝1.3，适应速度常数η＝1，组织沉积非线性常数γ＝2，临界应力定为pcrit＝0.3，初始密度ρ0＝1g/cm³，最大密度ρmax＝1.74g/cm³，弹性模量与密度的相关系数ξ＝7900，n＝2，重建系数β＝1，死区边界阈值s＝0.1，稳态值k＝0.004j/g。

本实施方式中，选取四个正常志愿者，每个志愿者选取大拇脚指远端3mm处的趾甲，制样前将趾甲放在零下20度条件下保存，实验前采用自凝塑料将趾甲样品镶嵌于内径20毫米的金相切片软胶模具中，将四个试样2个水平镶嵌，2个竖直镶嵌，依次采用p600、800、1000、1200的不同粗糙度砂纸间歇打磨至暴露趾甲，然后用金刚石抛光膏进行抛光，试样制备后放入封闭的盒子保存，采用纳米压痕仪对试样进行压痕实验，每个试样取5个测点，每次测量设定载荷为20mn，压痕间距为50μm，完整的压痕过程包括加载和卸载两个阶段，卸载过程中材料都有一定的弹性恢复，通过对卸载曲线分析，结合经验公式可以得到材料的弹性模量，测得趾甲甲板的弹性模量为7.9gpa。

常见软组织如血管、心脏弹性模量在3mpa左右，足部肌肉等弹性模量在1-10mpa之间，通过对多种软组织弹性模量的对比，将甲皱的弹性模量和泊松比分别取为2mpa和0.4，计算出拉梅第一常数λ＝2.857mpa，第二常数μ＝0.714mpa，通过临床和其他软组织的研究，将最大肌节长度比值设为θmax＝1.3，适应速度常数η＝1，组织沉积非线性常数γ＝2。脚趾的受力一般在0.1mpa附近，将临界应力定为pcrit＝0.3。

甲板和骨骼在结构、力学性能有很多类似的特性，甲板材料的泊松比取和骨骼相同的0.37，甲板的弹性模量取实验测得的数据7.9gpa，初始密度ρ0＝1g/cm³，最大密度ρmax＝1.74g/cm³，弹性模量与密度的相关系数ξ＝7900，n＝2，重建系数β＝1，由于甲板和骨骼的生长方式不同，在重建参数上也存在许多差异，死区效应参数通过参考其他硬组织、参数识别进行确定，死区边界阈值s＝0.1，稳态值k＝0.004j/g。

本实施方式中，运用abaqus二次开发，通过umat子程序对所述材料模型进行实现，其中，甲板重建的流程为：

(1)增量步开始；

(2)调用umat；

(3)读取初始应力应变；

(4)计算应变是否在死区，若是，则执行步骤(9)；否则，则运用向前欧拉法，选取甲板密度重建迭代参数，对甲板密度进行更新，具体的更新迭代公式为：

式中，取δt＝0.02，设置最大迭代次数为30，当时停止迭代。

(5)由步骤(4)计算的结构更新密度和弹性模量。

(6)n-r迭代收敛，若是，则执行步骤(8)，若否，则执行步骤(9)。

(7)计算修正量，并更新应力应变。

(8)计算雅克比矩阵，存储状态变量矩阵。

(9)umat退出，增量步结束。

甲皱的有限生长过程为：

(1)读取初始生长参数θn。

(2)由变形梯度张量f计算fe、be。

(3)计算出kirchhoff应力τ。

(4)检验应力：θ＝θn，如果|trτ|＞0.3，则计算出k(θ)，否则更新状态变量，返回abaqus主程序。

(5)牛顿迭代：计算残差rθ，用向后欧拉法选取甲皱生长迭代参数获取更新后的伸长率，残差计算公式如下：

r(θⁱ)＝-θⁱ+θn+k(θⁱ)max{|trτ|-pcrit，0}.δt

r(θⁱ⁺¹)＝r(θⁱ+δθ)＝r(θⁱ)+r′(θⁱ)δθ＝0

迭代公式如下：

θⁱ⁺¹＝θⁱ+δθ

其中，δt取0.01，设置最大迭代次数为30。

(6)当rθ小于设定的容差时，退出迭代，更新θ及柯西应力。

(7)计算更新雅克比矩阵更新状态变量，返回abaqus主程序。

本实施方式中，甲皱的长度为10mm，最大宽度为8mm，甲沟宽度为0.2mm，甲板和甲皱都采用solid固体单元，通过拉伸和剪切等操作生成几何模型，采用用户自定义材料模型设置材料参数，甲板弹性模量为7.9gpa，甲板泊松比0.37，甲皱拉梅常数λ为2.857mpa，μ为0.714mpa；不考虑甲板的向前生长，采用abaqus简化模型，在constraint菜单中选择tie连接；设置甲沟处甲板与甲皱的相互作用的接触对，在定义接触属性时将切向属性摩擦系数设置为0.3，法向属性选择为硬接触；建立合适的参考点和参考面对甲板甲皱进行剖分，设置甲板甲皱单元类型，甲皱可选择杂交单元，将甲板的缩减积分单元取消，选择六面体网格进行网格划分，控制单元的大小为0.4，将甲皱划分为,16588个单元，甲板划分为1305个单元；不考虑甲下趾骨增生以及底部作用力对嵌甲形成的影响，主要研究侧面的挤压对甲板甲皱生长及应力分布的影响，选择初始分析步将甲下趾骨用趾甲底部的固定端约束代替，甲板和甲皱近心端也采用固定端约束，在第一个分析步将载荷加在侧面的甲皱上。

步骤三，将待预测趾甲的参数输入到所述趾甲有限元模型中，所述趾甲有限元模型基于接收到的参数数据对趾甲的生长过程进行模拟，进而实现对嵌甲的预测。

本实施方式中，基于所述趾甲有限元模型对外载荷、临界应力、湿度等影响趾甲生长的过程进行分析，以对嵌甲的形成机制进行分析，具体地：基于有限生长理论和适应性生长理论开发的甲皱甲板材料子程序以及abaqus建立的几何模型，分别计算载荷为0.4mpa、0.3mpa以及0.2mpa下甲板和甲皱的应力分布，探索甲板甲皱重建的演化过程。

先天畸形、过度修剪或处在特殊外环境下的趾甲容易发生甲皱损伤，炎症，这些因素会使甲皱生长的临界应力发生改变，为了研究临界应力对嵌甲形成的影响规律，分别对载荷为0.3mpa，临界应力为0.25、0.26、0.27、0.28、0.29、0.3mpa的趾甲模型进行计算。

甲沟处湿度发生变化是嵌甲患者的一个明显特征，为了探索湿度对趾甲的影响规律，对载荷为0.3mpa时，临界应力为0.25mpa的湿度影响模型进行计算。湿度影响区甲板的等效密度用新的状态变量进行显示。

首先，为了探索趾甲的影响因素，先在甲板和甲皱上分别划分一个病变区域，根据前面的计算结果将该区域甲皱的生长参数设置为1.1，甲板密度设置为0.5建立病变模型；同时建立一个对病变模型进行甲板部分切除的模型，研究两者在0.1mpa小载荷时的应力分布差异。

随后对切口的影响进行研究，建立了切除甲板和切除甲皱恢复后的模型(甲皱生长参数和甲板密度都恢复正常但切口还在)，分析两者在0.1mpa小载荷时与正常趾甲的应力分布差异。

本发明提供的基于生物力学的嵌甲预测方法，所述嵌甲预测方法从生物力学角度出发，通过借鉴血管和心脏建模的软组织生长理论和骨骼建模的自适应生长理论，结合理论和实验建立了趾甲的生物力学模型，采用有限元软件abaqus的umat实现了该模型的仿真计算分析，对嵌甲形成机制及手术方法进行了探索。

本领域的技术人员容易理解，以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。