© M. Abellan et al., Published by EDP Sciences, 2023

Introduction

The human skin is a living organ constantly renewed by new cells. Its physiology and its evolutions during aging influence its responses to coupled chemico-biological-thermo-hydro-mechanical influences. Moreover the response can be different depending whether the load is applied on the outer surface, in the volume through the surface or directly in the tissues. Understanding these mechanisms is of prime importance for surgical and medical purposes and in dermatology for both curative and care needs. Nowadays, bio-robots and devices are being developed and used to improve patient-designed therapies. They must rely on robust, objective, qualitative and quantitative parameters.

In vivo human skin has a multi-layered structure well adapted to its barrier functions. It is a stratified medium with four main layers: the stratum corneum (SC), the viable epidermis (VE), the dermis (DE) and the hypodermis (HY) [1–3]. The stratum corneum is made of hard dead corneocytes embedded in a lipid mortar [4, 5]. On the contrary, the viable epidermis contains evolving soft cells bathed in a physiological fluid. Its cohesion stems from the large amount of desmosomes that bridge the cells together. The dermis has dense fiber networks of collagen and elastin surrounded by an amorphous ground substance [6]. The adsorbed part of this ground substance provides cohesion between the fibers [7]. The hypodermis is composed of lipid cells separated by fibrous walls [8]. Together the fibers and the adsorbed ground substance initiate a structural process which gives skin its mechanical resistance and structural integrity. The strong desmosome bridges, the adsorbed part of the lipid mortar, and the fibrous walls also contribute to this structural process. These particularities make human skin a highly complex material (non-homogeneous, anisotropic, non-linear viscoelastic material subjected to pre-stress in vivo [9–12]) difficult to model.

In this context, attempts have been made in the past to simulate the mechanical properties of skin. Models, numerical simulations and/or innovative experiments have been developed. In vivo experimental tests include tensile testing [13], suction methods [14–16], torsion tests [7, 17], and indentation experiments [18, 19]. These works were confronted by the fact that the device had to be attached to the skin throughout the tests and thus disturbed its natural state. To overcome this drawback, the present study is based on contact-free creep tests performed with a contact-free indenter: the Tonoderm^® device [20, 21].

Constitutive models are based on simplifying assumptions that consider skin as either homogeneous isotropic linear elastic [15, 22] or nonlinear isotropic hyper-elastic [14, 19]. Bischoff et al. [23] coupled a hyperelastic model for fibrous tissue with a viscoelastic reptation model for soft materials. Azzez et al. [24] proposed using a homogeneous linear viscoelastic model to simulate indentation tests at two ages. To capture the in vivo behavior of the skin of the arm, Flynn et al. [25] used an Ogden strain energy function coupled with viscoelasticity and anisotropic pre-stress. Abellan et al. [26] proposed comparing four viscoelastic models to simulate indentation tests and to describe the coupled transient water flow and ion transport through skin. But these studies were still based on experimental contact tests.

This paper goes further by proposing a coupled theoretical-experimental-numerical simulation approach to study the influence of the dwell time of numerical contact-free creep tests on the viscoelastic responses of human skin in vivo at two ages. This study is intended to improve understanding of the state (healthy or diseased) of the skin before and after treatment in order to provide practitioners with clues for adapting the protocol of contact-free creep tests to the special needs of their patients, based on the knowledge of their state to design the appropriate curative/care therapy.

The experimental contact-free creep tests were performed on the volar forearm of two volunteers: a young adult and an elderly subject. An inverse approach was taken to determine the viscoelastic characteristics of the skin of both subjects. To ensure admissible physical consistency between the experimental data and the numerical simulations, the numerical tests were based on these experimental data. Five numerical contact-free creep tests were formulated with dwell times of 0.1 s, 0.5 s, 1 s, 2 s and 100 s, respectively.

The 3-D numerical simulations of the five numerical contact-free creep tests were conducted using the SYSTUS^® software. The skin is seen as a stratified medium described by a viscoelastic law with a large displacement formulation. This model allows a highly nonlinear description of the skin’s behavior. The main advantage is to gain insights into the mechanical properties of each layer. However, this is not the sole goal. It appears that a simple observation of the response of the skin surface to an air flow provides knowledge of what occurs in the volume. The mechanical theoretical model considered in conjunction with the constitutive viscoelastic law presented in this paper makes it possible to analyze the evolutions of each layer on its own and the relative influences of each layer on the others. The aim is to deduce the contribution of each layer, or more generally of each material point of the human skin, to the overall response of the skin to the different numerical contact-free creep tests, and to find out if certain people respond differently as a function of dwell time and/or age.

In the following, the materials and methods are presented. In particular the contact-free indenter is described and the numerical contact-free creep tests used to study the possible influence of the dwell time chosen are defined. The theoretical model, the numerical simulations and the viscoelastic parameters estimation procedure used in this work are presented. Then the results of the numerical simulations are given. A discussion analyses the main features of the results. The conclusion will end this paper.

Materials and methods

Contact-free creep tests in vivo

Experimental device

Air flow as an external load applied has been used for a long time by ophtalmologists in order to measure eye pressure. Based on this principle, the contact-free indenters (Tonoderm^® and Waveskin^®) at the Laboratoire de Tribologie et Dynamique des Systèmes (LTDS, Lyon, France), has been developed by Prof. H. Zahouani’s team [12, 20]. This device allows performing very difficult tests even on damaged or painful skins in vivo. Indeed, in such cases, it is compulsory to avoid any contact with the skin tissue in order to prevent contamination or destruction of the default studied. Moreover, in the case of contact-free indentations, the load applied can be removed instantaneously. It is then possible to observe and measure the undisturbed and free evolutions of the skin in vivo coming back to its natural stress state.

The experimental set-up of the Tonoderm^® and the sketch of its network are given in Figures 1a and 1b, respectively. A compressor (A) provides the source of compressed air necessary to the system. The regulation of this compressed air is firstly done by two expansion valves (not shown on Fig. 1) connected in series. They convert the pressure of 6 bars provided by the compressor at the working pressure of the regulator which is 4.5 bars. The choice of two instead of one regulator is based on the fact that such assembly can effectively minimize the change in the output pressure due to a variation of the inlet pressure. The air flow is then controlled by a mass flow controller (B) with a flow range from 2 to 100 NL/min (Normalized litters per minute) expressed in normal conditions of pressure P_N (Pa) and temperature T_N (K) for a standard density ρ_N (kg/Nm³) which are for the air: T_N = 273.15 °K, P_N = 101,325 Pa and ρ_N = 1.293 kg/Nm³. At this regulator’s exit are connected two solenoid valves, the first one (C) aims to empty the air network before a test. The second one (D) is placed right before the exit of the air flow. Therefore this air flow can be regulated (diameter 2 mm, distance exit-skin 10 mm). The air flow reaches the skin surface with an inverted-bell shape (Fig. 1c). It is equivalent to a spherical indenter with a tip surface of 6.28E-6 m² (Fig. 1d). The deflection is measured by a laser beam (E) (LK-G82, Keyence). The experimental set-up is controlled with a computer (F) and the experimental data are recorded with an analogical/numerical card (16 bits, 6221 M, National Instruments).

Figure 1 (a) Principle of the Tonoderm® [20]. (b) Forearm set-up [31]. (c) Photo of the pressured area of the skin surface under the air flow [20]. (d) Inverted-bell shape of the air flow and typical contact profile which may develop.

Protocol and experimental data

The experimental tests of this study are based on the following protocol:

• Before the tests, the subjects were resting seated in a medical chair during 10–15 min in a room with controlled temperature (22 °C) and humidity (RH 50%).

• The imposed force tests were performed on the inner forearm of the volunteers.

• The tested area was 6 cm from the elbow.

• For the tests, the system applies an air flow onto the outer surface of the forearm until the pressure equivalent to a flux of 20 NL/min is obtained (loading branch of 9.1 s). Then the flux is kept constant during 1 s. Finally the load is removed (the flux is stopped) and the system discharges. The free-return is monitored for 14 s. The mean of 3 measurements has been done.

The experimental data recorded during the indentation tests give profiles of the load (Fig. 2) and of the deflection with respect to time.

Figure 2 Experimental loading curve as a function of time (Exp).

Numerical contact-free creep tests

The possible influence of the dwell time chosen when performing the in vivo contact-free creep tests described before is analyzed based on five numerical contact-free creep tests numerically designed and presented here. This is done in order to ensure admissible physical consistency between the experimental records and the numerical findings.

The protocols of the five numerical contact-free creep tests are: a loading branch of 9.1 s; dwell times of respectively 0.1 s, 0.5 s, 1 s, 2 s and 100 s; then the applied pressure simulating the air flow at the upper surface of the REV is stopped and the free-return is monitored for 3600 s for each test.

The choice of monitoring the free-return for 3600 s lies on the observation of the numerical results obtained. After 3600 s, the human skin has completed its free-return. There is no more evolution in the responses of the skin. It will be shown, discussed later in the paper, that the numerical results exhibit close hysteresis cycles after 3600 s both for the young subject and for the elderly subject.

The loading curves of the numerical tests are given in Figure 3.

Figure 3 Loading curves of the five numerical contact-free creep tests with dwell times of 0.1 s, 0.5 s, 1 s, 2 s and 100 s.

Theoretical model

Within the framework of the thermodynamics of irreversible processes [27] and following Bathe [28], in an updated-Lagrangian formulation, when neglecting inertia forces, the total balance of momentum reads:

$$\underset{\Omega }{\overset{}{\int }}\mathrm{div} \mathit{S}\left(\mathit{x}, t\right)\mathrm{d}\Omega =\mathbf{0},$$(1)where S(x, t) is the second stress tensor of Piola–Kirchhoff of the skin in vivo in Pa, Ω is the configuration of the REV gathering all the material points M(t) defining the REV for the updated-Lagrangian description and x are the spatial coordinates of the material point M(t). Following Latorre and Montans [29], the constitutive equations are derived here in terms of the strain measure: the Eulerian rate of deformation tensor D and the corresponding constitutive stress measure: the (objective) Green–Naghdi stress rate σ^△ defined by:

$${\mathit{\sigma }}^{△}=\dot{\mathit{\sigma }}+\mathit{\sigma }.\dot{\mathit{R}}.{\mathit{R}}^T-\dot{\mathit{R}}.{\mathit{R}}^T.\mathit{\sigma }.$$(2)R being the rotation tensor coming from the polar decomposition of the deformation gradient.

Skin is supposed to follow an isotropic behavior and a Zener rheological model is chosen for modelling the skin behavior (Fig. 4). The Zener model is constituted by an elastic branch in parallel with a viscoelastic Maxwell’s branch. Therefore, the Cauchy stress is defined as the sum of an elastic part σ^E and a viscoelastic Maxwell’s part σ^M. It reads: σ = σ^E + σ^M.

Figure 4 Rheological Zener model.

In the elastic branch, the generalized Hooke’s law (Eq. (2)) is written in a rate form:

$${{\sigma }^E}_{\mathrm{ij}}^{△}=\left[\frac{E_{\infty } \nu }{\left(1+\nu \right)\left(1-2\nu \right)} {\delta }_{\mathrm{ij}} {\delta }_{\mathrm{rs}}+ \frac{E_{\infty }}{2\left(1+\nu \right)} \left({\delta }_{\mathrm{ir}} {\delta }_{\mathrm{js}}+{\delta }_{\mathrm{is}} {\delta }_{\mathrm{jr}}\right)\right]D_{\mathrm{rs}},$$(3)where ν is the Poisson’s ratio, E_∞ is the long term modulus in Pa, δ is the Kronecker delta.

The long term modulus E_∞ is defined by:

$$\frac{1}{E_{\infty }}= \frac{1}{E_{\mathrm{SC}}}+\frac{1}{E_{\mathrm{VE}}}+\frac{1}{E_{\mathrm{DE}}}+\frac{1}{E_{\mathrm{HY}}},$$(4)where E_i is the Young modulus of layer i, i standing for layers SC, VE, DE and HY, respectively.

In the viscoelastic Maxwell’s branch, for an isotropic material, we suppose an additive decomposition of the Eulerian rate of deformation tensor:

$$\mathit{D}={\mathit{D}}^e+ {\mathit{D}}^v,$$(5)

The stress-strain rate relations then separate volumetric and deviatoric responses of the skin in vivo. Toward this end, the Maxwell part σ^M of the Cauchy stress tensor and the Eulerian rate of deformation tensor D are decomposed into a spheric and a deviatoric part:

$${\sigma }_{\mathrm{ij}}^M={\sigma }_{\mathrm{ij}}^{\mathrm{\text{'}M}}+{\sigma }_m^M {\delta }_{\mathrm{ij}}\hspace{1em}\mathrm{with}\hspace{1em}{\sigma }_m^M=\frac{1}{3}{\sigma }_{\mathrm{kk}}^M,$$(6)

$$D_{\mathrm{ij}}=D_{\mathrm{ij}}^{\text{'}}+D_m {\delta }_{\mathrm{ij}}\hspace{1em}\mathrm{with}\hspace{1em}{D}_m=\frac{1}{3}{D}_{\mathrm{kk}},$$(7)where σ′^M and ${\sigma }_m^M$ are respectively the deviatoric part of the stress tensor and the mean stress in the Maxwell’s branch, D′ is the deviatoric part of the eulerian rate of deformation tensor and D_m, the volumetric part.

The volumetric and deviatoric responses of the material are then expressed by:

$${\sigma }_m^M+ \frac{{\eta }_v}{3K^M} {\sigma }_m^{M^{△}}={\eta }_v D_m,$$(8)

$${\sigma }_{\mathrm{ij}}^{\mathrm{\text{'}M}}+ \frac{{\eta }_d}{2{ G}^M} {\sigma }_{\mathrm{ij}}^{{\mathrm{\text{'}M}}^{△}}={\eta }_d D_{\mathrm{ij}}^{\text{'}},$$(9)with η_v the volumetric viscosity in Pa s, ${ K}^M=\frac{E_1}{3\left(1-2\nu \right)}$, the compressibility modulus in Pa, η_d the deviatoric viscosity in Pa s, G^M = E₁/2(1 + ν) the shear modulus in Pa and E₁ the modulus in the Maxwell’s branch in Pa.

The balance equation (1) and the behavior law (2), (8) and (9) are complemented by physically admissible boundary conditions set on additional parts of the boundary ∂Ω of the REV:

$$\mathit{\sigma }\left(\mathit{x}, t\right)\bullet \mathit{n}= {\mathit{t}}_p\hspace{1em}\mathrm{and}\hspace{1em}\mathit{u}\left(\mathit{x}, t\right)= {\mathit{u}}_p,$$(10)with ∂Ω = ∂Ω_t ∪ ∂Ω_u and ∂Ω_t ∩ ∂Ω_u = ∅, σ(x, t) the Cauchy stress tensor, n the outward unit normal of the boundary ∂Ω_t, t_p the prescribed effort on the boundary ∂Ω_t in Pa, u(x, t) the displacement vector of the material point M(t) in meter and u_p the prescribed displacement on the boundary ∂Ω_u in meter.

Numerical simulations

The theoretical model is applied to the numerical simulations of the experimental contact-free creep tests in vivo and also to the numerical simulations of the numerical contact-free creep tests. The boundary value problem is solved, 6 times for the young subject and 6 times for the elderly subject (one time for each protocol for each subject), using the finite element method with the SYSTUS^® software.

When defining the representative elementary volume (REV), the structural architecture spread in all layers is of importance for defining a user-defined REV and work at the soft cutaneous tissue scale. The thicknesses of the forearm were extracted from echograph images available at the Laboratory for a panel of 28 volunteers (14 women 23.2 ± 1.6 years-old and 14 women 60.4 ± 2.4 years-old) [20] from which were extracted the records of the young and the elderly volunteers of this study. It has been found for the thicknesses of the three layers SC + VE + DE, a mean of 1.094 mm ± 0.128 mm for the thickness of the young volunteers and 1.033 mm ± 0.108 mm for the thickness of the elderly volunteers. An analysis of variance conducted on these results has led to a statistically non-significant difference between the two groups. These results are similar to the measurements reported by Agache [7]. Hence, the user-designed REV considered here (Fig. 5) is a stratified volume composed of four layers: SC, VE, DE and HY with dimensions 30 mm × 30 mm × 11.11 mm where the thicknesses of each layer is taken equal to: 0.01 mm for SC, 0.1 mm for VE, 1 mm for DE and 10 mm for HY. The same REV appears to be admissible for both subjects.

Figure 5 User-designed REV deduced from the echograph images and 3-D mesh.

A large displacement analysis with an updated Lagrangian formulation allows handling the geometrical non-linearities as well as the non-linearities resulting from the behavior law. For the time increment, the choice was made of 0.1 s for the loading branch, 0.001 s during the dwell time and 1 s for monitoring the free-return. These choices for the formulations and for the time increments ensured the convergence of the computations.

For symmetry reasons, all the simulations were carried out on ¼ of the REV (Fig. 5) with dimensions 15 mm × 15 mm × 11.11 mm (Fig. 6). Furthermore, it is important to work with a mesh adapted to the mechanical settings of our model. In order to do so, SC, VE and DE are meshed with prisms while HY is meshed with tetrahedrons. The numerical results have not shown any mesh dependence on the elements chosen. The mesh is refined in and close to the loading zone (blue front corner in Fig. 6) for reasons of convergence and precision. The mesh consists of 9689 nodes spread into 25,181 elements.

Figure 6 3-D mesh with the boundary conditions.

The boundary conditions prescribed are: a Gaussian function P₀ is used to model the pressure distribution given in Figures 2 and 3. It is applied on the top corner (blue zone) of the mesh (Fig. 6); a zero displacement vector for all the nodes of the bottom surface; zero horizontal components in both directions for the displacement vector for all the nodes of the vertical surfaces for reasons of symmetry (Fig. 6).

The mechanical characteristics should be experimentally determined. For the Poisson ratio, Zahouani et al. [30] have reported of contact indentation and static friction tests performed on volunteers aged 55–70 years. They found a Poisson ratio between 0.44 and 0.48. Delalleau et al. [18] performed numerical simulations of contact indentation tests which lead to a Poisson ratio of 0.48. In a first theoretical model [21], we have worked with a Poisson ratio of 0.3 for obtaining convergence of the numerical process. In 2017, Ayadh et al. [31] proposed an extension of the first model with a Poisson ratio of 0.45 which presented no problem of convergence. Based on that information, the choice is done here to take a Poisson ratio of 0.45 for all the layers both for the young subject and for the elderly subject. The Poisson ratio should be complemented with viscoelastic parameters experimentally determined using the inverse procedure presented here after.

Inverse procedure

The input file supplies the finite element program with initial values for all the parameters. The relative error e_rel between the experimental and the numerical deflection is calculated using Equation (11):

$$e_{\mathrm{rel}}=\frac{\sum u_{\exp }- \sum u_{\mathrm{num}}}{\sum u_{\exp }},$$(11)where u_exp and u_num are respectively the experimental and the numerical displacements in meter. If this error e_rel is less than a given tolerance, the procedure stops: the parameters used in the calculation are the parameters sought. If the error e_rel is greater than the tolerance, the procedure continues with new parameters. The procedure stops when the calculations converge at an estimate of the parameters that minimize the error e_rel. Validation of this approach with a sensitivity analysis has already been documented [32].

Results

Experimental contact-free creep records in vivo

One volunteer is chosen in each group, i.e. a young healthy adult 23 years-old and an elderly subject 62 years-old. Figure 7 presents the deflection of the outer skin surface of the forearms in the test zone recorded for the young and the elderly subjects.

Figure 7 Profiles of the experimental deflection as a function of time for the young subject (blue line) and for the elderly subject (red line) for a flux of 20 NL/min.

Determination of the mechanical characteristics for the young and the elderly subject: inverse procedure

The estimation of the mechanical parameters for the young and the elderly subjects was carried out using the inverse procedure presented in section “Materials and methods”. Figure 8 shows an example of comparisons between the experimental and the numerical displacements of the central point under the air flow for three different sets of parameters. The relative errors are also calculated and drawn in Figure 8 (Bottom) for the different sets of Figure 8 (Top). It can be noticed the quite good fit between the experimental and the numerical curve 5 for the young subject and between the experimental and the numerical curve 111 for the elderly subject. Table 1 gives the sets of mechanical parameters associated with the numerical results of Figure 8. The set 5 for the young subject and the set 111 for the elderly subject are showing the best fit with the minimum error. They will provide the viscoelastic parameters used in the numerical simulations of the numerical contact-free creep tests for the young and the elderly subjects presented here after.

Figure 8 (Top) Comparison of the experimental and numerical curves of the displacement for three sets of parameters (Bottom) Curves of the relative error for the same sets of parameters.

Numerical results of the numerical simulations of the numerical contact-free creep tests

Taking into account the mechanical parameters estimated in Table 1 for each volunteer, ten boundary value problems (five for the young subject and five for the elderly subject) are now solved for the boundary conditions given in Figures 2 and 3. A key point appeared when studying the numerical results. For a given layer, each material point of the layer exhibits the same kind of stress state as the point of the layer belonging to the central vertical under the air flux but with lower amplitudes. For this reason, the obtained numerical results are presented and analyzed in terms of evolutions of layers exemplified with the numerical results of seven nodes of the central vertical under the air flux: top node of SC (6966), middle node of SC (5832), boundary between SC and VE (4688), boundary between VE and DE (4116), middle node of DE (3515), boundary between DE and HY (2920) and middle node of HY (4251) (Fig. 9).

Figure 9 Sketch giving the position of the seven nodes of the central vertical under the air flux used in the analysis.

Figure 10 presents the numerical isovalues of the vertical component of the displacement of the central point under the air flux at the end of the dwell times 0.1 s, 0.5 s and 1 s for the young subject. It gives an example of what happens in the volume and more specifically inside the layers for the young subject. The top surface sinks under the applied air flow. The large blue zone faces the maximum applied load and its response gives the highest values for the displacements. Then for a given layer, the amplitudes of the displacements decrease slowly for material points becoming increasingly distant from the central vertical. At the end of the different dwell times 0.1 s, 0.5 s, 1 s, 2 s and 100 s, the isovalues have the same shape as the ones proposed in Figures 10a–10c for the three dwell times 0.1 s, 0.5 s and 1 s. Moreover the isovalues exhibit equal amplitudes with a maximum deflection of the central point under the air flow of 1.038 mm.

Figure 10 Isovalues of the vertical displacements at the end of the dwell times: (a) 0.1 s, (b) 0.5 s and (c) 1 s for the young subject. (d) Contour lines of the mean stress for the dwell time of 0.5 s for the elderly subject.

Equivalently for the elderly subject, the most highly stressed area is under the air flow. The influence of the external mechanical applied load decreases with increasing distance (Fig. 10d). However depending on the dwell time the maximum deflection of the central point under the air flow is different. It equals 1.246 mm for 0.1 s, 1.251 mm for 0.5 s, 1.255 mm for 1 s, 1.264 mm for 2 s and 1.470 mm for 100 s.

Figure 11 presents the numerical normal component of the stress tensor σ_zz drawn as a function of the numerical deflection for the central point of the top surface of SC under the air flux. It should be noticed the peculiar shape of the graphs of these numerical results available, in Civil Engineering, they are called “hysteresis cycle”. These hysteresis cycles are given for the young subject and for the elderly subject for the dwell-times 0.1 s, 0.5 s, 1 s, 2 s and 100 s These hysteresis cycles give information on the exchanges of energy taking place between the REV and the outside when the studied medium is subjected to an external mechanical load, i.e. when the material points of the REV are provided with external mechanical energy. They are available here because the theoretical model considers that the responses of the cutaneous soft tissues of the different layers of the REV can be modelled reasonably well with the rheological Zener model. As said before, this viscoelastic model treated with a large displacement analysis with an updated Lagrangian formulation allows to follow the non-linear evolutions of the responses and in particular the hysteresis cycles.

Figure 11 Numerical normal stress σzz as a function of the numerical deflection of the central point under the air flux for the young subject and for the elderly subject for the dwell times of 0.1 s, 0.5 s, 1 s, 2 s and 100 s.

Figure 12 gives the comparison between the numerical normal component of the stress tensor σ_zz as a function of the numerical deflection for dwell times 1 s and 100 s for the seven nodes of the central vertical under the air flux. It displays the results obtained for the young subject (Fig. 12a) and for the elderly subject (Fig. 12b).

Figure 12 Numerical normal stress σzz as a function of the numerical deflection for the dwell times 1 s and 100 s for seven nodes of the central vertical under the air flux: (a) for the young subject, (b) for the elderly subject.

During the numerical contact-free creep tests, the pressure distribution (simulating an air flow) applied at the top corner of the mesh provides mechanical energy to the material points of the REV for both subjects. For the dwell times 1 s and 100 s, the amounts of energy exchanged are extracted from the hysteresis curves of Figure 12 for the young subject and for the elderly subject. These amounts are given in Table 2 for the seven nodes of the central vertical under the air flow for the two subjects. The order of amplitude of the values given in Table 2 should be analyzed keeping in mind the dimensions of the REV (15 mm × 15 mm × 11.11 mm) and the amplitude of the external mechanical load (from 0 Pa to 7277 Pa). It should be noticed that for the elderly subject, the amount of mechanical energy of the dwell time of 100 s is 24.17 times higher than the amount for the dwell time of 1 s for the top of the SC, 21.06 times higher for the middle of SC, 20.51 times higher for the boundary SC-VE, 21.51 times higher for the boundary VE-DE, 20.14 times higher for the middle of DE, 17.62 times higher for the boundary DE-HY and 18.5 times higher for the middle of HY.

Discussion

As explained when introducing the numerical results, they are analyzed in terms of evolutions of layers exemplified with the numerical results of seven nodes of the central vertical under the air flux.

In civil engineering, hysteresis cycles (Figs. 11 and 12) are studied to quantify the amount of mechanical energy supplied to the material during loading and the amount of mechanical energy given back or dissipated by the material during unloading.

Here the hysteresis cycles of the young subject for the central point of the top surface of SC under the air flux (Fig. 11) overlap with each other whatever the dwell time during which the load was applied. Each test supplied the same amount of mechanical energy (equal to 1.01 mm²) to the top surface of SC. Further there is no storage of energy during the creep tests i.e. the starting and ending points of the cycles are the same. These mechanical findings mean that the evolutions of the material points of the top surface of SC are reversible and unique. In other words, let us consider a material point of the top surface of SC. Before one numerical creep test, it is in a given unknown initial state. During the test, the mechanical energy supplied by the load changes this state. After the test, the free evolution of this material point causes it to return to what appears to be its initial state (Fig. 11). This observation can be done for all material points of the top surface of SC and for all dwell-times and whatever the amount of energy provided.

At the opposite, the hysteresis cycles (Fig. 11) of the elderly subject appear to be dwell time dependent. The starting and ending points are still the same for the five tests. However the surface of the cycles of the top node of SC increases when the dwell time increases. As the amount of mechanical energy provided by the test is equal to the surface of the cycle, increasing the surfaces is equivalent to increasing these amounts: 0.77 mm² for 0.1 s, 0.86 mm² for 0.5 s, 0.97 mm² for 1 s, 1.20 mm² for 2 s and 23.45 mm² for 100 s. But the material points of the top surface of SC return to their initial state, i.e. the dissipated mechanical energy is equal to the mechanical energy supplied. Hence, each material point of the top surface of SC of the elderly subject uses all the supplied energy to return to its initial state after the load has been removed.

These findings are coherent with the qualitative observations done during the tests for both subjects i.e. the outer surface of the skin seems to come back to its initial position more or less quickly after the tests. The plus of the numerical simulations is to quantify the amount of energy exchanged during the tests and to be able to confirm that the skin outer surface comes really back to its mechanical initial state after the tests whatever can be the amount of energy provided. The numerical results have also shown that the evolutions of the outer surface of the skin are not influenced by the dwell time for the young subject while they are dwell time dependent for the elderly subject.

These questions remain open for the evolutions in the volume and will be discussed here after based on the other available numerical results.

It can be observed (Fig. 12) that the starting and ending points of a given hysteresis cycle are the same for both subjects. However each sub-layer of the skin appears to have its own particular evolution as emphasized by the different hysteresis cycles drawn in Figure 12. They all contribute differently to the overall response of the skin. In a sense, mechanically speaking, each sub-layer of the skin starts from an unknown initial state and returns to it.

For the young subject (Fig. 12a) each layer follows overlapping paths whatever the dwell time. This observation is strengthened for the two dwell times (1 s and 100 s) by the amounts of energy exchanged in each case by the seven nodes of the central vertical (Table 2). It should be noticed that either the amounts are equal or when they differ, the difference is not significant. Again the numerical results indicate that in the volume also the responses are not influenced by dwell times for the young subject. In other words, for the chosen REV, the response of the skin of the young subject is not influenced by the dwell time in all its volume when subjected to contact-free creep tests.

Observations are different for the elderly subject (Fig. 12b). It appears that each sub-layer has a dwell time dependent evolution. When the dwell time is increased, the surfaces of the hysteresis cycles of the different sub-layers increase (Table 2). These surfaces give a way to quantify the amount of mechanical energy supplied by the load. For each sub-layer, the amount of the mechanical energy of the dwell time of 100 s is 20.50 ± 2.13 times higher than the one of the dwell time of 1 s. Recalling the main features of the physiology of the different sub-layers it appears that the comparison of the different curves shows a link between the dwell time of the creep test and the viscous behavior of the skin. When the dwell time is increased, the viscous contribution of a given sub-layer has time to develop and drastically change the response and the surfaces of its hysteresis cycles (Table 2). Again the surfaces give a way to quantify the amount of dissipated mechanical energy used by the sub-layer to return to its unknown initial state after the load has been removed. The top surface of the SC faces the external load first and it is provided with the highest amount of supplied external mechanical energy. It is also, the sub-layer which needs the highest amount of energy for coming back to its unknown initial state after the load is removed because it dissipates all the energy supplied. Contrary, in a Civil engineering sense, HY is the last sub-layers reached by the charge transfer process. This sub-layer receives the smallest amount of supplied external mechanical energy. When the load is removed, HY dissipates all this energy and comes back to its unknown initial state without needing extra energy for completing this evolution. In between these two responses, the other sub-layers receive from the charge transfer process different amounts of external mechanical energy. Each sub-layer uses entirely this amount for coming back when the load is removed. More, each sub-layer satisfies itself with the amount available and complete its evolution taking it back to its unknown initial state without needing extra-energy. Hence, for the chosen REV, every layer of the elderly subject exhibits its own evolutions that are strongly dependent on the dwell time. These evolutions are triggered by the viscosity coefficient.

This work is a first contribution for establishing the causal relationship existing between mechanical fields and physiology associating the age of a healthy subject, his physiology, his mechanical viscoelastic characteristics and his numerical stress fields in the volume obtained through numerical simulations of in vivo contact-free creep tests.

Conclusion

In this paper we presented 3-D finite element simulations of five numerical contact-free creep tests with dwell times of respectively 0.1 s, 0.5 s, 1 s, 2 s and 100 s . These numerical tests were based on experimental data recorded during contact-free creep tests performed on the volar forearm of a young adult subject and an elderly adult subject.

The numerical results showed different behaviors for the young subject and for the elderly subject.

For the young subject, the upper surface of the skin faces the external load with the same viscoelastic response whatever the dwell time of the numerical creep test. This evolution was also retrieved in the volume of the skin. The choice of dwell time appeared irrelevant and the response obtained was equivalent to that of an indentation test (dwell time of 0.1 s).

On the contrary, the skin of the elderly subject exhibited responses highly dependent on the amount of energy supplied during the tests. The sensorial feelings associated with the tests can be very different also. Therefore the optimal choice may be to work with a patient’s specific experimental protocol with a given dwell time and keep to it throughout the study.

Therefore, for a practitioner, the contact-free indenter is an easy way of characterizing the state of their patient’s skin before and after treatments.

In the context of medical survey, a practitioner needs an easy and robust way to characterize the state of his patient. He can perform an in vivo contact-free creep test with a short dwell time. A 1 s dwell time looks to be a good compromise between the needs of clinical tests and to be sure to provide sufficient amount of information to the coupled experimental-theoretical-numerical approach.

Conflict of interest

The authors declare that they have no conflicts of interest in relation to this article.

Ethical statement

Each volunteer has provided his informed consent before participating.

References

All Tables

All Figures

	Figure 1 (a) Principle of the Tonoderm® [20]. (b) Forearm set-up [31]. (c) Photo of the pressured area of the skin surface under the air flow [20]. (d) Inverted-bell shape of the air flow and typical contact profile which may develop.
In the text

	Figure 2 Experimental loading curve as a function of time (Exp).
In the text

	Figure 3 Loading curves of the five numerical contact-free creep tests with dwell times of 0.1 s, 0.5 s, 1 s, 2 s and 100 s.
In the text

	Figure 4 Rheological Zener model.
In the text

	Figure 5 User-designed REV deduced from the echograph images and 3-D mesh.
In the text

	Figure 6 3-D mesh with the boundary conditions.
In the text

	Figure 7 Profiles of the experimental deflection as a function of time for the young subject (blue line) and for the elderly subject (red line) for a flux of 20 NL/min.
In the text

	Figure 8 (Top) Comparison of the experimental and numerical curves of the displacement for three sets of parameters (Bottom) Curves of the relative error for the same sets of parameters.
In the text

	Figure 9 Sketch giving the position of the seven nodes of the central vertical under the air flux used in the analysis.
In the text

	Figure 10 Isovalues of the vertical displacements at the end of the dwell times: (a) 0.1 s, (b) 0.5 s and (c) 1 s for the young subject. (d) Contour lines of the mean stress for the dwell time of 0.5 s for the elderly subject.
In the text

	Figure 11 Numerical normal stress σzz as a function of the numerical deflection of the central point under the air flux for the young subject and for the elderly subject for the dwell times of 0.1 s, 0.5 s, 1 s, 2 s and 100 s.
In the text

	Figure 12 Numerical normal stress σzz as a function of the numerical deflection for the dwell times 1 s and 100 s for seven nodes of the central vertical under the air flux: (a) for the young subject, (b) for the elderly subject.
In the text