Guaranteed upper and lower bounds on the uniform load of contact problems in elasticity

Two mathematical models are developed within the theoretical framework of large strain elasticity for the determination of upper and lower bounds on the total strain energy of a finitely deformed hyperelastic body in unilateral contact with a rigid surface or with an elastic substrate. The model problems take the form of two continuous optimization problems with inequality constraints, the solutions of which are used to provide an enclosure on the uniform external load acting on the body's surface away from the contact zone.


Introduction
In materials analysis and design, contact problems in elasticity are central to the modelling and investigation of many structural systems. In particular, the analysis of soft tissue biomechanics and the design of bioinspired synthetic structures involve non-linear hyperelastic models for which the mathematical and numerical treatment poses many physical, theoretical, and computational challenges [2,17,29,[31][32][33]39].
Hyperelastic materials are the class of material models described by a strain energy density function [14,15,35,40]. For these materials, boundary value problems are often equivalent to variational problems, which provide powerful methods for obtaining approximate solutions. They can also be used to generate finite element methods for which the numerical analysis stands on the shoulder of the mathematical analysis for the elastic model [6,23,30,34].
The complexity of contact problems modelling structural systems is generally associated with the detection of contacts and openings and the resolution of non-linear equations for contact. For example, in systems formed from linear elastic bodies (e.g. ceramics, metals), surface roughness can impede active contact when a surface is pressed against another, and contact forces cannot be transmitted where surface separation occurs [16,18,19]. By contrast, non-linear elastic bodies (e.g. rubber, soft tissue) are more pliable and thus capable of attaining more active support through which contact forces can be transmitted effectively [3,7,8,20]. These problems can be formulated and solved in the framework of variational inequalities [27], which originated in the paper by Fichera [11] on the existence and uniqueness of solution to the celebrated Signorini problem with ambiguous boundary conditions [38]. While Signorini problem consisted in finding the equilibrium state of an elastic body resting on a rigid frictionless surface, the influence of variational inequalities went beyond the fields of mathematical analysis, optimization, and mechanics [10,13,41], leading also to new models in game theory, economics, and finance.
In large strain elasticity, the study of finite dimensional structural models is motivated by the increasing need for effective computational techniques required in practical applications, but little attention has been paid to date to the underpinning infinite dimensional problems which are fundamental in the development of new methods. For example, while the reciprocal theorems of stationary potential and complementary energy are well known in the infinitesimal theory of elasticity and in the theory of structures [9,28], for large strain deformations, the stress-strain relation is non-linear, and its inverse is not uniquely defined in general, either locally or globally. For systems with a finite number of degrees of freedom, extensions of the reciprocal principles have been proposed in [21,26,36]. For an elastic continuum, assuming that the stress-strain relation is invertible, a complementary energy principle in terms of stresses was formulated via a Legendre transformation on the strain energy function in [24]. In order to avoid the difficulty of inverting the constitutive relation, in [22,37], trial functions for the deformation gradient instead of the stress tensor are used to obtain a variational principle of the complementary energy type. Then, under appropriate conditions, the stationary potential energy and complementary energy principles become extremum principles which can be used to provide lower and upper bounds on physical quantities of interest.
In the present study, the variational approach of [22] is extended first to the problem of a non-linear hyperelastic body in unilateral contact with a rigid surface, then to the case of two hyperelastic bodies in mutual non-penetrative contact. For each of these problems, the corresponding variational models take the form of two continuous optimization problems with inequality constraints, the solutions of which can be used to provide upper and lower bounds for the uniform external loading acting away from the contact zone. To illustrate the theory, analytical upper and lower bounds for the external load of a system formed from two elastic bodies in mutual unilateral contact and subject to large compression or bending, which can be maintained in every homogeneous isotropic incompressible hyperelastic material in the absence of body forces, are obtained.
In practice, contact problems with large stresses and strains at the adjoining material surfaces arise, for example, when creases are formed and self-contact takes place in soft materials and structures [42], or in biological systems where the attachment between cells are sufficiently weak so that cells separate, and failure through the appearance of gaps between cells occurs [25]. Although cell debonding is a spontaneous mechanism for crack initiation in many natural structures, it has been less investigated to date [4,5,12].
The remainder of this paper is organised as follows: in Section 2, the equilibrium problem is formulated for a non-linear elastic body in non-penetrative contact with itself or with a rigid obstacle, in the absence of friction forces; in Sections 3 and 4, variational problems are derived and analysed for the upper and lower bounds on the total strain energy of a finitely deformed elastic body made from a compressible or an incompressible material, respectively; in Section 5, the variational approach is extended to the case of two elastic bodies in mutual non-penetrative contact, and examples where the elastic bodies are made from a neo-Hookean material and subject to large compression or bending are presented.

Elastostatic Equilibrium with Unilateral Contact
A continuous material body occupies a compact domainΩ of the three-dimensional Euclidean space R 3 , such that the interior of the body is an open, bounded, connected set Ω ⊂ R 3 , and its boundary Γ = ∂Ω =Ω \ Ω is Lipschitz continuous (in particular, we assume that a unit normal vector n exists almost everywhere on Γ). The body is subject to a finite elastic deformation defined by the one-to-one, orientation preserving transformation: such that J = det (Gradχ) > 0 on Ω and χ is injective on Ω (see Figure 1). The injectivity condition on Ω guarantees that interpenetration of the matter is avoided. However, since selfcontact is permitted, this transformation does not need to be injective onΩ.
Let the spatial point x = χ(X) correspond to the place occupied by the particle X in the deformation χ. For the deformed body, the equilibrium state in the absence of a body load is described in terms of the Cauchy stress by the Eulerian field equation: div σ(x) = 0. The above governing equation is completed by a constitutive law for σ, depending on material properties, and supplemented by boundary conditions.
Since the domain occupied by the body after deformation is usually unknown, we rewrite the above equilibrium problem as an equivalent problem in the reference configuration where the independent variables are X ∈ Ω. The corresponding Lagrangian equation of non-linear elastostatics is: where P is the first Piola-Kirchhoff stress tensor. For a homogeneous compressible hyperelastic material described by the strain energy function W(F), the first Piola-Kirchhoff stress tensor is: where F = Grad χ is the deformation gradient.
Then the corresponding Cauchy stress tensor can be expressed as follows σ = J −1 PF T , and P = σcofF.  If the body is in non-penetrative contact with itself or with a rigid obstacle and subject to external dead loading conditions (see Figure 2), then the general boundary value problem is to find the displacement field u(X) = x − X, for all X ∈ Ω, such that the equilibrium equation • On Γ D , the Dirichlet (displacement) conditions : 3) • On Γ N , the Neumann (traction) conditions : P(X)N = g N (X), (2.4) where N is the outward unit normal vector to Γ N , and g N dA = τ da, where τ = σn is the surface traction measured per unit area of the deformed state.
• On Γ C , the conditions for unilateral (non-penetrative) contact either with itself or with a rigid obstacle, in the absence of friction forces: For the contact with a rigid obstacle, η : R 3 → R is the function describing the relative position of the body to the surface of the obstacle and N is the outward unit normal vector to this surface (oriented towards the obstacle), d ≥ 0 is the relative distance which cannot be exceeded between potential contact points, and g ≥ 0 is the cohesion parameter. For self- denotes the jump in the displacement at the points X = X ′ where self-contact may occur and N and N ′ represent the outward unit normal vectors to the boundary at those points, respectively. Then (2.5) sets the permitted relative distance between potential contact points; (2.6) gives the allowed normal force acting at a contact point; (2.7) is the complementarity condition that the maximum relative distance or the maximum normal force is attained at each point; and (2.8) states that, at the points of self-contact, the action and reaction principle holds.
For the study of the above contact boundary value problem (CBVP), the following result guarantees injectivity of the deformation mapping, so that self-penetration of the elastic body cannot occur [3].
satisfies the following condition: then χ is injective on Ω.

Unconstrained Materials
We now cast the CBVP (2.1)-(2.8) in variational forms similar to those described in [22]. However, in our case, we also take into account the kinematic constraints on the displacement field and the associated reaction force (2.5)-(2.6) at the potential contact zone Γ C . Throughout the analysis, we assume that the conditions (2.9)-(2.10) are satisfied. For simplicity of notation, we also set d = 0 and g = 0 in (2.5)-(2.7), which corresponds to the physical case of non-penetrative cohesionless contact. The case when these parameters are non-zero can be treated by analogy. For the given problem, the existence of a solution implies the mathematical correctness of the problem, while the (local) uniqueness of the solution indicates that the problem can be used to model of a physical phenomenon. We denote the closed convex set of kinematically admissible displacement fields by: for some s > 3/2 (see e.g. [23]). The closed convex set of statically admissible gradient fields is denoted by: for some q > 1 (see [23]).
Let u ′ , u ′′ ∈ K and F ′′ ∈ S, such that F ′′ = I + Grad u ′′ . By the divergence theorem, we obtain: where F ′ = I + Grad u ′ . From this identity we derive two important inequalities as follows. If, in (3.1), u ′ = u ′′ = u ∈ K and F ′′ = F ∈ S are the displacement and gradient fields, respectively, satisfying the CBVP, then: 1. First, in (3.1), we take u ′ ∈ K and set u ′′ = u ∈ K and F ′′ = F ∈ S to satisfy the CBVP. Subtracting (3.2) from (3.1) implies: where the second integral is non-negative.
The primal (displacement) variational problem is to find u ∈ K satisfying: for all admissible fields u ′ ∈ K.
2. Next, in (3.1), we take F ′′ ∈ S and set u ′ = u ∈ K and F ′ = F ∈ S to satisfy the CBVP. Subtracting (3.2) from (3.1) yields: where the second integral is again non-negative.
The complementary-type variational problem is to find F ∈ S such that: for all admissible fields F ′ ∈ S.
Assuming that the contact problem CBVP has a statically admissible solution F ∈ S and a corresponding kinematically admissible displacement u ∈ K, our objective is to determine upper and lower bounds on the total strain energy of the deformed body.

Upper and Lower Bounds on the Strain Energy
The potential energy of a kinematically admissible vector field u ′ ∈ K is defined by: Note that the contact constraint appears in the set of admissible fields K, but not in the expression of the energy functional E p . The principle of stationary potential energy states that the solution of the CBVP is characterised by the variational formulation: for all u ′ ∈ K.
Proof: Let u ′ = u + ǫ ∈ K, where u ∈ K is the solution to (3.3) and ǫ i ≪ 1, i = 1, 2, 3, such that: Then, the value of the potential energy at u ′ can be approximated to the second order in ǫ as follows: By (3.3): Hence: Assuming that the following condition holds for all non-zero ǫ satisfying (3 .7): we obtain (3.6), i.e. u ∈ K is a local minimum for the functional E p defined by (3.5). Note that condition (3.8) is also necessary for the potential energy to attain a local minimum at u ∈ K. ✷ Next, we define the complementary-type strain energy: Then the complementary energy of a stationary admissible field F ′ ∈ S is: In this case also, the contact constraints are present in the set of admissible fields S, but not in the expression of the energy functional E c . The principle of stationary complementary energy states that the solution of the CBVP is characterised by the variational formulation:
Proof: Let F ′ = F + Σ ∈ S, where F ∈ S is the solution to (3.4) and Σ ij ≪ 1, i, j = 1, 2, 3, satisfy the following conditions: We approximate the value of the complementary energy at F ′ to second order in Σ as follows: By (3.4), the above approximation implies: Assuming that the following condition is valid for all non-zero Σ satisfying (3.11)- (3.12): we obtain (3.10), i.e. F ∈ S is a local maximum for the functional E c defined by (3.9). The condition (3.13) is also necessary for the complementary energy to attain a local maximum at F ∈ S. ✷ Theorem 3.5 For unconstrained materials, the values of the potential energy functional E p defined by (3.5) and of the complementary-type energy functional E c defined by (3.9) represent an upper and a lower bound, respectively, for the minimum potential energy.
Proof: When u ∈ K and F ∈ S satisfy the CBVP, by (3.2): (3.14) Then, by (3.6), (3.10), and (3.14): for all u ′ ∈ K and F ′ ∈ S. ✷ The relation (3.15) guarantees that the values of the potential and complementary energy provide an upper and a lower bound, respectively, for the (local) minimum potential energy, under the assumption that such a minimum exists.

Materials with Internal Constraints
In this section, we extend the variational problems formulated in Section 3 to the case where the possible deformations of the body are restricted to those for which: where γ is a scalar-valued function representing the material internal constraint and F(X) is the deformation gradient at X. We denote the closed convex set of kinematically admissible displacement fields by: We denote the closed convex set of statically admissible gradient fields by: where λ ′ is the Lagrange multiplier associated with the material constraint (4.1). Let u ′ , u ′′ ∈ K γ and (F ′′ , λ ′′ ) ∈ S γ , such that F ′′ = I + Grad u ′′ . The divergence theorem implies: where F ′ = I + Grad u ′ . From this identity also, we derive two inequalities. If, in (4.2), u ′ = u ′′ = u ∈ K γ and (F ′′ , λ ′′ ) = (F, λ) ∈ S γ are the kinematically and statically admissible fields, respectively, satisfying the CBVP with the material constraint (4.1), then: 1. First, in (4.2), let u ′′ = u ∈ K γ and (F ′′ , λ ′′ ) = (F, λ) ∈ S γ satisfy the CBVP with the material constraint (4.1), and u ′ ∈ K γ . Subtracting (4.3) from (4.2) gives: In this case, the primal variational problem is to find u ∈ K γ that satisfies: for all kinematically admissible fields u ′ ∈ K γ .

Upper and Lower Bounds on the Strain Energy
In this case, the potential energy of a kinematically admissible field u ′ ∈ K γ is given by: (4.6)

Lemma 4.3
The energy functional E p defined by (4.6) has a local minimum at the solution u ∈ K γ of (4.4), i.e.
for all u ′ ∈ K γ .
Proof: Let u ′ = u + ǫ ∈ K γ , where u ∈ K γ is the solution to (4.4) and ǫ i ≪ 1, i = 1, 2, 3, such that: We approximate the value of the potential energy at u ′ to the second order in ǫ as follows: By (4.4): Hence: Assuming that the following condition holds for all non-zero ǫ satisfying (4.8)-(4.9): we obtain (4.7), i.e. u ∈ K γ is a local minimum for the functional E p defined by (4.6). ✷ We also define the complementary-type energy density as follows: The corresponding complementary energy of a statically admissible field (F ′ , λ ′ ) ∈ S γ is: Lemma 4.4 The energy functional E c defined by (4.12) has a local maximum at the solution (F, λ) ∈ S γ of (4.5), i.e.
Proof: Let (F ′ , λ ′ ) ∈ S γ , such that F ′ = F + Σ and λ ′ = λ + ǫ, where (F, λ) ∈ S γ is the solution to (4.5), and Σ ij ≪ 1, i, j = 1, 2, 3, and ǫ ≪ 1 satisfy the following conditions: ∂γ ∂F ij Σ ij = 0 in Ω, (4.14) We approximate the value of the complementary energy at F ′ to second order in Σ and ǫ as follows: By (4.5), the above approximation implies: Assuming that the following condition is valid for all non-zero Σ and ǫ satisfying (4.14)-(4.16): we obtain (4.13), i.e. (F, λ) ∈ S γ represents a local maximum for the functional E c defined by (4.12). ✷ Theorem 4.5 For materials with the internal constraint (4.1), the values of the potential energy E p defined by (4.6) and of the complementary type energy E c defined (4.12) provide an upper and a lower bound, respectively, for the minimum potential energy.

Extension to Two Elastic Bodies in Unilateral Contact
We now consider a system formed from two elastic bodies made from possibly different homogeneous hyperelastic materials described by the strain energy function W i , i = 1, 2, respectively, which are in mutual non-penetrative contact on part of their boundary, and denote by Ω i , i = 1, 2, the two open, bounded, connected, with Lipschitz continuous boundary, distinct domains occupied by the two bodies, respectively, such that Ω 1 ∩ Ω 2 = ∅ (see Figure 3).

S1
S2 external load obstacle Figure 3: Schematic representation of a system of two elastic bodies in mutual non-penetrative contact and subject to external loading.
Each body is subject to a finite elastic deformation: such that χ i satisfies (2.9)-(2.10) on Ω i .
The corresponding Lagrangian equation is: where F i = Grad χ i , i = 1, 2, p = 0 for compressible materials, and p is the hydrostatic pressure for incompressible materials. The formulation of the corresponding variational problems and their analysis are then parallel to that presented in Sections 3-4. We illustrate this with the following examples for large elastic deformations which can be maintained in every homogeneous isotropic incompressible material in the absence of body forces.
Example 5.1 Two thin elastic bodies which occupy the domains Ω 1 = (0, 1/2) × (0, d) × (0, d) and Ω 2 = (1/2, 1)×(0, d)×(0, d), respectively, where 1 ≪ d < +∞, are made from different neo-Hookean materials described by the strain energy function: where C i , i = 1, 2, are the corresponding material constants and I 1 is the first principal strain invariant. Each body occupying the domain Ω i , i = 1, 2, is subject to a triaxial stretch of the form:  For the deformation (5.2), the x, y and z-directions are principal directions, and the principal stretches are respectively: We first account for the unilateral contact constraints at the interface Γ C , where the relative normal displacement must satisfy the non-penetration condition (2.5): Since, by the action-reaction law the contact stresses at Γ C are equal in magnitude, if the hydrostatic pressure for the body occupying the domain Ω i is equal to −p i , i = 1, 2, then by the normal contact stresses condition (2.6): and by the complementarity condition (2.7): When the equality holds in (5.3), the conditions (4.10) and (4.17) are satisfied simultaneously if and only if: such that (5.4) is valid. Equivalently, the external load τ in the x-direction per unit area of the deformed configuration satisfies: and the action-reaction equality in (5.4) holds. If the strict inequality is satisfied in (5.3), then τ = 0.
Example 5.2 For the two elastic bodies described in the previous example and deformed by (5.2), we now consider the case where cohesive contact occurs at Γ C . Then the relation (5.4) is replaced by: where g > 0 is the given cohesion parameter, and the complementarity condition (5.5) becomes: In this case, if the equality is valid in (5.3), then by (5.6) and (5.7), the external load τ in the x-direction per unit area of the deformed configuration satisfies: , and the action-reaction equality in (5.7) holds. If the strict inequality is valid in (5.3), then τ = g.
Example 5. 3 We consider again two thin bodies made from different neo-Hookean materials as described by (5.1) and occupying the domains Ω 1 = (0, 1/2) × (0, d) × (0, d) and Ω 2 = (1/2, 1) × (0, d) × (0, d), respectively, where 1 ≪ d < +∞. Each body occupying the domain Ω i , i = 1, 2, is now subject to the a combined stretch and bending of the form: where A and a i , b i , i = 1, 2 are positive constants, (X, Y, Z) are the Cartesian coordinates for the reference configuration and (r, θ, z) are the cylindrical polar coordinates for the deformed configuration, such that r ∈ (r i−1 , r i ), i = 1, 2. During the deformation, non-penetrative contact between the two bodies at the interface Γ C = {1/2} × (0, d) × (0, d) is assumed (see Figure 5). This deformation is attained by uniformly loading the surface Γ N = {0} × (0, d) × (0, d) of the first body, while on the surface Γ D = {1} × (0, d) × (0, d) of the second body the deformation is prescribed. In this case also, since the dimension of the bodies in the X-direction is much smaller than in the other two directions, changes in the area of the contact surface during deformation may be neglected. For the deformation (5.9), the r, θ and z-directions are principal directions, and the principal stretches are: At the contact interface Γ C , the relative normal displacement satisfies the non-penetration condition: such that (5.11) is valid. Equivalently, the external load τ in the r-direction per unit area of the deformed configuration satisfies: − min C 1 max{a 1 /r 0 , Ar 1 / √ a 1 , 1/(A √ a 1 )} − C 1 a 2 1 r 2 1 , C 2 max{a 2 /r 1 , Ar 2 / √ a 2 , 1/(A √ a 2 )} − C 2 a 2 2 r 2 1 < τ < 0, and the action-reaction equality in (5.11) holds. If the strict inequality is satisfied in (5.10), then τ = 0.

Conclusion
In large deformations, elastic bodies are very likely to enter into contact with neighbouring obstacles or with themselves, and therefore both the mathematical models and their analysis must take these phenomena into account. In this paper, for an elastic body made from an homogeneous isotropic non-linear hyperleastic material in unilateral contact with a rigid obstacle or with another elastic body, trial functions for the deformation gradient were used to obtain a variational principle of stationary complementary energy type. Then the stationary potential energy and complementary energy principles were used to provide an enclosure on the total strain energy of the finitely deformed body. To illustrate the theory, this variational framework was applied to determine upper and lower bounds on the external load for a system made from two elastic bodies of neo-Hookean material under uniform compression or under combined stretch and bending. This investigation addresses the need for a better understanding of contact problems in natural and industrial systems where mechanical models that take into account the large stresses and strains at adjoining material surfaces are crucial.