Skin Research and Technology 2001; 7: 152–158

Copyright C Munksgaard 2001

Printed in Denmark. All rights reserved

Skin Research and Technology

ISSN 0909-752X

Nonlinear model of skin mechanical behaviour analysis

with finite element method

V. Retel, P. Vescovo, E. Jacquet, F. Trivaudey, D. Varchon and A. Burtheret

Laboratoire de Me´canique Applique´e R. Chaleat, Universite´ de Franche-Comte, Besanc¸on, France

Background/aims: A problem commonly encountered in plastic
and reconstructive surgery is the prediction of the stress put on
the skin when various types of skin flaps are used.
Methods/results: We developed a mathematical model based
on the finite element method, in order to determine the stress
field, by simulating the mechanical behaviour of human skin dur-
ing wound closure. We chose to take into account the low com-
pressive strength by modifying the mechanical parameters of the
model at each step of the calculation. The model has been de-
veloped and tested on a diamond-shaped incision and applied
on a Limberg incision.

S

URGICAL

techniques (1, 2) using skin flaps are in

constant evolution, mainly due to the surgeons

experimental know-how. The mathematical modelling
of the mechanical behaviour of the skin can help im-
prove incisions. The esthetic of and respect for the
functionality of the scar are among the main criteria
determining the incision shape. In some cases there
are complications (e.g., nodules, expansive scars and
necrosis). Adding to these complications, from a
mechanical point of view, are the large residual
stresses in the skin once the incision is closed. The
stress field depends on material mechanical character-
istics and boundary conditions, constrained by
wound closure. Some studies (3–8) show that it is
worth developing a mechanical behaviour model of
human skin to forecast the strain and stress fields in
the area around an incision. The study presented in
this paper is close to the one realized by Manios et
al. (9) in that both deal with an incision closure finite
element model. To account for the behaviour of hu-
man skin, we developed an isotropic, elastic, non-
linear model that gives a compression stiffness lower
than the traction stiffness. In the following text, we
call this model nonlinear. Skin flaps comply with the
laws of continuum mechanics that are fundamental to
the development of the finite element method. So we
introduced the nonlinear model in finite element soft-

152

Conclusion: The results presented are the Von Mises stress in
the area of skin around the scar and, when possible, the result-
ant closure force. They show the relevance of the method.

Key words: mechanical properties – human skin – finite element
method – nonlinear model

c

Munksgaard, 2001

Accepted for publication 26 January 2001

ware. After implementation, we obtained the stress
field. This model is compared with a diamond-shaped
incision with an isotropic and Hooke’s law elastic lin-
ear (10) model in which compression stiffness is the
same as traction stiffness. Then it is applied to a more
complex incision form.

Material and Methods

Human skin is an anisotropic, nonlinear, viscous ma-
terial (11, 12). In traction loading, it can sustain large
reversible strains [when the load is suppressed, after
a certain time, the material goes back to its initial
state – 25% (13)]. Because, of its structure, skin starts
to crease under very low compression loading. Our
aim is to forecast, as well as possible with a two-di-
mensional (2-D) model, the strain distribution in the
area of the skin close to the incision. Thus, we have
to consider the particular mechanical behaviour of
skin in compression. This was realized by considering
a special material, instead of a complex structure
whose behaviour needs a time-consuming, expensive
study. The nonlinear model is based on the finite ele-
ment software ANSYS54 (Anco Ingenieue, Venissi-
euse, France). The finite element techniqrie is a mod-
ern numerical method that permits the modelling of
complex structures by analysing them as an aggregate

Nonlinear model of skin mechanical behaviour

of small elements (14). A simple approach enabled the
validation of the method. Then we used it to predict
the mechanical behaviour of a complex skin flap.

Geometrical description of the model

We modelled the skin (epidermis and dermis) as a rec-
tangular plate of constant thickness (1 mm). A dia-
mond-shaped skin defect was chosen for the first
study (Fig. 1). The x and y directions are presented in
Fig. 1. We chose large dimensions for the plate, rela-
tive to the hole, in order to minimize the influence of
the boundary conditions on the area close to the in-
cision. It is a necessary condition for validating the
calculation. Because of the symmetry of the problem,
computational efficiency was improved by perform-
ing the calculations only on a quarter of the plate. The
model was based on the plane stress hypothesis: the
loading was applied in the plane of the skin and the
thickness of the plate plane was very low with respect
to the dimensions in the plane. This allowed us to
assume that the stress was homogeneous in the thick-
ness of the skin. A finite element meshing of 525
nodes and 158 quadrangular elements (PLANE 82)

Fig. 1. Geometry and mesh of the diamond-shaped skin defect.

153

was mapped on the representative quadrant. We
chose to study large strains and large displacements.

Boundary conditions

The nodes that belong to the edges of the diamond-
shaped hole were constrained to successive displace-
ments along the x axis until the complete closure of
the incision. So, the x direction is the direction of the
closing process. The total displacement was obtained
after 10 elementary displacements.

With respect to the symmetry of the problem, the

nodes belonging to the symmetry axes AAø and BBø
(Fig. 1) must stay on those axes and cannot be moved
perpendicular to the axis. The nodes of the external
boundaries were constrained (fixed in space).

Choice of material

In compression loading, the mechanical behaviour of
human skin is characterized by a low-resistance free
draping (8) and the compression stiffness is lower
than the traction stiffness. This particular behaviour
resembles that of a tissue and can be modelled by in-
troducing a buckling threshold [three-dimensional (3-
D) phenomenon] in the constitutive law. However, the
present study was based on 2-D modelling and the
plane stress assumption. To account for the effect of
the structure, an original material was created and as-
signed finite elements. The compressive behaviour of
this material differs from its behaviour in traction. But
this material does not exist in ANSYS software. There-
fore, a specific technique was developed to account
for this behaviour.

It has been estimated (15) that the modulus of elas-

ticity values of human skin vary from EΩ0.06 MPa to
EΩ4 MPa. Young’s modulus (E) represents the slope
of the stress-strain curve in the range defined for a
traction test. For orthotropic materials, the shear
modulus G

xy

is defined as:

G

xy

Ω

E

x

2(1πn

xy

)

with

n

xy

E

x

Ω

n

yx

E

y

where E

x

and E

y

are the components of Young’s

modulus in the x and y directions, respectively, and
where n

xy

and n

yx

represent the transverse Poisson’s

ratio when the traction is along the x and y directions,
respectively, as follows:

n

xy

Ω

ªe

y

e

x

where e

x

is the strain along the x direction and e

y

is

the strain along the y direction.

Based on these data, E

1

Ω1 MPa was eventually

Retel et al.

Fig. 2. Diamond-shaped defect, north-east quadrant. Compressed area along the x and y directions at steps 1 and 10 are shown.

chosen for the traction stiffness. The Poisson ratio was
chosen at nΩ0.4 (15). The compression stiffness was
chosen at E

2

Ω0.3 MPa. This value is the smallest one

that permitted the implementation.

We introduced four materials corresponding to the

different cases: Material 1 was assigned to the ele-
ments loaded in traction along the x and y directions.
Material 2 was assigned to the elements loaded in
compression along the x direction and in traction
along the y direction. Material 2

⬜ was assigned to the

elements loaded in compression along the y direction
and in traction along the x direction. Material 3 was
assigned to the elements loaded in compression along
the two directions. Materials 1 and 3 are isotropic and
linear. Materials 2 and 2

⬜ are orthotropic and linear.

Their characteristics are as follows:

154

Material 1:

E

x

ΩE

y

ΩE

1

n

xy

Ω0.4Ωn

yx

G

xy

Ω0.36 MPa

Material 2:

E

x

ΩE

1

, E

y

ΩE

2

n

xy

Ω0.4 n

yx

Ω0.12 G

xy

Ω0.3 MPa

Material 2

⬜: E

x

ΩE

2

, E

y

ΩE

1

n

xy

Ω0.12 n

yx

Ω0.4 G

xy

Ω0.3 MPa

Material 3:

E

x

ΩE

y

ΩE

2

n

xy

Ω0.12Ωn

yx

G

xy

Ω0.11 MPa

The method of assigning one material to one element
was based on a criterion defined below, permitting the
selection of elements in a systematic approach.

Criterion definition

During wound closure, the strains along the x and y
directions can be positive or negative, depending on
the nodes. We considered that an element was in com-
pression along the ‘‘i’’ axis if, after the first step of the
calculation, the average strain along the i axis (calcu-
lated by taking every vertex node of the element into

Nonlinear model of skin mechanical behaviour

account) had a negative sign and was 10% bigger than
the maximal positive strain of the nodes of the ele-
ment – that is to say:

冏

e

1

i

πe

2

i

πe

3

i

πe

4

i

4

冏

±10% ¡ e

max

where e

j

i

is the strain of the node j (jΩ1,2,3,4) along

the direction i and e

max

is the maximal strain among

the strains of every nodes, in both directions of the
plane after the first step of the calculation.

The choice of the criterion value was done in order

to give a realistic selection of elements in compression
in a diamond-shaped incision. The criterion made ap-
pear neither as isolated elements in compression nor
as elements in compression at the plate boundaries.
So we think that the corrections were efficient.

We took the following approach in four stages, for

each step:
1. do the calculation for this step;
2. pick out the elements that are in compression with

the criterion (Fig. 2);

3. do the calculation again at the same step after

modifying the material of the elements in com-
pression; and

4. go to the next step and resume the process at stage 1.
Comments: in a given step, the calculation at stage 1 is
a preliminary one and the calculation at stage 3 is the
final one. These two calculations, however, are
launched from the results of the final calculation of
the previous step. Before the preliminary calculation
of the first step, the same material (material 1) is
chosen for every element.

Actually, the global calculation is the succession of

the 10 final calculations corresponding to the 10 steps.
The 10 preliminary calculations are just used for pre-
dicting the elements in compression.

Results and Validation of the Method

Figure 3 shows the results of the calculations around
the diamond-shaped incision with the large strain and
large displacement assumptions. The resultant force
at the edge of the hole quadrant is the resultant force
needed to put the flaps edge to edge. It is called the
‘‘closure’’ force and is denoted by F. F

x

and F

y

are its

components along the x and y directions. F

x

Ωª3.12 N

and F

y

Ωª0.13 N.

Human skin is usually modelled in a finite element

calculation, by an isotropic elastic (10) material and a
thin plane shell. In this simple formulation, the com-
pression stiffness is equal to the traction stiffness and
the mechanical behaviour is the same in every direc-
tion in the plane of the skin.

155

Fig. 3. Diamond-shaped defect, north-east quadrant. The resulting
mesh with nonlinear elastic law is shown.

The calculation results presented above were com-

pared to the results obtained with the same meshing,
the same loading, the same boundary conditions and
a linear isotropic constitutive law for each element
(Young’s modulus EΩ1 MPa and Poisson’s coefficient
nΩ0.4). In the latter calculation, the ‘‘closure’’ force
was F

0

(F

0x

, F

0y

). F

0x

Ωª2.76 N and F

0y

Ωª0.18 N. The

discrepancy was low between the classical linear
model and our new method. The displacement fields
were almost the same. We actually expected to have
almost the same displacement fields because of the
small compression areas. This allowed us to say that
this original method led to consistent results and was
able to account for the specific mechanical behaviour
of human skin.

Application of the Method

Our original method was validated on a very simple
incision closure (diamond-shaped incision). The sec-
ond step consisted in applying it to a Limberg skin
flap (1). We chose an incision whose diamond-shaped
hole had the same dimensions as the hole of the
simple incision. The Limberg incision closure tech-
nique consists of moving the flap defined by lines 8,
9 and 10 (Fig. 4) upon the hole defined by lines, 7, 6

Retel et al.

Fig. 4. a) Geometry and mesh of Limberg defect. b) Evolution of Linberg flap during wound closure.

and 5 and then moving edge 11 to edge 12. Figure 4a
shows the geometry of Limberg incision, its position
in the plate specimen, the numbers of the lines defin-
ing the mesh geometry and the x and y directions.

The geometric model was developed with eight no-

dal quadrangular elements (PLANE82) and six nodal
triangular elements (TRIANGLE2). The model con-
tains 529 elements and 1163 nodes. The lack of sym-
metry in this complex incision prevented the simpli-
fication made in the previous case. So for this appli-
cation, the whole area around the incision was necess-
ary for the calculation.

Boundary conditions were applied at the edge of

the plate so that: the nodes belonging to lines 1 and

156

3 could not move along the y direction. The nodes
belonging to lines 2 and 4 could not move along the x
direction. The closure was calculated in 11 elementary
steps: First (5 steps), the nodes of the flap defined by
lines 8, 9 and 10 were moved in the x direction until
line 8 was superimposed on line 7. It came down to a
translation of the nodes of line 9 until line 8 reached
line 7 and line 10 was aligned with line 5 (Fig. 4b).
Next (5 steps), the nodes of lines 9 and 10 were moved
along line 5 until lines 9 and 10, respectively, were
superimposed on lines 6 and 5, respectively (Fig. 4b).
Then each node of lines 8, 9 and 10 was coupled with
the corresponding node on lines 7, 6 and 5 and, finally
(last step), the position constraints of the nodes be-

Nonlinear model of skin mechanical behaviour

longing to lines 8 and 7, 9 and 6, and 10 and 5 were
relaxed one after the other.

In this application, the behavioural law was based

on a method that had been established previously and
that took into account the low stiffness in compression
loading. The criterion was calculated with the average
of the strains at vertex nodes. The nonlinearity of the
problem made a great number of substeps necessary
for each increment in the calculation. In that case, the
resultant force could not be obtained, because the re-
sultant force was calculated from the constrained
nodes, and here no node was constrained at the end of
the last step. Therefore we chose the Von Mises stress
representation, which takes into account the stress in
the x and y directions at once and with the same ratio.
Engineers usually use Von Mises stress in order to
compare it to the plasticity threshold value (engineer-
ing material data), and in this way they know if there
are permanent strains (10). Von Mises stress results
appear in Fig. 5. The Von Mises stress maximal value
(S

Mx

) was 0.7 MPa.

We compared our study with one where the plain

linear elastic law led to a Von Mises stress maximal
value noticeably greater (S

Mx

Ω0.96 MPa, Fig. 6). The

present study showed that our method can be applied
to a real complex incision. In this particular case, the
area where the elements were in compression was

Fig. 5. Limberg defect, whole area. Von Mises stress field with non-
linear elastic law.

157

Fig. 6. Limberg defect, whole area. Von Mises stress field with linear
elastic law.

bigger than in the simple incision, which is why the
difference in the calculation results between the two
models was more significant.

The calculations for the two different skin flap types

showed that in the Limberg case the maximal stress
magnitude was lower than in the simple incision and
for the same plugged surface. With the nonlinear
model, we obtained a Von Mises stress maximal value
of 1.04 MPa for the diamond-shaped incision and 0.7
MPa for the Limberg incision. These results confirmed
the practitioners experimental know-how for the
choice of the incision shape.

Conclusions

Our method is an original approach for describing the
behaviour of skin when loading creates areas in com-
pression. The law used is based on a nonlinear elastic
approach; the nonlinearity is introduced directly in
the numerical model on elements in compression. The
operation is made manually and is based on a com-
pression criterion to select areas where the linear
model is not convenient. The first study of a simple
incision closure case proves the coherence of the
method. Its application to a Limberg incision shows
the importance of the new technique. We have made
possible the simulation of the specific mechanical be-

Retel et al.

haviour of skin during a wound closure. By imposing
hole dimensions, different incision forms can be
studied and compared from a mechanical point of
view. Furthermore, for each incision type, we can
change the dimensions and then optimize them. How-
ever, the modelling could be improved by introducing
the viscous behaviour of the skin and the mechanical
anisotropy in the plane of the skin.

Acknowledgments

This research was performed in collaboration with
Centre Hospitalier Universitaire. I wish to thank M.
Henry Paul Corral for his help.

References

1. Koss N. The mathematics of flaps. Pedicle flaps. In: Krizek

Tj, ed. Symposium on basic science and plastic surgery. St
Louis: Mosby, 1976.

2. Ohsumi N. A new skin suture technique for multiple Z-

Plasty. Plast Reconst Surg 1995: 96: 1713–1714.

3. Danielson DA. Human skin as an elastic membrane. J Bio-

mech 1973: 6: 539–546.

4. Chretien-Marquet B, Caillou V, Hartl Brasnu D, Bennaceur

S, Buisson T. Description of cutaneous excision and suture
using a mathematical model. Plast Reconstr Surg 1999: 103:
145–150.

5. Kawabata H, Kawai H, Masada K, Ono K. Computer-aided

analysis of Z-plasties. Plast Reconstr Surg 1989: 83: 319–325.

6. Wayne F, Larrabee JR. A finite element model of skin de-

formation. I. Biomechanics of skin and soft tissue: a review.
Laryngoscope 1986: 96: 399–405.

158

7. Wayne F, Larrabee JR, Wight Sutton D. A finite element

model of skin deformation. II. An experimental model of
skin deformation. Laryngoscope 1986: 96: 406–412.

8. Wayne F, Larrabee JR, Galt JA. A finite element model of

skin deformation. III. The finite element model. Laryngo-
scope 1986: 96: 413–419.

9. Manios A, Katsantonis J, Tosca A, Skulakis CH, Tsiftsis D.

The finite element method as a research and teaching tool
in the analysis of local skin flaps. Dermatol Surg 1996: 22:
1029–1033.

10. Lemaitre J, Chaboche JL. Mecanique des mate´riaux solides

Paris: Bordas, 1988.

11. Wijn PFF. The nonlinear viscoelastic properties of human

skin in vivo for small deformations. PhD Thesis, University
of the Netherlands, 1980.

12. Wijn PFF, Brakkee AJM, Stienen GJH, Vendrik AJH. Mechan-

ical properties of the human skin in vivo for small deforma-
tions; a comparison of uniaxial strain and torsion measure-
ments. London: Macmillan, 1976.

13. Agache P, Aubin F, Baron R, et al. Physiologie de la peau et

explorations fonctionnelles cutane´es. Cachan: Editions
Me´dicales Internationales, 2000: 408.

14. Zienkiewicz OC, Taylor RL. The finite element method, Vol.

1, 4th London: McGraw-Hill, 1989.

15. Agache P, Aubin F, Baron R, et al. Physiologie de la peau et

explorations fonctionnelles cutane´es. Cachan: Editions
Me´dicales Internationales, 2000: 664.

Address:
E. Jacquet
Laboratoire de Me´canique Applique´e R. Chale´at
Unite´ mixte de recherche CNRS-UFC 6604
24 Rue de l’Epitaphe
25030 Besanc¸on Cedex
France
Tel: [phone]
Fax: [phone]
e-mail: emmanuelle.jacquet/univ-fcomte.fr