Indentation of a circular membrane on an incompressible liquid

A model of the indentation of a circular elastic membrane (neo-Hookean material) on an incompressible liquid layer and an effective numerical solution method are presented. For the interface between the indenter and the membrane, stick or slip contact conditions are considered. The solution procedure identifies the liquid pressure and, for slip contact, the radius in the reference configuration displaced to the edge of the indenter. Assuming a range of material and test performance parameters regarding deep tonometry of soft subcutaneous tissue, the predictions for the distributions of stretches, line tensions and profiles of the indented membrane are analyzed. Among the key results of this work are the dependencies of the total indentation force and the liquid pressure on the Young’s modulus of the membrane’s material that appeared to be approximately linear for a range of Young’s modulus values. The profiles of the membrane for different indentation depths are close to parabolic and can be used to relate the longitudinal line tension at the edge of the indenter to the total indentation force.


Introduction
The indentation method is one of the well-developed testing techniques of materials. The usefulness of the indentation method is especially appreciated in situations where its noninvasive nature is important, for example in studies of biomaterials and biological materials and in clinical trials. Among these applications, the indentation method is frequently used to measure properties of soft materials, such as cells [18,24], gels [7], cartilage [22], skin and subcutaneous tissue [12,16], ocular tissues [25], and tympanic membranes [1]. The identified properties may refer to cells, tissues, or organs and express elastic or viscoelastic parameters as well as permeability. The indentation depth can be of the order of nanometers [18] or micrometers [16,25], and millimeters [7,22] up to one centimeter [26]. In some of the solutions, forces are controlled, and in others, penetration and speed are controlled. The sizes and shapes of indenters are also varied.
In the case of medical applications involving the indentation of soft tissues, the material frequently cannot be treated as a homogeneous half-space or a layer; nevertheless, its true structure should somehow be incorporated. The interpretation of the test results for such structures obtained via the indentation method requires the use of additional data delivered by other techniques that allow the separation of the impacts of individual layers. Such complementary techniques are the ultrasonic method and the magnetic resonance technique [8,23], which are also noninvasive and can provide information on the thickness of individual layers or their shear stiffness. However, these methods require complex and expensive measuring equipment, thereby limiting their applications. Another approach to study inhomogeneous materials is the use of different sizes of indenters and indentation depths as well as different relevant tissue models [2,5,16] to assess the structures and properties of the tested materials.
Among the inhomogeneous soft biomaterials or tissues are materials in which a very soft medium is covered by a layer with a much higher shear stiffness. For example, some cells [4,18,24], microcapsules [3], and soft tissues for which the subcutaneous tissue has a relatively low shear stiffness in relation to the skin [10] qualify as such materials. The last example refers to cases such as early lymphedema or lipedema. From the perspective of indentation tests, such materials can be modeled as a liquid compartment enclosed by a membrane [17][18][19]21].
A clinical test, to which the considered in this paper model may particularly refer to, is the deep tonometry applied to a very soft edematous tissue of a limb. Then, a relatively thin elastic skin (resembling membrane) covers a layer of soft subcutaneous tissue and both wrap around the core composed of muscles and bones. The stiffness of the muscles surrounded by fascia and of the bones is generally much higher than the stiffness of subcutaneous tissue. As a result for physiologically acceptable loads applied in tonometry deformations of muscles and bones caused by an indenter compared to deformations of skin and subcutaneous tissue are small and in the first approximation can be neglected. In the simplest version of the deep tonometry (see [26]), a dynamometer equipped with a flat cylindrical tip, with a surface area of 1 cm 2 , is pressed against the tissue at a constant rate of approximately 1 mm/s to the depth of 1 cm. The maximum force is recorded and is used as a descriptor of a state of the tissue. Since the period when the maximum force is registered is short undrained conditions and incompressibility of subcutaneous tissue can be assumed. The finite element simulations based on the two-phase hyperporoelastic model for subcutaneous tissue and the hyperelastic model of skin proved that for low stiffness of subcutaneous tissue (early lymphedema or lipedema) and the typical stiffness of skin the role of the latter one is significant (see [15]). However, a number of questions, important from the practical point of view of deep tonometry, are not yet answered, for example: what is the dependence of the force on the indenter, which is the contribution of skin on the skin's stiffness, how important is the interaction of the deformed skin with subcutaneous tissue and indirect influence of the interaction on the response of tissue on the indenter, and what is the role of friction between the indenter and the skin for acquired data. Taking into account the pinpointed features of the deep tonometry and the clinical observation that a few centimeters from the indenter no visible tissue deformations can be seen, answers to some of these questions can be found based on studies of the indentation of a membrane on incompressible liquid layer.
In this paper, a model, a method to obtain solutions to the model, and selected numerical results are presented for an indentation test of a system composed of a circular elastic membrane on a liquid, the volume of which does not change. Because the contact between the indenter and membrane is poorly controlled, two extreme cases are considered: sticking friction (no-slip) or frictionless (slip) interfaces. The membrane is modeled as a neo-Hookean material, and the liquid is assumed to be incompressible. Taking into account the fact that the pressure in liquid is unknown and changes with indentation depth, an iterative numerical technique is proposed to determine the pressure from the condition of a constant liquid volume. In the case of a slip interface, the membrane radius in the reference configuration, which corresponds to the material particles displaced to the edge of the indenter, must also be found. The results of the simulations are illustrated by showing the distributions of the stretches and the line tensions for each selected indentation depth. The dependence of the indentation force and the pressure in the liquid induced by the membrane deformation on the elasticity of the membrane material and the indentation depth is studied, and the role of the boundary conditions is determined. The interaction of the membrane and the liquid is also analyzed indirectly by determining the ratio of the net force of the pressure and the total indentation force.

Governing equations
We consider a circular membrane rigidly clamped at its rim on a layer of incompressible liquid indented at its center by a cylindrical punch of radius a. Because the volume of liquid is assumed as constant, the initial shape of the liquid reservoir does not play any role in this case. The contact between the membrane and the indenter is modeled by sticking friction or frictionless (slip) interfaces, and the indenter remains rigid. Figure 1 shows a half-cross section of such a system. The undeformed membrane occupies the interior of a circle lying on the plane z = 0 with radius R e . The membrane is clamped along the edge at r = R e . The membrane's deformation is constrained by the underlying volume of liquid, which must be constant (liquid is assumed to be incompressible). As a result, the liquid interacts with the membrane via pressure p = P 0 − ρgz, where P 0 is an unknown constant, ρ is the liquid density, and g denotes the gravitational acceleration.
To model the deformed material points on the membrane, a cylindrical coordinate system (r, ϕ, z) is adopted. Because the deformation is axisymmetric, attention is restricted to any cross section of the membrane, e.g., The set of equations governing the large deformation of the weightless membrane under pressure p are [11,13] where λ s and λ ϕ are the longitudinal (along the cross section curve in the r-z plane) and latitudinal (along the direction normal to the r-z plane) principal stretches, respectively: and T ϕ and T s are the longitudinal and latitudinal line tensions in the deformed configuration, respectively. For the incompressible neo-Hookean material (the simplest hyperelastic model) of the membrane, the tensions are where E is the small strain Young's modulus, and h 0 denotes the thickness of the undeformed membrane.

Boundary conditions
A part of the membrane that, in the deformed configuration, is located between the edge of the indenter (r (R 0 ) = a) and the clamping ring (r = R e ) is considered. The radial coordinate of the left boundary of the considered part of the membrane in the reference configuration is R = R 0 . Assuming a given indentation depth z 0 , we have the condition The other boundary condition for the same radial coordinate is derived assuming stick or frictionless (slip) contact of the indenter with the membrane. In the stick contact case, because the indenter is assumed to be rigid, we have a lack of deformation of the portion of the membrane for which R ≤ a. This results in the condition Alternatively, for the frictionless contact between the indenter and the membrane, we incorporate the equilibrium condition for longitudinal line tensions, the continuity of latitudinal stretches, and the homogeneity of the stretches in the plane part of the membrane under the indenter. Thus, the boundary conditions for R = R 0 become (see the "Appendix") where the stretches λ ϕ , λ s are defined for R = R + 0 . For the clamped edge of the membrane (R = R e ), the boundary conditions are

Method of solution
The numerical solutions of the boundary value problem (1), (4), (5 or 6), (7) and (8) with unknown P 0 for the stick contact and P 0 and R 0 for the frictionless contact are obtained using the iterative process shown schematically in Fig. 2. The code was written in the MATLAB environment (version 7.11). The calculations start with initial guesses for P 0 and R 0 . The integration of (1) using the bvp5c function (MATLAB) requires also initial guesses for dependent variables λ s , α, λ ϕ and z as functions of R. The initial values of stretches are assumed constant and slightly above 1. (Using values equal to 1 has generated convergence problems.) For the coordinate z(R) parabolic approximation is adopted, i.e., where c 1 , c 2 and c 3 are constants that are determined from Eqs. (4) and (7) and the requirement that the volume of liquid under membrane is constant: The initial guess for the angle α(R) is calculated from the derivative of z p defined by (9). The integration of (1) is followed by checking condition (10) with the assumed tolerance limit ε v . If the condition is not satisfied, then the value of P 0 is updated following the built-in algorithm of the function solve in MATLAB. When the volume change is smaller than ε v , the loop ends. Next, the actual position of r (R 0 ) is calculated from λ ϕ (R 0 ) and (2) 2 . The value of |r (R 0 ) − a| is compared with the assumed tolerance limit ε r . Updating R 0 is realized by the bisection method, assuming as the lower and the upper limits for searching R 0 the values a/2 and a, respectively, which guarantees a sufficiently fast and convergent algorithm. The results are visualized when the external loop ends.

Numerical results
The proposed model may roughly correspond to the indentation test of a very soft subcutaneous tissue covered by skin. Therefore, the results of the simulations are illustrated assuming a range of material parameters (thickness and stiffness of skin) and performance parameters (indentation depth and indenter radius) appropriate for deep tonometry [26]. In all the cases considered, the following parameters are used: the radius of the indenter a = 5.5 mm, the thickness of the undeformed membrane h 0 = 1 mm, and the membrane's outer radius R e = 0.05 m. Selection of this radius is motivated by clinical observations, which show that a few centimeters from the indenter no visible tissue deformations can be seen. The value of the indentation depth has an important influence on the stretches. The difference due to the type of boundary condition increases closer to the indenter. For both boundary conditions, the longitudinal stretches λ s (R) are greater than one. For the latitudinal stretches, λ ϕ (R) are always less than one for the stick contact and cross the line λ = 1 for slip contact.
The values of line tensions T ϕ and T s change more rapidly along the radial coordinate than the stretches and are also strongly sensitive to the indentation depth. The maximum values for each case are at the edge of the indenter. The difference between solutions for the stick and slip conditions increases closer to the indenter, and both T ϕ and T s may be as high as 20%. The latitudinal tension has a minimum between the edge of the indenter and the membrane's rim. The presence of the minimum is the result of particular distribution of stretches, and for more compliant material of the membrane the model may predict compressive latitudinal stress, which induces wrinkling of the membrane (see, e.g., [20]).
The examples of the deformed membrane profiles are plotted in Fig. 5. The profiles cross the line corresponding to the undeformed membrane almost at the same point, irrespective of the indentation depth (parabolic approximations meet at the same point for R ≈ 0.025 m). The role of the boundary conditions for the membrane's shape is less visible than for tensions or stretches.
The profiles of the indented membrane are close to parabolic functions of radius R; the differences between the model prediction and the parabolic shape defined by (9) normalized by the indentation depth z 0 , i.e., (z − z p )/z 0 , for three indentation depths of z 0 = − 2, − 6 and − 10 mm, are shown in Fig. 6. The differences depend on the type of contact between the indenter and the membrane as well as the radius R. For larger z 0 , the maximum values of the deviations are less than 2% for the whole range of R. The largest normalized deviations appear for the smallest z 0 and at certain distances from the indenter, potentially amounting to 5%.
For a given indentation depth z 0 , the total force F acting on the indenter is composed of the two contributions, namely the net longitudinal force in the membrane F m for r (R 0 ) = a and the net force of pressure F p acting on the surface r ≤ a:  Thus, taking into account Eq. (3) 1 and the dependence of pressure in a liquid p = P 0 − ρgz, (11) can be rewritten as where λ s0 = λ s (R 0 ), λ ϕ0 = λ ϕ (R 0 ), and α 0 = α(R 0 ). Taking into account the usefulness of using the indentation test as the method for the identification of the material properties of the membrane, the dependence of the total indentation force F and pressure P 0 on the elastic constant E in the range roughly corresponding to the stiffness of skin tissue, i.e., between 0.01 and 10 MPa [9], and the indentation depth z 0 applied for deep tonometry [26] are presented in Figs. 7 and 8. Note that the total indentation force F is almost linearly dependent on Young's modulus E for the majority of the considered range. Only for values of E less than 0.1 MPa do the deviation from linearity become visible, particularly for lower indentation depths. The differences between the solutions for the two considered boundary conditions for total force and pressure in the log-log scale are almost constant. The distributions for liquid pressure resemble those for the total force.
The differences between the results for the total force and the liquid pressure for stick contact and for slip contact of the membrane with the indenter increase with both the indentation depth and elastic constant, as shown in a linear plot in Fig. 8.
Because both the stress in the membrane and the liquid pressure contribute to the total indentation force, the ratio between the net force of pressure F p and the total force F for the studied ranges of elastic modulus and indentation depth is determined, see Fig. 9. The ratio assumes larger values for small indentation depths and soft materials and, for the studied range of parameters, does not exceed 35%. For a higher stiffness, the ratio tends toward a limit at the level of 3% independent of the indentation depth.

Summary and discussion
The large deformation model of the indentation of a circular elastic membrane on an incompressible liquid layer was developed. The membrane's material was modeled as a neo-Hookean continuum. An effective numerical solution methodology was proposed that provides the distributions of the stretches, line tensions, and profiles of the indented membrane for two extreme conditions regarding the interface between the indenter and membrane: stick and slip contacts. For each indentation depth, the solution procedure identifies the parameter that determines the liquid pressure and the radius in the reference configuration corresponding to the particles being displaced to the edge of the indenter (only for the case of slip contact).
The results were discussed assuming the range of material and test performance parameters appropriate for deep tonometry applied to very soft subcutaneous tissue (lymphedema or lipedema) and skin. The difference between the numerical predictions for the stick and slip contacts of the indenter and membrane is most pronounced for line tensions close to the indenter and may be as high as 20%. The difference is smaller but also noticeable for the indentation force, especially for high values of Young's modulus and large indentation depths. The assumed size of the indenter is relatively large, as it is used for deep tonometry [26]; in other cases, when the indenter is smaller, the difference between solutions should decrease.
The key results of this work concern the dependencies of the total indentation force and the liquid pressure as well as the contribution of the liquid pressure to the total indentation force on the Young's modulus of the membrane's material. The first two quantities were found to be approximately linear functions of Young's modulus when the modulus spans the range between 0.1 and 10 MPa, i.e., if the considered model is valid and membrane's thickness is known, then the indentation force is directly related to Young's modulus. The contribution of liquid pressure to the total indentation force (expressed by the ratio F p /F) is significant only for small indentation depths and low Young's modulus values. For a large indentation depth (on the order of 10 mm), the contribution is less than a few percent. This knowledge should be important for diagnostic applications of the considered system because liquid pressure expresses the interaction between the membrane and a soft material beneath the membrane. If the interaction between the two compartments can be disregarded, the total indentation force can be approximated by the sum of the resistance of membrane and resistance of the material under the membrane taken separately. In this case, the presented model can also be used to approximate the contribution of skin in the indentation test if the skin's properties (thickness and stiffness) are known from other methods (e.g., ultrasound and extensometer). Taking into account the fact that tonometry tests can be performed not necessarily vertically the effect of water gravity was assessed comparing model's predictions for the membrane's shape and stretches with and without the contribution of gravity to liquid pressure. In general the differences have appeared insignificant and only for latitudinal stretches they can be noticed.
The simulations showed that the profiles of the membrane for different indentation depths are close to parabolic in shape, and the coefficients of the parabolic function can be found from boundary conditions (4) and (7) and the constraint for the liquid volume under the membrane (10). The approximation for the angle   α(R 0 , z 0 ) can be then determined from Eq. (9). Thus, assuming that in Eq. (11), the net force of pressure F p is negligible, the longitudinal line tension at the edge of the indenter T s (R 0 , z 0 ) is directly related to the total indentation force F(z 0 ), i.e., T s (R 0 , z 0 ) ∼ = F(z 0 )/2πa sin α(R 0 , z 0 ). The analysis presented in this work assumes that the reference configuration is stress free. When in a flat state, the membrane is in equi-biaxial tension and the solution must be modified; the methodology that uses the transformation of stretches to a stress-free state can be adopted [14]. The model and solution method supplemented with an inversion algorithm can be useful for the identification of the stiffness of the membrane's material, including for some soft polymeric materials or systems such as the solid-on-liquid (SOL) nanostructures used in opto-electronic applications [6].