張竹筠 彭妙娟 黃能 程玉民

摘要：建立人體皮膚力學(xué)分析的多層皮膚有限元模型，對(duì)不同皮膚層分別采用超彈性和黏彈性與超彈性耦合的本構(gòu)關(guān)系，利用Abaqus對(duì)兩層和三層皮膚模型的非線性大變形問題進(jìn)行分析，模擬皮膚起皺，并與實(shí)驗(yàn)結(jié)果進(jìn)行比較。對(duì)比結(jié)果表明：與兩層模型相比，三層模型能更真實(shí)地模擬皺紋;在相同外力下，越緊致的皮膚產(chǎn)生的皺紋深度越淺，皮下組織越厚的皮膚產(chǎn)生的皺紋深度越深。數(shù)值分析結(jié)果符合實(shí)際情況。

關(guān)鍵詞：皮膚;皺紋;超彈性;黏彈性;變形

中圖分類號(hào)：R714.149;R322.991

文獻(xiàn)標(biāo)志碼：B

文章編號(hào)：1006-0871（2019）01-0062-08

0 引 言

皺紋是皮膚老化的結(jié)果。隨著年齡增長，人體皮膚細(xì)胞的分化會(huì)越來越慢，真皮漸漸變薄，真皮下的脂肪細(xì)胞開始萎縮，皮膚組織變得越來越松弛。當(dāng)皮膚失去彈性時(shí)，遭受擠壓后就不能再復(fù)原，從而造成皮膚下垂，產(chǎn)生皺紋。皺紋直接影響面部容貌，尤其是眼角的魚尾紋最能表現(xiàn)一個(gè)人的衰老。因此，了解皺紋的形成原因、減少皺紋的產(chǎn)生很有必要。

目前，國內(nèi)外關(guān)于人體皮膚的有限元研究已有很多。YOSHIDA等[1]建立人體皮膚傷口的有限元模型，并對(duì)傷口愈合過程進(jìn)行數(shù)值分析。GROVES等[2]使用老鼠皮膚和人體皮膚的力學(xué)參數(shù)，建立超彈性、各向異性的有限元模型。BISCHOFF等[3-4]認(rèn)為，在人體皮膚中，膠原網(wǎng)絡(luò)對(duì)各向異性非線性性質(zhì)有關(guān)鍵作用，以此作為出發(fā)點(diǎn)建立有限元模型，并重新建立各向異性、非線性的黏彈性材料基礎(chǔ)模型。RETEL等[5]建立臨床手術(shù)中菱形切口的有限元模型，并且分析這種切口術(shù)后生長過程產(chǎn)生的最大應(yīng)力。SHOEMAKER等[6]建立皮膚單層模型。MAGNENAT-THALMANN等[7]提出包含角質(zhì)層和真皮層的皮膚模型，并在簡單的二維模型中與真實(shí)皮膚行為進(jìn)行定性比較。CORMAC等[8]使用Abaqus模擬人體前臂的皮膚，并且找到志愿者進(jìn)行真實(shí)的皮膚擠壓實(shí)驗(yàn)，以此對(duì)比真實(shí)和模擬2種情況下皺紋的擬合程度。FLYNN等[9]對(duì)不同壓力下的同一塊皮膚產(chǎn)生的皺紋進(jìn)行實(shí)驗(yàn)，記錄卸載后的皮膚松弛狀況。KIRBY等[10]使用豬皮的泊松比和彈性模量，建立三維有限元模型。盧天健等[11]使用拉、壓等各種實(shí)驗(yàn)方法測(cè)量人體皮膚和不同的動(dòng)物皮膚的離體和活體力學(xué)參數(shù)，包括泊松比和彈性模量等?？紫榍宓萚12]使用Neo-Hooken本構(gòu)模型在Abaqus中完成人體皮膚建模，并且模擬蚊子口器刺入皮膚過程。孟慶領(lǐng)[13]使用五參數(shù)法重新定義真皮層，精確度高達(dá)96%，并且利用新的皮膚模型完成針刺皮膚的數(shù)值模擬。

雖然這些模型能夠模擬皮膚的幾種復(fù)雜力學(xué)特性，如各向異性和黏彈性，但上述學(xué)者通常認(rèn)為皮膚是均質(zhì)材料，忽略角質(zhì)層、真皮層和皮下組織各自的不同貢獻(xiàn)，還有一些多層模型忽略中間的真皮層或底層的皮下組織層。此外，在大多數(shù)模型中，活體皮膚的自然張力也沒有明確指出。

模擬皮膚起皺時(shí)，皮膚各層在模型中均應(yīng)是必不可少的。單層皮膚模型壓縮不會(huì)形成皺紋。為模擬皺紋，至少需要一層表示角質(zhì)層和一層表示真皮層的模型。SHOEMAKER等[6]和BISCHOFF等[4]的單層模型模擬皮膚雙向拉伸和應(yīng)力松弛，但無法模擬更復(fù)雜的變形，如皮膚皺紋。MAGNENAT-THALMANN等[7]的模型中包含角質(zhì)層和真皮層，但忽略皮下組織，皺紋的最大范圍不符合皮膚擠壓實(shí)驗(yàn)的趨勢(shì)，因此真皮層下還需要額外的一層。CORMAC等[8]的三層皮膚模型在考慮真皮層和皮下組織的層間關(guān)系時(shí)沒有加入相互作用。

本文提出三層皮膚有限元模型，在Abaqus中模擬壓縮過程中皮膚的起皺現(xiàn)象，比較在人體手臂實(shí)驗(yàn)中測(cè)量的皺紋與模擬中形成的皺紋，檢驗(yàn)多層模型。對(duì)有限元模型計(jì)算結(jié)果與實(shí)驗(yàn)結(jié)果進(jìn)行定量比較，對(duì)MAGNENAT-THALMANN等[7]的簡單二維模型與真實(shí)皮膚行為的定性比較進(jìn)行改進(jìn)。此外，模擬過程考慮層間的相互作用，研究影響人體皮膚皺紋深度的因素。分析數(shù)值結(jié)果可知：在相同外力作用下，越緊致的皮膚所產(chǎn)生的皺紋深度越淺，越胖的人產(chǎn)生的皺紋深度越深，符合實(shí)際情況。

1 皮膚各層本構(gòu)模型

人類皮膚由3層結(jié)構(gòu)構(gòu)成：最外層是角質(zhì)層，中間層是真皮層，最底層是皮下組織層。皮膚是復(fù)雜的非線性、超彈性和黏彈性材料，在活體內(nèi)處于張力狀態(tài)，并且不同的皮膚層力學(xué)性質(zhì)不同。為更真實(shí)地模擬皮膚起皺，任何模型都應(yīng)該考慮這些因素。

1.1 角質(zhì)層

根據(jù)KOUTROUPI等[14]的實(shí)驗(yàn)結(jié)論，角質(zhì)層可以看作是各向同性的超彈性材料，并且沒有黏彈性效應(yīng)。對(duì)于各向同性超彈性材料，可選擇Neo-Hooken的應(yīng)變能公式建立本構(gòu)方程，即

在相對(duì)濕度為96%時(shí)，角質(zhì)層的彈性模量E=24 MPa，C10=4 MPa，D1=0.010 13 MPa。[8]

1.2 真皮層

HENDRIKS等[15]研究認(rèn)為，真皮層是控制皮膚對(duì)外界載荷響應(yīng)的主要力學(xué)元素。假設(shè)真皮層為Mooney-Rivlin模型，模型的應(yīng)變能函數(shù)為

HENDRIKS等[15]結(jié)合在體抽吸實(shí)驗(yàn)和數(shù)據(jù)擬合分析，得到的C10和C01分別為9.4 kPa和82 kPa。

1.3 皮下組織

皮下組織又稱“皮下脂肪組織”，即脂肪層，具有復(fù)雜的復(fù)合材料性質(zhì)，即非線性、超彈性和黏彈性。皮下組織超彈性選用Yeoh模型建立本構(gòu)方程，即

根據(jù)FUNG[16]提出的準(zhǔn)線性黏彈性理論，皮下組織的黏彈性部分選用Prony級(jí)數(shù)表示的縮減松弛方程作為本構(gòu)模型，即

2 皮膚三層有限元模型

根據(jù)上述超彈性和黏彈性與超彈性耦合的本構(gòu)關(guān)系，在Abaqus中對(duì)2種模型在不同壓縮位移情況下的非線性大變形力學(xué)特性進(jìn)行計(jì)算。模型樣本的實(shí)際尺寸為90 mm×20 mm，見圖1。由于邊界條件對(duì)稱，數(shù)值模擬采用長度為45.0 mm、寬度為10.0 mm的皮膚模型，模型底部在z方向固定，在x和y的負(fù)方向施加2 kPa的拉力，并假定皮膚完全光滑。2種不同的有限元模型是角質(zhì)層、真皮層和皮下組織的不同組合，見表1。在模型中，角質(zhì)層厚度為20 μm，真皮層厚度為1.2 mm，皮下組織厚度為1.5 mm。此外，為研究皮膚皺紋深度與人體胖瘦的關(guān)系，采用相同厚度的角質(zhì)層和真皮層，分別與2.4和3.6 mm厚的皮下組織耦合，建立模型。

由于材料的特殊性質(zhì)，當(dāng)把角質(zhì)層、真皮層和皮下組織這3種材料約束在一起時(shí)，擠壓其表面會(huì)出現(xiàn)褶皺，這與實(shí)際情況下的人體皮膚一致。施加拉力的大小即人體皮膚的張力，此拉力是模擬人體皮膚緊致程度的重要考查標(biāo)準(zhǔn)。

第一個(gè)分析步：將拉力施加到光滑的皮膚上。在皮膚邊緣，即真皮層和皮下組織層，沿x和y的負(fù)方向施加2 kPa的拉力。為進(jìn)一步分析隨著年齡增長，皮膚老化對(duì)人體皺紋的影響，分析拉力大小為1和3 kPa的情況。

第二個(gè)分析步：模擬皮膚的面內(nèi)壓縮。將模型中左半部的上表面在z軸固定，然后對(duì)這部分的x方向施加載荷，模擬壓縮位移。在產(chǎn)生皺紋后，記錄z方向的最大應(yīng)變和最小應(yīng)變，即可知道皺紋的波峰和波谷的高度差，即皺紋最大范圍，也就是皺紋的深度。

在網(wǎng)格劃分中，角質(zhì)層采用4節(jié)點(diǎn)曲面薄殼單元S4R劃分單元，真皮層和皮下組織層采用8節(jié)點(diǎn)線性六面體單元C3D8R劃分單元。圖1模型共劃分為28 476個(gè)單元，接觸區(qū)域網(wǎng)格足夠細(xì)化，以確保高應(yīng)變、應(yīng)力區(qū)域計(jì)算結(jié)果的精確度。

3 計(jì)算結(jié)果

3.1 角質(zhì)層和真皮層皮膚兩層模型

在真皮層四周施加2 kPa的拉力，模擬皮膚的不同緊致程度，然后在模型左半部上表面沿x方向施加15～20 mm的壓縮位移，不同壓縮位移下皮膚兩層模型的位移云圖見圖2，皺紋最大范圍隨壓縮位移的變化曲線見圖3。由此可以看出：當(dāng)水平壓縮位移小于19 mm時(shí)，皺紋最大范圍隨壓縮位移的增加而增大，但當(dāng)水平壓縮位移大于19 mm時(shí)，皺紋最大壓縮范圍隨壓縮位移的增加而減小。

3.2 皮膚三層模型

3.2.1 皮下組織厚度為1.5 mm

在模型四周施加不同的拉力，模擬皮膚的不同緊致程度，然后在模型左半部上表面沿x方向施加15～20 mm的壓縮位移。由于材料的特殊物理性質(zhì)和層間關(guān)系，施加壓縮位移后會(huì)產(chǎn)生皺紋。

3.2.1.1 施加1 kPa的拉力

皮下組織厚度為1.5 mm、施加1 kPa拉力時(shí)，不同壓縮位移下皮膚三層模型的位移云圖見圖4，皺紋最大范圍隨壓縮位移的變化曲線見圖5。

3.2.1.2 施加2 kPa的拉力

皮下組織厚度為1.5 mm、施加2 kPa拉力時(shí)，不同壓縮位移下皮膚三層模型的位移云圖見圖6，皺紋最大范圍隨壓縮位移的變化曲線見圖7。皺紋產(chǎn)生以后，隨著壓縮位移的增加，皺紋最大范圍逐漸增大，且增長率基本保持不變。

3.2.1.3 施加3 kPa的拉力

皮下組織厚度為1.5 mm、施加3 kPa拉力時(shí)，不同壓縮位移下皮膚三層模型的位移云圖見圖8，皺紋最大范圍隨壓縮位移的變化曲線見圖9。皮下組織厚度為1.5 mm時(shí)，施加不同拉力的皮膚三層模型皺紋的最大范圍變化曲線見圖10。

由此可以看出：在壓縮位移相同的情況下，隨著皮膚模型中拉力的增大，生成皺紋的最大范圍減小。皮膚張力的增加可減少皮膚壓縮產(chǎn)生皺紋的最大范圍，即延緩皺紋的產(chǎn)生。由此表明，在皮膚模型中包含自然張力很重要。皮膚皺紋隨年齡的增加而增加，其原因之一就是皮膚自然張力的喪失。[17]考慮皮膚拉力提出的皮膚三層模型，能更好地模擬皮膚老化造成自然張力的損失。

3.2.2 皮下組織厚度為2.4 mm

皮下組織厚度為2.4 mm、施加2 kPa拉力時(shí)，不同壓縮位移下皮膚三層模型的位移云圖見圖11，皺紋最大范圍隨壓縮位移的變化曲線見圖12。壓縮位移小于16 mm時(shí)，皺紋最大范圍曲線斜率很大;壓縮位移為16mm時(shí)斜率放緩;隨后，隨著壓縮位移的增加，曲線斜率基本保持定值。

3.2.3 皮下組織厚度為3.6 mm

皮下組織厚度為3.6 mm、施加2 kPa拉力時(shí)，不同壓縮位移下皮膚三層模型的位移云圖見圖13，皺紋最大范圍隨壓縮位移的變化曲線見圖14。曲線變化情況與皮下組織厚度為2.4 mm時(shí)皮膚三層模型的模擬結(jié)果大致相同。相同厚度角質(zhì)層和真皮層分別與1.5 、2.4和3.6 mm皮下組織耦合的模擬結(jié)果見圖15。隨著皮下組織厚度增加，皺紋最大范圍增大。在皮下組織厚度為1.5 mm的模型中，產(chǎn)生皺紋以后，隨著壓縮位移的增加，皺紋最大范圍逐漸增大，且增長率基本保持不變。在皮下組織厚度為2.4和3.6 mm的模型中，開始時(shí)2條皺紋最大范圍曲線斜率很大，當(dāng)壓縮位移大于16 mm時(shí)斜率放緩，然后隨著壓縮位移的增加，曲線斜率基本保持定值。因此，皮膚的皮下組織越厚，或者說人越肥胖，越容易產(chǎn)生皺紋，符合實(shí)際情況。

皮膚兩層模型、三層模型的模擬結(jié)果與文獻(xiàn)[8]的前臂皮膚表面壓縮起皺實(shí)驗(yàn)結(jié)果對(duì)比見圖16。皮膚兩層模型的皺紋最大范圍與實(shí)驗(yàn)的皺紋最大范圍相差20.68%，三層模型的皺紋最大范圍與實(shí)驗(yàn)結(jié)果相差38.58%。雖然兩層模型與實(shí)驗(yàn)結(jié)果相差的百分比小于三層模型與實(shí)驗(yàn)的差異，但兩層模型模擬結(jié)果的趨勢(shì)與實(shí)驗(yàn)趨勢(shì)不一致。當(dāng)水平壓縮位移小于17 mm時(shí)，兩層和三層模型皺紋最大范圍結(jié)果相差不大，但施加更大的水平壓縮位移時(shí)，兩層模型模擬結(jié)果的曲線斜率減小至負(fù)值，而三層模型模擬結(jié)果的曲線斜率與實(shí)驗(yàn)測(cè)量一致，由此可定性說明本文有限元模型的有效性。

4 結(jié) 論

提出皮膚三層有限元模型，數(shù)值模擬皮膚壓縮起皺，并將模擬結(jié)果與實(shí)驗(yàn)結(jié)果進(jìn)行比較。該模型可比現(xiàn)有的單層和兩層模型更精確地模擬皮膚起皺，數(shù)值結(jié)果分析表明：

（1）本文所建立的角質(zhì)層、真皮層和皮下組織耦合的皮膚三層模型，在施加位移載荷后出現(xiàn)的物理變化特點(diǎn)符合實(shí)際情況，證明本文模型的準(zhǔn)確性。

（2）模型中施加皮膚張力是必要的：同樣厚度的皮膚，越緊致越難產(chǎn)生皺紋;在相同張力下，皮膚越厚，越容易產(chǎn)生皺紋。

（3）模擬皮膚皺紋時(shí)，皮膚各層以及層間的相互作用必不可少，準(zhǔn)確表示皮膚各層的力學(xué)特征也很有必要。更精確的本構(gòu)方程可以進(jìn)一步提高皮膚三層有限元模擬結(jié)果與實(shí)驗(yàn)結(jié)果的一致性。

此外，本文提出的模型假定角質(zhì)層表面光滑，使得初始皺紋深度與實(shí)驗(yàn)結(jié)果不一致;角質(zhì)層在實(shí)際狀態(tài)中有大量的毛孔，受到擠壓時(shí)會(huì)有應(yīng)力集中，這些也應(yīng)是建?？紤]的細(xì)節(jié)。

參考文獻(xiàn)：

[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] 盧天健， 徐峰. 皮膚的力學(xué)性能概述[J]. 力學(xué)進(jìn)展， 2008， 38（4）：393-426.

[12] 孔祥清， 吳承偉. 蚊子口針刺破人體皮膚過程的數(shù)值模擬[J]. 力學(xué)與實(shí)踐， 2010， 32（2）：90-94.

[13] 孟慶領(lǐng). 微針刺入皮膚的力學(xué)行為研究[D]. 大連：大連理工大學(xué)， 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.

（編輯 武曉英）