pp. X–XX NUMERICAL SOLUTION OF A VARIATIONAL PROBLEM ARISING IN STRESS ANALYSIS: THE VECTOR CASE

In this article, we discuss the numerical solution of a constrained minimization problem arising from the stress analysis of elasto-plastic bodies. 
This minimization problem has the flavor of a generalized non-smooth eigenvalue problem, with the smallest eigenvalue corresponding to the load capacity ratio of the elastic body under consideration. 
 
An augmented Lagrangian method, together with finite element approximations, is proposed for the computation of the optimum of the non-smooth objective function, and the corresponding minimizer. 
 
The augmented Lagrangian approach allows the decoupling of some of the nonlinearities and of the differential operators. 
Similarly an appropriate Lagrangian functional, and associated Uzawa algorithm with projection, are introduced to treat non-smooth equality constraints. 
Numerical results validate the proposed methodology for various two-dimensional geometries.

1. Introduction and Motivations.The stress analysis of elasto-plastic bodies is a crucial problem for engineers who want to estimate maximal stresses and fractures for a given body and geometry [28,31].The ultimate goal is to compute or estimate the threshold sustained by the material before it fractures.Stresses are induced by external forces, acting on the body and/or its boundary, such as tension, torsion, bending, etc.
In nonlinear elastic systems, the computation of stresses in a given body under external forces gives rise to non-smooth, nonlinear variational problems, see e.g.[9,20,40].Finite element methods are commonly used for the numerical solution of such problems, see e.g.[10,44], and numerous works exist in the literature concerning the numerical approximation (via finite elements, mixed finite elements, or discontinuous finite element methods) of stresses, fractures in elastic or plastic materials [22,29,37].On the other hand, the modeling of cracks can be tackled with level sets methods for instance, see e.g.[23].
One related problem is the so-called limit analysis problem, that corresponds to the computation of the limit (maximal) load the elasto-plastic body can sustain.Upper and lower bounds for the magnitude of the stresses have been obtained in [33,34,43].Mixed finite element methods have been discussed in [42] for the The numerical solution of a non-smooth minimization problem arising in limit analysis is investigated in this article.This problem is based on the limit analysis problem introduced in [40] and aims at computing the load capacity ratio as in [32] for given geometries.In this article, we consider an external traction on the boundary only (see [30,40] for a theoretical investigation of the addition of body forces).
We assume that the boundary of the bounded domain Λ is sufficiently smooth and can be decomposed into ∂Λ = γ 0 ∪ γ 1 ∪ γ 2 , such that γ i ∩ γ j = ∅, i, j = 0, 1, 2, i = j.An external surface traction acts on the boundary γ 1 , while the body is fixed on the boundary γ 0 .Consider the case when the domain Λ is an infinite cylinder oriented along one direction of space: when applying an external traction field parallel to the invariant direction, one obtains a scalar-valued problem, that has been addressed in [6].When the external surface traction field is oriented perpendicularly to the invariant direction of Λ, the limit analysis problem becomes a vector-valued two-dimensional problem (in a cross-section of the domain Λ perpendicular to the invariant direction) that is addressed here.
The model problems read as follows: let Ω be the two-dimensional domain obtained by cutting Λ perpendicularly to its invariant direction, and Γ i = γ i ∩ Ω, i = 0, 1, 2. We look for approximate solutions to the following minimization problems where D(v) = 1/2 ∇v + (∇v) T and Σ1 = v ∈ BD(Ω), v| Γ0 = 0, where BD(Ω) is the space of the bounded deformation vector-valued functions [41].Similar questions have been addressed with Newton methods in [18,27] in the field of elasto-plastic limit analysis.Such non-smooth optimization problems arise also in calculus of variations, and correspond to the computation of generalized, eigenvalues of non-smooth operators [4,7,26].Appropriate solution methods are required, see e.g.[11,20,21].In this article, an augmented Lagrangian approach is advocated, that does not introduce any regularization or convexification parameters [15,24,26], keeping in mind that augmented Lagrangian methods are equivalent to alternating direction implicit (ADI) methods, as shown in, e.g., [20,21].
In Section 2, the derivation of the model problem from the limit analysis problem is presented and the formulation as a vector-valued generalized eigenvalue problem is justified.In Section 3, we discuss the solvability of the problems considered in Section 2. Section 4 details an augmented Lagrangian algorithm for the solution of the non-smooth optimization problems, that decouples differential operators and nonlinearities in the objective functions and constraints.A finite element discretization method is presented in Section 5.The result of numerical experiments are presented in Sections 6 and 7 for various geometries, and the convergence of the finite element approximations is discussed in the last section.
The authors think that this article is very appropriate to a volume dedicated to Professor Roger Temam, since Professor Temam has contributed in some essential ways to many of the issues addressed in this article.

2.
From the Limit Analysis Problem to a Model Problem.We consider an elasto-plastic body that occupies a domain Λ ⊂ R 3 , whose smooth boundary is partitioned into ∂Λ = γ 0 ∪γ 1 ∪γ 2 such that γ i ∩γ j = ∅, i = j.Body forces f ∈ L 2 (Λ) 3  and boundary forces g ∈ L 2 (γ 1 ) 3 are applied to the elastic body.Classical results, see e.g.[10,21,40], allow to compute the elastic displacement as the solution of a variational problem, namely where D(v) = 1/2 ∇v + (∇v) T is the deformation tensor (also called rate of strain tensor by some practitioners), ψ(•) is a proper, lower semi continuous convex function that characterizes the material properties, and 3 , v| γ0 = 0 .The function ψ(•) can be expressed for instance in terms of the Lamé coefficients in linear elasticity models.Problem (1) has a finite value objective function when the elastic body sustains external forces.Following [40], the following problem gives an equivalent criterion to predict if the body fractures.Consider: where 2) is taken on the set of functions with divergence free constraint.Closely related to (2), we introduce: where W 0,+ = {v ∈ W 0 , L(v) > 0}.This particular optimization problem is not convex but, because of the homogeneity of the L 1 -norm, it is equivalent to the so-called limit analysis problem: where 4) is a non-smooth convex optimization problem.Following [40], we have the following result: Result.The infimum of ( 1) is finite, meaning that the elastic material does not fracture, if and only if the value of the objective function in ( 4) is greater than or equal to one.
Following the above statement, it makes perfect sense to look for a numerical approximation of the objective function in (4), and, later on, of δ 1 and δ 0 , for given three-dimensional elastic bodies Λ, in order to determine a priori the presence of fractures in the elastic body.
We consider here the case without body forces f = 0, and with normalized tractions g ∈ L 2 (γ 1 ) 2 on the surface γ 1 that satisfy |g| ≤ 1.In that case, by noticing that |v| = g • v for some g ∈ L 2 (γ 1 ) 2 , |g| ≤ 1, we consider the particular model problem (that corresponds to the load capacity problem proposed in [32]): where We focus in the sequel on infinite cylindrical three-dimensional domains, that is on infinite cylinders Λ, with a bounded two-dimensional cross-section Ω.Let us consider an infinite cylinder Λ ⊂ R 3 , oriented for instance along the Ox 3 axis, and traction fields g ∈ L 2 (γ 1 ) 2 , |g| ≤ 1 on the boundary γ 1 that are orthogonal to the Ox 3 axis, as illustrated in Figure 1.The elastic displacement v ∈ H 1 (Λ) 3 can be written therefore as v(x 1 , x 2 , x 3 ) = (v 1 (x 1 , x 2 , x 3 ) , v 2 (x 1 , x 2 , x 3 ) , 0) T , and v k (x 1 , x 2 , x 3 ) = v k (x 1 , x 2 ), k = 1, 2, by invariance along the Ox 3 axis.Problem ( 5) is written therefore as a two-dimensional vector-valued problem as follows: let Ω ⊂ R 2 be the bounded domain (in the Ox 1 x 2 -plane) obtained by cutting Λ perpendicularly to Ox 3 , with a smooth boundary ∂Ω = Γ 0 ∪ Γ 1 ∪ Γ 2 , Γ i = γ i ∩ Ω, j = 0, 1, 2. This notation is illustrated in Figure 2 in the case of a rectangular two-dimensional bar clamped into a vertical wall.
. Two-dimensional domain Ω: sketch and notation.The boundary ∂Ω is split into ∂Ω = Γ 0 ∪ Γ 1 ∪ Γ 2 ; the body is assumed to be fixed on Γ 0 (for instance when the bar is stuck into a wall); a traction is imposed on Γ 1 .The boundary Γ 2 may be empty.Within this framework, the solution of (5) corresponds to solving the following non-smooth problem (with dx = dx 1 dx 2 and v = (v 1 , v 2 ) T ): where If δ 1 ≥ 1, the infinite cylinder of cross-section Ω sustains the surface traction.Otherwise, the material presents a fracture.
Remark 1.In [6], the corresponding scalar case is obtained when applying a surface traction that is parallel to the Ox 3 axis.The corresponding problem in a crosssection Ω involves a scalar-valued function.
Remark 2. The quantity δ 1 defined by ( 6) is equivalent to the inverse of the load capacity ratio defined as in [32].
In parallel, we investigate also the related problem that does not contain the divergence free constraint, that is where The quantities δ 0 and δ 1 correspond to generalized eigenvalues of non-smooth operators.In the following sections, we present a numerical method to compute finite element approximations of these eigenvalues, together with approximations of the corresponding generalized eigenfunctions (i.e. the minimizers of the finite dimensional approximations of ( 6) and ( 7)).
Remark 3. It is worth noticing that in ( 6) and (7), we look for infima and not for minima 3. On the Solvability of the Constrained Optimization Problems.The problems ( 6) and (7) are formulated in H 1 (Ω) 2 for reasons that are related to the numerical method used in the following.Indeed, in Section 4, a quadratic regularization term that is typically related to the L 2 -norm of the gradient is added to the objective function by means of an augmented Lagrangian formulation.A minimizing sequence of approximations of the eigenfunctions lying in a sub-space of H 1 (Ω) 2 is shown to converge to a function that is not in H 1 (Ω) 2 .
However, both problems are actually not well-posed in H 1 (Ω) 2 , as they do not admit a solution in that space in the general case.Numerical results in Section 6 and 7 actually show that fractures arise in the plastic body, emphasizing that the eigenfunction is indeed not in H 1 (Ω) 2 .
With these definitions, one can show that (6) (7) have minimizers in BD(Ω); the proof can be found in [40], and the application to load capacity ratio is detailed in [32].The uniqueness of the solution depends on the boundary conditions that actually define the constraints of the optimization problem.Therefore, uniqueness is not guaranteed in the general case (see also [5,7] for other eigenproblems with a similar behavior).
Remark 4. Similar results hold for the least gradient problem [7,25,35,36], where the constrained minimization of Ω |∇v| dx (in the case of scalar-valued function) leads to 'natural' solutions that are in BV (Ω) the space of functions of bounded variations on Ω.
For the reasons specified earlier, we formulate from now on ( 6) and ( 7) in H 1 (Ω) 2 , keeping in mind that (i) these problems do not have minimizers in H 1 (Ω) 2 , but we can still exhibit a minimizing sequence that will converge to an eigenfunction in BD(Ω); (ii) we are mainly interested in the computation of the value of the objective function (the generalized eigenvalue) but not so much in the computation of the corresponding eigenfunction.
4. An Augmented Lagrangian Algorithm.In the sequel, we address the solution of ( 6) and ( 7) by an augmented Lagrangian method (keeping in mind the equivalence with ADI methods).We define q = D(v) ∈ L 2 (Ω) 2×2 and the following space: Problems (6)(7) are equivalent to where We define the following scalar product where Let r be a positive parameter.The augmented Lagrangian method discussed here is inspired from [20,21], and consists in searching for a saddle point of the following augmented Lagrangian functional: Namely, we are looking for {u, p, λ} for all {v, q, µ} ∈ An Uzawa-type algorithm for the solution of ( 12) reads as follows: let u −1 ∈ V i , λ 0 ∈ L 2 (Ω) 2×2 be arbitrary given functions.Then, for n ≥ 0: (a) Solve where (13) does not involve any derivatives and can be solved point-wise a.e. in Ω.It leads to the following closed form solution for p n [7, 12, 13]: More details are given in the next section, when discussing the finite element discretization.
(b) Solve The solution of problem ( 14) is detailed in the sequel.(c) Update the multipliers λ n+1 ∈ L 2 (Ω) 2×2 as follows: where ε is a given tolerance.The augmented Lagrangian algorithm produces a sequence of iterates {u n } n≥0 that eventually converges to the function realizing the infimum of (6) (7).
Remark 5. We denote by BV (Ω) the space of functions of bounded variation over Ω (see, e.g., [39,Chapter 37] for a complete definition).In [6], one has seen that the (scalar-valued) sequence {u n } n≥0 may converge to a limit in BV (Ω)\H 1 (Ω).Numerical results in Sections 6 and 7 show a similar behavior for the sequence {u n } n≥0 in the vector-valued case, the limit being in BD(Ω)\H 1 (Ω) 2 .4.1.On the Solution of the Constrained Elliptic Problems.In order to solve the constrained minimization problem (14), we first note that, for all and µ satisfies the fixed-point relation [14]: where ρ > 0 is given and P Λ is the orthogonal projector from L 2 (Γ 1 ) 2 onto Λ, and Actually, the orthogonal projection can be explicitly written as Based on this result, an embedded Uzawa algorithm for the solution of ( 14), that deals with the non-smooth constraint, is advocated.Let ξ 0 ∈ Λ be given, then, for k ≥ 0: (i) Compute: where For i = 0 (resp.i = 1), the solution of problem ( 18) is addressed in Sections 4.2 (resp., Section 4.3).(ii) Update the multiplier ξ k via (16), that is where ρ > 0 is small enough.until convergence is achieved for the variables ξ and u n .Typically the stopping criterion is u n,k − u n,k+1 L 2 (Ω) 2 < ε ′ , where ε ′ is a given tolerance.When convergence is achieved, say at iteration k, we set u n = u n, k and proceed to the next step.
4.2.Uzawa Algorithm for the Case without Incompressibility Constraint.Details about the solution of (18) are provided here when i = 0 and v ∈ V 0 .Namely, for ξ k ∈ Λ given, we want to solve min where troduced to handle the (linear) scalar equality constraint.The corresponding Lagrangian functional reads: The first order optimality conditions read as follows: find u n,k ∈ V 0 and χ k ∈ R such that, for all v ∈ V 0 : Equation ( 21) is linear with respect to the multiplier χ k .It therefore suffices to solve (21) for two values of χ k , e.g.χ k = 0 and χ k = 1, with corresponding solutions being u n,k 0 and u n,k 1 , and then consider the linear combination of u n,k 0 and u n,k 1 that satisfies (22), namely: 4.3.Uzawa Algorithm for the Case with Incompressibility Constraint. where A Lagrange multiplier χ k ∈ R is introduced to handle the scalar equality constraint Γ1 ξ k • vds = 1, while a second multiplier p n,k ∈ L 2 (Ω) is introduced to account for the incompressibility constraint ∇ • u n,k = 0 (that corresponds to a pressure in a generalized Stokes problem formulation).The corresponding Lagrangian functional reads: Remark 6.For well-posedness reasons, the multiplier p n,k should be in and no forces are applied on the boundary.This case of little interest is not included in the sequel, and p n,k ∈ L 2 (Ω).
The first order optimality conditions corresponding to (24) read as follows: find u n,k ∈ V 0 , p n,k ∈ L 2 (Ω) and χ k ∈ R such that, for all v ∈ V 0 and for all q ∈ L 2 (Ω): The Stokes-type system (25)(26) being linear with respect to the multiplier χ k , we can proceed as in Section 4.2.We solve therefore (25) (26) for χ k = 0 and χ k = 1, the corresponding solutions being (u n,k 0 , p n,k 0 ) and (u n 1 , p n,k 1 ), and then consider the linear combination of (u n,k 0 , p n,k 0 ) and (u n 1 , p n,k 1 ) that satisfies (27), namely: where f n,k is given by (23).
5.1.Generalities.Finite element techniques are used for the numerical implementation of algorithm ( 13)- (15).Let h > 0 be a space discretization step.A family {Ω h } h of polygonal approximations of the domain Ω is introduced such that lim h→0 Ω h = Ω, together with lim h→0 Γ 0,h = Γ 0 , and lim h→0 Γ 1,h = Γ 1 .We consider a triangulation T h of the domain Ω h satisfying the usual compatibility conditions between triangles (see for instance [19] for a precise definition).Let N e denote the number of elements of T h , N n the number of vertices of T h in Ω h \Γ 0,h , and N nt the total number of vertices of T h in Ω h .Let K denote a generic element (triangle) of T h , and P j , j = 1, . . ., N nt be the vertices of T h Let P k be the space of polynomials of degree k.We consider the functional spaces defined by Let ϕ j , j = 1, . . ., N n be the finite element basis functions of V 1 0,h , associated with the triangulation T h .Problem ( 6) is therefore approximated by min where will detail in the sequel the specific method to impose the divergence free constraint at the discrete level (which is closely related to the choice of the finite element method for the solution of the generalized Stokes problem).Problem ( 7) is approximated by min where In both cases, and in the whole article, the integrals over Γ 1,h are computed via quadrature formula (trapezoidal rule).
We also define As in the continuous case, we define a discrete scalar product: We are looking for a saddle-point To compute the above saddle-point, we use the following discrete version of the algorithm ( 13)-( 15): Let u −1 h ∈ V i,h and λ 0 h ∈ (V 0 h ) 2×2 be arbitrary given functions.Then, for n ≥ 0: (a) Solve ) + λ n h .We denote by q i (resp.X n i,h ) the values of q (resp.X n h ) on the ith element K. Locally on each element K, this problem corresponds to min The minimum occurs when q i = α i X n i,h , α i ∈ R, α i ≥ 0. Solving the first order optimality conditions leads to: The computation of u n h will be discussed in Section 5.2.(c) Update the multipliers λ n+1 h ∈ (V 0 h ) 2×2 via: until convergence is reached.

5.2.
On the Solution of the Discrete Constrained Elliptic Problems.In order to solve the constrained minimization problem (32), we employ the discrete analogue of algorithm ( 18)-( 19) described below.With the Uzawa algorithm we employ for the solution of (32) reads as follows: Let ξ 0 h ∈ Λ h be given, then, for k ≥ 0: (i) Compute: where The solution of ( 35) is discussed for i = 0 (resp., i = 1) in Section 5.3 (resp., Section 5.4).(ii) Update ξ k h via where ρ > 0 is small enough, until convergence is achieved.It follows from ( 17) that ( 36) can be explicitly written as: , for all P vertex of T h belonging to Γ 1,h .

On the Solution of Problem (35) :
The Case without Incompressibility Constraint.The solution of ( 35) when i = 0 and v ∈ V 0,h is detailed here.A Lagrange multiplier χ k h ∈ R is introduced to handle the scalar equality constraint The first order optimality conditions read as follows: find u n,k h ∈ V 0,h and χ k h ∈ R such that, for all v ∈ V 0,h : Equation ( 37) is linear with respect to the multiplier χ k h .It suffices therefore to solve (37) for two values of χ k h , e.g.χ k h = 0 and χ k h = 1, the corresponding solutions being u n,k 0,h and u n,k 1,h , and then consider the linear combination that satisfies (38), namely, for i = 1, . . ., N n : The boundary integrals on Γ 1,h are again approximated by numerical quadrature (trapezoidal formula).

On the Solution of Problem (35) :
The Case with Incompressibility Constraint.If i = 1 and v ∈ V 1,h , the solution of (35) involves a discrete Stokes problem.A Lagrange multiplier χ k h ∈ R is introduced to handle the scalar equality constraint Γ 1,h ξ k h • u n,k h ds = 1, while another multiplier p n,k h ∈ V 1 h (a pressure-like multiplier) is introduced to take into account, in a weak sense, the divergence free constraint.
The solution of the linear constrained generalized Stokes problem relies on stabilized P 1 − P 1 finite elements [16,17], in order to use, as mentioned above, T h associated with piecewise linear approximations for both velocity and pressure.The problem therefore reads as follows: find where τ K = αh 2 K /r, α > 0 is a constant, and h K is the diameter of the triangle K [16,17].
The discrete stabilized Stokes-type system (40)( 41) is linear with respect to the multiplier χ k h .It suffices therefore to solve (40) for two values of χ k h , e.g.χ k h = 0 and χ k h = 1, the corresponding solutions being (u n,k 0,h , p n,k 0,h ) and (u n 1,h , p n,k 1,h ), and then consider the linear combination that satisfies (42), namely, for i = 1, . . ., N n : where f n,k h is the coefficient given by ( 39).
Remark 7. Two generalized Stokes problems are solved at each iteration of the Uzawa algorithm, and at each step of the augmented Lagrangian algorithm.The discretization with stabilized finite P 1 − P 1 elements introduces an additional error that affects the accuracy of the approximation of the generalized eigenvalue and eigenvectors [16,17].
Numerical results are given in the next sections for various domains/bodies, with and without the divergence free constraint.These results are generalizing those obtained in [6] for the scalar-valued case.

Numerical Results without the Divergence Free Constraint.
Load Capacity in Bars.We consider first the horizontal bar Ω = (0, 5) × (0, 1).Let a ∈ (0, 5) be a parameter.As in [6], we consider two cases for the boundary conditions, namely: (a) The bar is fixed on the left side, while a traction is applied at the extreme right side of the bar only.We have thus Γ 2 ⊂ ∂Ω, Γ 2 = ∅, Γ 2 being load and imposed displacement free; more precisely: The bar is fixed on the left side, while a traction is applied on the complement of the boundary part where essential boundary conditions are imposed.We have thus: In the scalar case studied in [6], these two different sets of boundary conditions lead to two different behaviors.In particular, the presence of Γ 2 = ∅ prevents the material from fracturing.For the vector-valued case, Figures 3 and 4 visualize the results for the first and second set of boundary conditions respectively, for various values of the parameter a.
In Figure 3, one eigenfunction u h related to the first eigenvalue δ 0,h is represented.The eigenfunction is zero in the region of the domain Ω delimited by the essential zero boundary conditions, and shows a constant twist on the other part of the domain.Results shown in Figure 4 are identical and show that the choice of the boundary conditions does not influence the shape of the minimizer, but only the magnitude of the displacement as shown in Figure 5.   Table 1 illustrates the value of the approximations δ 0,h of the smallest eigenvalue δ 0 for various values of the parameter a.We can see that the first eigenvalue is larger when Γ 2 = ∅, and eventually converges to the same value when the boundary conditions become identical (i.e. when a → 5).
Table 1.Value δ 0,h of the objective function in (29) for various sets of boundary conditions and h = 0.05.First row: case (a) As expected, the magnitude of u h is an increasing function of a in both cases and, for the same value of a (that is the same Γ 0 ) the magnitudes in case (a) are larger than those in case (b).This implies (see Table 1) that for the same Γ 0 , δ 0,h is larger if Γ 2 = ∅; this last property is not surprising, since in case (a) the surface traction is concentrated on a part Γ 1 of the boundary which, for the same a (that is Γ 0 ), is the subset of the Γ 1 in case (b).Various Two-Dimensional Domains.Various domains are used to illustrate the shape of the approximations u h of the eigenfunctions.Figure 6 shows the approximation u h of the first eigenfunction on the unit disk centered at the origin, for various choices of boundary conditions; Γ 0,h is chosen as various sets of segments or points on the boundary, and Γ 1,h = ∂Ω h \Γ 0,h .The value of the parameter is r = 10 to guarantee the convergence of the augmented Lagrangian method.The initial guess for the Uzawa algorithm is given by u .Concerning now the convergence of algorithm (35)(36) at each step of the augmented Lagrangian algorithm, we use u where ε ′ = 10 −4 is a given tolerance.The results in Figure 6 are obtained after 100 iterations of the augmented Lagrangian algorithm, with a maximum of 50 iterations for the embedded Uzawa algorithm; they provide approximations that satisfy (a) δ 0,h = 0.586365802 (b) δ 0,h = 0.0933761933 (c) δ 0,h = 0.675285667 Figure 7 shows the approximation u h of the first generalized eigenfunction associated with the (non-simply connected) ring-shaped domain The boundaries are defined as Γ The value of the parameter is r = 10 to guarantee the convergence of the augmented Lagrangian method.Figure 7 (left) visualizes the global optimum obtained when considering u −1 h (x, y) = (−y, x), while Figure 7 (right) visualizes the local optimum obtained when u −1 h (x, y) = (0, 0), and λ 0 h = 0 in both cases.The results obtained after 100 iterations of the augmented Lagrangian algorithm, with a maximum of 50 iterations for the embedded Uzawa algorithm, provide approximations that satisfy u n h − u n−1 h 2 < 10 −2 for the local minimum and u n h − u n−1 h 2 < 10 −5 for the global minimum.Finally, we consider the (non-convex) L-shaped domain Ω L = (0, 2) × (0, 1) ∪ (0, 1) × (0, 2). Figure 8 illustrates the approximation u h of the first generalized eigenfunction for two sets of boundary conditions (both satisfying Γ 0 ∪ Γ 1 = ∂Ω).The value of the parameter is r = 10.After 100 iterations of the augmented Lagrangian algorithm, the stopping criterion satisfies u n h − u n−1 h 2 < 10 −3 for both cases.The numerical results for the eigenfunction exhibit a rotating displacement that is similar to the δ 0,h = 0.196780193 δ 0,h = 0.386533885 solution for the rectangular domain, and the values of the objective function indicate that the material is fracturing for both sets of boundary conditions.In all cases, the non-smooth equality constraint is satisfied accurately.
Load Capacity in Bars.We consider now the numerical solution of problem (6), that is we assume from now on that the vector-valued function v in (6) verifies ∇ • v = 0. We consider again the rectangular domain shown in Figures 3 and 4.
Concerning the choice for Γ 0 , Γ 1 , and of the corresponding boundary conditions, we consider only the case (b) described in Section 6, that is For the case with incompressibility constraint, the augmented Lagrangian method identifies a bifurcation phenomenon between local and global minima.
Table 2 shows the variation, as a function of a, of the minimal value δ 1,h (global and local) of the objective function in (28).Table 2 shows that the bar does fracture, since the smallest eigenvalue satisfies δ 1,h < 1 in all cases.28) for various configurations of Γ 0 (h = 0.05).Case when Γ 0 = {(0, x 2 ) : x 2 ∈ (0, 1)} ∪ {(x 1 , 0) : x 1 ∈ (0, a)} ∪ {(x 1 , 1) : x 1 ∈ (0, a)} and Γ 1 = ∂Ω\Γ 0 , with a = 0, 1, 2, 3, 4, 5 (left to right, top to bottom).Various Two-Dimensional Domains.As in Section 6, domains of various shapes are considered, in order to investigate the pattern of the minimizers u h of the objective function in (28), and the formation of fractures.Following Section 6, we consider first the unit disk centered at the origin.With Γ 0 , Γ 1 , and the choice of the associated boundary conditions, identical to those in Figure 6, we have visualized on Figure 11 the corresponding minimizers and provided the associated value of δ 1,h .
We have taken r = 10 to guarantee the convergence of the augmented Lagrangian method, and used u −1 h = 0 and λ 0 h = 0 as initial guess.After 200 iterations, the augmented Lagrangian algorithm verifies u n h − u n−1 h 2 < 10 −2 , typically, assuming that the number of iterations of the embedded Uzawa algorithm is limited to 50.By comparing with Figure 6, one sees that the divergence free constraint enforces vortices.Note also that, as expected and in all cases, δ 1,h > δ 0,h .
Assuming that Ω is the ring-shaped domain S 0.5,1 , already considered in Section 6, and that Γ 0 and Γ 1 , and the associated boundary conditions are identical to those in Section 6, we visualize on Figure 12   Finally, as in Section 6, we consider the L-shaped domain Ω L = (0, 2) × (0, 1) ∪ (0, 1) × (0, 2).The sets Γ 0 , Γ 1 and the associated boundary conditions are identical to those in Section 6.On Figure 13, we have visualized the computed minimizers u h in (28) and the corresponding values of δ 1,h .The values of δ 1,h indicate that in both cases the material fractures.8. On the Convergence of the Approximated Generalized Eigenvalues δ 0,h and δ 1,h .Finally, we investigate, computationally, the order of convergence with respect to h of the generalized eigenvalues δ 0,h and δ 1,h (we have to keep in mind that δ 0 and δ 1 are not known, complicating the comparison with δ 0,h and δ 1,h ).We consider the horizontal bar Ω = (0, 5) × (0, 1).In order to avoid any mesh effects, the sets Γ 0 and Γ 1 are defined as follows: Γ ]}, and Γ 1 = ∂Ω\Γ 0 (as shown in Figure 15, top left).In that case, the discontinuity (break) line is not aligned with the mesh and mesh effects are avoided.A typical triangulation of Ω is shown in Figure 14; it corresponds to h = 0.05.On Table 3, we have reported the computed values of δ 0,h and δ 1,h , obtained for h = 1/10, 1/20 and 1/40; they clearly suggest convergence as h → 0 (observe also that, as expected, δ 0,h < δ 1,h (rather significantly, indeed)).Table 3. Convergence investigation of the approximate load capacity ratios δ 0,h and δ 1,h for the case of an horizontal bar with an oblique fracture line.h 0.10 0.05 0.025 δ 1,h 0.060845 0.059568 0.057253 δ 0,h 0.040369 0.038460 0.037720 On Figure 15 (resp., Figure 16), we have visualized, for h = 1/10, 1/20 and 1/40, the minimizers u h of the objective function in (29) (resp.( 28)).They correspond to the values of δ 0,h and δ 1,h reported in Table 3.Note that, even though the minimizers corresponding to δ 0,h and δ 1,h look alike and the equality constraints are satisfied, the value of the objective function satisfy (as expected), δ 0,h < δ 1,h .Table 3 and Figures 15 and 16 show that u h , δ 0,h and δ 1,h converge to some limit.However, to the best of our knowledge, the exact values for δ i , i = 0, 1 are neither known, nor easy to guess.Remark 8.In order to estimate the order of convergence of the errors |δ i − δ i,h |, i = 0, 1, one can estimate first the exact values δ 0 and δ 1 , from the values at h = 0 of a polynomial approximation δi (h) of δ i,h obtained by interpolation.Interpolating the data of Table 3 provides such polynomial approximations (we have thus δi ∈ P 2 , i = 0, 1).When using δ0 (0), resp.δ1 (0) as an exact value for δ 0 , resp.δ 1 , the orders of convergence for the errors |δ 1,h − δ 1 | and |δ 0,h − δ 0 | are between h and √ h, as shown by Figure 17.These results are consistent with those reported in [5,7] for related problems with known exact solutions.3).9. Conclusions.A non-smooth constrained optimization problem arising in stress analysis and related to the limit analysis of elastic materials has been investigated.The value of this infimum corresponds to the smallest generalized eigenvalue of a non-smooth operator.It determines if an elastic body with a given geometry fractures under external surface forces.
An augmented Lagrangian method, together with a piecewise linear finite element discretization, has been presented for the approximation of the minimizers.It relies on Uzawa-type iterative algorithms that take into account the various nonlinearities of the problem and decouple them from the treatment of the differential operators.The methodology discussed in the preceding sections has been validated by the results of numerical experiments associated with various two-dimensional geometries.
Figure 5 visualizes, when a = 2, the magnitude |u h | 2 of the generalized eigenfunction u h (|•| 2 being the canonical Euclidean norm in R 2 ) for the two cases (a) and (b).As expected, the magnitude of u h is an increasing function of a in both cases and, for the same value of a (that is the same Γ 0 ) the magnitudes in case (a) are larger than those in case (b).This implies (see Table1) that for the same Γ 0 , δ 0,h is larger if Γ 2 = ∅; this last property is not surprising, since in case (a) the surface traction is concentrated on a part Γ 1 of the boundary which, for the same a (that is Γ 0 ), is the subset of the Γ 1 in case (b).

Figure 5 .
Figure 5. Magnitude of the elastic displacement u h minimizer of (29) for the cases (a) and (b) (h = 0.05).

1 h 2 < 1 h 2 <
the global and local minimizers u h , and the corresponding values of δ 1,h .Similarly to the case without divergence free constraint, one observes a bifurcation between local and global optimum for the same set of initial conditions; the corresponding generalized eigenvalues and eigenfunctions are actually identical to the case without the incompressibility constraint.The value of the parameter is r = 10.The results obtained after 100 iterations of the augmented Lagrangian algorithm, with a maximum of 50 iterations for the embedded Uzawa algorithm, satisfy u n h − u n−10 −2 for the local minimum and u n h − u n−10 −5 for the global minimum.

Figure 15 .
Figure 15.Visualization of the elastic displacement u h , minimizer of(29), for an horizontal bar with an oblique discontinuity of the boundary conditions.Visualization of Γ 0 and Γ 1 (top left).Visualization of u h for h = 1/10, 1/20 and 1/40 (left to right, top to bottom).

Figure 16 .
Figure16.Visualization of the elastic displacement u h , minimizer of(28), for an horizontal bar with an oblique discontinuity of the boundary conditions.Visualization of Γ 0 and Γ 1 (top left).Visualization of u h for h = 1/10, 5/80 and 1/20 (left to right, top to bottom).

1 − 1 Slope 1 / 2 Figure 17 .
Figure 17.Order of convergence (log-log scales) for the errors |δ i,h − δ i |, i = 0, 1, for the case of an horizontal bar with an oblique break line (the exact value of δ i is obtained by extrapolation at h = 0 of the computed values reported in Table3).

Table 2 .
(28)e δ 1,h of the objective function in(28)for various sets of boundary conditions and h = 0.05.First column: local minimum; Second column: global minimum.