First-Order Shape Derivative of the Energy for Elastic Plates with Rigid Inclusions and Interfacial Cracks

Within the framework of Kirchhoff–Love plate theory, we analyze a variational model for elastic plates with rigid inclusions and interfacial cracks. The main feature of the model is a fully coupled nonpenetration condition that involves both the normal component of the longitudinal displacements and the normal derivative of the transverse deflection of the crack faces. Without making any artificial assumptions on the crack geometry and shape variation, we prove that the first-order shape derivative of the potential deformation energy is well defined and provide an explicit representation for it. The result is applied to derive the Griffith formula for the energy release rate associated with crack extension.


Introduction
Recent and ongoing advances in engineering and material science have increased the need for mathematical tools in order to design and optimize in an efficient way three-dimensional highly inhomogeneous thin structures.
In this paper, we confine our attention to the model of a composite Kirchhoff-Love plate proposed in [17]. More precisely, we consider an elastic plate containing a partially debonded rigid inclusion of regular shape. The displacement field of the B Viktor Shcherbakov shcherbakov@mathematik.uni-kassel.de Evgeny Rudoy rem@hydro.nsc.ru 1 Lavrentyev Institute of Hydrodynamics, Lavrentyev Ave. 15, Novosibirsk, Russia 630090 2 Institute of Mathematics, University of Kassel, Heinrich-Plett-Str. 40, 34132 Kassel, Germany inclusion is specified by an equality constraint and determined completely by six unknown constants. In order to eliminate nonphysical interpenetration of the crack faces, an inequality constraint that involves both the normal component of the longitudinal displacements and the normal derivative of the transverse deflection of the crack faces is imposed. This implies a full coupling between the longitudinal displacements and transverse deflection of the plate. Our interest is in investigating the differentiability of the potential deformation energy with respect to the shape variations of the reference configuration and in finding an explicit form of the corresponding first-order shape derivative. The motivation for such a study lies in the fact that this derivative can be used to calculate efficiently the energy release rates associated with variation of defects, which are utmost of importance to predict crack propagation [4,5,12,20].
The principal difficulty in the problem at hand comes from the equality and inequality constraints, which do not preserve their structure under domain variations. To overcome this difficulty, we develop an approach that is purely based on the minimizing properties of variational solutions. The most important idea behind our technique originates from [40], where the Griffith formula for the energy release rate associated with crack extension was deduced in a rigorous way for two-dimensional elastic bodies with curvilinear cracks subjected to the Signorini conditions. Following this idea, in the case of Signorini-type constraints imposed on crack faces, general results on the shape differentiability of the potential deformation energy supported by explicit formulae for the first-order shape derivative were derived in [41] for linear elastic materials, in [32] for a Mindlin-Timoshenko elastic plate model, and in [42] for a Kirchhoff-Love elastic plate model. Another approach to the shape differentiability of the potential deformation energy based on primal-dual Lagrange formulations of the nonlinear crack problems goes back to [25,26] and the subsequent [15]. Recent developments of the primal-dual shape sensitivity analysis in smooth domains and some applications to fluid mechanics problems with divergence-free constraints were given in [7,13,28,47]. We finally refer to [1,24] for a derivation of the Griffith formula for elastic materials with cracks by using the Piola transform. We stress here that all of the cited works related to solid mechanics, except [42], tackle inequality constraints on the crack faces that do not involve the gradient of displacement fields, while in the model under consideration the nonpenetration condition depends on the normal derivative of the transverse deflection. This additional difficulty essentially affects the computation of the shape derivative of the potential deformation energy.
In the context of Kirchhoff-Love plate theory, the nonpenetration condition for vertical cracks was put forward in [16]. A dimension reduction from 3D to 2D in the Kirchhoff-Love energy scaling regime for rate-independent fracture processes has been the subject of research in two past decades. To our best knowledge, all of the results are restricted to the use of the notion of global energetic solutions [34]; here we just mention a couple of them. With the help of evolutionary -convergence [35], a two-dimensional model for the propagation of a single crack in a plate from a threedimensional linearly elastic model involving a Barenblatt-like cohesive crack surface energy was obtained in [8]. A similar scenario for deriving two-dimensional models in the cases of Griffith-type delamination and an adhesive contact was performed in [9]. The advantages of the energetic formulation are that the initial problems and their two-dimensional counterparts have a derivative-free form and hence do not require any formulae for the energy release rate associated with crack extension; moreover, the crack path is not prescribed in advance.
In the presence of a preexisting interfacial crack between strongly dissimilar materials, it is natural to assume that the interface is the weakest fracture path and to apply the Griffith criterion [12] in order to describe propagation of the crack along the prescribed path. We recall that the Griffith criterion states: a crack is stationary if the energy release rate associated with crack extension is less than the fracture toughness of a material, and otherwise the crack grows. The main result of the present study gives a closed formula for the first-order shape derivative of the potential deformation energy (Theorem 2.2), and its application consists in deducing the Griffith formula for the energy release rate associated with extension of the crack along the interface, under minimal assumptions on the regularity of the crack path (Remark 2.11). Recent results in the same spirit of our work are [29,39,43,45].
Fracture behavior depends on many factors, among them the mutual dislocation and interaction of inclusions and cracks. In a series of papers [23,30,31,44], this issue was treated by an optimization technique based on regular shape variations. A more delicate shape-topological analysis of the nonlinear crack problems aimed at retarding or avoiding the crack propagation process was carried out in [22,27,33,49]. Inverse problems for the nondestructive determination of elastic and rigid inclusions embedded in Kirchhoff-Love bending elastic plates were investigated in [3,36,37].
The layout of the paper is as follows. In Sect. 1, we discuss the mechanical model and present the functional framework used to assert the existence of the variational solutions. After the necessary preliminaries regarding the velocity method, which we employ to describe shape variations of the reference configuration, we formulate and prove in Sect. 2 our main result. The paper closes with a short discussion of the case when the elastic properties of the plate are inhomogeneous and an application of the result to fracture mechanics: we derive the Griffith formula for the energy release rate associated with crack extension.

Formulation of the Problem
We consider a composite plate, which is a planar thin-walled structure, within the confines of linear Kirchhoff-Love theory. The 3D undamaged plate in the reference configuration is the right cylinder where ⊂ R 2 defines the midsurface of the plate and 2h denotes the thickness, which we assume to be constant. The midsurface is a bounded simply connected domain with boundary ∂ of class C 1,1 , and let ω 0 be a compact simply connected subdomain of (i.e., ω 0 ), with boundary ∂ω 0 of class C 1,1 . The boundary ∂ω 0 is the union of two disjoint simple curves γ 0 and ∂ω 0 \ γ 0 , where γ 0 is a relatively open set. The outward pointing unit normal to ∂ω 0 as well as to ∂ is denoted by ν 0 = (ν 01 , ν 02 ) T . Figure 1 provides a schematic visualization of the setup. We make also the following assumption on the geometry of the non-Lipschitz domain 0 = \ γ 0 : G1 There exist two domains 1 and 2 with Lipschitz boundaries such that 1 ∪ where H 1 stands for the one-dimensional Hausdorff measure.
The region ( \ ω 0 ) 3D is occupied by a linearly elastic homogeneous material whose properties are characterized by a fourth-order tensor C 3D that is supposed to be positive, that is, there exists c C 3D > 0 such that and monoclinic symmetric with respect to the (x 1 , x 2 )-plane, which implies We refer the reader to Remark 2.10 for some comments on the choice of the elasticity tensor. The region ω 3D 0 corresponds to a rigid inclusion, while the cylindrical surface γ 3D 0 is an interfacial crack. As usual in the theory of thin-walled structures, we deal with quantities that refer to the midsurface of the undeformed plate 0 . Let W = (w 1 , w 2 ) T be the longitudinal displacements of the midsurface points in the (x 1 , x 2 )-plane, and w is the vertical deflection of the midsurface points along the x 3 -axis. In what follows, we use a variational approach to equilibrium and thus look for (W , w) belonging to the Cartesian product of Sobolev spaces H 1 ( 0 ) 2 × H 2 ( 0 ).
Let C be the fourth-order tensor with components
Hence the resultant membrane stress tensor σ can be expressed in terms of the membrane strain tensor ε(W ) and the bending moment tensor m is related to the bending strain tensor −∇ 2 x w through the constitutive equations The restriction of a displacement field (W , w) on the domain ω 0 is an element of the are the spaces of infinitesimal rigid displacement fields on ω 0 that represent ker ε and ker ∇ 2 x , respectively: We observe that the following identities hold true [ The regularity assumptions on the domains ω 0 and \ ω 0 provide that the trace operators where the superscript k * is equal to 2 if k = 1 and 1 if k = 2, are well defined and continuous. Moreover, these operators have right continuous inverses (lifting operators) [14, Theorem 1.5.1.2], which we denote by Tr −1 k,ω 0 and Tr −1 k, \ω 0 , respectively. According to the positive and negative directions of the unit normal ν 0 on ∂ω 0 , there is a positive face ∂ω + 0 and a negative face ∂ω − 0 . The double brackets v = v + − v − stand for the jump of a function v across ∂ω 0 . For convenience, the same notation serves as well for the jump of functions across another smooth curves, it always clear from the context. In order to eliminate interpenetration of the crack faces, we impose the inequality restriction which is understood in the sense of traces. We emphasize that, from a mechanical point of view, the nonpenetration condition (4) describes the interaction of the opposite crack faces, admitting their frictionless contact (the equality case) or the absence of contact (the inequality case), the contact surface being unknown a priori. We finally suppose that the plate is clamped along its outer edge, which leads to the homogeneous boundary conditions The potential deformation energy of the system associated with a displacement field (W , w) is the stored elastic energy minus the work of external loadings The variational principle of minimum potential deformation energy forces that the composite plate 3D 0 clamped along the outer edge, with elastic part ( \ ω 0 ) 3D , rigid part ω 3D 0 , interfacial crack γ 3D 0 , and potential deformation energy given by (5), is in equilibrium if the displacement field of its midsurface (W 0 , w 0 ) is a solution of the minimization problem inf where the set of kinematically admissible displacement fields reads as follows Since E is convex and Gâteux-differentiable, the minimization problem (6) is equivalent to the variational inequality Thanks to G1, the Korn inequality and the Friedrichs inequality are valid in i , i = 1, 2, which implies the coercivity of E. Moreover, the energy functional E is weakly lower semicontinuous. The minimization problem (6) (and hence the variational inequality (7)) therefore admits a solution (W 0 , w 0 ) ∈ K 0 ( 0 ). It is straightforward to check that the solution to (7) is unique.
We are now ready to give the classical formulation of the equilibrium problem. We seek in the domain \ ω 0 satysfying the following equations and boundary conditions Here the traction σ + ν 0 , transverse force t + ν 0 , and bending moment m + ν 0 acting on ∂ω 0 are and the traction σ + ν 0 is decomposed into normal and tangential components given by A repeated Roman subscript signifies summation from 1 to 2. Equations (8) and (10) are equilibrium equations for the elastic part of the plate with respect to the (x 1 , x 2 )-plane and x 3 -axis, respectively. Conditions (13)- (15) are the complementarity conditions for contact of the crack faces. The nonlocal boundary conditions (18) and (19) ensure that the rigid part of the plate in equilibrium, that is, the resultant force and moment acting on it vanish. It is a routine matter to verify that our minimizer (W 0 , w 0 ) ∈ K 0 ( 0 ) is indeed a weak solution of the Euler-Lagrange system (8)- (19), see [17]. Remark 1. 1 We elucidate briefly in what sense the natural boundary conditions (14)- (19) for the minimizer (W 0 , w 0 ) are satisfied. For each k = 1, 2, we introduce the Lions-Magenes space is a Banach space. Since ∂ω 0 ∈ C 1,1 , the following equivalence holds [18, Proposition Denote by ·, · 00 k−1/2,γ 0 the duality pairing between the space H k−1/2 00 (γ 0 ) and its and the boundary conditions (14) and (15) read as follows [17] σ + τ 0 , v whereas the interpretation of the nonlocal conditions (18) and (19) is provided by Here ·, · k−1/2,∂ω 0 is the duality pairing between the spaces H k−1/2 (∂ω 0 ) and H −k+1/2 (∂ω 0 ). We exploit some of these conditions to derive an explicit formula for the first-order shape derivative of the potential deformation energy, see Theorem 2.2.

Shape Derivative of the Energy
We begin this section by recalling the velocity (speed) method [46, Section 2.9], which we employ to describe shape variations of the reference domain 0 . Let V ∈ W 2,∞ loc (R 2 ) 2 be a given shape velocity field. We assume further that V is compactly supported in the domain . This assumption does not limit the generality of our analysis, but it simplifies the presentation. The following proposition is a core of the velocity method and allows us to construct a family of diffeomorphisms of R 2 generated by V .
loc (R 2 ) 2 be a given shape velocity field satisfying the property suppV . Then there exists δ 0 > 0 such that: (a) the Cauchy problem Proof The proof is carried out in [21, Section 2] for W 1,∞ vector fields. Without any changes, the arguments are also applicable to W 2,∞ vector fields V .
Before proceeding further, we record here a consequence of [2, Theorem 2.2] that will be useful throughout the paper. For each k = 1, 2, the space W k,∞ ( ) is algebraically and topologically isomorphic to C k−1,1 ( ).
Let us fix δ ∈ (−δ 0 , δ 0 ) and set ω δ = T δ (ω 0 ) and γ δ = T δ (γ 0 ). Since T δ is an orientation-preserving C 1,1 -diffeomorphism between and T δ ( ), coinciding with the identity on ∂ , we have T δ ( ) = , see [6,. The domain ω δ is a simply connected one, and its boundary ∂ω δ = T δ (∂ω 0 ) is of class C 1,1 , with the outward pointing unit normal ν δ . Moreover, the boundary ∂ω δ is the union of two disjoint simple curves γ δ and ∂ω δ \ γ δ , and γ δ is a relatively open set. So we can put δ = \ γ δ . Against this background, we conclude that the mapping T δ induces a C 1,1 -diffeomorphism between the non-Lipschitz domains 0 and δ . To complete our considerations on the geometry of the perturbed domain δ , let us note that from the condition suppV follows the validity of G1 for δ . As above, the equilibrium of the composite plate 3D δ clamped along the outer edge, with elastic part ( \ ω δ ) 3D , rigid part ω 3D δ , interfacial crack γ 3D δ , and potential deformation energy is achieved when the displacement field of its midsurface (W δ , w δ ) is a solution of the minimization problem the set of kinematically admissible displacement fields being The minimization problem (20) possesses a unique solution satisfying the variational inequality By definition, the first-order shape derivative of the potential deformation energy E at 0 in the direction of V , if it exists, is given by the limit The questions we are now interested in are whether (22) does exist and how this shape derivative can be expressed in an explicit form.
To state our main result, let us fix some notation. For U , W ∈ H 1 ( 0 ) 2 , we define the bilinear form where the second-order tensor E reads as where We finally set Theorem 2.2 Under the assumptions of Proposition 2.1, the limit in (22) exists and is equal to (25).
The remainder of the paper is devoted to the quite long and technical proof of Theorem 2.2. The key idea behind the proof consists in applying the mapping T δ : 0 → δ to transform the potential deformation energy E( δ ; W δ , w δ ) to the reference domain 0 and calculating the desired limit as δ → 0 in 0 . The particular steps are presented below.
We begin with the study of properties of the mapping T δ : 0 → δ . In view of [38, Chapter 2, Lemma 3.4], it generates an isomorphism between the spaces H 1 Then the derivatives of first and second orders of such transformed functions can be calculated from the relations On the other hand, for (U δ , Let K δ ( 0 ) be the image of the set of the admissible displacement fields K δ ( δ ) under the mapping T δ . In order to identify K δ ( 0 ), we first note that any element (W , w) ∈ K δ ( 0 ) satisfies the inequality constraint The transformed outward pointing unit normal ν δ = ν δ • T δ is related to ν 0 through the equation [6, Section 1.7] where · 2 stands for the Euclidean norm in R 2 . Inequality (26) then reads Additionally, the restriction (W , w) on the domain ω 0 belongs to the Cartesian product Thus, K δ ( 0 ) may be defined as follows and so T δ is a bijection between the sets K δ ( δ ) and K δ ( 0 ). The change of variables y = T δ (x) in the variational inequality (7) implies that (W δ , w δ ) ∈ K δ ( 0 ) is a unique solution of the variational inequality Here, (F δ , f δ )(x) = (F, f )(T δ (x)).
As a preliminary step towards the verification of the uniform boundedness of (W δ , w δ ) in H 1 ∂ ( 0 ) 2 × H 2 ∂ ( 0 ), we establish the next proposition.

Proposition 2.3
The following Taylor formulae as δ → 0 are valid: Proof This is an immediate consequence of Proposition 2.1.

Proposition 2.5
with c > 0 independent of u and v.
We label the part ∂ ∪ ∂ω 0 \ γ 0 of the outer boundary of the domain 0 \ ω 0 as . For each k = 1, 2, we introduce the extension operator Ext k : It is clearly that Ext k are well defined, linear, and continuous operators. With these preliminaries behind us, we are able to state the key assertion of our analysis that guarantees the existence of certain test elements.

and uniform in δ bounds
are valid.
Proof For each i = 1, 2, we construct the required tuple (W i δ , w i δ ) as sums W i δ = P i δ + Q i δ and w i δ = p i δ + q i δ , where the set {P i δ , p i δ } provides the appropriate inequality constraint for ( W i δ , w i δ ) on γ 0 , and the set {Q i δ , q i δ } is responsible for the correct structure of ( W i δ , w i δ ) in ω 0 , which is specified by the equality constraint. The case δ = 0 is trivial, so we assume that δ = 0.
These considerations allow us to establish convergence properties of the transformed solutions (W δ , w δ ) as δ → 0.

Lemma 2.7
There exists a constant c > 0 that is independent of δ such that it holds for every δ ∈ [−δ 2 , δ 2 ] : Proof As in the proof of Lemma 2.4, we change the integration domain \ω 0 by 0 in the variational inequality (7). The resulting variational inequality is referred to as (7) 0 . Inserting ( W 1 δ , w 1 δ ) and ( W 2 δ , w 2 δ ) as test functions for (27) 0 and (7) 0 , respectively, and adding both variational inequalities, we invoke the asymptotic expansions (33)- (34) and (37)- (38) to obtain We finally apply the Korn inequality and the Friedrichs inequality to complete the proof.
The strong convergences ∇ Hence (51)-(52) imply the remaining convergences C δ → C 0 in R 2 and a δ → a 0 in R, so we are done.
We intend now to identify the limit of the sequence (W i δ , w i δ ) for each i = 1, 2.
Proof We first note that, by construction, ; furthermore, P 0 ≡ (0, 0) T and p 0 ≡ 0 in ω 0 . Since all the functions P i δ and p i δ vanish in ω 0 , it is enough to investigate their convergence in \ ω 0 .
Now we have enough tools at hand to prove Theorem 2.2.
Proof of Theorem 2.2 For the calculation of limit (22), we first transform the potential deformation energy E( δ ; ·, ·) to the reference configuration 0 , which leads to . Utilizing the same machinery as in the proof of Lemma 2.4 yields where Since T δ is a bijection between the sets K δ ( δ ) and K δ ( 0 ), it follows that Let δ > 0. Using ( W 1 δ , w 1 δ ) ∈ K δ ( 0 ) as a competitor in (55), we discover which implies (57) On the other hand, considering ( W 2 δ , w 2 δ ) ∈ K 0 ( 0 ) as a competitor in (6), we arrive at so that (59) Invoking the asymptotic expansion of the transformed potential deformation energy (54), we can calculate the limit superior on the right-hand side of (57) and the limit inferior on the left-hand side of (59). Direct computations based on Lemmas 2.7 and 2.9 reveal that they are finite and equal. Therefore, Taking into account that p + 0 = 0 and (∇ x p 0 ) + · ν 0 = 0 on ∂ω 0 \ γ 0 , we thereby obtain Furthermore, the equilibrium equation (10) and Remark 1.1 yield Combining (60)-(62) completes the proof.

Remark 2.10
The structural assumption (1) arises in dimension reductions from 3D to 2D within the Kirchhoff-Love energy scaling regime [8][9][10]. In the case of an undamaged elastic plate, this assumption provides in particular that the limit stored elastic energy is a sum of uncoupled bending and stretching energies. From a mathematical standpoint, we have no difficulty in treating a more general case in which the elastic properties of the plate depend on the spartial variable x as in [17]. Indeed, let C p , C v ∈ C 1 ( ) 2×2×2×2 be fourth-order tensors. For each n = p, v, the tensor C n is assumed to be positive, i.e., there exists c C n > 0 such that C n (x)X : X ≥ c C n |X | 2 for all X ∈ R 2×2 sym and all x ∈ , and symmetric C n i jkl = C n i jlk = C n kli j , i, j, k, l = 1, 2.
We denote by C n,V the fourth-order tensor with components C n,V i jkl = ∇ x C n i jkl · V , i, j, k, l = 1, 2, and replace the constitutive relations (2),(3) (see also (9),(11)) by Exploiting the asymptotic representations C n • T δ = C n + δC n,V + r n6 δ , r n6 we can easily adapt the arguments above to show that all our results remain valid.