Plastic Deformation Instabilities : Lambert Solutions of Mecking-Lücke Equation with Delay

The aim of this paper is the study of instabilities during plastic deformation at constant cross-head velocity. The deformation is supposed to be controlled by the emission of dislocation loops. Under some hypothesis analogous to the Mecking-Lücke relation, we derive a linear delay differential-difference equation. The “retarded” time term appears as the phase shift between the time of loop nucleation and the time at which the mean strain is recorded. We show the existence of the solution of strain equation. We give an analytic approach of solution using Lambert functions. The stability is also investigated close to the stable solution using a linearization of the number of nucleated loops functions.


Introduction
The morphological change of solids has been largely developed these last years in mechanical engineering and materials science.Localization of plastic deformation in homogeneous materials can be associated with instabilities of the stress-strain curves.This phenomenon can have very different aspects: Portevin-Le Chatelier PLC effect, Lüders bands, twinning, thermomechanical effect, avalanches of dislocations.Some criteria of plastic localization are proposed in [1].
The PLC effect was first related to a negative strain rate sensitivity by Penning [2].The physical origin of this phenomenon is the dynamic strain ageing associated with the interaction between mobile dislocations and diffusing solute atoms [3].In [4], localization of plastic deformation associated to the PLC effect was investigated by Kovács et al. in an Al−Mg alloy subjected to tensile deformation at constant stress rates and room temperature.In [5], a new model and numerical results are presented for a physically-consistent description of plastic material instabilities referred to as PLC effect, namely the oscillatory plastic flow that may be observed in metallic alloys subjected to load-or displacementcontrolled plastic deformation in a certain range of strain, strain rate, and temperature.Demirski and Komnik present in [6] an investigation of the relation between the dynamic (inertial) properties of the loading system and the plastic deformation kinetics at the jump and formulate the condition of the jumplike character of the deformation: the strain rate sensitivity of the flow stress must be below critical, which depends on the effective mass and the stiffness of the crystal-machine system.Graff et al. [7,8] propose finite element simulations and experimental observations of PLC effect and Lüders bands propagation in notched and compact tensile specimens of aluminum using the macroscopic PLC constitutive model.Stüwe and T óth [9] analyze the role of crystal orientation in the stability of tensile testing when plastic deformation and Lüders bands are produced by crystallographic slip.Yang and Tong [10] consider the interaction between solute atoms and mobile dislocations during plastic deformation in an aluminum alloy and explain the formation of the coarse slip bands in terms of dynamic strain aging under both single-and multiple-slip conditions.Louchet and Brechet [11] present the different types of dislocations patterning during uniaxial deformation as a function of significant physical parameters such as crystalline structure; they have shown that there is a competition between dislocation production and rearrangements and they prove that this phenomenon is controlled by strain rate and temperature.Miguel et al. [12] propose a simplified numerical model to study the deformation of ice single crystals; they analyze the avalanche-like rearrangements of dislocations during the dynamic evolution and characterize the viscoplastic deformation of ice or similar crystalline materials in the from of nonequilibrium statistical mechanics.
In this paper, following the work of Grilhé et al. [13] restricted to graphically analysis of the stability of solution, we give a complete mathematical study of the problem of plastic deformation instability.Under some assumptions and by using a linear analysis, we deduce a differential equation with delay.We show the theoretical existence and stability of the solution according to the characteristics of material and the time lag.We use the Lambert W-functions to give an analytical solution of our problem.This will enable us to validate our theoretical result of stability.The paper is presented as follows: in Section 2, we present the modeling of the physical problem of plastic deformation instability.In Section 3, we show the existence of the solution and we use the Lambert W-functions to give an analytical approach of the solution.In Section 4, we theoretically prove the asymptotic stability of the solution around the initial stress and we complete this work by a numerical validation of stability.

Presentation of the problem
2.1.Mecking-Lücke equation.We consider the elastic-plastic deformation.In this case, the strain rate ε is the sum of the plastic strain rate εp of the specimen and of the elastic strain rate εe = σ/M of the combined sample and loading system (with a stiffness M) ε(t) = εp (t) + εe (t). (2.1) Saïd Hilout et al. 3 For instance, we are interested in an expression of plastic strain rate.We can write where V is the volume of sample which is supposed to remain constant, b is the magnitude of Burgers vector along the tensile axis, and Σ(t) is the area in which the dislocation appears.The relation (2.2) is depending of the dislocation mechanism which is operative in the crystal.In the case where the dislocation density ρ(t) and mean velocity vary slowly and in a monotonous way and we suppose that plastic deformation is governed by Peierl's forces, (2.2) becomes εp (t) = bρ(t)v(t). (2. 3) The plastic deformation may also be controlled by the emission of dislocation loops from Frank-Read type sources model.We denote by n(t) the number of loops arising at time t in the unit volume and during unit time and by S the mean area swept by the loops supposed constant during times which are long enough compared with the period of instabilities.Equation (2.2) becomes εp (t) = bn(t)S. ( The area S depends on the instantaneous density of the forest and thus on the previous strain history of the sample.We suppose that S varies slowly.If we denotes by L the mean free-path of a dislocation, then we can write the rate of creation of dislocation in the form Replacing in (2.4) the term n(t)S by L ρ(t) obtained in (2.5), we get the relation given by Mecking and Lücke [14]: The relations (2.4) and (2.6) are established assuming that the area S is instantaneously swept by each dislocation as soon as it is emitted.

Derivation of the evolution equation [13].
Delay differential equation arises in many areas of mathematical modelling; for example: population dynamic, chemical kinetics, biosciences problems, and more general control problems.In this study, the system is governed by a principle of causality: the future state of the system is independent of the past states and is determined solely by the present.With a finite dislocation velocity, the area swept by a loop nucleated at a time t = 0 is a function S(t) which is depends on the mechanism considered (slip, twinning, etc.) and on the state of the crystal.After the flight-time τ , the mobile dislocation gets pinned or reaches the free surface of the sample having covered a constant area S(τ ) = S since it was emitted.Then only loops generated at a time t = t , with 0 < t < τ , will contribute to the deformation at a time t.We can write (2.4) in the form The number of loops nucleated at a time t is a function of time through the applied stress.Thus, the strain rate is given by the relation In a tensile test with constant strain rate ε, a stationary solution of (2.8) is σ = constant = σ 0 and we have ε(t) = bn σ 0 S(τ ). (2.9) n(σ 0 ) and σ 0 can be considered as remaining constants only during periods shorter than the duration of the tensile test.To simplify the problem, we suppose that [13] where S is a constant equal to S(τ ), δ(t − τ) is a Dirac distribution, and τ is the mean delay time given by (2.11) The approximation (2.10) amounts to replacing S(t) by a step function.The time lag given by relation (2.11) can be interpreted by the phase displacement between the time of loop nucleation and the time at which the main strain is recorded.Under the assumption (2.10), we can rewrite (2.8) in the form (2.12) We remark that the solution given in (2.9) is always valid.To investigate the stability of system strain-stress curves, we linearize the function n(σ) close to the value σ 0 Then (2.12) becomes where Saïd Hilout et al. 5 Finally, using (2.9) we derive the Mecking-Lücke equation with delay describing the temporal evolution of stress σ(t) = −aσ(t − τ) + aσ 0 . (2.16)

Resolution of the evolution equation
A linear constant-coefficient ordinary differential equation with a constant delay-time can be solved by Laplace transform techniques [15].In this section, we show the existence of the solution of problem (2.16).We also give an asymptotic approach to complete solution for system (2.16) based on the concept of Lambert W-functions.The stability of stress function close to a constant stress σ 0 is presented using an explicit solution of (2.16) in the form of an infinite series of modes written in terms of Lambert W-functions.Equation (2.16) is a linear retarded differential difference equation with delay time τ.The solution of either equation is determined uniquely when initial data ϕ defined on an initial interval is prescribed (ϕ is not necessary differentiable).To define a function σ in (2.16) for t ≥ 0, we impose an initial data on the interval [−τ,0] (e.g., we consider ϕ ≡ 1 in [−τ,0]).In fact, let ϕ be a given continuous function on [−τ,0] (ϕ is called preshape function) and we consider the problem (2.16) with initial data ϕ: The following result proves the existence of the solution of problem (3.1).

3.1.
Lambert W-functions.Some physical problems and modern engineering problem use the Lambert W-functions [17,18].Corless et al. [19,20] present some properties and applications in pure and applied mathematics of the Lambert functions.They have also developed the asymptotic expansions of the branches of W. The Lambert W-function is defined to be the inverse of the function ω → ωe ω .This function W(z) satisfies the following: W is a multivalued function from C to C. If z is real and z < −1/e then W(z) is multivalued complex.If z, real and −1/e ≤ z < 0, there are two possible real values of W(z): the branch satisfying −1 ≤ W(z) is denoted by W 0 (z) and called the principal branch of Lambert W-functions and the other branch satisfying W(z) ≤ −1 is denoted by W −1 (z).If z is real and z ≥ 0, there is a single real value for W(z) which also belongs to the principal branch W 0 (z).We can write the principal branch W 0 (z) in the following: Calculation of other branches of Lambert W-function, for k = 0;±1;±2;..., are given by where the coefficient C pm = (−1) p S(p + m, p + 1)/m!, S is a nonnegative Stirling number of the first kind [21], computable via the generating function

Analytic approach solution.
In this section we give an analytic approach to the solution of (3.1) using the Lambert W-functions.For this we consider the characteristic equation associated to (2.16) σ(t) = −aσ(t − τ).The analytic approach solution of characteristic (3.7) can be written as C k e (Wk(−aτ)/τ)t , (3.13) where where Λ(τ,N) is defined as a matrix with the functions ζ k (t) = e (Wk(−aτ)/τ)t as its elements [17], and (•) k represents the kth element of the corresponding vector (see the appendix).
We can choose the preshape ϕ = 1, then if we denote by 1 2N+1 the (2N + 1)-vector of coefficients 1, then we have Remark 3.2.We can write the solution of (3.7) with initial data ϕ in the form [22] where [z] denotes the integer part of z.Now, we consider the problem (3.1), the solution of (3.1) with initial data ϕ ≡ 1 can then be written as k ds e (Wk(−aτ)/τ)t , (3.17) where C k is defined by (3.15) and is the (2N + 1)-vector of coefficients aσ 0 , that is,

Stability analysis
Many techniques has been used to prove the stability of solution of differential-difference equations.In this section, we make use of the method of Liapunov functionals [16] that generalize the second method of Liapunov for ordinary differential equations.The exact region of asymptotic stability of solution of system (3.1) is obtained by the roots of a characteristic equation when they are in the left half-plane.We specify that the asymptotic stability of the solution of (3.1) around σ 0 amounts to study the asymptotic stability of the solution of the homogeneous equation (3.7) near of the origin, since any solution of (3.1) can be written in the form σ = σ 1 + σ 2 , where σ 1 is the solution of the homogeneous equation (3.7) having the initial data ϕ on [−τ,0], and where σ 2 is the particular solution of (2.16) having zero initial values on [−τ,0].Thus we see that if the zero solution of homogeneous equation (3.7) is asymptotically stable, then all solutions of (3.1) are asymptotically stable around a constant stress σ 0 .

Theoretical stability.
We give sufficient conditions for the stability and instability of the solution of (3.1).It is well known that the solution σ of (3.1) is asymptotically stable for every preshape continuous function ϕ on [−τ,0] if (3.8) has no zeros in C + where C + = {λ ∈ C/ Re(λ) ≥ 0}.We show the asymptotic stability using the following lemma (see [16, page 416]). where Now we present the result of asymptotic stability of solution of system (3.1).
Theorem 4.2.For every preshape continuous function ϕ on [−τ,0], the solution of (3.1) is asymptotically stable if and only if 0 < aτ < π/2 and the solution is unstable if and only if Proof of Theorem 4.2.We can write (3.8) in the form The solution of (3.1) is asymptotically stable if and only if the roots of (4.2) are in the left half-plane.Using Lemma 4.1 for (4.2) with α = 0 and β = aτ, the roots (4.2) have negative real parts if and only if ζ = π/2 and 0 < aτ < π/2.The proof of Theorem 4.2 is complete.

Numerical validation of stability.
We use the analytical approach solution in the form (3.17) of problem (3.1) to make the asymptotic analysis of stability of solution.
The following numerical results do not give the exact solution of (3.1), but they show asymptotic stability and instability of solution of (3.1) according to the parameter aτ.We use only the first terms (N = 0 and N = 1) in the expression (3.17) to make the asymptotic analysis of stability.Various calculations are made by using the MAPLE software and [23].These numerical results validate the theoretical result obtained in Theorem 4.2.Physical discussions.The results of this paragraph show that if aτ (a = MbS(∂n/∂σ)(σ 0 )) is lower than π/2, then the stress-strain curve is a horizontal σ = σ 0 .On the other hand, if aτ becomes higher than π/2, then periodic instabilities must appear.Several physical factors (M, bS, (∂n/∂σ)(σ 0 ), and τ) play a role in the strain-stress curve stability.
Saïd Hilout et al. 9 (1) First, the stiffness machine value M.This has been checked during deformation of Cu−Al alloys by Coujou and Vergnol [24]: with a hard stiffness machine, serrated stress-strain curves are observed and these curves become smooth with a soft stiffness machine.(2) Sb is the amplitude of an elementary step of deformation.In the case of twinning, these elementary steps are microtwins [25][26][27] so that Sb is large and aτ is higher than π/2.This must explain the observed twinning instabilities [25][26][27].(3) In the case of PLC effect, dislocations are pinned by impurities and are unlocked when stress becomes large.In this case, instabilities can be attributed to great values of (∂n/∂σ)(σ 0 ) [7,8].

Conclusion
In this article, we established a differential-difference equation with delay allowing to describe the plasticity of a solid becoming deformed by loops of dislocations or microtwinning.The linearized problem is used for study of deformation: we showed the existence of solution.The analytic approach of solution via the Lambert W-functions is presented.We could state the criterion of stability and describe the beginning of the deformation in the stable and unstable regions.For a long time, it is necessary to use the autonomous nonlinear equation where (5.2) Equation (5.1) derives from (2.12) by replacing n(σ(t)) by the second-degree Taylor polynomial expansion of n at σ 0 .
where L k is the kth Lambert coefficient and ζ k (t) = e (Wk(−aτ)/τ)t .The real R can go to infinity.
To find the values of the coefficients L k , we assume that the most dominant modes are the first N modes, where N is a large number, we can write h(t) in the form h(t) ≈ We denote by The vector L represents an approximation for the coefficient L k for large values of N. We assume that the matrix Λ(R,N) is invertible, then we can write the system (A.3) in the form L = Λ(R,N) −1 Θ. (A.5)

Lemma 4 . 1 .
All roots of the equation (z + α)e z + β = 0, where α and β are real, have negative real parts if and only if

Figures 4 . 1 (
a) and 4.2(a) show the asymptotic stability of the solution of (3.1) near to σ 0 .The beginning of phase instability of the solution of (3.1) is shown in Figure 4.3 and Figures 4.1(b) and 4.2(b).