NON-ISOTHERMAL CYCLIC FATIGUE IN AN OSCILLATING ELASTOPLASTIC BEAM

We propose a temperature dependent model for fatigue accumulation in an oscillating elastoplastic beam. The full system consists of the momentum and energy balance equations, and an evolution equation for the fatigue rate. The main modeling hypothesis is that the fatigue accumulation rate is proportional to the dissipation rate. In nontrivial cases, the process develops a thermal singularity in finite time. The main result consists in proving the existence and uniqueness of a strong solution in a time interval depending on the size of the data. Introduction. Elastoplastic materials subject to cyclic loading exhibit increasing fatigue, which is manifested by material softening , and material failure in finite time with strong heat release. The analysis of the so-called rainflow method for cyclic fatigue accumulation in uniaxial processes [2] has discovered a qualitative and quantitative relationship between accumulated fatigue and dissipated energy. Indeed, the rainflow algorithm counts closed hysteresis loops in the loading history, and with each closed loop associates a number depending on its amplitude – the contribution of the loop to total damage – taken from an experimental diagram (the Wöhler line). This corresponds to the mechanism of energy dissipation: The number associated with a closed loop is its area in this case. In multiaxial loading processes, the concept of closed loop is meaningless, and no counterpart of the rainflow algorithm is known to the authors. On the other hand, the notion of energy dissipation is independent of the experimental setting. In [5], we have proposed a model for fatigue accumulation based on the hypothesis that 2000 Mathematics Subject Classification. Primary: 74R20; Secondary: 74K10, 74H35, 74N30.

Introduction. Elastoplastic materials subject to cyclic loading exhibit increasing fatigue, which is manifested by material softening, and material failure in finite time with strong heat release. The analysis of the so-called rainflow method for cyclic fatigue accumulation in uniaxial processes [2] has discovered a qualitative and quantitative relationship between accumulated fatigue and dissipated energy. Indeed, the rainflow algorithm counts closed hysteresis loops in the loading history, and with each closed loop associates a number depending on its amplitude -the contribution of the loop to total damage -taken from an experimental diagram (the Wöhler line). This corresponds to the mechanism of energy dissipation: The number associated with a closed loop is its area in this case.
In multiaxial loading processes, the concept of closed loop is meaningless, and no counterpart of the rainflow algorithm is known to the authors. On the other hand, the notion of energy dissipation is independent of the experimental setting. In [5], we have proposed a model for fatigue accumulation based on the hypothesis that there exists a proportionality between fatigue and dissipated energy in the multiaxial case as well, and the model was demonstrated on the example of a transversally oscillating thermoelastoplastic plate. This idea is perhaps not completely unrealistic. Plastic deformations are manifested on the microscale by moving dislocations which, at the same time, damage the material cohesion, reduce the stiffness, and dissipate mechanical energy into heat. Note that in engineering practice, e.g. in aircraft industry, temperature tests for fatigue estimation are commonly used, and similar considerations appear also in theoretical engineering literature, going back e.g. to [11].
In the whole generality, a rigorous existence and uniqueness theory is still missing. Here, we make a first step in this direction and investigate the well-posedness of the system describing non-isothermal fatigue accumulation in a transversally oscillating elastoplastic beam. The equation for transversal oscillation of an Euler-Bernoulli beam with the single yield von Mises plasticity law was shown in [15] to give rise to a 1D beam equation with a multiyield Prandtl-Ishlinskii hysteresis operator as a result of dimensional reduction. Here, we allow the hardening modulus to depend on the accumulated fatigue, couple the momentum balance equation with the full energy balance involving temperature, and prove that the resulting PDE system admits a unique solution in the whole existence time interval until a thermal shock occurs. As the primary goal of the model is to describe material failure in finite time, we cannot expect to have solutions that are global in time. For some details about modeling issues, we refer the interested reader to [5]. More generally, the analysis of nonisothermal problems in mathematical modelling of advanced materials is gaining more and more importance in the recent years. We refer for instance, without aiming at completeness, to [1,7,8,18,19] in the setting of elastic materials, plastic materials, shape memory alloys and liquid crystals.
Some extensions of our analysis to a higher-dimensional setting and further generalizations are in progress. A thermodynamic model for fatigue accumulation in an oscillating elastoplastic Kirchhoff plate based on the same hypothesis, namely that the fatigue accumulation rate is proportional to the dissipation rate, is considered in [6]. The existence of a unique solution in the whole time interval before a singularity (material failure) occurs is proved there under the simplifying hypothesis that the temperature history is a priori given.
In Section 1, we derive the governing system of equations from the thermodynamic principles. The main existence and uniqueness results are stated in Section 2. Section 3 is devoted to the analysis of a spatially semi-discrete scheme and estimates for the approximate solution, and the proof of the main theorem is carried out in Section 4. 1. The model.

1.1.
Transversal oscillations of a beam. In [15], it was shown that the Euler-Bernoulli dimensional reduction, applied to the 3D momentum balance for transversal oscillations of a thin elastoplastic body with a single yield von Mises plasticity law, leads to the system for x ∈ (0, ) and t ∈ [0, T ], where the subscripts x and t denote partial derivatives with respect to x and t, the function w is the transversal displacement, 2h is the thickness of the beam, B > 0 is the hardening modulus, P is the multiyield elastoplastic Prandtl-Ishlinskii operator (1.12) with a specific form of the density function γ, is the mass density and f is the time integral of the external load. In fact, (1.1) results from the strain-stress law, as ε = w xx is, up to the factor −z, the only nontrivial component ε 11 of the strain tensor at distance z from the midsurface, and σ = u t is a weighted average over the thickness of the σ 11 stress component. Note that the multiyield character of the Prandtl-Ishlinskii constitutive law in the beam or plate bending problem, cf. also [9], reflects the fact that eccentric layers are subject to larger deformations than the central ones, so that plastic yielding propagates gradually from the outer surface towards the midsurface as on Figure 1. The hysteresis memory state of the Prandtl-Ishlinskii operator P keeps the information about the time dependent stress distribution across the profile after elimination of the transversal variable z, see Subsection 1.2 below. In the present paper, we want to include fatigue and temperature effects by extending the constitutive law into where m > 0 is a scalar time and space dependent fatigue parameter, ν is the viscosity, β represents thermal expansion, θ > 0 is the absolute temperature, and θ ref is a fixed referential temperature. The function B(m) accounts for material softening by the action of fatigue. In fact, physically, ε = w xx is the curvature and σ in (1.3) is the bending moment. The hypothesis that one component of the bending moment is proportional to the temperature with factor β can be justified in case that the beam has a layered structure with mass densities, elasticity moduli, and thermal expansion coefficients depending on the transversal variable z (e.g. bimetal, say, with two different values of the material constants below (z < 0) and above (z > 0) the midsurface). A discussion about the form of the Prandtl-Ishlinskii operator in the case of layered materials can be found in [15,Remark 2.2].
With the constitutive law (1.3), we associate the free energy operator where V is the Prandtl-Ishlinskii potential (1.15), and c (the specific heat capacity) is a given constant. The entropy operator S and internal energy operator U then read We consider the first and the second principles of thermodynamics in the form This is certainly satisfied provided the fatigue rate m t is positive, i.e. fatigue is increasing in time, and B (m) ≤ 0, i.e. the elasticity modulus decreases with increasing fatigue in agreement with experimental observations. The system has to be complemented with an evolution equation for m. Motivated by the rainflow method of cyclic fatigue evaluation, we assume, as mentioned in the Introduction, that the fatigue rate is proportional to the elastoplastic dissipation rate (1.10) We assume that the fatigue rate at a point x is proportional to the bulk elastoplastic dissipation rate in a neighborhood of the point x. This means, for a nonnegative function λ with (small) compact support we set The modeling requirement that the system should develop a thermal singularity at the moment of failure, is clearly satisfied in (1.11). In nontrivial cases, if the external load is large, ε 2 (x, t) also becomes large on a large set, and m t blows up to +∞. Consequently, also the heat supply in the energy balance (1.9) becomes infinite, producing thus a thermal shock in the system.

The Prandtl-Ishlinskii operator.
Assuming that a nonnegative function γ ∈ L 1 (0, ∞) is given, we define the Prandtl-Ishlinskii operator P : where s r is the stop operator with threshold r. The modeling idea goes back to [10,17]. In the literature, e.g. [16], it is also called generalized Saint-Venant model . Let us first recall the definition of the stop operator.
Definition 1.1. Let ε ∈ W 1,1 (0, T ) and r > 0 be given. The variational inequality where Q r is the projection of R onto the interval [−r, r], defines the stop and play operators s r and p r by the formula, see Figure 2, (1.14) The stop and play operators were introduced in [12]. The parameter r is the memory variable, and for each given time t 0 , the functions r → p r [ε](t 0 ), r → s r [ε](t 0 ) represent the memory state of the system. They keep the information about important local minima and maxima of the input ε in the time interval [0, t 0 ], see [3, Sections 2.3, 2.4]. Let us illustrate the Prandtl-Ishlinskii memory structure on the example corresponding to an input history of ε with an absolute minimum ε( The memory curves are represented in the left part of Figure 3. The distribution of the stress σ at time t 0 across the beam profile along the transversal variable z ∈ [0, h] from the midsurface z = 0 to the upper boundary z = h in case that the elasticity modulus E is constant in (0, h), is given by the formula [15,Eq. (44)], that is, where R is the von Mises yield limit and E is the Young elasticity modulus, see the right part of Figure 3. The interval [R/Ea, h] corresponds to the plasticized zone near the outer boundary z = h. Let us list here some basic properties of the play and stop operators. The proofs are elementary and can be found e.g. in [3,13].
It is easy to see that the variational inequality (1.13) can be equivalently written in the formε representing the energy balance. Indeed,ε(t)σ(t) is the power supplied to the system, which is partly used for the increase of the potential 1 2 σ 2 (t), and the rest r|ξ(t)| is dissipated. This enables us to establish the energy balance for the Prandtl-Ishlinskii operator (1.12). Indeed, if we define the Prandtl-Ishlinskii potential 15) and the dissipation operator we can write the Prandtl-Ishlinskii energy balance in the forṁ As a consequence of Proposition 1.2 (iv), we have We see that Hypothesis 2.1 (i) below enables us to estimate the dissipation from above in terms of the input velocity. Furthermore, for ε 1 , ε 2 ∈ W 1,1 (0, T ), we estimate the difference of the dissipation rates using (1.16) and Proposition 1.2 (ii) as follows: Statement of the problem. The mathematical analysis of the problem derived in Subsection 1.1 is independent of the exact values of the material constants. For any T > 0, we denote Ω T := (0, 1) × (0, T ). In an a priori unknown domain Ω T * , we consider the system for unknown functions u, w, θ, m, with initial and boundary conditions under the following hypothesis: and D is its associated dissipation operator (1.16). We assume that its distribution function γ ∈ L 1 (0, ∞) is such that γ ≥ 0 a.e., and The main result reads as follows.
Theorem 2.2. Let Hypothesis 2.1 hold. Then there exists T * such that system (2.1)-(2.6) has a unique solution in Ω T * with the regularity Remark 2.4. We restrict ourselves to zero initial conditions for u, w, and m. The motivation for this choice is twofold. First, in practice, it is impossible to determine the initial fatigue if we do not know the full loading history. Hence, the only realistic hypothesis for the problem of estimating the lifetime of a material subject to cyclic fatigue, is to assume that it is at rest, and completely undamaged at the beginning. The second reason is that general speculative initial conditions would require complicated compatibility conditions involving derivatives of hysteresis operators, which do not bring any added value, while the resulting additional technical difficulties in the proof would make the presentation less transparent.
3. Approximation. We choose a cutoff parameter R > 1 and an integer n ∈ N, and consider the truncated space discrete approximations of (2.1)-(2.4) for k = 1, . . . n − 1 and f k (t) = n k/n We prescribe initial conditions for k = 1, . . . , n − 1 and "boundary conditions" This is a system of ODEs for u k , w k , θ k with locally Lipschitz continuous right hand sides, and one integrodifferential equation for m (n) . We proceed as follows: We first check that the system admits a local solution on a time interval [0, T n ), and then establish a lower bound for 0 < T R ≤ T n independent of n. On the interval [0, T R ), we derive estimates for the approximate solutions which enable us, on the one hand, to show that the truncation in (3.5) never becomes active if R is sufficiently large, and can be removed. On the other hand, we select a convergent subsequence indexed by n and pass to the limit as n → ∞ to obtain a solution to (2.1)-(2.6).
3.1. Existence of local solutions to (3.1)-(3.6). For given functions a, r ∈ L 1 (Ω T ), such that We define an auxiliary function µ as the maximal solution to the differential equatioṅ that is, We see thatμ(t) blows up to +∞ for t 1/(2LAR). We plan to show that the solution of (2.1)-(2.6) can be continued until the blowup of µ occurs. To this end, we choose an arbitrarily small b ∈ (0, 1) that we keep fixed throughout the paper, and set We define the set   , if (a 1 , r 1 ), (a 2 , r 2 ) are two sets of data satisfying (3.10), then the corresponding solutions m 1 , m 2 satisfy the inequality (3.17) Here and in the sequel, we denote by | · | p for p ∈ [1, ∞] the norm in L p (0, 1).
Proof. We choosem 1 ,m 2 ∈ M , and definem 1 ,m 2 ∈ M by the formulâ a.e., hencem i ∈ M . We further have For a 1 = a 2 and r 1 = r 2 we obtain This result enables us to conclude that (3.1)-(3.6) admits a local W 1,∞ solution in an interval [0, T n ] for every n and R. First, denoting by w the vector (w 1 , . . . , w n−1 ), and ε = (ε 1 , . . . , ε n−1 ), we have −ε = Sw with a positive definite matrix S, which has the form Hence, the left hand side of (3.2) reads (I + S)ẇ. Furthermore, the dissipation rates |D[ε k ] t | depend Lipschitz continuously onε by virtue of (1.18). Hence, by Lemma 3.1, Eq. (3.6) defines a locally Lipschitz continuous mapping that with ε andε associates the solution m (n) and m (n) t . By (3.2),ε itself is a Lipschitz continuous mapping of u = (u 1 , . . . , u n−1 ). We see that (3.1)-(3.4) can be considered as an ODE system in u k , w k , θ k , with a right hand side, which is locally Lipschitz in L 1 (0, t) for every t ∈ [0, T ], and the existence of a local absolutely continuous solution in an interval [0, T n ] follows from the standard theory of ODEs. Consequently, the right hand side is bounded, and we conclude that the solution belongs to W 1,∞ (0, T n ).
In the sequel, we will systematically use the "summation by parts formula" (3.22) At the end of this subsection, we prove that θ k remain positive in the whole range of existence. As a first step, we test (3.4) by −θ − k (the negative part of θ k ). We obtain from (3.22) that d dt and Gronwall's argument yields θ − k (t) = 0 for all k and t, hence θ k (t) ≥ 0 for all k and t. Furthermore, we havė Let p be the solution of the differential equatioṅ for all k and all t ∈ [0, T n ].
In the remaining four subsections of Section 3, we derive three types of estimates for the approximate solutions of (2.1)-(2.4) defined by (3.1)-(3.6). The first estimate in Subsection 3.2 is the energy estimate. Its formal "continuous" counterpart consists in testing Eq. (2.1) by w xxt , differentiating Eq. (2.2) in t and testing by w t , testing (2.3) by 1, and summing up. This yields estimates for θ in L ∞ (0, T n ; L 1 (0, 1)), w t , w xx in L ∞ (0, T n ; L 2 (0, 1)), and w xt in L 2 (0, T n ; L 2 (0, 1)). Furthermore, it implies that T n are bounded from below by some T R > 0. This is indeed related to the energy balance (1.7). The equations of motion (2.1)-(2.2) can be formally written as w tt − w xxtt = −σ xx + f t , and the term ε t σ on the right hand side of (1.7), when integrated by parts, yields w t (−w tt + w xxtt + f t ). Denoting by K(w) = 1 2 (w 2 t + w 2 xt ), we obtain from (1.7) the equation d dt which has the usual interpretation that the time increment of the total (potential plus kinetic) energy equals the energy supply. The second estimate in Subsection 3.4 is the so-called Dafermos estimate. It consists in testing (2.3) by θ −1/3 and getting upper bounds for w xxt in L 2 (0, T n ; L 2 (0, 1)) and θ in L 8/3 (0, T n ; L 8/3 (0, 1)).
Higher order estimates are derived in Subsection 3.5. We first formally differentiate (2.1) in x and square both sides of the resulting equation to obtain In the right hand side, we have integrated by parts and replaced u t from (2.1) and u xx from (2.2). This enables us to estimate the L 2 -norm of w xxxt by the L 2 -norms of θ t and θ x , which are in turn estimated by testing the heat equation (2.3) by θ t and using the Gagliardo-Nirenberg inequality. Finally, we differentiate (2.1) in t and test by w xxtt , differentiate (2.2) twice in t and test by w tt , sum up the results, eliminating the terms in u tt by integrating by parts, and obtain an L 2 -bound for w xxtt . We now carry out this programme more rigorously.
By hypothesis, B(m k ) is bounded from below by a constant. Furthermore, with the convention that here and in the sequel, C denotes any constant independent of n and R. Here, C depends only on the lower bound for B(m k ), and on f t L 1 (0,T ;L 2 (0,1)) .
We see that (3.6) is of the form (3.10)-(3.11), with a(x, t) = − 1 2 (ε (n) ) 2 (x, t) and with 1 0 a(y, t) dy = 1 2 We conclude that the solution to (3.1)-(3.6) exists in the whole interval [0, 3.3. Discrete Gagliardo-Nirenberg inequality. Let p, q, s ∈ [1, ∞] be such that q > s, and let | · | p denote the norm in L p (0, 1). The Gagliardo-Nirenberg inequality states that there exists a constant C > 0 such that for every v ∈ W 1,p (0, 1) we have In fact, (3.30) is straightforward: If we introduce an auxiliary parameter r = 1 + s(1 − 1 p ) and use the chain rule d dx |v(x)| r ≤ r|v(x)| r−1 |v (x)| a.e., we obtain from Hölder's inequality the estimate Combined with the obvious interpolation inequality |v| h ≤ |v| The discrete counterpart of (3.30) reads . The two dissipation terms on the right hand side of (3.4) corresponding to (1.10) are nonnegative and can be omitted. To deal with the heat flux term θ k+1 − 2θ k + θ k−1 , we use (3.22) with ξ k = θ −1/3 k , η k = θ k , and the elementary inequality We obtain for all t ∈ [0, T R ], after summing up from k = 1, ..., n − 1 and integrating in time, that (3.33) The last term on the right hand side is bounded by virtue of (3.28). By Hölder's inequality,  3.5. Higher order estimates. We define ε 0 , ε n as solutions to the differential equations with initial conditions ε 0 (0) = ε n (0) = 0. Then (3.1) holds for all k = 0, . . . , n, and we havė This yields d dt We estimate the terms on the right hand side of (3.43) as follows: where we have used (3.15), (3.28), (3.37), and (3.38). To estimate in (3.43) the "fatigue" term (B(m k )ε k − B(m k−1 )ε k−1 ) 2 , we notice that We have Moreover, by Proposition 1.2, From the above considerations and from (3.43) we conclude that Gronwall's argument then yields 1 2n We now test (3.4) by θ k and obtain d dt where, by (3.37), (3.50) Using the vector notation (3.31), we have by (3.39) and (3.37) that and we rewrite (3.50) as (3.51) We estimate the right hand side of (3.51) using (3.32) as follows: In (3.6), we thus have in particular We thus may fix R sufficiently large depending only on the data of the problem such that for all x and t, so that the truncation in (3.6) is never active, we may set T * = T R , and all the above estimates are thus independent of R. In particular, there exist a constant C > 0 such that for t ∈ [0, T * ]. By comparison, we also have t 0 n 3 and similarly for u k . Finally, we differentiate (3.1) once in t and test byε k , (3.2) twice in t and test byẅ k , and sum the two equations up. Using (3.55)-(3.56), we immediately get the estimate Similarly, We check in the same way that there exist u, w, θ ∈ C(Ω T * ) such that, selecting again a subsequence if necessary, Finally, by Lemma 3.1, using also the estimate |m (n) x (x, t)| ≤ C analogous to (3.45) to control the differencem (n) − m (n) , we obtain m (n) t → m t strongly in L 1 (0, T * ; L ∞ (0, 1)), where m satisfies (2.4). This enables us to pass to the limit in (3.6) and (4.9)-(4.12) and conclude that (u, w, θ) is a strong solution to (2.1)-(2.4) with the regularity indicated in Theorem 2.2 and satisfying the initial conditions (2.5). It remains to check that the boundary conditions (2.6) hold. We have w n (t) = 0, hence and similarly for w(0, t), u(1, t), u(0, t). To complete the existence proof, we only have to verify the homogeneous Neumann boundary condition for θ. In other words, we have to check that for every ψ ∈ C 1 (Ω T * ) we have T * 0 1 0 (θ x ψ x + θ xx ψ)(x, t) dx dt = 0 . x ψ x +θ (n) xx ψ)(x, t) dx dt = 0 and (4.14) follows.

4.2.
Uniqueness. Let (u, w, θ), (ũ,w,θ) be two solutions of (2.1)-(2.6). We integrate the difference of (2.3) for θ andθ in time and test by θ −θ. Taking into account the estimates derived in the previous section and Proposition 1.2, we obtain We estimate the terms on the right hand side as follows. Using Lemma 3.1, Proposition 1.2, and the fact that w xx belongs to L ∞ (Ω T * ), we have For the last two integrals in (4.16), we obtain from Hölder's inequality In the next step, we test the difference of the time derivatives of (2.2) for w andw by w t −w t , the difference of (2.1) for u andũ by w xxt −w xxt , and sum up. Arguing as above, we obtain |w xxt −w xxt | 2 + |θ −θ| 2 (x, t) dx Gronwall's argument now yields that w =w, θ =θ, and the proof of Theorem 2.2 is complete.