Existence and asymptotic results for an intrinsic model of small-strain incompatible elasticity

A general model of incompatible small-strain elasticity is presented and analyzed, based on the linearized strain and its associated incompatibility tensor field. Strain incompatibility accounts for the presence of dislocations, whose motion is ultimately responsible for the plastic behaviour of solids. The specific functional setting is built up, on which existence results are proved. Our solution strategy is essentially based on the projection of the governing equations on appropriate subspaces in the spirit of the Leray decomposition of solenoidal square-integrable velocity fields in hydrodynamics. It is also strongly related with the Beltrami decomposition of symmetric tensor fields in the wake of previous works by the authors. Moreover a novel model parameter is introduced, the incompatibility modulus, that measures the resistance of the elastic material to incompatible deformations. An important result of our study is that classical linearized elasticity is recovered as the limit case when the incompatibility modulus goes to infinity. Several examples are provided to illustrate this property and the physical meaning of the incompatibility modulus in connection with the dissipative nature of the processes under consideration.


Introduction.
1.1. The intrinsic approach to elasticity. An intrinsic approach to elasticity simply means that the main and primal variable is the strain, together with its derivatives, and that the displacement and rotation fields are possibly recovered in a second step, in case they are needed. This approach is most probably the first historically, being the strain indeed considered to measure deformation, that is, variation in length and in mutual orientation of infinitesimal fibers within a solid body. As a matter of fact, for the geometer the strain is a metric from which all other geometric concepts (such as curvature, torsion and other sophisticated tensors, see [4,11]) are retrieved. Specifically, given a smooth strain tensor field ε, the classical Kirchhoff-Saint Venant construction (for a historical review, see [4,28,36]) in linearized elasticity basically consists in 3770 SAMUEL AMSTUTZ AND NICOLAS VAN GOETHEM • introducing the Frank tensor 1 F = Curl t ε, where Curl t ε stands for the transpose of the curl of the symmetric tensor ε (computed row-wise); • defining the rotation field as ω(x) = ω(x 0 ) + x x0 F(ξ)dl(ξ), on a smooth curve joining the endpoints x 0 and x; • defining the displacement field as u(x) = u(x 0 ) + x x0 (ε − (ω)) (ξ)dl(ξ), where (ω) stands for the skew-symmetric rotation tensor constructed from ω, namely (ω) il := ilk ω k , with ilk standing for the Levi-Civita symbol. Obviously such definitions are a priori path-dependent. In order for u and ω to be well-defined, i.e., to be path-independent, it is immediately seen that a sufficient and necessary condition in a simply connected domain be that inc ε := Curl Curl t ε = Curl F = 0, where inc ε denotes the strain incompatibility tensor defined index-wise by ( inc ε) ij = ikm jln ∂ k ∂ l ε mn , that is easily seen to be symmetric. When inc ε = 0 one retrieves the well-known expressions ∇ S u = ε and ∇u = ε − (ω), where ∇ S u = (∇u+∇ t u)/2 stands for the symmetrized gradient of the displacement field u. On the contrary, if the tensor inc ε is not vanishing, this means exactly that the rotation and/or displacement fields exhibit a jump around what is classically called a Burgers circuit (to this respect an important role is played by the choice of the origin x 0 as shown in [42]). So far, it appears clear that the important geometric quantities are ε, Curl t ε, and inc ε. This is the seminal motivation for our model which is precisely designed from these variables. In particular, this choice makes the model of gradient-type. Note that it appears natural to consider the curl instead of the full gradient, and the inc instead of the full Hessian.
Let us also stress that this approach, despite rarely seen today, has a long history: we date the origin of the intrinsic view to Riemann with ground-breaking applications in general relativity and later in mechanics (in particular see the Hodge and Prager approach in perfect plasticity [34]). In general, as explained in [4], Riemann's view is in contrast with Gauss' standpoint of immersions, that is a displacement or velocity-based formulation. Furthermore, as far as dislocations are involved, this geometric approach was very much developed and enhanced by the physicist E. Kröner in the second half of last century [27]. It should however be mentioned that in finite as well as in linearized elasticity the intrinsic approach was recently considered and developed during the last decades in a systematic way by Ciarlet and co-authors [2,3,[12][13][14][15], and Geymonat and co-authors [20][21][22] (see also [47] for a geometric approach). In particular, their aim was to write, for the elasticity system, the homogeneous boundary condition on the displacement in terms of the elastic strain only.
Although the incompatibility operator has been used in the engineering literature for a long time, the mathematical study of spaces of square integrable tensor-valued functions with square integrable incompatibility was not yet considered and thus our first step was to dedicate a paper to the subject [5].

1.2.
A model of incompatible elasticity. The approach we propose was introduced in [6] from a physical standpoint. It aims at accounting for the macroscopic effect of the motion of dislocations, since it is known from the works of E. Kröner that elastic strain incompatibility is related to the dislocation density tensor [27]. The proposed model is expected to ultimately provide an original framework to the modeling of elasto-plastic behaviors. At the current stage of development, it can be termed generalized elasticity or incompatible elasticity, since one important feature is that classical linearized elasticity is recovered as a limit case. The main point of our approach is that the strain is a symmetric tensor, but not necessarily a symmetric gradient. Moreover, our theory relies on the following rationales.
1. Strain rate is preferred to strain and is given the following, primordial definition. The medium is considered as a collection of infinitesimal cells that individually deform smoothly, so that within each cell one can identify and follow fibers. Denote by a 1 , a 2 , a 3 three such fibers, which at time t originate from point x and are oriented along the axes of a Cartesian coordinate system and scaled to be of unit lengths. Then the strain rate is defined at x as (see, e.g., [18]) (1.1) Having fixed an initial time t 0 = 0, the time integral of the objective tensor d, called the strain or deformation tensor, reads ε(t) = t 0 d(s)ds. Note that this latter expression is only valid in the small strain setting, while (1.1) is general. 2. The strain rate defined in this way admits a mathematical decomposition: it is the sum of a compatible part and an incompatible part, that is given by a structure theorem called Beltrami decomposition [28]: This decomposition is unique once boundary conditions for v, which can be seen as a velocity field, are prescribed. Therefore, while d is an objective field by definition, neither ∇ S v nor E 0 are objective, see the discussion in [6]. For this reason the model will be constructed upon the full strain rate d (or ε if small strain is considered) and its space derivatives. 3. The governing equations should generalize those of classical linear elasticity in the sense that they must take into account the possible strain incompatibility. The idea behind this is to represent the macroscopic effect of dislocations in the micro-structure, as strain incompatibility is related to the density of dislocations [27,39,[41][42][43]. Our model can be briefly described as follows in the simplified case of a homogeneous material (see [6] for details). One considers linearized gradient elasticity in the sense of Mindlin [31]. One assumes that the virtual strain rated and its gradient produce intrinsic work, and by the virtual power principle we write where σ, τ are the Cauchy stress and hyperstress tensors, respectively, and K is a tensor representing external efforts. Assuming first a natural initial configuration (at t = 0), constitutive relations are taken as σ = Cε and τ = D∇ε, where C, D are the isotropic Lamé and Mindlin tensors, respectively (see [31]). Next, we require that the intrinsic power induced by the hyperstress Ω τ · ∇d dx vanishes as soon as the deformation is compatible, i.e., that it is only due to micro-structural defects in the form of dislocations. Then it was shown in [6] (see also [4] for different arguments) that the components of D are related through a scalar , called incompatibility modulus, which eventually yields that − div τ = inc ε. Therefore the virtual power principle leads to the weak form where E is the set of admissible virtual strain rates. To see that this equation generalizes linearized elasticity, taked = ∇ Sv witĥ v = 0 on Γ D ⊂ ∂Ω and take K such that − div K = f in Ω and KN = g on Γ N := ∂Ω \ Γ D . Then, plugging this into (1.3) immediately yields which is exactly the system of linearized elasticity in case of compatible strain, i.e. with ε = ∇ S u, u = u 0 on Γ D , since for such strains inc ε = 0. More generally, we believe that our model of incompatible linearized elasticity is able to represent inhomogeneous material properties and finite deformations through an incremental formulation. Indeed, nonlinear problems in continuum mechanics are classically solved through the finite element method used in conjunction with an incremental solution procedure. In this way, nonlinear problems are reduced to a sequence of iterations consisting of linearized problems. Therefore, in our model we will rather consider the strain increment (later denoted by E) in place of ε, together with the generalized tangent parameters (C, ). We emphasize that the procedure we suggest is Eulerian by essence, with all coordinates related to the deformed configuration, and is not to be confounded with Lagrangian incremental methods (as described e.g. in [10]). In case coordinates in a fixed reference configuration are needed for practical purposes, standard transformation rules might be applied. Of course, the evolution of the tangent moduli between increments should be driven by constitutive laws in order to account, for instance, for hardening phenomena. A possible thermodynamic approach is to relate changes in these coefficients with dissipation. In order to reach this aim, as a first step, a sensitivity analysis of the dissipation functional with respect to a variation of within a small inclusion was conducted in [6] for a simplified model. The extension to the full model and the numerical implementation, for which the existence results of the present work constitute a mandatory preliminary step, is an ongoing work. Our approach and model have been put in a historical perspective of intrinsic views in geometry in the survey paper [4]. For philosophical thoughts about modeling in physics, we refer to [29].
1.3. Summary of our results. Set E = L 2 (Ω, S 3 ). Subsequently the model variable will be denoted by E, and can represent either a strain, a strain increment or a strain rate by time derivation. The main purpose of this work is to prove that (1.3), or equivalently the associated strong form CE + inc E = K in Ω, has a unique solution in the space of square integrable functions with square integrable incompatibility, with the additional condition on the dislocation flux at the boundary inc E N = h on ∂Ω. The main ingredients to achieve the proof are (i) an orthogonal decomposition of L 2 (Ω, S 3 ) related to the Beltrami decomposition, and (ii) Fredholm's alternative. It is also to be stressed that our model has no internal variational structure in the sense that the solution is not seen as the minimizer of some energy. Moreover we analyze the limit case | | → ∞, showing that our model reduces to classical linearized elasticity. We conclude by three explicit computations to illustrate our approach.
2. The incompatibility operator: Generalities and preliminary results.
Let Ω be a sufficiently regular bounded domain of R 3 . We denote by ∂Ω its boundary and by N its outward unit normal. For simplicity we will assume that ∂Ω is C ∞ , but weaker assumptions could be considered for each specific result, depending on the traces and liftings involved. We recall the definition of the incompatibility of a symmetric second order tensor E: with the operator Curl intended row-wise, and with Curl t denoting its transpose. In a Cartesian frame, the incompatibility of a symmetric tensor is obtained by taking its curl column-wise then row-wise, or vice-versa by symmetry. Note that some authors define the operator "curl" as our Curl t , then what they call "curl-curl" coincides with our inc operator, since by symmetry one has inc E = ( inc E) t = Curl t Curl t E.
2.1. The curvilinear frame. For all x ∈ ∂Ω, the system (τ A (x), τ B (x)) is an orthonormal basis of the tangent plane to ∂Ω, that can be naturally extended along N (x) in a tubular neighborhood W of ∂Ω (see [5]). The curvatures along τ A and τ B are denoted by κ A and κ B , respectively. Define the normal derivative as ∂ N := N ·∇ and the tangential derivatives as ∂ R := τ R · ∇, for R ∈ {A, B}. We will also use the notation The following results are proved in [5].
Theorem 2.1. There exist smooth scalar fields ξ, γ A , γ B in W such that ) are oriented along the principal directions of curvature then ξ(x) = 0. These spaces, endowed with the norms defined by E 2 the corresponding inner products are obviously Hilbert spaces. Also, by classical regularization arguments (see e.g. [9,16,37] is dense in each of these spaces. We also define H inc 0 (Ω, S 3 ) = the closure of D(Ω, S 3 ) in H inc (Ω, S 3 ), where the notation D stands for compactly supported C ∞ functions, as well as the trace spaceH where the subscript T stands for the tangential part given by the components (G T ) RR = Gτ R · τ R , R, R ∈ {A, B}. In addition, such a lifting can be obtained through a linear continuous operator Define the subset of C ∞ (∂Ω, S 3 )  Here and in the sequel, for the sake of readability, duality pairings are denoted by integrals.

Green formula and applications.
Recall that the Green formula for the divergence allows us to define, for any E ∈ H div (Ω, S 3 ), its normal trace withφ ∈ H 1 (Ω, R 3 ) an arbitrary lifting of ϕ, see e.g. [23,37]. For the incompatibility operator one has the following counterpart.
Remark 2.1. We have defined T 1 (E) against test functions which admit divergencefree liftings, because spaces of divergence-free tensors arise naturally in problems involving the incompatibility, see the Beltrami decomposition and its consequences in the next sections. But we could also have defined T 1 (E) ∈ H −3/2 (∂Ω, S 3 ) by using a classical lifting in H 2 (Ω, S 3 ). Upon adopting the convention that representatives in H −3/2 (∂Ω, S 3 )/G satisfying the gauge condition (2.10) are chosen, the two definitions are equivalent.
3. Properties of trace operators in H inc (Ω, S 3 ). In this section, homogeneous displacement-like boundary conditions are analyzed in terms of traces of the symmetric strain. These results should be put in perspective with previous results about this problem obtained by Ciarlet and co-authors by means of change-of-metric and change-of-curvature tensors (see [14,15]). Though our characterization is different, the objective of expressing Dirichlet boundary conditions in terms of intrinsic quantities is the same.
To begin with, as particular cases of the two Green formulae recalled in the previous section, one readily obtains the following.
Hereafter, we consider an open and sufficiently regular (C ∞ for simplicity) subset ω ⊂⊂ Ω. If u is a vector or tensor field defined over Ω with well-defined traces on each side of ∂ω, we denote by u the jump of u across ∂ω with inner term counted positively.
Proof. Let Φ ∈ D(Ω, S 3 ). We have in the sense of distributions inc E, Φ = Ω E · inc Φdx, and the Green formula yields In the forward implication the last two integrals vanish by assumption, hence the distribution inc E ∈ D (Ω, S 3 ) is actually an L 2 function. In the converse implication the distribution inc E identifies with an L 2 function, whereby the last two integrals must vanish.
If E ∈ H inc (Ω, S 3 ) and T 0 (E) = T 1 (E) = 0 on ∂Ω, then extending E by 0 and applying Lemmas 3.3 and 3.2 yields inc E N = 0 on ∂Ω. If v ∈ H 1 0 (Ω, R 3 ), then by density of D(Ω, R 3 ) and continuity of the trace operators in H inc (Ω, S 3 ), it follows These two remarks admit the following local versions.
Considering a relatively open subset Γ of ∂Ω and given E ∈ H inc (Ω, S 3 ), we say that T 0 (E) = T 1 (E) = 0 on Γ if the corresponding distributions vanish on Γ, namely and similarly that inc By lifting, this holds true for any v ∈ C ∞ (∂Ω, R 3 ) with support in B. By linearity, covering and partition of unity this extends to any v ∈ C ∞ (∂Ω, R 3 ) with support in Γ.
Proof. Let z ∈ Γ and B be an open ball of center z such that We conclude as in Lemma 3.4.
Proof. On Γ it holds by Lemma 3.5.
We remark that the condition In particular one sees the role of the boundary condition expressed in terms of the Frank tensor Curl t E, namely E = 0 and Curl t E ×N = 0 is equivalent to E = 0 and T 1 (E) = 0.
Lemma 3.8. We have the characterization . By local charts, shifting and convolution with mollifiers, we can define through a standard construction E n ∈ D(R 3 ,

Compatibility conditions and Beltrami decomposition.
In this section we first recall the Beltrami decomposition of symmetric tensor fields, stated here in an L p version for the sake of generality. A specific proof of the L 2 version, which in fact is our main concern, can be found in, e.g., [21,22]. This structure theorem is named after Eugenio Beltrami (1835-1900), an Italian physicist and mathematician known in particular for his works on elasticity, in particular by stating the equilibrium equations of a body in terms of the stress in place of the strain [8] 2 , but also in non-Euclidean geometries in the wake of Gauss and Riemann 3 . We need to introduce first the so-called Saint-Venant-Beltrami condition, originally considered by Saint-Venant in [7], then extended by Donati [17], Ting [38], Moreau [32], Ciarlet and Ciarlet in [13], Geymonat and Krasucki [20], and eventually by Amrouche at al. [3]. Below we give the version found in [28] (originally from [3]).

Theorem 4.1 (Saint-Venant-Beltrami compatibility conditions).
Assume that Ω is simply-connected. Let p ∈ (1, +∞) be a real number and let E ∈ L p (Ω, S 3 ). Then, Let us also refer to [25] and [33] for more details and references on this topic. The following decomposition will show crucial in our model. Pioneer version of this result can be found in [22] for p = 2.
Theorem 4.2 (Beltrami decomposition [28]). Assume that Ω is simply-connected. Let p ∈ (1, +∞) be a real number and let E ∈ L p (Ω, S 3 ). Then, for We call v and F the velocity and incompatibility fields, respectively, associated with E. The following result is the dual counterpart of Saint-Venant's conditions.

Corollary 4.3 (Representation of solenoidal symmetric tensors). Assume that
We now specialize Saint-Venant's decomposition in the case of homogeneous boundary conditions.

SAMUEL AMSTUTZ AND NICOLAS VAN GOETHEM
By the Green formula and the assumptions it holds We arrive at We can now state a converse to Lemma 3.5.
Hence there exists a rigid displacement field r such that v = w + r. On ∂Ω this reduces to v = r.

5.
Orthogonal decompositions of symmetric tensors in L 2 . We assume in this section that Ω is simply-connected.

5.1.
Orthogonal decomposition of L 2 (Ω, S 3 ). In this section we obtain a decomposition of L 2 (Ω, S 3 ) into orthogonal subspaces, in the same spirit as in [22], but to account for more general boundary conditions. We define the spaces and, given a subset Γ of ∂Ω, This is usually stronger than vanishing in the sense of distributions, see e.g. [24] for density and extension properties in fractional Sobolev spaces.
2) Moreover, the velocity field v in (5.1) is unique up to a rigid displacement field. The incompatibility field F in (5.2) is unique.
Remark 5.2. If |Γ| > 0 then the velocity field v in the definition of V 00 Γ is unique.
We have the orthogonal decomposition We infer This completes the proof. We have the following additional property.
Proof. By the Green formula, we obtain Writing K = ∇ S w and applying the Green formula for the divergence yields However, we have w = 0 on Γ 1 while incF N = 0 on Γ 2 , achieving the proof.
We now gather some properties of the spaces Z and Z 0 .
Proof. Let Clearly, X and Y are Banach spaces and Φ is continuous. If E ∈ X and Φ(E) = 0 then inc E = 0 and E is a symmetric gradient by Theorem 4.1. From div E = 0 and EN = 0 on ∂Ω, one obtains E = 0. Hence Φ is injective. By Corollary 4.3, Φ is also surjective. The open mapping theorem entails that Φ −1 is continuous. Hence there exists c > 0 such that and the result follows.
The following result is proved in [23,Theorem 3.8.], see [26,35,45] for L p versions, generalizations and extensions to non simply connected domains.
Theorem 5.5. There exists a constant c > 0 such that for all u ∈ L 2 (Ω, R 3 ) such that div u ∈ L 2 , Curl u ∈ L 2 and u · N = 0 on ∂Ω.
Proposition 5.6 (Poincaré's inequality in Z). There exists c P > 0 such that for Proof. Let E ∈ Z. By Proposition 5.4 we already have Then Theorem 5.5 yields for some other constant c. This completes the proof.
We infer in particular that Z is imbedded in H 1 (Ω, S 3 ) and compactly imbedded in L 2 (Ω, S 3 ).
Proposition 5.7. We have the representation Lemma 5.8. Given a symmetric uniformly positive definite fourth order tensor field B (i.e. B(x)T · T > α|T | 2 ∀T ∈ S 3 for some α > 0 independent of x) such that |B| ∈ L ∞ (Ω), define the linear map L B : Z → Z by Then L B is an isomorphism from Z into Z .
Proof. By Proposition 5.4, L B E, E defines a norm in Z equivalent to the H incnorm. Let T ∈ Z . By the Riesz representation theorem, there exists T ∈ Z such that T , Φ = L B T, Φ for all Φ ∈ Z. Therefore L B is an isomorphism.
We define the inverse map L −1 B : Z → Z, that is continuous by Banach's continuous inverse theorem. Since Z ⊂ L 2 (Ω, S 3 ) ⊂ Z , the restriction L −1 B : L 2 (Ω, S 3 ) → L 2 (Ω, S 3 ) is also well-defined. Proof. The compactness stems from the compact embedding Z → L 2 (Ω, S 3 ), consequence of Proposition 5.6. One has for all E, F ∈ L 2 (Ω, S 3 ) It follows that L −1 B is self-adjoint and positive definite, achieving the proof.
6. Two elliptic boundary value problems for the incompatibility. Lemmas 5.8 and 5.9 yield the following proposition. Similarly we have the following.
Proposition 6.2 (Weak form in Z 0 ). Let K ∈ L 2 (Ω, S 3 ) and B a symmetric uniformly positive definite fourth order tensor field. There exists a unique E ∈ Z 0 such that Moreover, the solution map Φ 0 : K ∈ L 2 (Ω, S 3 ) → E ∈ L 2 (Ω, S 3 ) is linear and compact.

Proposition 6.3 (Strong form in Z).
Let K be such that div K = 0 in Ω and KN = 0 on ∂Ω. Then, the strong form of (6.1) reads Proof. Eq. (6.1) holds actually true for allÊ ∈ Z + V = H inc (Ω, S 3 ), in particular forÊ ∈ D(Ω, S 3 ) and forÊ with arbitrary traces T 0 (∂ NÊ ) andÊ on ∂Ω, by Theorem 2.3. Then the Green formula provides the strong form, which is seen to be equivalent to the weak form.   The second and third order tensor fields σ and τ are called the (Cauchy) stress and hyperstress tensors, respectively. Moreover, the material is supposed to be linear, homogeneous and isotropic within each Ω p , which is represented by the constitutive laws where E is the strain, C p is standard Hooke's tensor and D p is Mindlin's tensor [31]. These constitutive laws read componentwise σij = λδijE kk + 2µEij, (7.2) where λ, µ, c 1 , ..., c 5 are constants assigned to each Ω p (index p is dropped for readability). Assumption 1 yields Σ = σ − div τ ∈ L 2 (Ω, S 3 ).
Assumption 4. The hyperstress τ does not produce any virtual intrinsic power as soon as the strain E is compatible. This means or equivalently inc E = 0 ⇒ − div τ = 0 in Ω. From expression (7.3) we derive the existence within each Ω p of a constant p such that c 1 + c 5 = − p , c 2 = p , c 3 = − p /2, c 4 = p /2, leading to − div τ = p inc E (see details in [6]).
Conclusion. We denote = p χ Ωp and C = C p χ Ωp , whereby σ = CE and − div τ = inc E in Ω. The expression of the internal virtual power is The scalar field is called incompatibility modulus, as it expresses the resistance of the material against incompatible deformations. Subsequently we will extend the model to the case where is a sufficiently regular function of x ∈ Ω.

7.2.
Power of external efforts. The power of external efforts is assumed to be a linear functional on L 2 (Ω, S 3 ). By Riesz representation, there exists K ∈ L 2 (Ω, S 3 ) such that We emphasize that the power of external efforts may be at first expressed in terms of the non-objective fieldsv andF of the Beltrami decomposition ofÊ. However, provided attention is paid to the uniqueness of the decomposition, these fields are themselves linear functions ofÊ. This will specified in Section 10.1.

Virtual power principle. The virtual power principle in the absence of inertia reads
for allÊ ∈ L 2 (Ω, S 3 ) satisfying possible kinematical constraints. In the absence of kinematical constraints, (7.4) is obviously equivalent to CE + inc E = K. Within an incremental formulation, C and are generalized elastic tangent moduli. They need to be updated at each increment as soon as nonlinear phenomena occur. The stress-strain relation is therefore piecewise linear. Typically, in a region with plastic deformations, the Lamé coefficients and the incompatibility modulus are expected to be less than in purely elastic regions. The way these coefficients evolve is driven by nonlinear constitutive laws that substitute to flow rules and hardening models.
8. Solution of elasto-plasticity equations with natural boundary condition. The main problem we address is the following: given K ∈ L 2 (Ω, S 3 ), find E solution of (7.4).
8.1. Kinematical setting. We will limit ourselves to the case where no kinematical constraint is assumed on the virtual strainÊ. Therefore, the problem reduces to (7.5). However, we will see that the absence of constraint on E leads to nonunique solutions, and that uniqueness is obtained by prescribing the incompatibility flux inc E N on ∂Ω. The homogeneous case inc E N = 0 is studied first. However prescribing a given value, either 0 or not, for the incompatibility flux may seem artificial from a modeling point of view. This issue is related to the characterization of the behavior of dislocations at interfaces, whose difficulty is emphasized in [30] in the case of grain boundaries for polycrystals. An attempt to determine the incompatibility flux through a domain extension technique is proposed in Section 10.3.
Remark 8.1. A particular kinematical setting is to require K ∈ V, and a very special case occurs when K = ∇ S v with div v = tr K constant. Then for C constant a solution to CE + inc E = K is E = C −1 K. Indeed by the structure of C −1 one has E proportional to K plus a constant tensor hence inc E = 0.

8.2.
Well-posedness with vanishing incompatibility flux. For the subsequent mathematical analysis it is not required that C be an isotropic Hooke tensor. We denote by S 3×3 the set of symmetric fourth order tensors acting on 3 × 3 matrices. Our main result is the following.
. Let c P be the Poincaré constant of Proposition 5.6. If C is uniformly positive definite and either > c P |C| a.e. or < −c P |C| a.e., then there exists one and only one E ∈ F such that Moreover we have the a priori estimate Proof. We write the problem as with B := C −1 and H := C −1 K. Note that B is uniformly positive definite if > 0 and uniformly negative definite if < 0. We will first prove uniqueness and then existence.
Step 1. Uniqueness. Let E ∈ F be such that Then By incF N = 0 on ∂Ω the first integral vanishes. Specifically, takeF = E i . We obtain By Proposition 5.6 we obtain L ∞ < c −1 P we infer inc E i = 0 then E i = 0, by Proposition 5.6. Thus (8.4) yields E c = 0, and eventually E = 0.
itself, by Proposition 5.7, equivalent to The operator L −1 B M : L 2 (Ω, S 3 ) → L 2 (Ω, S 3 ) is compact, since it is continuous from L 2 (Ω, S 3 ) to Z 0 and Z 0 is compactly embedded in L 2 (Ω, S 3 ) as consequence of Proposition 5.6. Furthermore, under the condition B −1 L ∞ < c P −1 , the operator is injective due to the uniqueness claim. Thus, Fredholm's alternative provides the existence of E i ∈ L 2 (Ω, S 3 ) solution of (8.9).
Let us turn to (8.7)(a). We have to find This is a standard linear elasticity problem.
Proposition 5.6 yields from which we arrive at (8.1).
Remark 8.2. The solution space F encompasses the transmission conditions stated in Lemma 3.3. In particular, no tangential slip along internal surfaces can occur. This is in contrast with classical formulations in perfect plasticity where spaces of bounded deformations are involved, see e.g. [19]. Moreover, the continuity of the incompatibility flux across internal surfaces stated in Lemma 3.2 shows that inc E N = 0 at the interface between a region with incompatible strain and a purely compatible region. Prescribing inc E N = 0 on ∂Ω can be interpreted as considering the exterior of Ω as filled with a purely compatible phase, with Ω and its exterior forming a continuum. Other assumptions, such as discussed in Sections 8.3 and 10.3, are needed in order to enhance incompatibilities near the boundary. implies that G := inc E converges weakly in L 2 (Ω, S 3 ) to some G as → ±∞, up to a subsequence. More precise limiting results will be given in the next section.
On the other side, the condition | | > c P |C| suggests that some degeneracy occurs when | | goes to 0, unless C also tends to 0. When → 0 the incompatibility is no longer controlled, in particular the transmission conditions recalled in Remark 8.2 do not hold any more. However one can let C go to 0 keeping constant, in order to represent a nearly void material. This will be considered in Section 10.3.
Let c P be the Poincaré constant of Proposition 5.6. If C is uniformly positive definite and and either > c P |C| a.e. or < −c P |C| a.e. then there exists one and only one E ∈ H inc (Ω, S 3 ) such that Moreover there exist constants c 1 and c 2 such that Existence and uniqueness follow from Theorem 8.1. The a priori estimate of Theorem 8.1 combined with Proposition 5.6 and standard elliptic regularity provide (8.11).
Obviously the same holds for a sequence k ∈ L ∞ (Ω, R * − ) with inf Ω | k | → +∞. Proposition 9.2. If is constant, K ∈ L 2 (Ω, S 3 ), and E ∈ F satisfies Note that the above relation is similar to a linear elasticity system, nevertheless E may not be a symmetric gradient.
Proof. TakeÊ ∈ V and observe that due to the assumptions, one has Ω inc E ·Êdx = 0. Theorem 9.3 (Elastic limit: homogeneous flux). Assume that C ∈ L ∞ (Ω, S 3×3 ) uniformly positive definite and K ∈ L 2 (Ω, S 3×3 ) are fixed and that = 0 is constant. Let E ∈ F be the unique solution of CE + inc E = K in Ω and let E ∞ ∈ V be the unique solution of the linear elasticity problem Proof. Existence and uniqueness for (9.1) is a consequence of the Riesz representation theorem in the Hilbert space V for the inner product (E,Ê) → Ω CE ·Êdx.
Substracting (9.1), one has By Propositions 9.1 and 5.6 we have Hence, as | | → +∞, one retrieves the standard linear elasticity problem with Neumann boundary conditions. In case of non-vanishing incompatibility flux the following holds.
Proof. Existence and uniqueness for (9.2) follows from the Riesz representation theorem. For allv ∈ H 1 (Ω, R 3 ) one has Hence Consider the decomposition E = E c + E i ∈ V ⊕ Z, see Proposition 5.3. By (8.11), one infers inc E L 2 → 0, and subsequently, by Proposition 5.6, E i L 2 → 0. Finally, 10. Interpretation of the kinematical framework and external efforts.
10.1. External efforts. Consider a virtual strainÊ ∈ L 2 (Ω, S 3 ) decomposed aŝ The work of the external efforts againstÊ reads The fields − div K and KN are recognized as classical body and contact forces, while inc K and (T 0 (K), T 1 (K)) are body and contact forces working against the divergence-free part of the virtual strain. The above fields are in principle known in the first place. The issue is then how and under which conditions it is possible to construct a corresponding tensor field K. Formally the boundary forces KN , T 0 (K) and T 1 (K) exhibit some coupling, as stressed in [6]. To address these points one must specify a kinematical framework ensuring the uniqueness of the decomposition (10.1).

10.2.
Kinematical framework. TakeÊ = ∇ Sv + incF with ∇ Sv ∈ V 00 Γ1 and incF ∈ W 0 Γ2 for some partition Γ 1 ∪ Γ 2 of ∂Ω, with Γ 1 relatively open in ∂Ω. As said above, f := − div K is identified with the body force, and g := KN is identified with a surface load on Γ 2 . Now, if K ∈ V 00 Γ1 the last two integrals of (10.2) vanish by virtue of Lemma 5.2. Then (10.2) rewrites in the classical form More generally we have the following existence result.

SAMUEL AMSTUTZ AND NICOLAS VAN GOETHEM
The standard strong form reads Setting K = inc F +∇ S v fulfills all the required conditions (recall Corollary 3.6).
Observe that one cannot prescribe T 0 (K) and T 1 (K) on Γ 2 , although they might yield a contribution in the last integral of (10.2). In fact, if G = 0, one can write K = ∇ S w for some w ∈ H 1 (Ω, R 3 ), and as in the proof of Lemma 5.2 one has for allF ∈ H inc (Ω, S 3 ) In the framework under consideration where incF ∈ W 0 Γ2 , the integral at the righthand side of (10.4) is localized on Γ 2 , and specifically w appears as work-conjugate to the virtual incompatibility flux on Γ 2 .
Remark 10.1. If one prescribes T 0 (K) = T 1 (K) = 0 on Γ 1 , then there is no virtual work associated with this boundary condition even if the virtual strain is allowed to have tangential components. Indeed, in the present model, both for finite and in the elastic limit, the values of T 0 (E) and T 1 (E) are unconstrained, hence allowing for slip at the boundary. But as | | → ∞, one may want to retrieve a Dirichlet boundary condition on Γ 1 , including T 0 (E ∞ ) = T 1 (E ∞ ) = 0, a property which is not verified, see Theorem 9.3. To obtain such a condition in the elastic limit, one may either prescribe kinematical constraints in the space of solutions or incorporate some generalized force working against boundary slip. Let us recall that tangential slip along Dirichlet boundaries are present in classical models of plasticity, see e.g. [19] where it is related to concentrated (in a measure sense) plastic deformation. For instance, setting a boundary condition of the type T i (K) = αT i (E) (i = 0, 1) (i.e., taking P = αE in Proposition 10.1), for some constant α and with E the sought solution itself, we are led to solve a coupled system for (K, E), for which the existence of solutions for α small enough can be inferred from the Banach fixed point theorem or the Neumann series. To see that this formulation is likely to impose a Dirichlet boundary condition in the elastic limit, consider a case where ∂Ω = Γ 1 , G = 0 and E = ∇ S v for some v ∈ H 1 (Ω, R 3 ). Then the condition is reached if and only if T 0 (E) = T 1 (E) = 0 on ∂Ω, by Theorem 2.3, which is itself equivalent to v being a rigid displacement on ∂Ω, by Proposition 4.5. This is expected if one lets α tend to infinity with in the governing system. However, the precise asymptotic analysis, together with the existence issue, are beyond the scope of the present work.
10.3. Alternative to the vanishing incompatibility flux condition. The condition inc E N = h on ∂Ω is related to the flux of dislocations at the boundary of Ω. This can be specified using Kröner's formula inc E = Curl κ, with κ the contortion tensor related to the dislocation density tensor. Prescribing an a priori given h (possibly zero) may seem an artificial or ad-hoc condition. An alternative is to consider that the exterior of Ω is filled with a fictitious material that mimicks void, with transmission conditions representing the fact that the two phases constitute a continuum. Therefore we restrict ourselves to the pure Neumann type boundary condition, that is, a surface load without kinematical restriction. For a full space extension, K = ∇ S w is defined by with f extended by 0 outside Ω. The equilibrium equation CE + inc E = K is fulfilled over R 3 with (C, ) extended by (C ext , ext ) outside Ω. In order to approximate a Neumann condition one needs that C ext be chosen significantly smaller than C within Ω. Then on ∂Ω one has with K = ∇ S w as found above. It is reasonable to assume that is continuous across ∂Ω, in such a way that ∂Ω has a neutral effect on the transport of dislocations (transmission without reflection). Under this assumption we have Existence and uniqueness of a solution has been shown in Theorem 8.2. In addition we derive on ∂Ω which is obviously the standard Neumann condition on the Cauchy stress.
Let us now examine the limit case. In view of Theorem 9.4, we have This rewrites as The standard Neumann elasticity problem in Ω is retrieved.
11. Interpretation of the incompatibility modulus by numerical examples. The purpose of the following examples is to show how the incompatibility modulus can be related to dissipative processes. Of course, irreversible behaviors in the form of plasticity are most often localized within regions that typically grow when loading increases. Here, in order to obtain explicit expressions for all quantities and illustrate the limits stated in Theorems 9.3 and 9.4, we assume that is constant in space.
Since we have only set up linearized governing equations, we simply define the free energy of the solid body submitted to a given load as the external work that can be recovered by purely elastic unloading, keeping C invariant for simplicity. To be more precise, in the small strain framework, call K(t) the load at current time t, and E(t) the associated total strain. In our model, E(t) is typically obtained through CĖ + incĖ =K over [0, t], E(0) = 0 and a boundary condition on the incompatibility flux. For some T > t consider a virtual extension of K over [t, T ] by a smooth function such that K(T ) = 0. The free energy is then defined by

SAMUEL AMSTUTZ AND NICOLAS VAN GOETHEM
where E rev (s) = E(t) + ∇ S u rev (s), ∇ S u rev (s) = s t ∇ Su rev (τ )dτ , and for all τ ∈ [t, T ],u rev (τ ) ∈ H 1 (Ω) satisfies Equation (11.1) ensures thatΨ(t) = Ω K(t) ·Ė rev (t)dx, which is the external power expenditure in case of reversible transformation. Furthermore, on the one hand, we have Ψ(t) = − T t Ω K · ∇ Su rev dxds, which after integration by parts and using (11.2) yields On the other hand, integrating (11.2) in time leads to The two above equations give a practical way to compute the free energy, i.e. Ψ(t) = − 1 2 Ω K(t)·∇ S u rev (T )dx, and show that this latter is independent of the path along which K is driven to 0.
The dissipation rate is then classically defined as D := P −Ψ, where P(t) := Ω K(t) ·Ė(t)dx is the external power expenditure. We arrive at In the following examples we will consider finite time increments in which K is constant and compute the dissipated energy. As said above, we will limit ourselves to spatially constant incompatibility moduli. Such configurations will be thermodynamically consistent under the condition that the dissipation be positive, which we expect.
This ordinary differential equation leads us to assume that < 0, since in the other case the solutions do not decay when |z| → ∞. We denote We obtain: 1. For |z| < h, Observe that which is the classical elastic solution with uniaxial strain. Let ∇ S U +E 0 be the Beltrami decomposition of E in the domain Ω a := {|z| < a} such that E 0 ∈ W 0 ∂Ωa (see Theorem 5.1), i.e., div E 0 = 0 in Ω a and E 0 N = 0 on ∂Ω a . One has div ∇ S U = div E in Ω a , Denoting by u the z component of U we obtain u = ψ , u (a) = ψ(a), u (−a) = ψ(−a). Thus u = ψ and, setting u(0) = 0, We obtain in particular

SAMUEL AMSTUTZ AND NICOLAS VAN GOETHEM
The functions ϕ, ψ are plotted in Figure 1 for h = 1, Young modulus Y = 10 and Poisson ratio ν = 1/3. The value of u(h) as a function of is also depicted. As expected, we observe an increase of elongation as | | decreases. Therefore, if we consider a quasi-static small-strain experiment consisting of a first increment with finite and a further purely elastic increment with opposite load (unloading), a residual positive elongation remains. This shows that a portion of the external work has been dissipated. One of the main outcomes of this example is that physically acceptable solutions are obtained on the condition that < 0. Henceforth this condition will be assumed. 11.2. 2D case: Cylinder under uniform radial traction. We consider a twodimensional model of the variables (x, y) and we assume that One has in Cartesian coordinates Elementary algebra leads to Within regions where inc K = 0, div K = 0 and λ, µ, are constant, substituting the above relations in the last equation of (11.6) entails Then, Hooke's law together with E rr + E θθ = u + v = (1 − 2λh − ∆h)/(2(λ + µ)) yield and The condition E θθ = 0 rewrites as Next, from T 1 (E) = −he r ⊗ e r + (h − h )e θ ⊗ e θ we infer h = 0. Coming back to (11.7) one looks for withh solution of the homogeneous equation in B and R 2 \B. In view of (11.8) and recalling that < 0,h is spanned by the Bessel-type functions J 0 (k + r), Y 0 (k + r), I 0 (k − r) and K 0 (k − r) with Due to boundedness and decay at infinity, it remains h(r) = aJ 0 (k + r) + bI 0 (k − r) if r < 1 cK 0 (k − r) if r > 1.
The three transmission conditions fix a, b and c through a linear system. In the following computations we take a material inside B with Young modulus Y = 10 and Poisson ratio ν = 1/3, and a nearly void exterior phase such that C e = 10 −5 C i . The incompatibility modulus is taken constant over R 2 . The plots of the strain are given in Figure 2, and compared with the classical plane strain elastic strain. The external work Table 1. This shows that decreasing the incompatibility modulus (in absolute value) increases the external work for a given load, thus increases dissipation.
−1000 −100 −20 W 0.3006 0.4616 1.0620 Table 1. External work x ∈ R 3 , |x| < 1 . We assume a uniform unit radial traction on ∂Ω. We treat the two kinematical frameworks addressed in this paper, namely the case of vanishing incompatibility flux and the case of free incompatibility flux through domain extension described in Section 10.3. 11.3.1. Case 1: Vanishing incompatibility flux. We assume in this case the condition inc E N = 0 on ∂Ω. We have K = ∇ S w with − div ∇ S w = 0 in Ω, The solution is immediately found as w = re r and K = I. Considering the form E = ϕ(r)I + ψ(r)e r ⊗ e r we have CE = ((3λ + 2µ)ϕ + λψ) I + 2µψe r ⊗ e r , and (see [40]) inc E = ϕ + ϕ r − ψ r I + −ϕ + ϕ r + ψ r − 2ψ r 2 e r ⊗ e r .
The exterior deformation field is completely determined. We now compute the displacement U such that E = ∇ S U + E 0 , div E 0 = 0, E 0 N → 0 at infinity. Therefore U solves div ∇ S U = div E in R 3 ∇ S U N → EN at ∞. (11.12) For U = u(r)e r , the first equation reads when r = 1 u + 2 u r − 2 u r 2 = ϕ + ψ + 2 ψ r .
Substituting (11.11) leads to the equation 1 r 2 (r 2 u) = 2λ 3 (λ + 2µ) Using (11.10) we arrive at u(r) = −λ 3µ(3λ + 2µ) for some constant e. On ∂Ω, (11.12) amounts to ∇ S U N = EN , that is u + 2u = ψ , hence we have the jump relations u = 0 and u = ψ . This fixes the constants d and e through a linear system. A Taylor expansion provides i 1 (s) = s/3 + o(s 2 ) as s → 0. Then for → −∞ it is easily checked that the elastic solution is retrieved, namely ϕ ∞ (r) = 1 3λ + 2µ χ {r<1} + 1 3λ + 2µ The functions ϕ, ψ and u are plotted on Figure 3 for different values of . The curves of displacement show an increase of dilation due to inelastic deformation. This again highlights a dissipative process, since the work R 3 K · Edx is proportional to u(1). 12. Concluding remarks. The main goal of this work was to give a rigorous mathematical basis to the model presented in [6], as well as to complete and refine the analysis of the incompatibility operator conducted in [5]. Our study showed that a general model of small-strain generalized elasticity could be considered to account for the strain incompatibility and hence for the presence of dislocations at the micro-scale. Moreover, classical infinitesimal elasticity is recovered as a limit case as the incompatibility modulus tends to minus infinity. Further, it was shown on examples that this model was able to represent dissipative processes through residual deformations after unloading. Our next step is twofold: first to devise a computational algorithm, based on shape / topological sensitivity analysis, in the spirit of [1,46], to simulate the time evolution of nonlinear irreversible effects, and second to understand the relation between this model and accepted models of elasto-plasticity. This challenge will require theoretical as well as computational efforts, expected in future works.