张竹筠 彭妙娟 黄能 程玉民

摘要：建立人体皮肤力学分析的多层皮肤有限元模型，对不同皮肤层分别采用超弹性和黏弹性与超弹性耦合的本构关系，利用Abaqus对两层和三层皮肤模型的非线性大变形问题进行分析，模拟皮肤起皱，并与实验结果进行比较。对比结果表明：与两层模型相比，三层模型能更真实地模拟皱纹;在相同外力下，越紧致的皮肤产生的皱纹深度越浅，皮下组织越厚的皮肤产生的皱纹深度越深。数值分析结果符合实际情况。

关键词：皮肤;皱纹;超弹性;黏弹性;变形

中图分类号：R714.149;R322.991

文献标志码：B

文章编号：1006-0871（2019）01-0062-08

0 引 言

皱纹是皮肤老化的结果。随着年龄增长，人体皮肤细胞的分化会越来越慢，真皮渐渐变薄，真皮下的脂肪细胞开始萎缩，皮肤组织变得越来越松弛。当皮肤失去弹性时，遭受挤压后就不能再复原，从而造成皮肤下垂，产生皱纹。皱纹直接影响面部容貌，尤其是眼角的鱼尾纹最能表现一个人的衰老。因此，了解皱纹的形成原因、减少皱纹的产生很有必要。

目前，国内外关于人体皮肤的有限元研究已有很多。YOSHIDA等[1]建立人体皮肤伤口的有限元模型，并对伤口愈合过程进行数值分析。GROVES等[2]使用老鼠皮肤和人体皮肤的力学参数，建立超弹性、各向异性的有限元模型。BISCHOFF等[3-4]认为，在人体皮肤中，胶原网络对各向异性非线性性质有关键作用，以此作为出发点建立有限元模型，并重新建立各向异性、非线性的黏弹性材料基础模型。RETEL等[5]建立临床手术中菱形切口的有限元模型，并且分析这种切口术后生长过程产生的最大应力。SHOEMAKER等[6]建立皮肤单层模型。MAGNENAT-THALMANN等[7]提出包含角质层和真皮层的皮肤模型，并在简单的二维模型中与真实皮肤行为进行定性比较。CORMAC等[8]使用Abaqus模拟人体前臂的皮肤，并且找到志愿者进行真实的皮肤挤压实验，以此对比真实和模拟2种情况下皱纹的拟合程度。FLYNN等[9]对不同压力下的同一块皮肤产生的皱纹进行实验，记录卸载后的皮肤松弛状况。KIRBY等[10]使用猪皮的泊松比和弹性模量，建立三维有限元模型。卢天健等[11]使用拉、压等各种实验方法测量人体皮肤和不同的动物皮肤的离体和活体力学参数，包括泊松比和弹性模量等。孔祥清等[12]使用Neo-Hooken本构模型在Abaqus中完成人体皮肤建模，并且模拟蚊子口器刺入皮肤过程。孟庆领[13]使用五参数法重新定义真皮层，精确度高达96%，并且利用新的皮肤模型完成针刺皮肤的数值模拟。

虽然这些模型能够模拟皮肤的几种复杂力学特性，如各向异性和黏弹性，但上述学者通常认为皮肤是均质材料，忽略角质层、真皮层和皮下组织各自的不同贡献，还有一些多层模型忽略中间的真皮层或底层的皮下组织层。此外，在大多数模型中，活体皮肤的自然张力也没有明确指出。

模拟皮肤起皱时，皮肤各层在模型中均应是必不可少的。单层皮肤模型压缩不会形成皱纹。为模拟皱纹，至少需要一层表示角质层和一层表示真皮层的模型。SHOEMAKER等[6]和BISCHOFF等[4]的单层模型模拟皮肤双向拉伸和应力松弛，但无法模拟更复杂的变形，如皮肤皱纹。MAGNENAT-THALMANN等[7]的模型中包含角质层和真皮层，但忽略皮下组织，皱纹的最大范围不符合皮肤挤压实验的趋势，因此真皮层下还需要额外的一层。CORMAC等[8]的三层皮肤模型在考虑真皮层和皮下组织的层间关系时没有加入相互作用。

本文提出三层皮肤有限元模型，在Abaqus中模拟压缩过程中皮肤的起皱现象，比较在人体手臂实验中测量的皱纹与模拟中形成的皱纹，检验多层模型。对有限元模型计算结果与实验结果进行定量比较，对MAGNENAT-THALMANN等[7]的简单二维模型与真实皮肤行为的定性比较进行改进。此外，模拟过程考虑层间的相互作用，研究影响人体皮肤皱纹深度的因素。分析数值结果可知：在相同外力作用下，越紧致的皮肤所产生的皱纹深度越浅，越胖的人产生的皱纹深度越深，符合实际情况。

1 皮肤各层本构模型

人类皮肤由3层结构构成：最外层是角质层，中间层是真皮层，最底层是皮下组织层。皮肤是复杂的非线性、超弹性和黏弹性材料，在活体内处于张力状态，并且不同的皮肤层力学性质不同。为更真实地模拟皮肤起皱，任何模型都应该考虑这些因素。

1.1 角质层

根据KOUTROUPI等[14]的实验结论，角质层可以看作是各向同性的超弹性材料，并且没有黏弹性效应。对于各向同性超弹性材料，可选择Neo-Hooken的应变能公式建立本构方程，即

在相对湿度为96%时，角质层的弹性模量E=24 MPa，C10=4 MPa，D1=0.010 13 MPa。[8]

1.2 真皮层

HENDRIKS等[15]研究认为，真皮层是控制皮肤对外界载荷响应的主要力学元素。假设真皮层为Mooney-Rivlin模型，模型的应变能函数为

HENDRIKS等[15]结合在体抽吸实验和数据拟合分析，得到的C10和C01分别为9.4 kPa和82 kPa。

1.3 皮下组织

皮下组织又称“皮下脂肪组织”，即脂肪层，具有复杂的复合材料性质，即非线性、超弹性和黏弹性。皮下組织超弹性选用Yeoh模型建立本构方程，即

根据FUNG[16]提出的准线性黏弹性理论，皮下组织的黏弹性部分选用Prony级数表示的缩减松弛方程作为本构模型，即

2 皮肤三层有限元模型

根据上述超弹性和黏弹性与超弹性耦合的本构关系，在Abaqus中对2种模型在不同压缩位移情况下的非线性大变形力学特性进行计算。模型样本的实际尺寸为90 mm×20 mm，见图1。由于边界条件对称，数值模拟采用长度为45.0 mm、宽度为10.0 mm的皮肤模型，模型底部在z方向固定，在x和y的负方向施加2 kPa的拉力，并假定皮肤完全光滑。2种不同的有限元模型是角质层、真皮层和皮下组织的不同组合，见表1。在模型中，角质层厚度为20 μm，真皮层厚度为1.2 mm，皮下组织厚度为1.5 mm。此外，为研究皮肤皱纹深度与人体胖瘦的关系，采用相同厚度的角质层和真皮层，分别与2.4和3.6 mm厚的皮下组织耦合，建立模型。

由于材料的特殊性质，当把角质层、真皮层和皮下组织这3种材料约束在一起时，挤压其表面会出现褶皱，这与实际情况下的人体皮肤一致。施加拉力的大小即人体皮肤的张力，此拉力是模拟人体皮肤紧致程度的重要考查标准。

第一个分析步：将拉力施加到光滑的皮肤上。在皮肤边缘，即真皮层和皮下组织层，沿x和y的负方向施加2 kPa的拉力。为进一步分析随着年龄增长，皮肤老化对人体皱纹的影响，分析拉力大小为1和3 kPa的情况。

第二个分析步：模拟皮肤的面内压缩。将模型中左半部的上表面在z轴固定，然后对这部分的x方向施加载荷，模拟压缩位移。在产生皱纹后，记录z方向的最大应变和最小应变，即可知道皱纹的波峰和波谷的高度差，即皱纹最大范围，也就是皱纹的深度。

在网格划分中，角质层采用4节点曲面薄壳单元S4R划分单元，真皮层和皮下组织层采用8节点线性六面体单元C3D8R划分单元。图1模型共划分为28 476个单元，接触区域网格足够细化，以确保高应变、应力区域计算结果的精确度。

3 计算结果

3.1 角质层和真皮层皮肤两层模型

在真皮层四周施加2 kPa的拉力，模擬皮肤的不同紧致程度，然后在模型左半部上表面沿x方向施加15～20 mm的压缩位移，不同压缩位移下皮肤两层模型的位移云图见图2，皱纹最大范围随压缩位移的变化曲线见图3。由此可以看出：当水平压缩位移小于19 mm时，皱纹最大范围随压缩位移的增加而增大，但当水平压缩位移大于19 mm时，皱纹最大压缩范围随压缩位移的增加而减小。

3.2 皮肤三层模型

3.2.1 皮下组织厚度为1.5 mm

在模型四周施加不同的拉力，模拟皮肤的不同紧致程度，然后在模型左半部上表面沿x方向施加15～20 mm的压缩位移。由于材料的特殊物理性质和层间关系，施加压缩位移后会产生皱纹。

3.2.1.1 施加1 kPa的拉力

皮下组织厚度为1.5 mm、施加1 kPa拉力时，不同压缩位移下皮肤三层模型的位移云图见图4，皱纹最大范围随压缩位移的变化曲线见图5。

3.2.1.2 施加2 kPa的拉力

皮下组织厚度为1.5 mm、施加2 kPa拉力时，不同压缩位移下皮肤三层模型的位移云图见图6，皱纹最大范围随压缩位移的变化曲线见图7。皱纹产生以后，随着压缩位移的增加，皱纹最大范围逐渐增大，且增长率基本保持不变。

3.2.1.3 施加3 kPa的拉力

皮下组织厚度为1.5 mm、施加3 kPa拉力时，不同压缩位移下皮肤三层模型的位移云图见图8，皱纹最大范围随压缩位移的变化曲线见图9。皮下组织厚度为1.5 mm时，施加不同拉力的皮肤三层模型皱纹的最大范围变化曲线见图10。

由此可以看出：在压缩位移相同的情况下，随着皮肤模型中拉力的增大，生成皱纹的最大范围减小。皮肤张力的增加可减少皮肤压缩产生皱纹的最大范围，即延缓皱纹的产生。由此表明，在皮肤模型中包含自然张力很重要。皮肤皱纹随年龄的增加而增加，其原因之一就是皮肤自然张力的丧失。[17]考虑皮肤拉力提出的皮肤三层模型，能更好地模拟皮肤老化造成自然张力的损失。

3.2.2 皮下组织厚度为2.4 mm

皮下组织厚度为2.4 mm、施加2 kPa拉力时，不同压缩位移下皮肤三层模型的位移云图见图11，皱纹最大范围随压缩位移的变化曲线见图12。压缩位移小于16 mm时，皱纹最大范围曲线斜率很大;压缩位移为16mm时斜率放缓;随后，随着压缩位移的增加，曲线斜率基本保持定值。

3.2.3 皮下组织厚度为3.6 mm

皮下组织厚度为3.6 mm、施加2 kPa拉力时，不同压缩位移下皮肤三层模型的位移云图见图13，皱纹最大范围随压缩位移的变化曲线见图14。曲线变化情况与皮下组织厚度为2.4 mm时皮肤三层模型的模拟结果大致相同。相同厚度角质层和真皮层分别与1.5 、2.4和3.6 mm皮下组织耦合的模拟结果见图15。随着皮下组织厚度增加，皱纹最大范围增大。在皮下组织厚度为1.5 mm的模型中，产生皱纹以后，随着压缩位移的增加，皱纹最大范围逐渐增大，且增长率基本保持不变。在皮下组织厚度为2.4和3.6 mm的模型中，开始时2条皱纹最大范围曲线斜率很大，当压缩位移大于16 mm时斜率放缓，然后随着压缩位移的增加，曲线斜率基本保持定值。因此，皮肤的皮下组织越厚，或者说人越肥胖，越容易产生皱纹，符合实际情况。

皮肤两层模型、三层模型的模拟结果与文献[8]的前臂皮肤表面压缩起皱实验结果对比见图16。皮肤两层模型的皱纹最大范围与实验的皱纹最大范围相差20.68%，三层模型的皱纹最大范围与实验结果相差38.58%。虽然两层模型与实验结果相差的百分比小于三层模型与实验的差异，但两层模型模拟结果的趋势与实验趋势不一致。当水平压缩位移小于17 mm时，两层和三层模型皱纹最大范围结果相差不大，但施加更大的水平压缩位移时，两层模型模拟结果的曲线斜率减小至负值，而三层模型模拟结果的曲线斜率与实验测量一致，由此可定性说明本文有限元模型的有效性。

4 结 论

提出皮肤三层有限元模型，数值模拟皮肤压缩起皱，并将模拟结果与实验结果进行比较。该模型可比现有的单层和两层模型更精确地模拟皮肤起皱，数值结果分析表明：

（1）本文所建立的角质层、真皮层和皮下组织耦合的皮肤三层模型，在施加位移载荷后出现的物理变化特点符合实际情况，证明本文模型的准确性。

（2）模型中施加皮肤张力是必要的：同样厚度的皮肤，越紧致越难产生皱纹;在相同张力下，皮肤越厚，越容易产生皱纹。

（3）模拟皮肤皱纹时，皮肤各层以及层间的相互作用必不可少，准确表示皮肤各层的力学特征也很有必要。更精确的本构方程可以进一步提高皮肤三层有限元模拟结果与实验结果的一致性。

此外，本文提出的模型假定角质层表面光滑，使得初始皱纹深度与实验结果不一致;角质层在实际状态中有大量的毛孔，受到挤压时会有应力集中，这些也应是建模考虑的细节。

参考文献：

[1] YOSHIDA H， TSUTSUMI S， MIZUNUMA M， et al. A surgical simulation system of skin sutures using a three-dimensional finite element method[J]. Clinical Biomechanics， 2001， 16（7）：621-626. DOI：10.1016/S0268-0033（01）00028-6.

[2] GROVES R B， COULMAN S A， BIRCHALL J C， et al. An anisotropic， hyperelastic model for skin：Experimental measurements， finite element modelling and identification of parameters for human and murine skin[J]. Journal of Mechanical Behavior of Biomedical Materials， 2013， 18（2）： 167-180. DOI：10.1016/j.jmbbm.2012.10.021.

[3] BISCHOFF J E， ARRUDA E M， GROSH K. Finite element modeling of human skin using an isotropic， nonlinear elastic constitutive model[J]. Journal of Biomechanics， 2000， 33（6）：645-652. DOI：10.1016/S0021-9290（00）00018-X.

[4] BISCHOFF J E， ARRUDA E M， GROSH K. A rheological network model for continuum anisotropic and viscoelastic behavior of soft tissue[J]. Biomechanics and Modeling in Mechanobiology， 2004， 3（1）：56-65. DOI：10.1007/s10237-004-0049-4.

[5] RETEL V， VESCOVO P， JACQUET E， et al. Nonlinear model of skin mechanical behaviour analysis with finite element method[J]. Skin Research and Technology， 2001， 7（3）：152-158. DOI：10.1034/j.1600-0846.2001.70302.x.

[6] SHOEMAKER P A， SCHNEIDER D， LEE M C， et al. A constitutive model for two-dimensional soft tissues and its application to experimental data[J]. Journal of Biomechanics， 1986， 19（9）：695-702. DOI：10.1016/0021-9290（86）90193-4.

[7] MAGNENAT-THALMANN N， KALRA P， LVQUE L J， et al. A computational skin model：Fold and wrinkle formation[J]. IEEE Transactions on Information Technology in Biomedicine， 2002， 6（4）：317-323. DOI：10.1109/TITB.2002.806097.

[8] CORMAC F， MCCORMACK B A. Finite element modelling of forearm skin wrinkling[J]. Skin Research and Technology， 2008， 14（3）：261-269. DOI：10.1111/j.1600-0846.2008.00289.x.

[9] FLYNN C， MCCORMACK B. Simulating the wrinkling and aging of skin with a multi-layer finite element model [J]. Journal of Biomechanics， 2010， 43（3）：442-448. DOI：

[10] KIRBY S D， WANG B， TO C W， et al. Nonlinear， three-dimensional finite-element model of skin biomechanics[J]. Journal of Otolaryngology， 1998， 27（3）：153-160.

[11] 盧天健， 徐峰. 皮肤的力学性能概述[J]. 力学进展， 2008， 38（4）：393-426.

[12] 孔祥清， 吴承伟. 蚊子口针刺破人体皮肤过程的数值模拟[J]. 力学与实践， 2010， 32（2）：90-94.

[13] 孟庆领. 微针刺入皮肤的力学行为研究[D]. 大连：大连理工大学， 2013.

[14] KOUTROUPI K S， BARBENEL J C. Mechanical and failure behaviour of the stratum corneum[J]. Journal of Biomechanics， 1990， 23（3）：281-287. DOI：10.1016/0021-9290（90）90018-X.

[15] HENDRIKS F M， BROKKEN D， VAN EEMEREN J T， et al. A numerical-experimental method to characterize non-linear mechanical behaviour of human skin[J]. Skin Research and Technology， 2003， 9（3）：274-283. DOI：10.1034/j.1600-0846.2003.00019.x.

[16] FUNG Y C. Biomechanics：Mechanical properties of living tussues[M]. New York：Springer， 1993.

[17] PIERARD G E， UHODA I， PIERARD-FRANCHIMONT C. From skin microrelief to wrinkles：An area ripe for investigation[J]. Journal of Cosmetic Dermatology， 2003， 2（1）：21-28. DOI：10.1111/j.1473-2130.2003.00012.x.

（编辑 武晓英）