Terence Bayen

Olivier Cots

About the prior-saturation phenomenon for minimal time problems in the plane

We consider minimal time problems governed by control-affine-systems in the plane. We focus on the synthesis problem in presence of a singular locus that involves a saturation point for the singular control. We show that the minimum time synthesis can exhibit a prior-saturation point at the intersection of the singular locus and a switching curve. We also provide a set of non-linear equations to compute the prior-saturation point, and, at this point, we show a tangency property involving the switching curve.

I. INTRODUCTION

In this paper, we consider minimal time problems governed by control-affine-systems in the plane

ẋ = f (x) + u(t)g(x), |u(t)| ≤ 1,
where f, g : R 2 → R 2 are smooth vector fields. Syntheses for such problems have been investigated a lot in the literature [START_REF] Bonnard | Singular Trajectories and their role in Control Theory[END_REF], [START_REF] Boscain | Optimal Syntheses for Control Systems on 2-D Manifolds[END_REF], [START_REF] Piccoli | Classification of generic singularities for the planar timeoptimal synthesis[END_REF], [START_REF] Sussmann | The structure of time-optimal trajectories for singleinput systems in the plane : the general real analytic case[END_REF], [START_REF] Sussmann | The structure of time-optimal trajectories for singleinput systems in the plane: the C ∞ nonsingular case[END_REF], [START_REF] Sussmann | Regular synthesis for time-optimal control of singleinput real analytic systems in the plane[END_REF]. We focus here on the notion of singular arc which appears in the synthesis when the switching function (the scalar product between the adjoint vector and the controlled vector field g) vanishes over a time interval I. In that case, the corresponding singular control u s (which allows the corresponding trajectory to stay on the singular locus) can be expressed in feedback form x → u s (x). However, it may happen that u s becomes non admissible, i.e., x → |u s (x)| takes values above the maximal value for the control. Such a situation naturally appears in several application models, see, e.g., [START_REF] Bayen | Time-optimal control of concentrations changes in the chemostat with one single species[END_REF], [START_REF] Bayen | Analysis of an optimal control problem connected to bioprocesses involving a saturated singular arc[END_REF], [START_REF] Ledzewicz | Antiangiogenic therapy in cancer treatment as an optimal control problem[END_REF], [START_REF] Rapaport | Dynamical modelling and optimal control of landfills[END_REF]. In that case, we say that a saturation phenomenon occurs. This implies the following (non-intuitive) property that if a singular arc is optimal, then it should leave the singular locus at a so-called prior-saturation point before reaching the saturation point. This property has been studied in the literature in various situations such as for control-affine systems in dimension 2 and 4 (see, e.g., [START_REF] Schättler | A synthesis of time-optimal controls in the presence of saturated singular arcs[END_REF], [START_REF] Schättler | Geometric Optimal Control[END_REF], [START_REF] Bayen | Analysis of an optimal control problem connected to bioprocesses involving a saturated singular arc[END_REF] and references herein).

Our objective in this paper is to provide new qualitative properties on the minimum time synthesis in presence of a saturation point. More precisely, our objective is twofold:

• We first give a set of conditions on the system that ensure the occurence of the prior-saturation phenomenon, and we show that (under certain assumptions) the system leaves the singular arc at this point (before reaching the saturation point) with the maximal value for the control. This last arc is usually called bridge following the terminology as in [START_REF] Bonnard | Geometric optimal control of the contrast imaging problem in nuclear magnetic resonance[END_REF], [START_REF] Bonnard | Geometric and numerical methods in the contrast imaging problem in nuclear magnetic resonance[END_REF] (see also [START_REF] Bonnard | Time minimal synthesis for planar systems in the neighborhood of a terminal manifold of codimension one[END_REF], [START_REF] Bonnard | Singular Trajectories and their role in Control Theory[END_REF]).

• Second, we introduce a shooting function that allows an effective computation of the prior-saturation point.

In addition, we show that when the system exhibits a switching curve emanating from the prior-saturation point, then a tangency property occurs between this curve and the bridge. The tangency property has been pointed out in several application models (see, e.g., [START_REF] Bayen | Analysis of an optimal control problem connected to bioprocesses involving a saturated singular arc[END_REF], [START_REF] Bonnard | Geometric optimal control of the contrast imaging problem in nuclear magnetic resonance[END_REF]). To the best of our knowledge, this property has not been addressed previously in a general setting in the literature.

The paper is structured as follows: in Section II, we recall classical expressions and properties of singular controls for control-affine-systems in the plane introducing the saturation phenomenon. In Section III, we provide a set of conditions on the system that ensure occurrence of the priorsaturation phenomenon. In Section IV, we show a tangency property between the switching curve emanating from a prior-saturation point and the bridge, and we describe how to compute the prior-saturation point thanks to a shooting function constructed via the Hamiltonian lifts of f and g.

II. SATURATION PHENOMENON

The purpose of this section is to recall some facts about the minimum time control problems in the plane that will allow us to introduce the saturation phenomenon. Throughout the paper, | • | stands for the Euclidean norm in R 2 associated with the standard inner product written a • b for a, b ∈ R 2 , and a ⊥ denotes the vector a ⊥ := (-a 2 , a 1 ). The interior of a subset S ⊂ R n is denoted by Int(S).

A. Pontryagin's Principle

We start by applying the classical optimality conditions provided by the Pontryagin Maximum Principle (PMP), see [START_REF] Pontryagin | Mathematical theory of optimal processes[END_REF]. Let f, g : R 2 → R 2 be two vector fields of class C 2 , and consider the controlled dynamics:

ẋ = f (x) + u(t)g(x), (1)
with admissible controls in the set

U := {u : [0, +∞) → [-1, 1] ; u meas.}.
Given an initial point x 0 ∈ R 2 and a non-empty closed subset T ⊂ R 2 , we focus on the problem of driving (1) in minimal time from x 0 to the target set T :

inf u∈U T u s.t. x u (T u ) ∈ T , (2)
where x u (•) denotes the unique solution of (1) associated with the control u such that x u (0) = x 0 and T u ∈ [0, +∞] is the first entry time of x u (•) into the target set T . We suppose hereafter that optimal trajectories exist1 and we wish to apply the PMP on [START_REF] Bayen | Analysis of an optimal control problem connected to bioprocesses involving a saturated singular arc[END_REF]. The Hamiltonian associated with (2) is the function

H : R 2 × R 2 × R × R defined as H(x, p, p 0 , u) := p • f (x) + up • g(x) + p 0 .
If u is an optimal control and x u is the associated trajectory steering x 0 to the target set T in time T u ≥ 0, the following conditions are fulfilled:

• There exist p 0 ≤ 0 and an absolutely continuous function p : [0, T u ] → R2 satisfying the adjoint equation: almost everywhere over [0, T u ]:

ṗ(t) = -∇ x H(x u (t), p(t), p 0 , u(t)). (3)
• The pair (p 0 , p(•)) is non-zero.

• The optimal control u satisfies the Hamiltonian condition almost everywhere over [0, T u ]:

u(t) ∈ argmax ω∈[-1,1] H(x u (t), p(t), p 0 , ω). (4)
• At the terminal time, the transversality condition 1) and ( 3)-( 4) is abnormal whenever p 0 = 0 and normal whenever p 0 = 0 (in the latter case, we take p 0 = -1 and the corresponding extremal is denoted by (x u (•), p(•), u(•)) and we shall then write H(x, p, u) in place of H(x, p, p 0 , u)). Since T u is free and (1) is autonomous, the Hamiltonian H is zero along any extremal: for a.e. t ∈ [0, T u ],

2 p(T u ) ∈ -N T (x u (T u )) is fulfilled. Recall that an extremal (x u (•), p(•), p 0 , u(•)) satisfying (
H = p(t) • f (x u (t)) + u(t)p(t) • g(x u (t)) + p 0 = 0. (5)
The switching function φ is

φ(t) := p(t) • g(x u (t)), t ∈ [0, T u ], (6)
and it gives us the following control law:

   φ(t) > 0 ⇒ u(t) = +1, φ(t) < 0 ⇒ u(t) = -1, φ(t) = 0 ⇒ u(t) ∈ [-1, +1]. (7)
A switching time is an instant t c ∈ (0, T u ) such that the control u is discontinuous at time t c . We say that the corresponding extremal trajectory has a switching point at time t c . Of particular interest is the case when there is a time interval I := [t 1 , t 2 ] such that the switching function vanishes over I, i.e.,

φ(t) = p(t) • g(x u (t)) = 0, t ∈ I.
We then say that the extremal trajectory has a singular arc over I. Note that we shall suppose such an extremal to be normal, i.e., p 0 = 0. Indeed, recall from [

φ(t) = p(t) • [f, g](x u (t)), t ∈ [0, T u ],
where [f, g](x) is the Lie bracket of f and g at point x. The singular locus ∆ SA (in the state space) is then the (possibly empty) subset of R 2 defined as:

∆ SA := {x ∈ R 2 ; det(g(x), [f, g](x)) = 0}. (8)
For future reference, we set

δ SA (x) := det(g(x), [f, g](x)), x ∈ R 2 .
Note that if an extremal is singular over an interval [t

φ(t) = p(t) • [f, [f, g]](x u (t)) + u(t)p(t) • [g, [f, g]](x u (t)),
for a.e. t ∈ [0, T u ]. Therefore, u s becomes:

u s (t) := - p(t) • [f, [f, g]](x u (t)) p(t) • [g, [f, g]](x u (t)) , t ∈ [0, T u ], (9)
provided that p(t)

• [g, [f, g]](x u (t)) is non zero for t ∈ [t 1 , t 2 ]
. This expression of the singular control does not guarantee that u s is admissible, that is, u s (t) ∈ [-1, 1]:

• When u s (t) ∈ [-1, 1], the point x u (t) is said hyperbolic if p(t) • [g, [f, g]](x u (t)) > 0, and elliptic if p(t) • [g, [f, g]](
x u (t)) < 0 (see [START_REF] Bonnard | Time minimal synthesis for planar systems in the neighborhood of a terminal manifold of codimension one[END_REF]). • When |u s (t)| > 1 for some instant t, we say that a saturation phenomenon occurs and that the corresponding points of the singular locus are parabolic (see [START_REF] Bonnard | Time minimal synthesis for planar systems in the neighborhood of a terminal manifold of codimension one[END_REF]). Our purpose in what follows is precisely to investigate properties of the synthesis when saturation occurs.

B. Singular control and saturation phenomenon

In this part, we derive classical expressions of the singular control in terms of feedback that will allow us to introduce saturation points (in terms of the data defining the system). The collinearity set associated with (1) is the (possibly empty) subset of R 2 defined as

∆ 0 := {x ∈ R 2 ; det(f (x), g(x)) = 0}. ( 10
)
Consider now the functions δ 0 , ψ : R 2 → R defined respectively by δ 0 (x) := det(f (x), g(x)) and

ψ(x) := - det(g(x), [f, [f, g]](x)) det(g(x), [g, [f, g]](x)) . ( 11
) Lemma 2.1: Suppose that ∆ SA = ∅, that x → det(g(x), [g, [f, g]](x)
) is non-zero over ∆ SA , and consider a singular arc defined over an interval [t 1 , t 2 ]. Then, one has:

u s (t) = ψ(x(t)), t ∈ [t 1 , t 2 ]. (12)
where x(•) is the corresponding singular trajectory such that

x(t) ∈ ∆ SA for t ∈ [t 1 , t 2 ].
Proof: Along I, the adjoint vector can be expressed as p

(t) = -g ⊥ (x(t)) δ0(x(t)) , t ∈ [t 1 , t 2 ]. Since for t ∈ [t 1 , t 2 ], {f (x(t)), g(x(t))} is a basis of R 2 , we deduce (decomposing [f, [f, g]](x(t)) and [g, [f, g]](x(t)) on this basis): p(t) • [f, [f, g]](x(t)) = det(g(x(t)), [f, [f, g]](x(t)))Λ(x(t)), p(t) • [g, [f, g]](x(t)) = det(g(x(t)), [g, [f, g]](x(t)))Λ(x(t)), where Λ(x) := g ⊥ (x)•f (x) δ0(x) 2 , x / ∈ ∆ 0 . Moreover, g ⊥ (x(t)) • f (x(t)) = -δ 0 (x(t)
) and thus, this scalar product is nonzero because p 0 = 0. This ends the proof.

Remark 2.1: Steady-state singular points are defined as the points x ∈ ∆ SA ∩ ∆ 0 such that g(x ) = 0, see [START_REF] Boscain | Optimal Syntheses for Control Systems on 2-D Manifolds[END_REF] (if ∆ SA ∩ ∆ 0 = ∅). Such points are equilibria of ( 1) with u = ψ(x). A singular arc defined over a time interval [t 1 , t 2 ] does not contain such a point because f (x(t)) and g(x(t)) must be linearly independent over

[t 1 , t 2 ]. But, it can contain points x ∈ ∆ 0 ∩ ∆ SA such that g(x ) = 0.
To introduce the notion of saturation point, it is convenient to consider a parametrization of ∆ SA as follows. When ∆ SA ∩ ∆ 0 is non-empty, ∆ SA \∆ 0 may consist of several parts (or components), and we write this set as

∆ SA \∆ 0 = j∈K γ j , where K is an index set. Lemma 2.2: Suppose that ∆ SA is non-empty and that x → det(g(x), [g, [f, g]](x)) is non-zero over ∆ SA . Then, each component γ of ∆ SA \∆ 0 can be parametrized by a one-to-one parametrization ζ : J → γ, τ → ζ(τ ) of class C 1 , where J is an interval of R. Proof: For x ∈ ∆ SA \∆ 0 , one has span{f (x), g(x)} = R 2 . Hence, there exist α(x), β(x) ∈ R such that [f, g](x) = α(x)f (x) + β(x)g(x). (13)
By taking the determinant, we find that for x ∈ ∆ SA \∆ 0 ,

α(x) = - det(g(x), [f, g](x)) δ 0 (x) , β(x) = det(f (x), [f, g](x)) δ 0 (x) . Consider now a component γ of ∆ SA \∆ 0 and x ∈ γ. By computing [f, [f, g]](x) thanks to (13), we get det(g(x), [g, [f, g]](x)) = -δ 0 (x)∇α(x) • g(x), x ∈ γ. Since x → det(g(x), [g, [f, g]](x)
) is non-zero over ∆ SA , the preceding equality implies that the scalar product ∇α(x)• g(x) is non-zero. On the other hand, γ is defined by the implicit equation δ SA (x) = 0. Observe that for x / ∈ ∆ 0 , δ SA (x) = -α(x)δ 0 (x). By taking the derivative, we find that for x / ∈ ∆ 0 , one has ∇δ SA (x) = -δ 0 (x)∇α(x)α(x)∇δ 0 (x). Therefore, for x ∈ γ, we obtain ∇δ SA (x) = -δ 0 (x)∇α(x). We can conclude that for any point x ∈ γ, the derivative ∂ 1 α(x) or ∂ 2 α(x) is non-zero. We are then in a position to apply the implicit function theorem to δ SA locally at each point x ∈ γ, which then implies the desired property.

Under the assumptions of Lemma 2.2, given a component γ of ∆ SA , there is a parametrization ζ such that

γ := {ζ(τ ) ; τ ∈ J}, where ζ : J → R 2 is C 1 -mapping (injective) and J is an interval. Definition 2.1: A point x * := ζ(τ * ) with τ * ∈ Int(J) is called saturation point if ψ(x * ) = 1, ψ(ζ(τ )) ∈ (-1, 1) for any τ ∈ J such that τ < τ * , and ψ(ζ(τ )) > 1 for any τ ∈ J such that τ > τ * .
As well, we can define saturation points x such that ψ(x ) = -1, that is, when the lower bound of the admissible control set is saturated. Our next aim is to study the optimality of singular arcs in presence of a saturation point.

III. EXISTENCE OF A PRIOR-SATURATION POINT

In this section, we show that a prior-saturation phenomenon can occur whenever the system exhibits a saturation point. We start by introducing our main assumptions. Assumption 3.1: System (1) satisfies the following:

(i) One has ∆ 0 = ∅ and δ 0 (x) < 0 for all x ∈ R 2 . (ii)
The set ∆ SA is non-empty, simply connected, and has exactly one saturation point x * with ψ(x * ) = 1. (iii) Along the singular locus, the strict Legendre-Clebsch optimality condition is satisfied, that is, any singular extremal

(x u (•), p(•), u(•)) defined over [t 1 , t 2 ] satisfies: ∂ ∂u d 2 H u dt 2 (x u (t), p(t), u(t)) > 0, ∀t ∈ [t 1 , t 2 ]. (14) (iv) If Γ -
is the forward semi-orbit of (1) with u = -1 with the initial condition x * at time 0, then

T ∩ Γ -= ∅. (15)
(v) The target T is reachable from every point x 0 ∈ R 2 .

Remark 3.1: (i) The hypothesis ∆ 0 = ∅ is not restrictive since we could restrict our analysis to a component γ of ∆ SA in place of ∆ SA . (ii) By the previous computations, we can observe that ( 14) is equivalent to

det(g(x), [g, [f, g]](x)) > 0, ∀x ∈ ∆ SA .
Recall that, under the strict Legendre-Clebsch condition, the singular arc is a turnpike, i.e., it is time-minimizing in every neighborhood of a hyperbolic point of ∆ SA , [START_REF] Bonnard | Time minimal synthesis for planar systems in the neighborhood of a terminal manifold of codimension one[END_REF]. This property can be retrieved by the clock form argument [START_REF] Hermes | Functional analysis and time optimal control[END_REF].

Under Assumption 3.1, the singular locus ∆ SA is written ∆ SA := ζ(J) where J ⊂ R is an interval and ζ : J → ∆ SA is a C 1 -mapping. In addition, ∆ SA partitions the state space into two simply connected subsets ∆ ± SA :

∆ + SA := {x ∈ R 2 ; det(g(x), [f, g](x)) > 0}, ∆ - SA := {x ∈ R 2 ; det(g(x), [f, g](x)) < 0}. Given a normal extremal (x u (•), p(•), u(•)), the function t → γ u (t) := β(x u (t)) -α(x u (t))u(t), t ∈ [0, T u ],
is well-defined since ∆ 0 = ∅. Lemma 3.1: Along a normal extremal (x u (•), p(•), u(•)), the switching function φ satisfies the ODE

φ(t) = γ u (t)φ(t) + α(x u (t)) a.e. t ∈ [0, T u ].
(16) Proof: The proof follows using the expression of φ and the fact that the Hamiltonian H is constant equal to zero.

The next proposition shows that an extremal containing a singular arc until the point x is not optimal. Proposition 3.1: Suppose that Assumption 3.1 holds true, and consider an optimal trajectory steering x 0 to the target T in time T u . Then, the corresponding extremal

(x u (•), p(•), u(•)) does not contain a singular arc defined over a time interval [t 1 , t 2 ] such that x u (t 2 ) = x * .
Proof: Suppose by contradiction that there is a time interval [t 1 , t 2 ] such that the trajectory is singular over [t 1 , t 2 ] and such that x u (t 2 ) = x * . We claim that, at time t 2 , the vector f (x u (t 2 ))+g(x u (t 2 )) is tangent to ∆ SA . Indeed, it is enough to check that the vector f (x * ) + g(x * ) is orthogonal to ∇δ SA (x * ) = -δ 0 (x * )∇α(x * ). As we have seen in the proof of Lemma 2.2, one has for x ∈ ∆ SA :

det(g(x), [g, [f, g]](x)) = -δ 0 (x)∇α(x) • g(x), det(g(x), [f, [f, g]](x)) = -δ 0 (x)∇α(x) • f (x). (17)
These equalities imply that

-δ 0 (x * )∇α(x * ) • (f (x * ) + g(x * )) = det(g(x * ), [g, [f, g]](x * )) + det(g(x * ), [f, [f, g]](x * )).
Since ψ(x * ) = 1, the right member of the above equality is zero which shows the claim. Note also that for x ∈ ∆ SA , (17) also implies the equalities

∇δ SA (x) • (f (x) + g(x)) = det(g(x), [g, [f, g]](x))(1 -ψ(x)), ∇δ SA (x) • (f (x) -g(x)) = det(g(x), [g, [f, g]](x))(-1 -ψ(x)).
Consider now the unique solution x -of (1) with u = -1 starting from x * at time t 2 . This trajectory enters into the set ∆ - SA for t > t 2 , t close to t 2 , because one has ∇δ SA (x * ) • (f (x * ) -g(x * )) < 0. Going back to the optimal trajectory, there are now two possibilities for x u (•): for t > t 2 , t close to t 2 , either x u (•) enters into ∆ + SA or into ∆ - SA (because the singular control becomes non admissible).

Suppose first that x u (•) enters into ∆ + SA . Then, there is ε > 0 such that one has α(x u (t)) > 0 for t ∈ (t 2 , t 2 + ε]. It follows from ( 16) that one has u = +1 on this interval. But the velocity set being convex, we obtain a contradiction with the non-admissibility of the singular control at x * (because x -enters into ∆ - SA ). It follows that the optimal trajectory necessarily enters into the set ∆ - SA . But then, since α < 0 in ∆ - SA , [START_REF] Sussmann | The structure of time-optimal trajectories for singleinput systems in the plane : the general real analytic case[END_REF] implies that u = -1 in some time interval (t 2 , t 2 + ε].

From Assumption 3.1, the forward semi-orbit with u = -1 starting from x * does not reach the target set. Hence, x u (•) must have a switching point to u = +1 in ∆ - SA or it must reach ∆ SA with the control u = -1. We see from ( 16) that the first case is not possible because at a switching time t c such that x u (t c ) ∈ ∆ - SA , we would have φ(t c ) ≥ 0 in contradiction with α(x u (t c )) < 0.

Suppose now that x u (•) reaches ∆ SA at some point x := ζ(τ ) with τ < τ * . Then, we obtain ∇δ SA (x) • (f (x)g(x)) < 0 since ψ(x) > -1. But, this contradicts the fact that x u (•) reaches ∆ SA with u = -1 at point x (indeed, because at this point, the singular control is admissible, one must have ∇δ SA (x)•(f (x)-g(x)) ≥ 0). In the same way, the trajectory cannot reach a point x ∈ ∆ SA such that x = ζ(τ ) with τ > τ * .

We can conclude that for any time t ≥ t 2 , one has u(t) = -1, but then, the optimal trajectory cannot reach the target set which is a contradiction (Assumption 3.1 (iv)). This concludes the proof.

As an example, if x 0 belongs to the singular locus and is such that x 0 := ζ(τ 0 ) with τ 0 < τ * , and if in addition, the optimal trajectory starting from x 0 contains a singular arc, then the trajectory should leave the singular locus before reaching x * . Let us insist on the fact that this property of leaving the singular locus before reaching x * relies on the fact that the optimal trajectory should contain a singular arc. In the fed-batch model presented in [START_REF] Bayen | Analysis of an optimal control problem connected to bioprocesses involving a saturated singular arc[END_REF], this property can be easily verified.

We now introduce the following definition (in line with [START_REF] Ledzewicz | Antiangiogenic therapy in cancer treatment as an optimal control problem[END_REF], [START_REF] Schättler | A synthesis of time-optimal controls in the presence of saturated singular arcs[END_REF], [START_REF] Schättler | Geometric Optimal Control[END_REF]). Hereafter, the notation S [τ 0 ,τ0] denotes a singular arc passing through the points ζ(τ 0 ) and ζ(τ 0 ) with τ 0 ≤ τ 0 < τ * . Definition 3.1: Let τ 0 < τ * . A point x e := ζ(τ e ) ∈ ∆ SA with τ 0 < τ e < τ * is called a prior-saturation point if the singular arc S [τ0,τ ] ceases to be optimal for τ ≥ τ e . This definition makes sense only for initial conditions ζ(τ 0 ) with τ 0 < τ * because for τ 0 ≥ τ * , optimal controls are not singular (since the singular control is non-admissible). We highlight the dependency of x e w.r.t. the initial condition ζ(τ 0 ) ∈ ∆ SA as follows.

Moreover, for any τ 0 ∈ (τ e , τ * ] an optimal trajectory starting at ζ(τ 0 ) leaves the singular locus at ζ(τ 0 ). Proof: Let E := {τ ∈ J ; S [τ1,τ ] is optimal} and F := {τ ∈ J ; S [τ0,τ ] is optimal} where τ 0 ∈ [τ 1 , τ 2 ) is fixed. Take a point τ ∈ F . Then, from our assumption, S [τ1,τ ] is also optimal (by concatenation) which shows that τ ∈ E. On the other hand, if τ ∈ E, then S [τ0,τ ] is optimal (as a sub-arc). It follows that E = F and, in addition, since x e is defined as the point such that S [τ1,τ ] ceases to be optimal, we obtain [START_REF] Sussmann | Regular synthesis for time-optimal control of singleinput real analytic systems in the plane[END_REF].

Finally, for every τ 0 ∈ (τ e , τ * ), a singular arc S [τ0,τ 0 ] with τ 0 < τ 0 < τ * cannot be optimal, since otherwise, this would contradict the definition of τ e . It follows that for every τ 0 ∈ (τ e , τ * ], no singular arc occurs and thus x e = ζ(τ 0 ). This property implies in particular that for some initial conditions in ∆ SA (e.g., for x 0 := ζ(τ 1 )), then the optimal path has a singular arc until x e and a switching point at the prior-saturation point. Remark 3.2: In addition to Assumption 3.1 (in particular [START_REF] Schättler | Geometric Optimal Control[END_REF]), if T is not reachable with the constant control u = -1 from those points of ∆ SA located between x e and x * , then the maximal value for the control u = +1 is locally optimal from x e . In other words, the bridge (the last arc leaving ∆ SA ) corresponds to u = +1. This can be proved by using similar arguments as for proving Proposition 3.2. Since the singular arc is a turnpike, this additionnal hypothesis also implies the existence of a switching curve emanating from x e . Our next aim is precisely to investigate more into details geometric properties of optimal paths at the point x e .

IV. TANGENCY PROPERTY AND PRIOR-SATURATION

A. Introduction to the tangency property 1) Prior-saturation lift: To introduce this concept, let us start with an example. Let us consider a target T := {x f }, x f ∈ R 2 , with an optimal trajectory of the form γ -γ s γ + , where γ -, γ + and γ s are arcs, respectively, with control u = -1, u = +1 and u = u s , where u s is the singular control. The PMP gives first order optimality conditions satisfied by this extremal trajectory, that we can write as a system of nonlinear equations, the so-called shooting equations. To introduce this set of equations, we introduce some notation: we define H f (z) := p • f (x) and H g (z) := p • g(x), z := (x, p), the Hamiltonian lifts of f and g. Any other Hamiltonian lift is defined like this. Define also the Hamiltonians H ± := H f ± H g and H s := H f + u s H g , where u s is viewed as a function of z:

u s (z) := - p • [f, [f, g]](x) p • [g, [f, g]](x) = - H [f,[f,g]] (z) H [g,[f,g]] (z) .
For any Hamiltonian H we define the Hamiltonian system #-H := (∂ p H, -∂ x H), and finally, we introduce the exponential mapping, such that the solution at time t of any differential equation ż(s) = ϕ(z(s)), with initial condition z(0) = z 0 , is written e tϕ (z 0 ) or exp(tϕ)(z 0 ). The shooting equations are then defined by the equation

S(y) = 0, y := (p 0 , t 1 , t 2 , t f , z 1 , z 2 ) ∈ R n+3+(2n)×2 , n = 2, with S(y) :=         H g (z 1 ) H [f,g] (z 1 ) H + (exp((t f -t 2 ) # - H + )(z 2 )) + p 0 π x (exp((t f -t 2 ) # - H + )(z 2 )) -x f z 1 -exp(t 1 # - H -)(x 0 , p 0 ) z 2 -exp((t 2 -t 1 ) # - H s )(z 1 )        
, where π x (x, p) = x, and x 0 ∈ R 2 , p 0 ≤ 0 are given. The two first equations mean that the trajectory is entering the singular locus at z 1 . Hence, the second arc is a singular arc. The third equation takes into account the free terminal time.

The fourth equation implies that the last bang arc reaches the target T := {x f } at the final time t f , and the last two equations are the so-called matching conditions (which are not required but improve numerical stability). Now, let y * := (p 0 , t 1 , t 2 , t f , z 1 , z 2 ) be a solution to S(y) = 0. Then, assuming t 2 > 0, the point π x (z 2 ) is a prior-saturation point and we call z 2 a prior-saturation lift.

2) The prior-saturation lift is locally unique: Let us consider a smooth and local one-parameter family of initial conditions x 0 (α), α ∈ (-ε, ε), ε > 0. Let us assume that for any α, the optimal trajectory is of the form γ -γ s γ + , and we denote by y(α) the corresponding solution of S α (y), where S α is defined as in section IV-A.1, but with initial condition x 0 (α) in place of x 0 . In addition, suppose that the lengths t 1 , t 2 -t 1 and t f -t 2 are positive, that is, each arc is defined on an interval with non-empty interior. Under this setting, for any α, we have z 2 (0) = z 2 (α), that is, the prior-saturation lift is locally unique. Hence, we can define z e := z 2 (0) the prior-saturation lift. In relation with the property 3.2, we have π x (z e ) = x e , where x e is the prior-saturation point. See the following illustration in the state space:

γ -γ - γ s γ + x 0 (0)
x 0 (α)

x e ∆ SA
3) Computation of the prior-saturation lift: Actually, in the previous example, we can compute the prior-saturation lift with a smaller set of equations. The prior-saturation lift is the switching between the singular and positive bang arcs. What happens before is useless and we also do not need the matching equations. With these considerations we can set

F (t b , z b ) :=     H g (exp(-t b # - H + )(z b )) H [f,g] (exp(-t b # - H + )(z b )) H + (z b ) + p 0 π x (z b ) -x f     ,
where F : R 5 → R 5 , and where we use the notation t b , z b (b stands for bridge) in relation with the concept of bridge defined in Ref. [START_REF] Bonnard | Geometric optimal control of the contrast imaging problem in nuclear magnetic resonance[END_REF]. Note that the exponential mapping is here computed by backward integration, and with the notation of the two previous parts: we have z b = exp((t f -t 2 )

# -H + )(z 2 ) and t b = t f -t 2 . Hence, the prior-saturation lift is simply given by z e = exp(-t b # -H + )(z b ). From a general point of view, we can assume that the prior-saturation lift is given by solving a set of nonlinear equations of the following form:

F (t b , z b , λ) := H [f,g] (exp(-t b # - H + )(z b )) G(t b , z b , λ) , (19)
where λ ∈ R k , k ∈ N, represents some parameters, F : R 5+k → R 5+k , and G is defined by

G(t b , z b , λ) :=   H g (exp(-t b # - H + )(z b )) H + (z b ) + p 0 Ψ(z b , λ)   , (20)
with Ψ(z b , λ) ∈ R 2+k . It is important to notice that Ψ does not depend on t b . In the previous example, we would have Ψ(z b ) = π x (z b ) -x f , that is, k = 0. For a more complex structure of the form γ -γ s γ + γ -, the parameter λ would be the last switching time between the γ + and γ -arcs. In this case, Ψ would contain the additional switching condition H g = 0 at this time.

4) The tangency property: Let us come back to the simple example where the solution is of the form γ -γ s γ + . Let us consider a smooth and local one-parameter family of initial conditions x 0 (α), α ∈ (-ε, ε), ε > 0. Let us also assume that for α = 0, the optimal solution is of the form γ -γ s γ + , but with γ s reduced to a single point, that is, t 2 (0) -t 1 (0) = 0, with y(0) := (p 0 (0), t 1 (0), t 2 (0), t f (0), z 1 (0), z 2 (0)) the solution of the associated shooting equations. Assume also that for α > 0, we are in the case of Sections IV-A.1 -IV-A.2, that is, t 2 (α) -t 1 (α) > 0, with y(α) := (p 0 (α), t 1 (α), t 2 (α), t f (α), z 1 (α), z 2 (α)) the solution of the associated shooting equations. The prior-saturation lift is given by z e = z 2 (α) for α ∈ [0, ε). Let Γ + denotes the forward semi-orbit of ż = # -H + (z) starting from z e . The idea is now to consider the case where there is a bifurcation in the structure of optimal trajectories when α = 0. We thus assume that for α ∈ (-ε, 0), the solutions are of the form γ -γ + and we denote by z 1 (α) the switching point (in the cotangent space) between the two arcs. In this setting, there is a birth of a switching locus (still in the cotangent space) defined by Σ -:= {z 1 (α); α ∈ (-ε, 0)}.

The aim of the next section is to prove that the semiorbit Γ + is tangent to the switching curve Σ -at the priorsaturation lift z e . This is the tangency property. See the following illustration where Σ s : H g = H [f,g] = 0 is the singular locus and where we define z 0 (α) := (x 0 (α), p 0 (α)):

Γ + z 0 (0) z 1 (α) z 0 (α), α < 0 z e Σ s Σ -

B. Results

We start by a general definition of a prior saturation lift. Next, we introduce the following assumptions.

F (t * b , z * b , λ * ) =   -a * * -b ∂G ∂z b (t * b , z * b , λ * ) ∂G ∂λ (t * b , z * b , λ * )   , where b = (H [f,g] (z e ), 0, 0) = 0 since F (t * b , z * b , λ * ) = 0 and a = H [f,[f,g]] (z e ) + H [g,[f,g]] (z e ) = 0 since u s (z e ) < 1.
Note that the point z e is locally unique by the inverse function theorem. From Assumption 3.1 and by Proposition 3.2, there exists a prior-saturation point x e that is locally unique, and from Assumption 4.2, x e has a locally unique lift in the cotangent space, given by the solutions to F = 0.

Lemma 4.1: Suppose that Assumption 4.2 holds true. Then, there is a curve, solution to G = 0, given by the graph of an implicit function

t b → σ(t b ) := (z b (t b ), λ(t b )) ∈ R 4+k satisfying σ(t * b ) = (z * b , λ * ), defined over an interval of the form (t * b -ε, t * b + ε), ε > 0, and such that 3 σ (t * b ) = 0 R 4+k . Proof:
The existence of σ is a simple application of the implicit function theorem. Its derivative is given by:

σ (t b ) = - ∂G ∂z b [t b ] ∂G ∂λ [t b ] -1 • ∂G ∂t b [t b ].
where [t b ] stands for (t b , σ(t b )). Since # -H + (z) starting from z e . The curve Σ is called a switching curve since H g (ϕ(t b )) = 0. However, the switching curve is not necessarily optimal, that is, the optimal synthesis, with respect to the initial condition, may not contain Σ. Let us stratify Σ according to Σ = Σ -∪ Σ 0 ∪ Σ + , with A typical situation is the following: Σ -∪ Σ 0 is contained in the optimal synthesis while Σ + is not optimal for local and/or global optimality reasons. Our main result is the following.

Theorem 4.1: Suppose that Assumptions 3.1 and 4.2 hold true. Then, there exist a prior-saturation lift z e and a switching curve Σ which is tangent to the semi-orbit Γ + at z e .

Proof: From Assumption 3.1, there exists a priorsaturation point. Since the system F from eq. ( 19) is assumed to be obtained from the PMP, the prior-saturation lift z e := exp(-t * b # -H + )(z * b ) is well defined. From Assumption 4.2, z e is locally unique (since u s (z e ) < 1) and we can define Σ and Γ + according to Definition 4.2. To prove the tangency property, we have to show that ϕ (t * b ) is colinear to Proof: Since ξ is a submersion at z e , then Σ s is locally a regular submanifold of codimension two near z e . Its tangent space at z e is given by the kernel of ξ (z e ). But V. CONCLUSION In presence of a saturation phenomenon of the singular control, the tangency property is a useful information for the computation of an optimal synthesis and optimal paths near the prior-saturation point. It also appears in other settings such as in Lagrange control problems governed by onedimensional systems, see, e.g., [START_REF] Kalbussi | Optimal control of membrane filtration systems[END_REF]. Future works could then investigate the prior-saturation phenomenon and the tangency property in other frameworks or in dimension n ≥ 3.

Proposition 3 . 2 :

 32 Suppose that Assumption 3.1 holds true and that there are τ 1 < τ 2 < τ * such that any optimal trajectory starting from ζ(τ 0 ) with τ 0 ∈ [τ 1 , τ 2 ) contains a singular arc S [τ0,τ2] . Then, for any initial condition τ 0 ∈ [τ 1 , τ 2 ), one has x e = ζ(τ e ) with τ e := sup{τ ∈ J and S [τ1,τ ] is optimal}.

  Definition 4.1: Let (t * b , z * b , λ * ) ∈ R 5+k be a solution to the equation F = 0 (recall (19)) and define z e := exp(-t * b # -H + )(z * b ) ∈ Σ s . The point z e is called a priorsaturation lift if π x (z e ) is a prior-saturation point.

  b ] = (H [f,g] (z e ), 0 R 3+k ) = 0 R 4+k , the result follows. Definition 4.2: We define the switching curve Σ := {ϕ(t b ); t b ∈ (t * b -ε, t * b + ε)} with ϕ(t b ) := exp(-t b # -H + )(z b (t b )), and Γ + as the forward semi-orbit of ż =

Σ

  -:= {ϕ(t b ); t b ∈ (t * b -ε, 0)}, Σ 0 := {ϕ(t * b )} = {z e }, Σ + := {ϕ(t b ); t b ∈ (0, t * b + ε)}.

  # -H + (z e ).For anyt b ∈ (t * b -ε, t * b + ε), we have ϕ (t b ) = -# -H + (ϕ(t b )) + Φ(t b , z b (t b )) • z b (t b ),where Φ(t, z 0 ) is defined as the solution at time t of the Cauchy problem Ẋ(s) = A(s, z 0 )X(s), X(0) = I, withA(s, z 0 ) := -# -H + (exp(-s # -H + )(z 0 )). By lemma 4.1, z b (t * b ) = 0 whence the result. Let ξ(z) := (H g (z), H [f,g] (z)), then Σ s = ξ -1 ({0 R 2 }).Proposition 4.2: Suppose that ξ is a submersion at z e and that Assumptions 4.1 and 4.2 hold true. Then the switching curve Σ is transverse to the singular locus Σ s at z e .

ξ

  (z e ) • ϕ (t * b ) = -ξ (z e ) • # -H + (z e ) = -(H [f,g] (z e ), H [f,[f,g]] (z e ) + H [g,[f,g]] (z e )) = (0, 0),since z e is prior-saturation lift.

  1 , t 2 ], then x u (t) ∈ ∆ SA for any t ∈ [t 1 , t 2 ] because p(•) must be non-zero and orthogonal to span{g(x u (t)), [f, g](x u (t))} for t ∈ [t 1 , t 2 ]. The singular control u s is then the value of the control for which the trajectory stays on the singular locus ∆ SA . Differentiating φ w.r.t. t gives:

  Assumption 4.1: The point z e satisfies u s (z e ) < 1.

	Assumption 4.2: The function G from eq. (20) satisfies:
	∂G ∂z b	(t * b , z * b , λ * )	∂G ∂λ	(t * b , z * b , λ * )	invertible.
	Proposition 4.1: Suppose Assumptions 4.1 and 4.2 hold
	true, then F (t * b , z * b , λ * ) is invertible. Proof: The Jacobian of F at (t * b , z * b , λ

* ) is given by:

If the target can be reached from x 0 and if f, g have linear growth, then (2) admits an optimal solution, thanks to Filippov's Existence Theorem.

Here, N T (x) stands for the (Mordukovitch) limiting normal cone to T at point x ∈ T , which coincides with the normal cone in the sense of convex analysis when T is convex, see[START_REF] Vinter | Optimal Control, Systems and Control: Foundations and Applications[END_REF].

We assume that all the functions are smooth enough.

VI. ACKNOWLEDGMENTS

We are very grateful to E. Trélat for helpful discussions and suggestions about the tangency property at the priorsaturation point.


Lev Truskinovsky

Anna Vainchtein

Quasicontinuum modeling of short-wave instabilities in crystal lattices

Keywords: Crystal lattice, Short-wave instability, Quasicontinuum model

We propose a hybrid quasicontinuum model which captures both long and short-wave instabilities of crystal lattices and combines the advantages of weakly nonlocal (higher gradient) and strongly nonlocal (integral) continuum models. To illustrate the idea, we study the simplest onedimensional lattice exhibiting commensurate and incommensurate short-wave instabilities. We explicitly compute stability limits of the homogeneous states using both discrete and quasicontinuum models. The new quasicontinuum approximation is shown to be capable of reproducing a detailed structure of the discrete stability diagram.

Introduction

Lattice instabilities are responsible for various transformations of crystal structures. Macroscopic or long-wave instabilities give rise to martensitic phase transformations, which proceed through the formation of finite-size domains of the new phase [START_REF] Pitteri | Continuum theories for phase transitions and twinning in crystals[END_REF][START_REF] Toledano | Reconstructive phase transitions: in crystals and quasicrystals[END_REF]. Microscopic or short-wave instabilities lead to the formation of multi-lattices and modulated "tweed" patterns [START_REF] Toledano | Reconstructive phase transitions: in crystals and quasicrystals[END_REF][START_REF] Castan | [END_REF]4]. Both micro and macroinstabilities are detectable in the phonon dispersion spectra. The loss of stability takes place when when the minimum of the dispersion curve touches the zero frequency level [5,[START_REF] Wallace | Thermodynamics of crystals[END_REF][START_REF] Triantafyllidis | [END_REF]. Macroscopic instabilities correspond to infinitely long waves and can be linked to the softening of the appropriate combinations of macroscopic elastic moduli. Typical microscopic instabilities occur at finite wave lengths associated with special points in the Brillouin zone corresponding to modulations commensurate with the lattice. Less frequent microinstabilities with the generic wave vectors give rise to incommensurate phases [8].

Classical continuum description of lattice instabilities is provided by the Landau theory. When the unstable wave vector is equal to zero, the order parameter can be chosen to coincide with a component of the macroscopic strain (e.g. [9]). If the unstable wave vector is different from zero, the macroscopic order parameter can be identified with the amplitude of an unstable normal mode (e.g. [10]). The resulting coarse-grained description remains adequate until external length scales become comparable with the lattice spacing as in the cases of ultra-thin layers, boundaries with singularities or strongly interacting defects.

To preserve the fine structure of the elastic fields in the continuum setting we propose a new quasicontinuum model which captures both long and short-wave instabilities associated with the acoustic branch of the phonon spectrum. The model is obtained by a long-wave expansion of the lattice energy and it shares with the discrete model the property that strain is the only order parameter. To make the truncated gradient expansion well-posed, we extend the polynomial dispersion relation by zero outside the first Brillouin zone. The model therefore restricts modulations in the physical space to length scales larger than lattice spacing. The suppression of the gradient model outside the Brillouin zone allows one to avoid the unphysical small-scale instabilities; a similar hard cutoff procedure is routinely used in coarse-grained continuum field theories [START_REF] Chaikin | Principles of condensed matter physics[END_REF]. As we show, the cutoff makes the gradient model equivalent to a fully nonlocal model with a slowly decaying oscillatory kernel.

To illustrate the idea, we consider a one-dimensional lattice involving interactions of up to third nearest neighbors which can be either of ferromagnetic or antiferromagnetic type. The discrete model is shown to generate both commensurate and incommensurate short-wave instabilities (e.g. [START_REF] Houchmandzadeh | [END_REF]13]). We show that the quasicontinuum model captures all unstable modes exhibited by the discrete model and provides good quantitative bounds for the instability thresholds.

Discrete model

Consider an infinite chain of interacting particles with the total energy

W = ε ∞ n=-∞ q p=1 pφ p u n+p -u n pε . (1)
Here φ p (w) is the energy density of an effective spring with reference length pε (representing interaction of pth nearest neighbors) and u n (t) is the displacement of nth particle. The equilibria in this system satisfy the following infinite system of difference equations:

q p=1 φ p u n+p -u n pε -φ p u n -u n-p pε = 0. ( 2
)
To access stability of a homogeneous equilibrium state u 0 n = nεw, where w is the average strain, we need to introduce perturbations v n = u n -u 0 n and study the positive 2 definiteness of the quadratic part of the energy expansion. This leads to the following eigenvalue problem:

-ω 2 v n = q p=1 K p (v n+p -2v n + v n-p ), (3)
where K p = φ p (w)/p and ω 2 is the square of the characteristic frequency. By representing the normal modes (phonons) in the form v n = exp(ink), where k is a real wave number, we obtain the dispersion relation

ω 2 (k) = 4 q p=1 K p sin 2 pk 2 . ( 4
)
By symmetry and periodicity it suffices to consider the wave numbers in the interval 0 ≤ k ≤ π (Brillouin zone). A uniform deformation is then stable if and only if ω 2 (k) > 0 for all k ∈ (0, π]. This condition is the far reaching generalization of the Legendre-Hadamard condition of strong ellipticity in continuum elasticity [14].

One can obtain some immediate necessary conditions for stability by requiring that

d 2 ω 2 (0) dk 2 > 0 and ω 2 (π) > 0.
The first condition, written as

E = q p=1 p 2 K p > 0, (5)
means physically that the effective elastic modulus along the homogeneous branch of equilibria is positive; the corresponding eigenmode, v n = 1, is infinitely long wave. The second condition, q k=1, k odd

K p > 0, (6)
is less transparent: the corresponding unstable eigenmode v n = (-1) n is microscopic, commensurate and has the smallest possible wave length. To characterize the complete set of stability conditions (necessary and sufficient) explicitly, we need to specify the number of interactions. One can show that the first generic case is q = 3, when nearest (NN), next to nearest (NNN) and next to next to nearest (NNNN) neighbors interact. The corresponding dispersion relation reads

ω 2 (k) = 4 sin 2 k 2 K 1 + 4K 2 + 9K 3 -4(K 2 + 6K 3 ) sin 2 k 2 + 16K 3 sin 4 k 2 (7)
A straightforward analysis of the positive definiteness of the right hand side produces the following necessary and sufficient conditions of stability:

K 1 + 4K 2 + 9K 3 > 0 K 1 + K 3 > 0 K 2 2 < 4K 3 (K 1 + K 2 ) if -6 < K 2 K 3 < 2 (8)
The first two of these conditions have been already obtained as necessary (see ( 5), ( 6)) and can be identified with macroinstabilities (k = 0) and commensurate microinstabilities (k = π), respectively. The third condition (8) 3 is associated with the incommensurate mode 0 < k = 2 arcsin (K 2 + 6K 3 )/(8K 3 ) < π. The full stability diagram in the plane

-1 1 2 1 2 K 1 > 0 ¡ M N ¢ ¢ ¢ £ ¥¤ £ ¤ £ ¤ k k k (i) (ii) (iii)
γ 2 = γ 3 = 1 (path 1); γ 2 = 1, γ 3 = 0.4 (path 2); γ 2 = 0.4, γ 3 = 0.2 (path 3); γ 2 = γ 3 = 0.4 (path 4). In all cases µ 1 = µ 2 = µ 3 = γ 1 = 1.
of the nondimensional parameters is presented in figures 1 for K 1 > 0 and 2 for K 1 < 0.

α = K 2 K 1 , β = K 3 K 1 (9) ¡ £¢ k (i) ¡ £¢ k (iii) -2 -1 -2 -1 ¤ ¥ R -4 -3 3 A 3 K 1 < 0 stable unstable continuum theory
To relate the microscopic stability limits with the macroscopic stress-strain relation, consider a family of homogeneous configurations u n = wnε parametrized by the strain w. Assume that all three interactions (NN, NNN and NNNN) are governed by the Morse potentials

U p (r) = µ p ε 2γ 2 p 1 -exp -γ p r ε -1 2 , 1 ≤ p ≤ 3. ( 10
)
The energy density associated with pth interaction in the homogeneous state is

φ p (w) = U p (pε(w + 1))/(pε); (11)
the elastic moduli are given by:

K p (w) = µ p exp[-pγ p w](2 exp[-pγ p w] -1), 1 ≤ p ≤ 3. ( 12
)
Using ( 9) and ( 12), we can compute the strain dependence of the dimensionless parameters α(w) and β(w). Depending on the choice of the microscopic parameters µ p and γ p we obtain different paths in the α-β plane; each path starts in the stable region and eventually reaches the instability threshold. If the potentials are identical for NN, NNN and NNNN interactions (µ

1 = µ 2 = µ 3 , γ 1 = γ 2 = γ 3
), the resulting homogeneous equilibrium branch (path 1 in figure 1) crosses the macroinstability boundary M N , activating the instability mode v n = 1. The corresponding macroscopic stress-strain relation σ(w) is shown in figure 3. In this case the bifurcation point A 1 correspond to the failure of the Legendre-Hadamard conditions for the macroscopic energy (σ (w) = 0). Suppose next that the three considered interactions (NN, NNN and NNNN) are governed by different Morse potentials. For instance, we may assume that U 1 (r) = U 2 (r) but that the force associated with the NNNN interaction U 3 (r) decays at large r slower than U 1 (r). If γ 3 is sufficiently smaller than γ 1 = γ 2 , the corresponding equilibrium branch (path 2 in figure 1) crosses the incommensurate microinstability boundary. The stability is lost at point A 2 where σ (w) > 0, i.e. before the macroscopic threshold. Now suppose that U 2 (r) decays slower than U 1 (r) but faster than U 3 (r) (e.g. γ 1 = 1, γ 2 = 0.4, γ 3 = 0.2). The corresponding equilibrium path 3 starts in the stable region where all three elastic moduli are positive (see figure 1). As w increases, K 1 (w) approaches zero faster than the other two moduli, so that α and β both tend to infinity. When K 1 changes sign, parameters α and β become infinite and the path re-emerges in the stable region with K 1 < 0 (figure 2). The failure of the homogeneous state at A 3 takes place via incommensurate microinstability before the macroinstability threshold.

Finally, assume that NNN and NNNN interactions potentials are identical and that the corresponding forces decay slower than U 1 (r) (γ 1 = 1, γ 2 = γ 3 = 0.4). Then we obtain path 4 which crosses the commensurate microinstability boundary at point A 4 to the right of point N in figure 1. The corresponding microinstability again precedes the macroinstability (see figure 3). In the special case when the only two interactions are NN and NNN (q = 2, K 3 = 0), the stability conditions (8) reduce to K 1 + 4K 2 > 0, K 1 > 0 [15]; see the shaded region in figure 4a. In this case the microinstability is necessarily commensurate, and it takes place at K 1 = 0. As before, the path with identical Morse potentials for NN and NNN interactions (path 1) leads to macroinstability. If U 2 (r) decays slower than U 1 (r), the instability becomes microscopic (see figure 4b).

A 1 A 2 A 2 A 1 ¡ ¢ 2 k ¡ ¢ 2 k (ii) (i)

Quasicontinuum approximation

A higher gradient quasicontinuum approximation of a lattice model is obtained by replacing the discrete dispersion relation (4) by the first few terms of its Taylor expansion around k = 0 [16,17]. To capture all three types of instabilities (long-wave, commensurate and incommensurate) exhibited by the discrete model with q ≥ 3, it is necessary to consider a polynomial expansion of at least sixth order:

ω 2 (k) ≈ k 2 (E + A 1 k 2 + A 2 k 4 ), ( 13
)
where the Taylor coefficients are

E = K 1 + 4K 2 + 9K 3 A 1 = - K 1 + 16K 2 + 81K 3 12
A 2 = K 1 + 64K 2 + 729K 3 360 . ( 14
)
The quadratic part of the continuum energy function corresponding to (13) takes the form

W = 1 2 ∞ -∞ [Eu 2 x + A 1 ε 2 u 2 xx + A 2 ε 4 u 2 xxx ]dx. ( 15
)
To access stability of the homogeneous state in the continuum problem we need to solve the following eigenvalue problem:

-ω 2 v = Ev xx -A 1 ε 2 v xxxx + A 2 ε 4 v xxxxxx . ( 16
)
One can immediately see that if A 2 < 0 the energy is unbounded from below (ω 2 (k) < 0 for sufficiently large |k|). This short-wave instability is unphysical if the unstable wave length is shorter than the lattice spacing. To eliminate this possibility we can limit the class of perturbations by imposing a constraint |k| ≤ π. This is achieved by replacing (13) with ω 2 (k) = 0 outside the first Brillouin zone (for |k| > π). The resulting dispersion relation with a short wave cutoff is compared in figure 5 with the dispersion relation for the exact quasicontinuum model [START_REF] Kunin | Elastic Media with Microstructure I: One-Dimensional Models[END_REF] obtained by replacing discrete dispersion by zero at |k| > π.

Despite its appearance, the proposed higher gradient approximation with a cutoff is essentially a continuum model with long-range spatial memory. Indeed, if we compute inverse Fourier transform of the truncated dispersion relation (13), we obtain the integral model

W = ∞ -∞ ∞ -∞ Φ(x -ξ)u(x, t)u(ξ, t)dxdξ (17)
with the kernel

Φ(x) = 1 πx 7 [2πx((E + 2A 1 π 2 + 3A 2 π 4 )x 4 -12(A 1 + 5A 2 π 2 )x 2 + 360A 2 ) cos πx + (π 2 x 6 (E + A 1 π 2 + A 2 π 4 ) + 2(E + 6A 1 π 2 + 15A 2 π 4 )x 4 + 24(A 1 + 15A 2 π 2 )x 2 -720A 2 ) sin πx]. (18)
Meanwhile, the exact quasicontinuum model has the kernel (e.g. [START_REF] Kunin | Elastic Media with Microstructure I: One-Dimensional Models[END_REF])

Φ D (x) = 4 sin πx πx K 1 x 2 -1/2 x 2 -1 - 2K 2 x 2 -4 + K 3 x 2 -9/2 x 2 -9 (19)
K 1 = 1, K 2 = 0.5, K 3 = -0.02.
The two kernels are compared in figure 5b; one can see that the nonlocal interactions decay in both cases as a power of distance which is characteristic for the models with long-range interactions.

The proposed quasicontinuum model can be used to study stability of the homogeneous state. Following the same procedure as in the discrete case, we obtain the following necessary and sufficient conditions of stability:

E > 0 E + A 1 π 2 + A 2 π 4 > 0 4A 2 E -A 2 1 > 0 if 0 < - A 1 2A 2 < π 2 . (20)
The first condition, indicating macroinstability, coincides with (8) 1 . The second condition is the analog of (8) 2 , and the unstable mode is again k = π. The last inequality in ( 20) is analogous to (8) 3 ; its failure corresponds to incommensurate microinstability with the wave number

k = k * = - A 1 2A 2 . (21)
A combined stability diagram illustrating conditions 20 is presented in figure 6 where it is compared with the stability diagram for the discrete model. Observe first that in the discrete and quasicontinuum models at K 1 > 0 the stability boundaries are tangent at point M where the unstable wave is infinitely long (k = 0).

-1 Above point M the quasicontinuum model reproduces the main qualitative features of the discrete diagram. Between points M and Q = (-4(π 2 -3)/(13π 2 -30), (π 2 -6)/(3(13π 2 -30)) the instability mode is macroscopic in both models; this interval, however, is much shorter in the quasicontinuum model. Between points Q and P = (-(120 -80π 2 + 3π 4 )/(10(30 -26π 2 + 3π 4 )), -(30π 2 -π 4 -90)/(45(30 -26π 2 + 3π 4 )) the stability loss is via microscopic mode with k = π. Below P the instability mode is again incommensurate; such transition from commensurate to incommensurate instability is not observed in the discrete case. Overall, the quasicontinuum theory underestimates stability domains in the region {K 1 > 0, K 2 > 0, K 3 < 0}. This is due to the larger contribution of the oscillation-producing K 3 terms in the polynomial model compared to the discrete case.

1 2 3 -0.5 ¡ M Q P ¢ £ ¥¤ ¢ £ ¥¤ ¢ £ ¤ £ ¥¤ ¢ (a) K 1 > 0 1 2 4 B 2 B 1 B 4 (i) (ii) (iii) (iii) -2 -1 -2 -1 ¢ £ ¤ (iii) 3 (b
At K 1 < 0 the quasicontinuum approximation has the opposite effect: it overestimates the stability of the homogeneous state. This is again caused by the very nature of the polynomial approximation: for positive K 2 and K 3 (which at K 1 < 0 result in negative α and β) the smoothening quadratic and quartic terms in the dispersion relation easily dominate the unstable contribution due to K 1 . Notice also that while in the discrete case the boundary of the stable domain contains a segment corresponding to commensurate instability (which at point R becomes incommensurate), in the polynomial quasicontinuum model the whole instability boundary is due to incommensurate mode.

It is instructive to compare the instability mechanisms along the four deformation paths considered earlier for the discrete model with the predictions of the quasicontinuum models (see figure 6). Along path 1 the quasicontinuum approximation predicts a commensurate microinstability (at point B 1 ), whereas the discrete model predicts macroscopic instability. Both models predict incommensurate instability along paths 2 and 3, with stability failure somewhat delayed in the quasicontinuum model. Finally, path 4 is predicted to pass through microinstability boundary by both models; the mode of instability is, however, commensurate in the discrete case and incommensurate in the quasicontinuum approximation. The analysis of the relative position of the bifurcation points on the stress-strain curves, presented for both models in figure 6c, shows a reasonable quantitative agreement.

Conclusions

It has been long recognized that higher gradient approximations of lattice models generate in the continuum limit either unbounded or nonpositive definite operators leading to ill-posed mathematical problems (e.g. [START_REF] Christov | [END_REF]20]). To overcome this difficulty we propose to regularize such operators by restricting them to a finite sphere in the Fourier space. At short wave lengths the resulting quasicontinuum model effectively replaces partial differential equations by integral equations. The dual nature of such cutoff polynomial models may be used to design hybrid computational schemes filtering parasitic small-scale oscillations while taking full advantage of the availability of simple partial differential equations for slowly varying fields. As we showed, the new quasicontinuum model combines the analytical simplicity of the gradient models at long waves with the physically correct description at small scales provided by models with long spatial memory. In particular, the approximation was shown to be sufficient to capture the whole spectrum of short-wave instabilities exhibited by the prototypical lattice model.

Figure 1 :

 1 Figure 1: Stability diagram for the discrete model with q = 3 and K 1 > 0. Regions of stability are shown in gray. Dispersion relations corresponding to different modes of stability loss are shown in inserts. Stability boundaries: (i) commensurate microinstability, β = -1 (to the right of N ); (ii) macroinstability, 1 + 4α + 9β = 0 (between N and M ); (iii) incommensurate microinstability, α 2 = 4β(1+α) (above M and below R). Classical continuum stability boundary is given by a gray line, 1 + 4α + 9β = 0. Dashed lines with numbers indicate deformation paths: γ 2 = γ 3 = 1 (path 1); γ 2 = 1, γ 3 = 0.4 (path 2); γ 2 = 0.4, γ 3 = 0.2 (path 3); γ 2 = γ 3 = 0.4 (path 4). In all cases µ 1 = µ 2 = µ 3 = γ 1 = 1.

Figure 2 :

 2 Figure 2: Stability diagram for the discrete model with q = 3 and K 1 < 0. Stability boundaries: (i) commensurate microinstability, β = -1 (to the left of R); (iii) incommensurate microinstability, α 2 = 4β(1 + α) (below R). Dashed line indicates deformation path 3 which is also shown in figure 1.

Figure 3 :

 3 Figure 3: Stress-strain relations along the deformation paths shown in figures 1 and 2. Points A i mark the loss of stability of the homogeneous deformation.

Figure 4 :

 4 Figure 4: (a) Stability diagram for the discrete model with q = 2. The trivial solution is stable in the gray area. Stability boundaries: (i) commensurate microinstability, K 1 = 0; (ii) macroinstability, K 1 + 4K 2 = 0. In the classical continuum model the trivial solution is stable above the gray line K 1 + 4K 2 = 0. The dotted curves show deformation paths: γ 2 = 1 (path 1) and γ 2 = 0.4 (path 2). In both cases µ 1 = µ 2 = γ 1 = 1. (b) The corresponding stress-strain curves.

Figure 5 :

 5 Figure 5: (a) Dispersion relations ω(k) for the exact quasicontinuum model (thick solid line) and the cutoff polynomial approximation (dashed). (b) The kernels Φ(x) of the exact quasicontinuum (solid) and cutoff polynomial (dashed) nonlocal models in physical space. Parameters: K 1 = 1, K 2 = 0.5, K 3 = -0.02.

) K 1 Figure 6 :

 16 Figure 6: Stability diagram for the quasicontinuum polynomial model with the cutoff compared with the discrete stability diagram. Regions of stability are shown in gray, the corresponding graphs of ω 2 (k) (dispersion relations) are shown in the inserts. Dashed curves indicate stability bounds for the discrete model. The two main diagrams correspond to the cases (a) K 1 > 0 and (b) K 1 < 0. Stability boundaries: (i) commensurate microinstability, E + A 1 π 2 + A 2 π 4 = 0 (P Q); (ii) macroinstability, E = 0 (M Q); (iii) incommensurate microinstability, 4A 2 E -A 2 1 = 0 (above M , below P and entire boundary of K 1 < 0 stability region in (b)). Classical continuum model predicts stability above the gray line in (a) and below it in (b). In (c) we present stressstrain relations along the deformation paths (the same as in figure 1) shown by the dotted curves in (a) and (b), with bifurcation points for quasicontinuum model shown by gray circles and for the discrete model by black circles.

Acknowledgements. This work was supported by the NSF grants DMS-0102841 (L.T.) and DMS-0137634 (A.V.).


Marc Bonnet

Sub Ra

T A Mukherjee

IMPLICIT BEM FORMULATIONS FOR USUAL AND SENSITIVITY PROBLEMS IN ELASTO PLASTICITY USING THE CONSISTENT TANGENT OPERATOR CONCEPT

published or not. The documents may come

Implicit BEM formulations for usual and sensitivity problems in elasto-plasticity using the consistent tangent operator concept

l. INTRODUCTION

In this paper, we address two important topics within the context of boundary element method (BEM) analysis of (small strain) elasto-plastic problems. The fi rst is an implicit BEM formulation for usual elasto-plastic analysis and the second is a sensitivity formulation for such problems. Both these formulations involve the consistent tangent operator (CTO see [START_REF] Simo | Consistent tangent operators for rate-independent elastoplasticity[END_REF].

Sensitivity analysis of nonlinear (material and/or geometrical) problems in solid mech anics is an a<:tive research area at present. In this context, design sensitivity coefficients (DSCs) are rates of change of response quantities, such as stresses or displacements in a loaded body, with respect to design variables. These design variables could be shape parameters, sizing parameters, boundary conditions, material parameters etc. DSCs are useful in diverse applications, a very important one being optimal design using gradient based optimization algorithms. Such analyses can be applied, for example, to optimal design of certain manufacturing processes.

Currently, the direct differentiation approach (DDA) or the ad joint structure approach (ADA) are popular for accurate sensitivity analysis. Either of these can be applied in conjunction with general purpose numerical methods such as the fi nite element method (FEM) or the BEM. The FEM has been used for sensitivity analysis of nonlinear problems by, among others, Arora and his eo-workers [START_REF] Arora | Variational principle for shape sensitivity analysis[END_REF]Yao and Arora, 1992a, b), Choi and his eo-workers [START_REF] Choi | Variational methods for design sensitivity analysis in nonlinear structural systems; part I : theory[END_REF][START_REF] Santos | Sizing design sensitivity analysis of nonlinear structural systems, part II : numerical method[END_REF], Haber and his eo-workers [START_REF] Vidal | The consistent tangent operator for design sensitivity of history dependent response[END_REF]Vidal and Habet, 1993), Kleiber and his eo-workers [START_REF] Kleiber | Computational coupled non-associative thermoplasticity[END_REF][START_REF] Kleiber | Incremental finite-element sensitivity analysis for non-linear mechanics applications[END_REF][START_REF] Kleiber | Parameter sensitivity of elasto-plastic response[END_REF], [START_REF] Michaleris | Tangent operators and design sensitivity formulations for transient non-linear coupled problems with applications to elastoplasticity[END_REF] and [START_REF] Badrinarayanan | A sensitivity analysis for the optimal design of metal forming processes[END_REF]. Haber, Kleiber and their associates were the fi rst to point out that the 1 consistent (or algorithmic) tangent operator (CTO) (as opposed to the continuum tangent operator) plays a key role in nonlinear sensitivity analysis. The CTO was originally proposed by [START_REF] Simo | Consistent tangent operators for rate-independent elastoplasticity[END_REF] for FEM analysis of (usual) nonlinear problems. The sensitivity problem is always linear (even if the usual problem is not) and the global or system matrix related to the CTO is precisely the stiffness matrix for these problems. Use of the CTO, as pointed out by Kleiber, Haber and their associates, provides very accurate numerical results for sensitivities, while other approaches (e.g., using the continuum tangent) might lead to significant errors. These researchers present numerical results for materially nonlinear problems. [START_REF] Michaleris | Tangent operators and design sensitivity formulations for transient non-linear coupled problems with applications to elastoplasticity[END_REF], in a recent paper, present sensitivity formulations for general transient nonlinear coupled problems, together with an accurate numerical pro cedure for the calculation of the CTO. Badrinarayanan and Zabaras (1995) present a consistent scheme for sensitivity analysis of nonlinear (both material and geometric) prob lems in solid mechanics. This paper presents very accurate numerical results for sensitivities at the end of a large deformation (extrusion) process.

All the researchers cited above have employed the FEM in order to obtain their numerical results. Mukherjee and his eo-workers have been active in solving nonlinear (both material and geometric) sensitivity problems by the BEM. Examples of the work of this group, using the explicit BEM, are Zhang et al. (1992a, b) and [START_REF] Leu | Sensitivity analysis and shape optimization in nonlinear solid mechanics[END_REF] (see also the forthcoming book: [START_REF] Chandra | Boundary Element Methods in Manufacturing[END_REF]. [START_REF] Wei | Shape optimization in elasticity and elastoplasticity[END_REF] have used sensitivities to carry out shape optimal design of an elasto-plastic problem.

Most of the publications on BEM analysis of (usual) nonlinear problems in solid mechanics report on the use of the explicit approach for time integration of the appropriate rate equations. Banerjee and his eo-workers [START_REF] Banerjee | The Boundary Element Method in Engineering[END_REF] have presented variable stiffness formulations for such problems. Implicit BEM formulations have been presented by [START_REF] Jin | Boundary element formulation in finite deformation plasticity using implicit integration[END_REF] and Telles andCarrer (1991, 1994). Mukherjee and his eo-workers have been interested in implicit sensitivity calculations, using the BEM, during the last few years. Leu and Mukherjee (1994a, b) have presented implicit objective integration schemes for recovery of stress sensitivities at a material point. This work addresses large strain visco plastic problems but only considers integration of the algorithmic constitutive model (analogous to the radial return algorithm) at a material point. They have coupled this analysis with the BEM (Leu and Mukherjee, 1995) to solve general boundary value prob lems. The CTO, however, has not been employed in the work by Leu and Mukherjee cited above. It is observed (Leu and Mukherjee, 1995) that stress sensitivities at some material points, at the end of a large deformation process, can exhibit significant numerical errors.

The remedy appears to be an implicit BEM formulation that employs the consistent tangent operator. Within the context of the BEM, this paper presents, for the first time, an implicit scheme that explicitly utilizes the CTO. Small-strain elasto-plastic problems, with isotropic and kinematic hardening behavior, are considered in this paper, but further generalizations present no conceptual difficulties. Next, the corresponding sensitivity for mulation is derived. It is shown that the (converged value of) the "global" CTO appears, as expected, as the stiffness matrix for the linear system of equations that govern the sensitivity of the strain increment over a time step. Initial numerical results for sensitivities, using the formulation presented in this paper, are very accurate.

THE GLOBAL CONSISTENT TANGENT OPERATOR USING FEM

Let n denote the spatial region occupied by the elastic-plastic solid under consideration. Attention is restricted to small-strain, quasistatic loading processes.

Constitutive law

Let u, a and s denote respectively the displacement, stress tensor and total strain tensor (with 2s =(V+ V�u, where V is the gradient operator). Following [START_REF] Simo | Consistent tangent operators for rate-independent elastoplasticity[END_REF], considering the evolution problem from a discrete incremental standpoint for a finite time step /'J.. t (as opposed to continuous time), the elastoplastic constitutive law reduces to giving a rule which outputs, an+ 1 consistent with the yield criterion, for any given strain increment /'J.. sn = 8n+ 1 -8n (input):

(1)

Here, the notation a symbolically denotes the action of the radial return algorithm (RRA)

of [START_REF] Simo | Consistent tangent operators for rate-independent elastoplasticity[END_REF] ; eP is the cumulated equivalent plastic strain:

where dP is the plastic strain rate, with tr(dP) = 0. Also, the subscript n in (tno ... ) refers to time (or pseudo-time) tno tr is the trace of the tensor dP and

Radial return algorithm (RRA)

The RRA, which we now briefl y summarize, was given in [START_REF] Simo | Consistent tangent operators for rate-independent elastoplasticity[END_REF], in the context of rate-independent plasticity with Von Mises pressure insensitive yield condition, atlid an associative flow rule. The yield condition is

(2)

where�= s--IX in terms of the deviatoric stress s = a-�(tra)l and the back stress IX. Here, 1 is the second order unit tensor. Also, eP-+ K(eP) is the hardening rule.

A trial deviatoric stress is introduced as

(3)

where e = t-• �(tr t)l and G is the shear modulus of the material. If /(��+1, 1Xno Kn) � 0, i.e., ��+ I is elastic, one has (4) (with K the bulk modulus and ® denoting the tensor product) i.e., the elastic constitutive equation in incremental form. On the other hand, if /(��+ 1, 1Xno Kn) > 0, it is given by the following equations, which constitute the RRA:

-p - where H'(eP) is the plastic modulus. Here, a superposed "prime" indicates differentiation with respect to the argument.

-p + � [ At] en+ I -en V 3 YLl

Equilibrium constraint

In a FEM framework, the supplementary constraint provided by the need to satisfy the equilibrium equation at any time t (in e.g., the virtual work form) is then used to determine which Aa exactly should be input in (1). The necessary condition for equilibrium at step n + 1 is, using (1) :

= 0 ('Vv E "// ) (10)

with 2Aan =(V+ VT) (u-un), and v is an admissible variation. Then, the displacement increment Aum such that Un+l = un+Aun solves the nonlinear eqn (10), is sought for iter atively using a Newton method: the additive correction bu� = Au� The global CTO is the stiffness operator associated with the bilinear form t5u,

V--+ Jn(en+ 1 : t5u: v) d V.
This is essentially the approach followed, in conjunction with the RRA, in [START_REF] Simo | Consistent tangent operators for rate-independent elastoplasticity[END_REF], where the expression of en+ 1 associated with the RRA (5)-( 9) is also given.

It is important to note that when ii(am am e�, Asn) is elastic one has

(12)
where e is the fourth-order tensor of elastic coefficients.

THE GLOBAL CONSISTENT TANGENT OPERATOR USING BEM

The virtual-work form (10) of the equilibrium equation is not suitable for boundary element-oriented formulations. In this section, it is shown how to build the global CTO when the equilibrium constraint is formulated in terms of boundary integral equations and internal representation formulas.

BEM formulation for elastic problems with initial strain

We consider the initial strain approach of elastic-plastic problems. For a given initial strain distribution sP, the elastic constitutive law reads:

J�

1n n

where x is any fixed source point on the boundary an, holds for any elastic state with initial strain eP but otherwise no domain forces; U�, P � denote the components of the elastic singular fundamental displacement and traction, i.e., those created in the infinite space IR 3 by a unit point force applied at x along the k-direction, p = a• n is the traction vector. The variable field point is denoted by z in the above equation.

The above BIE symbolically reads:

[H]{u}-[G]{p} = [Q]{C:sP} (15)
In the standard boundary element method, the above equation is discretized and then recast as

[A]{y} = {f}+[Q]{C:sP} ( 16
)
where {y} collects the boundary unknowns and { f} is the contribution of known boundary variables, i.e., values prescribed by the boundary conditions.

BEM representation at internal points. The displacement at any point x not on the boundary is then given in terms of the boundary variables u, p and initial strain eP by the representation formula:

Differentiation of the above formula with respect to x1 yields the representation formula for the displa.cement gradient. A suitable regularization procedure which has been discussed elsewhere [START_REF] Chandra | Boundary Element Methods in Manufacturing[END_REF]) is needed at this stage for the domain integral.

As a result, one has uk,,(x) = I U;(z)D�(x, z) dSZ-I p;(Z) uz,(x, z) dSZ Jao Jao

using the notation D� = Cijahn j U!,ht• The total strain at xis then readily obtained from the above equation. In symbolic form, one has, with {s}: "vector" of strains at all internal points:

{s} = [G']{p}-[H']{u} + [Q']{C: sP} = -[A']{y} + {f} + [Q']{C: sP}
Substituting for {y} from (1 6) into the above equation, we have:

where

{s} = {n}+[S]{C:sP} {n} = {f'}-[A'][A)-1{f} [S] = [Q'] -[A'][ A] -I [Q] ( 19
)
Note that { 111 } denotes the purely elastic solution, i.e., the one obtained for the same loading but in the absence of initial strain. Then, (13)2 is incorporated (in the form {C: sP} = {C:s}-{a}) into (19), giving:

{a} = {n} + [S]({C: s}-{a})
Finally, the total strain and the stress are related through:

[S]{ a -Cs} -{ n} + [I]{ B } = { 0} (20)
Following [START_REF] Telles | Implicit procedures for the solution of elastoplastic problems by the boundary element method[END_REF], the above formulae for elastic problems with "initial" strain are given in accumulated form (as opposed to the rate form).

Consistent tangent operator for the BEM elastic-plastic formulation

A new implicit BEM formulation, using the CTO, is presented in this section. With sP the plastic strain, we consider the evolution of the structure between time tn and tn+ 1• One has, from (20) and using the notation 8(

)n = ( )n+ 1 -( )n : (21)
which includes the equilibrium constraint.

On the other hand, the radial return algorithm (1) relates ii = O"n+I = O"n+L\an to L\sn.

Combining the constitutive and equilibrium equations in the form 1)), we obtain a nonlinear equation for L\sn of the form:

(where L\an(L\sm • •, • • ) comes from the BEM equation while ii(L\sm • •, • •) comes from the RRA (
The Newton method can also be applied in this case, and it is readily seen that the consistent tangent operator c+ I appears here as well. The additive correction c5a� = L\a� + I -L\s� to L\a� solves :

(23) The quantity ([S][C-C�+ d-[I]) = ([D�+ d-[I]
) is hereafter called the global CTO (see Kleiber et al., 1994 for the FEM version). Once the nonlinear eqn ( 22) is solved for L\sm all the variables at time tn+ 1 are readily computed. It is interesting to note that the Newton step (23) involves the difference [C-C�+ d between the elastic constitutive law and the local CTO, rather than the local CTO itself; this is entirely consistent with the fact that eqn (21) accounts for equilibrium as well as the elastic constitutive law, while, for the FEM, eqn (10) from which the Newton step (1 1) stems accounts for equilibrium only.

Also, it is important to note that the elastic constitutive law and the local CTO differ only at points (referred to as 'incrementally plastic') where the current strain increment has a non-zero plastic component. Hence, the Newton step (23) admits the following block decomposition :

(24) (25) Equations ( 24), ( 25) use the fact that (26)

The subscripts E, P indicate vectors and matrices restricted to the currently elastic (E) or plastic (P) internal nodes and collocation points. Thus, only the restriction to currently plastic nodes of the global CTO [D�+ 1-I] pp need to be factored. This shows that the global CTO has to be set up and factored only at currently plastic nodes, the currently elastic part { t:5s� } E being given explicitly by ( 25) once ( 24) is solved for { t:5s� } P

SENSITIVITY ANALYSIS

Now we are interested in computing the sensitivities of the mechanical variables associated with an infinitesimal perturbation of some design parameter b of an unspecified nature. As a matter of fact, such sensitivities are expected to be undefi ned at some particular locations, e.g., on the transition interface between the current elastic and plastic zones. However, in the interest of simplicity, we elect to disregard this issue, so that sensitivities of mechanical variables are treated as continuous in the derivations that follow. The reader is referred to a very recent paper, [START_REF] Lee | A computational method for design sensitivity analysis of elastoplastic structures[END_REF], for a detailed discussion of this issue.

Since sensitivity in the present context means a comparison of two history-dependent mechanical processes, the sensitivity computation must also proceed in an incremental way, resulting in an accumulation of sensitivity increments, each of which is a solution to a linear problem. Th1� adjoint structure approach (ASA), very powerful for steady-state situations, leads to regre:ssive computations when applied to evolution problems. Applied to nonlinear, history-dependent problems, the ASA would require the storage of the whole mechanical history and of the converged tangent operators for all time steps, thus becoming rather impractical. An excellent reference for this area of research is a recent paper by Michaleris et al. (1994). In contrast, the direct differentiation approach (DDA) fits perfectly into the progressive tlime-stepping computation scheme. The DDA is thus preferable for nonlinear evolution problems; it is used in the present work.

The DDA: principle Differentiation of the boundary integral eqn ( 16) and internal strain representation (19) with respect to b for given sP is a well understood process [START_REF] Chandra | Boundary Element Methods in Manufacturing[END_REF]. One n�adily sees that the evaluation of the derivative tP plays a key role. We now concentrate on this particular task. (33) where the new parameter p (as in [START_REF] Simo | Consistent tangent operators for rate-independent elastoplasticity[END_REF] is: Also, P-� is called y in Simo and Taylor (1985). Combining the above equation with the sensitivity version of (5), we obtain the sensitivity of the radial return algorithm as � = {P(I-� 1 ® 1) -(P-�)n ® n}: (t +K(l ® t):n * +{Kt® 1+2GP(I-� 1 ® 1)-2G(P-�)n ® n} : .1sn +0(1-�g�+ � ��+ {(1-P)I+(P-�)n ® n} :�n (34) Note that �n lis updated using: (35) which results from a differentiation of (9) w.r.t. b. It is very important to note that the * factor multiplying .1sn in (34) is equal to the converged value of the consistent tangent operator en+ I as given in [START_REF] Simo | Consistent tangent operators for rate-independent elastoplasticity[END_REF].

As mentioned at the start of this section, sensitivity equations such as (34) are not strictly valid throughout a history dependent elasto-plastic process. For example, as expected, strain sensitivities suffer jump discontinuities at the onset of plasticity. This is clear from tht: numerical results presented in Section 6 of this paper. Usually, however, this is not a major problem in practical optimization applications because possible local errors are not expected to significantly affect sensitivities of global functionals (usually defi ned over a part of a body or its boundary).

Sensitivity analysis

* Finally, substitution of (34) into (27) gives a linear equation for {.1sn} of the form (36)

The function {F} is completely known at this stage. It depends on the converged values of the variables at states n and n + 1 and known sensitivities-those of the stress etc. at state n and the known loading sensitivity {�n+ I}. The matrix multiplying Xsn is the converged value of the "global (or system) consistent tangent matrix" (see eqn (23)). Here it is also advisable to use the block decomposition, so that the sensitivity increment computation at currently plastic and elastic nodes is uncoupled:

ALGORITHM

Let us assume that a standard discretization has been adopted, using boundary elements and domain integration cells. The latter are associated to nodal DOFs for the total strains and stress a, using e.g. fi nite element type interpolations. Note, however, that the domain di:;cretization is restricted to the potentially plastic part ofQ. Then [A], [S], etc. and {s}, etc. fl!duce to ordinary matrices and DOF vectors, respectively.

The following algorithm is then proposed for solving the incremental elastic-plastic problem, from initial time t0 to final time tNr• The initial time t0 is assumed to correspond to the first yield load.

ForO� n �(Nr-l):

1. Compute { Lln n } (purely elastic internal strain). 2. Initialize {Ll£�} (e.g., to the elastic value).

Iterative solution of ( 22 (c) Update: { t +d = { t }+{ A Bn}, { � �+d = { � n+ v1{[yAt]} (using eqn (30)), { ; n+d ={�} (using eqn (34)).

Continue.

EXAMPLE-A SPHERICAL CAVITY IN AN INFINITE DOMAIN

A one-dimensional numerical example is discussed in some detail in this section. The example chosen here has closed form solutions, for both the usual mechanics problem as well as the sensitivity problem, for the case of linear work hardening materials. The purposes of this example are to illustrate the main features of the present general formulation and to check the accuracy of (some of) the numerical solutions by comparing them with the available analytical ones. Also, the efficiency of obtaining numerical solutions, as well as their robustness, are examined for this particular example.

Let n be the infinite three-dimensional domain exterior to a spherical cavity (radius a, center 0). A time-dependent uniform pressure p(t) is applied on the cavity surface: p = -p(t)n. Usual and sensitivity elastic-plastic problems will be considered.

We denote by (r, (), r/1) a system of spherical coordinates centered at 0, so that n = [a,+ oo [ x Y and an = { r = a} x Y, Y being the unit sphere associated with the spherical angles ((), r/1). Due to the spherical symmetry, the example is essentially one dimensional: the only nonzero mechanical variables, namely Un e m e ee = e q,q,, a m G ee = (} q,q,, depend only on the radial coordinate r. The strain and stress deviatoric tensors e and s then have the form (39) where k0, k1 and m are material parameters.

For this problem, fi rst yielding occurs at r =a at the yield pressure py = (4G/3)k0• A nondimensional load parameter A. is defi ned as (40)

The exact solution of the problem, for the linear hardening case (m= I) is given in the Appendix: for completeness.

Integral equations and representations

It is instructive to observe the one-dimensional counterparts of the key equations that are given in general form in the earlier sections of this paper. Some of these are given below and the rest are given in the Appendix.

Let x(r, () = 0, <P = 0), with r � a, denote a fixed point of Q. Due to the spherical symmetry of the problem and assuming the 'initial' strain sP to be purely deviatoric, the integral representation ( 17) reduces to u,(r)+u p(a) I P ;(x,z(a,B,l/J)) a2dw-rrpp(a) r U�(x,z(a,B,l/J))a2dw

Jg Jg

with dw = sin cP d() d</J. Here U'(x, z), P'(x, z) are the Kelvin fundamental solutions for isotropic elasticity, with the point force applied in the fixed r-direction at x, and E'(x, z) is the strain tensor associated with the displacement U'. The one-dimensional counterparts of eqns ( 19) and ( 20) are given in the Appendix.

Local and global consistent tangent operators

One of the primary contributions of this paper is the proposal to combine the BEM eqn ( 21) with the RRA to get eqn ( 22), and then to set up a Newton scheme to solve eqn (22). We feel that it is very important to illustrate this key procedure, in some detail, for the spherical cavity problem at hand. This is done below.

Let us consider a loading process of the form (40), so that each load step is defined by a given increment LUm n � 0. Equation (A.5), taken at the beginning and the end of a load step, gives the following relation, which has the structure of eqn ( 21): The governing operator is the converged global CTO, which has already been set up and factored during the elastoplastic state solution step.

Internal discretization

The numerical treatment needs a discretization of the potential plastic zone [a, c] (c being arbitrarily chosen by the user) into integration cells. Here, a constant size L\r = (ca)/N is used, together with a piecewise linear and continuous interpolation for the scalar deviatoric strain e and stress s. Substitution of this interpolation into eqns (44), (48) collocated a:t the radii r = rk = a+kl\r leads to the discretized nonlinear equation and associated Newton step.

Since the operator [D�+ d yields nonzero results only when applied to a currently plastic { c5e�+ 1} , it is readily seen that the Newton step (48) takes the following matrix form, where the currently elastic and plastic nodal values are uncoupled: with

Numerical results: linear isotropic hardening case

In this ex:ample (referred to as 'example I'), the following material parameter values were used: m= I (linear hardening), k0 = k1 = 0.001 (G, a were used as scaling parameters in the implementation, so that their precise values have no influence on the results). A potential plastic radius cja = 4 has been used.

First, an increasing loading process, with �nal = 3 (i.e., so that P fi nai = 4py), is considered. Figure I shows the computed final values of the total strain e and accumulated plastic strain e', compared against the analytical solution, while Figs 2 and 3 depict the final sensitivity with respect to the hardening coefficient k1 of e and e', respectively. Two grid sizes were considered: L\rja = 0.1 (coarse mesh, N = 30) and L\rja = 0.05 (fi ner mesh, N = 60). Both the elastoplastic state and its k1-sensitivities compare well with the analytical solutions. The finer mesh gives very good sensitivity values, and catches accurately the sensitivity jump at the plastic-elastic transition radius r = b;:::; 2.70a. Figure 4, which shows the k1-sensitivities for 2.6a:::; ; r:::; ; 2.9a (zoom around r =b) and using two fine meshes (L\rja = 0.05 and L\rja = 0.01), highlights the convergence towards the analytical solution; the computed values for the fi nest mesh reproduce accurately the expected sensitivity jumps (see Fig. 4).

Similar results have been obtained with the same data and domain grids as above but for a smaller fi nal load ).final = I (pfinal = 2py). The accuracy for both state and sensitivity variables improve for). = I compared to). = 3. The results for both cases are acceptable.

It should be emphasized that in both cases the numerical results obtained proved to be virtually independent on the number of load steps used. Indeed, all previous results as shown were computed using just one single load increment : .:\). = ).final• This fact has been consistently observed in all numerical tests performed on this problem. For example, another run with ).final= 5 and L\rja = 0.05 also gave almost identical results; I load step reached conwrgence after 6 Newton iterations on the equilibrium eqn (52), while 5 and 10 equal load inc:rements consumed respectively 14 and 22 Newton iterations.

On the other hand, for a fixed grid size L\rja, the final accuracy deteriorates as the final load increases; this trend has also been observed consistently in all tests performed. Tables 1 and2 show clearly the convergence of the computed state and sensitivity variables towards the analytical solution with decreasing grid size !irja, for lfinat = 1 and lfinat = 3 respectively. Consistent with previously mentioned observations, the error levels 0.00 --------------------------------r-O---<)--<?-+-P--� �� ��H HK�I--€>�--4-+� with a fixed mesh size are higher for ).final = 3 than for ).final = 1 ; however the convergence is clear in both cases. Here also, one single load increment was used. The tables also show the total number of Newton iterations (from 2 to 5) needed for convergence and report the bracketing of the elastic-plastic radius b. Numerical results: non/inear isotropic hardening case This example (referred to as 'example 2') used the values k0 = k1 = 0.001 together with various values of the hardening exponent m. No exact solution is available for this example.

Except when indicated otherwise, the results were computed using one load increment.

Figure 5 shows the computed final values of the total strain e and accumulated plastic strain eP for several grid sizes (final load: ).. final = 3, potential plastic zone boundary: c = 3a); the hardening exponent value is m = 0.1. The results for a very fine mesh (500 As expected, the supplementary CPU time consumed by the sensitivity calculation is very small compared with that necessary for the elastoplastic state calculation. This com ment is true for all sensitivity calculations presented in this paper. For example, with N = 600 and 50 load steps, CPU times of about 153 and 147 seconds were observed, depending on whether the sensitivities were evaluated or not.

Finally, it is important to point out that even in the presence of nonlinear hardening the results for a given final load and grid size appear to be virtually independent of the number of load steps used. This is apparent in Fig. 8, in which the computed m-sensitivities are shown for a large final load (Afinat = 5) and a fixed grid size (Arja = 0.04): the final results obtained using 1, 5, 10 or 50 load steps (and respectively 8, 24, 40 and 147 Newton iterations) are almost identical.

All computations were carried out on an APOLLO/HP 735 workstation.

CONCLUSIONS

Within the context of the finite element method, the value of the consistent tangent operator (CTO) for obtaining implicit solutions of elastoplastic problems has been dem onstrated by, among others, [START_REF] Simo | Consistent tangent operators for rate-independent elastoplasticity[END_REF]. It has been shown that use of other operators, such as the continuum tangent, leads to loss of quadratic convergence of the associated Newton method that is used for obtaining stepwise iterative solutions of the nonlinear problem. The associated sensitivity problem is stepwise linear. In this case, use of the CTO leads to a stiffness matrix which is the converged value of the global CTO for that load step (see, for example, the work of Haber, Kleiber and their associates cited in the references). Now, use of operators other than the CTO can lead to significant numerical errors in the computed values of sensitivities. Thus, the consequences of using, for example, the continuUlm tangent operator, can be more serious for sensitivity analysis than for the usual analysis of the nonlinear mechanics problem.

This paper presents, for the first time, the use of the CTO for solving usual and sensitivity problems in elastoplasticity by the boundary element method. An implicit BEM scheme is proposed in which one proceeds from step n to step n + 1 with the stress increment obtained from the BEM equations and the stress at step n + 1 obtained from the radial return algorithm (RRA). As expected, the global CTO appears when the Newton method is applied to solve the usual problem, and its converged value appears in the corresponding linear sensitivity equation for that step. Numerical results are presented for pressurization of a spherical cavity in an infinite elastoplastic medium. The results for all the examples are uniformly excellent. The algorithm proposed here appears to be both robust and very powerful. Very large load steps (up to 6p Y , where p y is the pressure for first yielding in the body) can be accommodated by the iterative algorithm and accurate results are obtained within a few iterations. It is also interesting to observe that, for the numerical examples presented in this paper, the final converged results are almost independent of load step size. On the other hand, as expected, for a fi xed spatial grid size !lrja, the fi nal accuracy decreases as the final load increases. Thus, as is usual in numerical computations, the chosen mesh size must be fine enough for the problem being solved in order to assure an acceptable solution.

The approach presented in this paper appears to be extremely promising for both the usual as well as sensitivity analysis of elasto-plastic problems.

  [yAt] solves the nonlinear consistency equation ll ��+ l ll-v1K(e�+v1 [yAt]) -2G[yAt] -v1AHn = 0 with AHn = Hn+ 1 -Hno Hn = H(e�), Kn = K(e�). The evolution of the back-stress 1Xn is

+ 1 -

 1 Au� to Au� solves with 2Aa� =(V+ VT)Au�. The local consistent tangent operator (CTO) is the fourth-order tensor en+ I = oiifoAan. It depends on the particular algorithm Asn--+ O"n+ I chosen.

  The boundary integral equation (BIE) I [u;(Z)-u;(x)]P7(x,z)dSZ-I p;(z)U7(x,z)dSZ = i uux,z)Cijabe�b(z)dVZ(14)

  Let us differentiate eqn (22) w.r.t. b, so that: (an * above a variable denotes a derivative of that variable with respect to b) {u(Bn, Un, e;, Asn) -Un -CAsn } + [I]{Asn } * * +[S]{it(8n, Un,e�, A8n)-un-CAsn } -{An n } = {0 }(27) (assuming that b is such that the elastic coefficients do not depend on it). This derivation * aims at finding a governing equation for Asn in terms of known variables at times tn and tn+ 1 and sensitivities at time tn. Obtaining the derivatives [S] and {�n } usually requires quite lengthy derivations and is thus by no means trivial. However, this aspect has been discussed, at length, in previous work on sensitivity analysis of elastic and elasto-plastic BEM for mulations (e.g., in[START_REF] Chandra | Boundary Element Methods in Manufacturing[END_REF], where the DDA has been applied on the appropriate integral eqns (14, 17 and 18). The focus in this paper is on obtaining {a( sm u m e�, Asn) } from the appropriate constitutive algorithm, and group the terms mul- * tiplying Asn and the rest separately. Of course, one expects where Y'n denotes the state at time tn. Note that the derivative of (22) w.r.t. b must be taken for the converged state Y'n+ �>i.e., with the value of Asn such that G(sm um �. Asn) = 0. Sensitivity of the RRA Let us first differentiate the consistency eqn (8) w.r.t. band with the converged value of [yAt]. Noting the following formulae: (a n+2Gilen -�n) �n+i = ( � �+v1 [)�At])K�+l Xnn = ilH� � �+V1[ y At]H�+i (28) (where I is the fourth order symmetric unit tensor with components �(i5;ki5 j 1+i5;,.b j k)) one gets, after grouping terms, the following scalar linear equation for [ y At] ] = (l -15)n: (an+ 2Gilen -�n) -0(1-15)11�n+ 1lle� where the abbreviated notation has been used. b = K�+l +H�+ I K�+i +H�+I +3G Then, differentiation of (6) with respect to b gives The sensitivity of ii is fi (8) n ): ( ; n+2G!en) -(l-n (8) o) :�n} (32) ll�n+lll Using eqns (28) 3 , (28)4, (30) and (32) in eqn (31) and grouping terms appropriately, one gets t +l = [(1-P ) I+ ( P-15 )n®n] : : n + [p(I -� 1 ® 1)-(p -15)n ® n]: (;n + 2G ! en) + 0(1-<5)�!+ I��

  As a result, the practical global CTO for both the usual mechanics and sensitivity problems is ([ D�+ d-(1]1 ) PP•

  the residual { G(Ll£�)}. (c) Convergence test: if {G(Ll£�)} � EPS, GOTO 3. (d) i:=i+l . (e) Compute the local CTOs C�+ 1 at all internal nodes; determine the sets of currently elastic (E) and plastic (P) internal nodes. (f) Set up and factor the global CTO [S(C -C�+ 1 ) -I) pp, set up [S(C-c�+ 1 ) -I] EP • (g) Solve (24) for { &�} P and compute { b£�} E using (25). (h) Update: { Ll£�} := { Ll£�} + { b£�} .

  (i) GOTO (b) (start new iteration). 3. Update: {eP}n+l = {e'} n+ v1 {[yAt]}, {an+d = { a (Ll ll � +J )}, {lln+d = { lln } + {Ll lln } • 4. Sensitivity problem: (a) Set up the right-hand sides {F} p, {F} E• * * (b) Solve (37) for { Ll£"} P and compute { LlBn} E using (38).

  e = ea s = s a Wlt a= 1ag 3•-3•-3 in terms of the scalar deviatoric strain e = arr-e ee and stress s = arr-aee • The material has isotropic elastic properties and obeys a von Mises yield criterion with nonlinear isotropic hardening (H' (eP) = 0), which, since l l a l l = v1, reduces to (see eqn 2):

  Integrations w.r.t. (), ljJ can be carried out analytically, and as a result the above equation becomes u,(r) =(-u,(a)----' -' defined in the Appendix. Using the boundary condition rr"(a) = -p, one gets the explicit form {k f+oo dp} u,(a)=a 3°( l+A.) -

  )-L\en(r) = S(v), V= Xn+ldn+l {en+ I [(K,b)n+l +K�+Ie�] -sn} (51)

Fig. 1 .

 1 Fig. 1. Example 1: total (e) and accumulated plastic (eP) strains, using 30 or 60 integration cells l load step , ).final = 3 .

Fig. 3 .

 3 Fig. 3. Example 1 : k1-sensitivity of accumulated plastic strain (eP), using 30 or 60 integration cells-1 load step, Afinal = 3 .

Fig. 4 .

 4 Fig. 4. Exa.mple I : k1-sensitivity of total (e) and accumulated plastic (eP) strains: zoom around the plastic/elastic transition radius-I load step, }.final = 3 .

5 Fig. 6 .

 56 Fig. 5. Example 2: convergence of total (e) and accumulated plastic (e") strains-! load step, A.finaL = 3, m= 0.1.

Fig. 7 .

 7 Fig. 7. Example 2: m-sensitivities of accumulated plastic strain if' for various values of exponent m-total load increment: ).final = 3; �rja = 0.05.

Fig. 8 .

 8 Fig. 8. Example 2: independence of m-sensitivities of total (e) and accumulated plastic (i!P) strainswith respect to load step--total load increment : ).final= 5; 8r/a = 0.04.

Table I .

 I Elastoplasti c and sensi ti vi ty results : evoluti on of L2 relati ve error wi th decreasing integration cell size-I load step, A.finaL = I

	l'lr/a	0.1	0.05	0.01	0.005
	Total number of iterations bfa bracketing (bexact = l.39525a) L2 relative error (e) L 2 relative error ( IJ") L2 relative error (def dki) L2 relative error (diJ"/dk1)	3 [1.4, 1.5] 7.64E-03 1.19E-2 7.5E-02 5.9\E-02	2 [1.35, 1.4] 1.60E-03 2.55E-03 6.54E-03 4.2\E-03	3 [1.39, 1.4] 3.76E-05 6.2\E-05 5.95E-05 5.95E-05	3 [1.395, 1.4] 6.66E-7 l.l1E-6 2.88E-5 2.15E-5

Table 2 .

 2 Elastoplastic and sensitivity results : evolution of L 2 relative error with decreasing integration cell size-! load step, A.finaL = 3

	A.r/a	0.1	0.05	0.01	0.005
	Total number of iterations bfa bracketing (bexa" = 2.70296a) L2 relative error (e) L 2 relative error ( IJ") L2 relative error (defdk,) L2 relative error (diJ"/dk1)	4 [2.8,2.9] 3.34E-02 3.55E-02 1.32E-01 l.29E-Ol	4 [2.7, 2.75] 5.35E-03 5.70E-0 3 1.59E-02 1.56E-02	5 [2.7, 2.71] 9.3\E-05 9.96E-05 3.66E-04 3.57E-04	5 [2.7, 2.705] l.64E-05 l.76E-05 9.08E-06 8.77E-06

Acknowledgements-Most of the research presented was performed when S. Mukherjee was on sabbatical leave at the Laboratoire de Mecanique des Solides, Ecole Polytechnique, Palaiseau, France. S. Mukherjee also acknowl edges partial financial support from a research contract from the General Motors Research Laboratories, Warren, Michigan, USA, with Cornell University.

and the unknown Ae. is to be found using Newton iterations, i.e., by solving the following linear equation on the additive correction be� until convergence of the residual G:

(48)

Note that the integration range in the integral operators S, D is in practice limited to those points p E [a,+ oo [ at which the current s. +As� is plastic, i.e., X�+ 1 (p) = 1.

Governing equations for the sensitivity problem

Sensitivities with respect to material parameters have been considered in this paper.

Exact solution

In the linear hardening case (i.e., m= 1), the outer radius b = b(A.) of the plastic region rE (a, b] for a given load A. ;;. 0 is given by the unique solution of the transcendental equation with X = (bfa) and

The analytical solution in the plastic region a ::;; r ::;; aXis then given as follows :

with the constants

In the absence of hardening, the previous solution still holds with k, = 0; moreover the plastic radius is explicitly given by X= bfa = exp(A./3) Exact solutions for various sensitivities, e.g., with respect to the hardening parameter k" can fi nally be obtained in a straightforward way.

To the best of the authors' knowledge, no analytical solution is known for general power-law hardening, i.e., when m# 0, 1 in (39) .

In tegral representation

The one-dimensional counterparts of eqns