A model for shape memory alloys with the possibility of voids

The paper is devoted to the study of a mathematical model for the thermomechanical evolution of metallic shape memory alloys. The main novelty of our approach consists in the fact that we include the possibility for these materials to exhibit voids during the phase change process. Indeed, in the engineering paper has been recently proved that voids may appear when the mixture is produced by the aggregations of powder. Hence, the composition of the mixture varies (under either thermal or mechanical actions) in this way: the martensites and the austenite transform into one another whereas the voids volume fraction evolves. The first goal of this contribution is hence to state a PDE system capturing all these modelling aspects in order then to establish the well-posedness of the associated initial-boundary value problem.


Introduction
Shape memory alloys are mixtures of many martensites variants and of austenite. They exhibit an unusual behavior: even if they are permanently deformed, they can totally recover their initial shape just by thermal or mechanical means. There may be voids in the mixture, which may appear when the mixture is produced by the aggregations of powders, as it has been recently proved in the engineering paper [57]. The composition of the mixture varies: the martensites and the austenite transform into one another whereas the voids volume fraction evolves. These phase changes can be produced either by thermal actions or by mechanical actions. The striking properties of shape memory alloys result from interactions between mechanical and thermal actions (cf., e.g., [9,42]).
We assume that the phases can coexist at each point and we suppose that, besides austenite, only two martensitic variants are present. However, this choice provides a sufficiently good description of the phenomenon, as we want describe a macroscopic predictive theory which can be used for engineering purposes. The phase volume fractions, which are state quantities, are subjected to constraints. In particular, their sum must be lower than 1, but not necessarily equal to 1, due to the presence of voids (cf. also [40] and [41] where this property is introduced in order to treat solidliquid phase transitions with the possibility of voids and [1] for the corresponding numerical results). It is shown that most of the properties of shape memory alloys result from careful treatment of those internal constraints (cf., e.g., [36]- [39]). All these quoted references are related to a three dimensional model taking the temperature, the macroscopic deformation and the volumetric proportion of austenite and martensite as state variables. Moreover, let us note that there are not too many mathematical models describing phase transitions in which the interactions between different types of substances and the possibility of having voids is taken into account: we can quote only the two contributions [40,41].
It is beyond our purposes to give a complete description of the existing literature on models for SMA. However, restricting ourselves to the macroscopic description of these phenomena, we can refer to the main contributions [36]- [39], [4,5,34,35,55] and [3,11,12,18,19,47,54,62] (and references therein) describing full thermomechanical models and studying the resulting PDEs from mathematical viewpoint respectively. We shall instead focus here on a generalization of the Frémond model for SMA introduced in [36]- [39].
Let us then explain in detail which is the main aim of this contribution, compare it with the results already present in the literature, and show the main mathematical difficulties encountered. As already mentioned, in this paper we deal with a generalization of the model introduced in [36]- [39] and later on studied in many contributions starting from the pioneering paper [25], where an existence and uniqueness result has been proved for the solution of a simplified problem, where all the nonlinearities in the balance of energy are neglected and the momentum balance equation is considered in the quasi-stationary form and fourth order terms (related to the second gradient theory) are taken into account. In the case when the fourth-order term is omitted an existence result dealing with the linearized energy balance equation has been proved in [21], while [22] one can find the proof of the existence of solutions to the linearized problem by including an inertial term in the momentum balance. We can report also of some results when some or all the nonlinearities are kept in the energy balance. The full one-dimensional model is shown to admit a unique solution both in the quasi-stationary case in [30] and in the case of a hyperbolic momentum equation in [31,59]. Existence results have been proved also for the three-dimensional model (cf. [23,28,43]). Finally, let us mention the uniqueness result for the full quasi-static three-dimensional model proved in [17] and an updated and detailed presentation of the Frémond model and related system of equations and conditions, applying to the multidimensional case as well, which is provided in [11,12], [37,Chapter 13], and [39]. Let us also point out [11,12] for recent existence and uniqueness results in the threedimensional situation, where the various nonlinear terms arising in the derivation of the model are accounted. The large time behavior of solutions is investigated in [26] in connection with the convergence to steady-state solutions and in [27,24] where the authors characterize the large time behavior according to the theory of dissipative dynamical systems.
However, all these contributions were dealing with the case in which no voids can occur between phases. To model this possibility and to solve rigorously the results PDE system is just our aim here. First, in the next Section 2, we derive a model taking the possibility of having voids into account, introducing rigorously the pressure which has a paramount importance on the mechanical behaviour. In order to do that, we follow the ideas of [40] in which this was done in case of a two-phase transition phenomenon. Then, in Section 3, we give a rigorous formulation of initial boundary value problem associated with the resulting PDEs and we state our main results: existence, uniqueness, and continuous dependence of solutions from the data. The proofs are carried over in Section 4 and Section 5. The main mathematical difficulties are due to the nonlinear and singular coupling between the equations. In particular, in order to describe the evolution of the absolute temperature variable, we shall use the entropy balance equation (cf. [13]- [15] for a complete derivation and motivation of this equation). This equation turns out to be singular in ϑ but the main advantage of using it is that once one has proved that a solution component ϑ does exists then it turns automatically out to be positive and the proof of positivity of the absolute temperature is historically one of the main difficulties of these types of problems. The idea here is to approximate the nonlinearities with regular functions, to solve the regularized system by means of a Banach fixed point argument and then to use compactness and lower semicontinuity arguments in order to pass to the limit and obtain a solution of the original problem.

The derivation of the model
In this section we explicitly derive a macroscopic model describing the evolution od SMA with the possibility of voids. The model is obtained by properly choosing the state quantities, the balance laws and the constitutive relations in agreement with the principle of thermodynamics and with experimental evidence.

The State Quantities
We deal only with macroscopic phenomena and macroscopic quantities. To describe the deformations of the alloy, the macroscopic small deformation ε(u), (u being the small displacement) and the temperature ϑ are chosen as state quantities.
The properties of shape memory alloys result from martensite-austenite phase changes produced either by thermal actions (as usual) or by mechanical actions. On the macroscopic level, some quantities are needed to take those phase changes into account. For this purpose, the volume fractions β i of the martensites and austenite are chosen as state quantities. For simplicity, we assume that only two martensites exist together with austenite. The volume fractions of the martensites are β 1 and β 2 . The volume fraction of austenite is β 3 . These volume fractions are not independent: they satisfy the following internal constraints due to the definition of volume fractions. Since we assume that voids can appear in the martensite-austenite mixture, then the β 's must satisfy an other internal constraint being the voids volume fraction. This is the case when the alloy is produced by aggregating powders as shown in [57]. In case no voids are considered, this sum should be equal to 1 and this considerably simplifies the analysis (cf., e.g., [25]). We denote by β the vector of components β i (i = 1, 2, 3) and the set of the state quantities is while the quantities which describe the evolution and the thermal heterogeneity are The gradient of ∇β accounts for local interactions of the volume fractions at their neighborhood points.

The mass balance
Assuming the same constant density ρ (the reader can refer to [41] for a model in which different densities of the substances are taken into account in a general twophase change phenomenon) and the same velocity U = u t for each phase, the mass balance reads Within the small perturbation assumption, this equation gives i 's are the initial values of the β i . For the sake of simplicity, we assume and have ∂ t (β 1 + β 2 + β 3 ) + div U = 0, hence Mass balance is a relationship between the quantities of δE , indeed, its effects will be included in the cinematic relations (cf. (2.9) in the following subsections).

The equations of motion
They result from the principle of virtual power involving the power of the internal forces, (cf, e.g., [37]) where V and δ are virtual velocities, the actual velocities being U and β t . The internal forces are the stress σ , the phase change work vector B, and the phase change work flux tensor H . The equations of motion are where ρ is the density, U t the acceleration of the alloy which occupies the domain Ω, with boundary ∂Ω and outward normal vector n. The alloy is loaded by body forces f and by surface tractions g , and submitted to body sources of damages A and surfaces sources of damage a (in the following we will suppose, for simplicity, A = a = 0).

The free energy
As explained above, a shape memory alloy is considered as a mixture of the martensite and austenite phases with volume fractions β i . The volume free energy of the mixture we choose is where the Ψ i 's are the volume free energies of the i phases and h is a free energy describing interactions between the different phases. We have assumed that internal constraints are physical properties, hence, we decide to choose properly the two functions describing the material, i.e., the free energy Ψ and the pseudopotential of dissipation Φ, in order to take these constraints into account. Since, the pseudopotential describes the kinematic properties (i.e., properties which depend on the velocities) and the free energy describes the state properties, obviously the internal constraints (2.1) and (2.2) are to be taken into account with the choice of the free energy Ψ. For this purpose, we assume the Ψ i 's are defined over the whole linear space spanned by β i and the free energy is defined by We choose the very simple interaction free energy where I C is the indicator function of the convex set Moreover, and by (k/2)|∇β| 2 we mean the product of two tensors ∇β multiplied by the interfacial energy coefficient (k/2) > 0. The terms I C (β) + (k/2)|∇β| 2 may be seen as a mixture or interaction free-energy. The only effect of I C (β) is to guarantee that the proportions β 1 , β 2 and β 3 take admissible physical values, i.e. they satisfy constraints (2.1) and (2.2) (cf. also (2.8)). The interaction free energy term I C (β) is equal to zero when the mixture is physically possible (β ∈ C ) and to +∞ when the mixture is physically impossible (β / ∈ C ). Let us note even if the free energy of the voids phase is 0, the voids phase has physical properties due to the interaction free energy term (ν/2)|∇β| 2 which depends on the gradient of β . It is known that this gradient is related to the interfaces properties: ∇β 1 , ∇β 2 describes properties of the voids-martensites interfaces and ∇β 3 describes properties of the voids-austenite interface. In this setting, the voids have a role in the phase change and make it different from a phase change without voids. The model is simple and schematic but it may be upgraded by introducing sophisticated interaction free energy depending on β and on ∇β .
For the volume free energies, we choose where K i are the volume elastic tensors and C i the volume heat capacities of the phases. Stresses σ i (ϑ) depend on temperature ϑ and the quantity l a is the latent heat martensite-austenite volume phase change at temperature ϑ 0 (see Remark 2.1 below).
Because we want to describe the main basic properties of the shape memory alloys with voids, we assume that the elastic matrices K i and the heat capacities C i are the same for all of the phases: Always for the sake of simplicity, we assume that where I stands for the identity matrix. Concerning the stress τ (ϑ), it is known that at high temperature the alloy has a classical elastic behaviour. Thus τ (ϑ) = 0 at high temperature, and we choose the schematic simple expression with τ ≤ 0 and assume the temperature ϑ c is greater than ϑ 0 . With those assumptions, it results

The pseudo-potential of dissipation
The dissipative forces are defined via a pseudo-potential of dissipation Φ introduced by J.J. Moreau (it is a convex, positive function with value zero at the origin, [33], [49], [50]). As already remarked, the mass balance (2.4) is a relationship between velocities of δE . Thus we take it into account in order to define the pseudo-potential and introduce the indicator function I 0 of the origin of R as follows ¿From experiments, it is known that the behaviour of shape memory alloys depends on time, i.e., the behaviour is dissipative. We define a pseudopotential of dissipation where λ ≥ 0 represents the thermal conductivity and c ≥ 0, υ ≥ 0 stand for phase change viscosities.

The constitutive laws
The internal forces are split between non-dissipative forces σ nd , B nd and H nd depending on (E, x, t) and dissipative forces by The nondissipative forces are defined with the free energy and the dissipative forces are defined with the pseudo-potential of dissipation where the subdifferential of Φ is with respect to δE . Relationship (2.14) gives where p is the pressure in the mixture and it results The state laws (2.10)-(2.13), besides implying that the internal constraints are satisfied, give also the value of the reactions, during the evolution, to these internal constraints. Relationships (2.10)-(2.13) and (2.15)- (2.19) give the constitutive laws It can be proved that our choice is such that the internal constraints and the second law of thermodynamics are satisfied (cf., e.g., [39,37] and the next Subsection 2.7).

The entropy balance
By denoting the entropy balance is due to (2.19), ϑQ is the heat flux vector, Rϑ is the exterior volume rate of heat that is supplied to the alloy, ϑπ is the rate of heat that is supplied by contact action, ε(u t ) is the strain rate. The constitutive laws, within the small perturbation assumption and (2.3), become

The set of partial differential equations
We assume also quasi-static evolution and, using again the small perturbation assumption, we get the following set of partial differential equations coupling the equations of motion (2.5), the entropy balance (2.27) and constitutive laws (2.29)-(2.34) This set is completed by suitable initial conditions and the following boundary conditions: where ∂ n is the normal outward derivative to the surface ∂Ω, g is the exterior contact force applied to Γ 1 , where (Γ 0 , Γ 1 ) is a partition of ∂Ω and Γ 0 , Γ 1 have positive measures.

Remarks on the model
The evolution of a structure made of shape memory alloys, i.e., the computation of ) depending on the point x of the domain Ω occupied by the structure and on time t, can be performed by solving numerically the set of partial differential equations resulting from the equations of motion (2.5), (2.6) the energy balance (2.27), (2.28) and the constitutive laws (2.20)-(2.25), completed by convenient initial and boundary conditions (cf., e.g., [30], [17], [61], [53]). The model we have described here is able to account for the different features of the shape memory alloys: in particular, their macroscopic, mechanical and thermal properties. We have used schematic free energies and schematic pseudopotentials of dissipation.
There are still many possibilities to upgrade the basic choices we have made to take into account the practical properties of shape memory alloys. Let us, for instance, mention that the pseudopotential of dissipation can be modified in order to describe more precisely the hysteretical properties of the materials. There is no difficulty in having more than two martensites, for instance, to take care of 24 possible martensites! In the same way, it is possible to take into account of the different forms of a single martensite variant, as explained in [55].
Note that the physical quantities for characterizing an educated shape memory alloys are K ,C , l a , ϑ 0 , ϑ c , τ , the two martensite volume fractions (for the free energy) and c, k , λ (for the pseudopotential of dissipation). They are indeed not so many in order to have a complete multidimensional model which can be used for engineering purposes.
Other models and results may be found in [4,5,10,45,52,44]. Let us also note the very important role of internal constraints and of the reaction B ndr to those internal constraints which are responsible for many properties.

Main results
In order to give a precise formulation of our problem, let us denote by Ω a bounded, convex set in R n (n = 1, 2, 3) with Lipschitz boundary Γ, by T a positive final time, and by Q the space-time cylinder Ω × (0, T ). Let (Γ 0 , Γ 1 ) be a partition of ∂Ω into two measurable sets such that both Γ 0 and Γ 1 have positive surface measure. Finally, denote by Σ : and identify, as usual, H (which stands either for the space L 2 (Ω) or for (L 2 (Ω)) 2 or for (L 2 (Ω)) 3 ) with its dual space H ′ , so that V ֒→ H ֒→ V ′ with dense and continuous embeddings. Moreover, we denote by · X the norm in some space X and by ·, · the duality pairing between V and V ′ and by (·, ·) the scalar product in H . Set, for simplicity of notation and without any loss of generality Then, in order to write the variational formulation of our problem (2.35-2.39), we need to generalize the relationship B ndr ∈ ∂I C (β) stated in (2.38) (cf. also [8] for similar generalizations). Hence, we need to introduce the following ingredients j : R 3 → [0, +∞] a proper, convex, lower semicontinuous function such that j(0) = 0 and its subdifferential (3.2) Moreover, we consider the associate functionals with their subdifferentials (cf. [7, Chap. II, p. 52]) Now we denote by W the following space endowed with the usual norm. In addition, we introduce on W × W a bilinear symmetric continuous form a(·, ·) defined by Note here that (since Γ 0 has positive measure), thanks to Korn's inequality (cf., e.g., [20], [32, p. 110]), there exists a positive constant c such that Next, in order to rewrite the problem (2. 35-2.39) in an abstract framework, let us introduce the operators Moreover, we make the following assumptions on the data Then, we introduce the functions R ∈ L 2 (0, T ; V ′ ) and F ∈ W 1,1 (0, T, W ′ ) such that and take J V as in (3.7), then we are ready to introduce the variational formulation of our problem as follows.
Problem (P). Given ϑ 0 > 0 find (u, w, β 1 , β 2 , β 3 ) and (ξ 1 , ξ 2 , ξ 3 , p) with the regularities . The difficult point in the proof of our result will be in fact the passage to the limit in the two non-smooth nonlinearities in equation (3.31). By the contrary, we aim to remark that this viscous term in (3.31) gives more spatial regularity to ∂ t β which furnish more regularity to w in (3.30) (cf. the following Theorem 3.2) and consequently to ϑ = γ(w) in (3.31). This regularity is needed in order to prove well-posedness for our problem. However, we can also notice that, from the mechanical viewpoint, since we have introduced in the model the elastic (non-dissipative) local interaction term −∆β (in (2.36)), it seems also reasonable to include in the model the dissipative local interaction term −∆β t .
We are now ready to state our main result which is the following global existence and uniqueness theorem. Moreover, in the last Section 5, we will get a proof for the following continuous dependence result for Problem (P). Theorem 3.3. Let T be a positive final time, (u i 0 , w i 0 , β i 0 ) (i = 1, 2) be two sets of initial data satisfying conditions (3.16-3.18), (R i , F i ) (i = 1, 2) be two data of Problem (P) satisfying assumptions (3.21-3. 22) with (f i , g i , R i , Π i ) (i = 1, 2) as in (3.19-3.20). Let (u i , w i , β i 1 , β i 2 , β i 3 ) (i = 1, 2) be two solutions of Problem (P) corresponding to these data. Moreover, besides conditions (3.2-3.3), suppose that the following hypothesis holds. Then, there exists a positive constant M , depending on the data of the problem, such that the following continuous dependence estimate holds for any t ∈ (0, T ).
Remark 3.4. Let us note that in this paper we can treat the difficult coupling between the phase-equations (3.31) in which appears the temperature ϑ (= γ(w)) and the entropy balance equation (3.30) in which only the function log ϑ (= w ) plays some role, using the L ∞ -bound on the w -component of solution to Problem (P) (cf. (3.25)). Indeed it was just due to the lack of regularity of solutions that in [13] (where there was the same type of coupling without the ∆∂ t β -term in (3.31)) the authors did not obtain uniqueness of solutions (cf. also [13,Remark 5.2]). However, let us, finally, observe that the main advantage of taking the entropy balance equation instead of the internal energy balance equation is that once one has solved the problem in some sense and has found the temperature ϑ := γ(w), it is automatically positive because it stands in the image of the function γ (cf. (3.23)). Indeed in many cases it is difficult to deduce this fact only from the internal energy balance equation (cf., e.g., [29] in order to see one example of these difficulties). Let us note that within the small perturbations assumption the entropy balance and the classical heat equation are equivalent in mechanical terms (cf. [13,14,15]).

Proof of Theorem 3.2
The following section is devoted to the proof of Theorem 3.2. First we approximate our Problem (P) by a more regular Problem (P ε ), then (fixed ε > 0) we find well-posedness for the approximating problem using a iterated Banach contraction fixed-point argument and then we perform some a-priori estimates (independent of ε) on its solution, which allow us to pass to the limit in Problem (P ε ) as ε ց 0, recovering a solution to Problem (P).

The approximating problem
We take a small positive parameter ε > 0 and approximate ∂ V,V ′ J V in (3.31-3.32) , let us take the Lipschitz continuous Yosida-Moreau approximation α ε = ∂j ε = (j ε ) ′ (cf. [16,Prop. 2.11,p. 39]) of α and the associated functional J H,ε (v) = Ω j ε (v(x)) dx, whose differential (∂ H J H,ε ) is the Yosida-Moreau approximation of ∂ H J H (cf. [16,Prop. 2.16,p. 47]). Now we aim to recall some properties of this approximation which will be useful in order to pass to the limit as ε ց 0. Note that the proof of the following lemma is a consequence of holds true. Moreover, the following properties hold true (for ε ց 0) (4.5) Finally, for ε ց 0, it holds Then, let us call γ ε the following Lipschitz continuous approximation of the function γ(w) = exp(w), i.e. the function Moreover let δ ε be the inverse function of γ ε , i.e.
In this way, we have defined an operator T : X → X such that (β ε 1 , β ε 2 , β ε 3 ) =: T (β 1 ε ,β 2 ε ,β 3 ε ). What we have to do now is to prove that T is a contraction mapping on X for a sufficiently smallt ∈ [0, T ] and moreover, repeating the procedure step by step in time (this is possible thanks to the regularities properties of the solution listed above), we can prove well-posedness for the Problem (P ε ) on the whole time interval [0, T ] and conclude the proof of Theorem 4.3. In order to prove that T is contractive, let us proceed by steps and forget of the apices ε. 1, 2, 3). Then, writing two times (4.15) withβ 3 i (i = 1, 2) instead of ∂ t (β 3 ), making the difference, testing the resulting equation with (w 1 − w 2 ) t , and integrating on (0, t) with t ∈ [0, T ], we get the following inequality 20) for some positive constant C 0 independent of t. Hence, we get, for all t ∈ [0, T ], being C 1 := C 0 e T /2 . Second step. Let us take ϑ i = γ ε (w i ) and εp i := ∂ t (β i 1 + β i 2 + β i 3 ) + div (u i ) t (i = 1, 2) and write (4.13) with β i 1 and β i 2 , make the difference, test the resulting equation with (u 1 − u 2 ) t , integrate on (0, t) with t ∈ [0, T ], and use equation (4.14), getting the following inequality Third step. Write equation (4.16) for ϑ i and u i , make the difference between the two equations written for i = 1 and i = 2, and test the resulting vectorial equation by the vector (β 1 Moreover, using the definition (4.2) of γ ε and the assumption (3.4) on τ , we get the following inequality Fourth step. Summing up the two inequalities (4.22) and (4.23), two integrals cancels out, and using (4.20), we get where C ε does not depend on t. Hence, choosing t sufficiently small (this is ourt), we recover the contractive property of T . Moreover, applying the Banach fixed point theorem to T , we get a unique solution for the Problem (P ε ) on the time interval [0,t]. Due to this estimate it is easy to prove that there exists m ∈ N such that T m is a contraction on X . Hence we have a unique solution on the whole time interval [0, T ]. This concludes the proof of Theorem 4.3.

A priori estimates
In this subsection we perform a-priori estimates on Problem (P ε ) uniformly in ε which will lead us pass to the limit as ε ց 0 and recover a solution of Problem (P).
We denote by c all the positive constants (which may also differ from line to line) independent of ε and depending on the data of the problem. For simplicity, we omit the subscript ε when it is not necessary.
We want to apply this result in order to find the uniform (in ε) bound on p ε . First of all let us note that from comparison in (4.13), using also the bound (4.28) on u with the assumption (3.19) on F , we immediately deduce that Moreover, always by comparison in (4.13), we have that Following the idea of [40], we can choose v ⋆ ∈ W such that Note that, since Ω is regular (it suffices for Ω to be a Lipschitz domain), we can always find a v ⋆ ∈ W such that (4.36) is satisfied, because, if we take B ε (x) the ball in R 3 centered in x ∈ Γ 1 with radius ε such that B ε (x) ∩ Γ 0 = ∅ and consider the parametrization of Γ 1 ∩ B ε (x) through the Lipschitz function (x 1 , x 2 ) → (x 1 , x 2 , ϕ(x 1 , x 2 )), then the normal unit vector associated is Then, if we take v ⋆ = (0, 0, ζ) with then v ⋆ ∈ W and moreover we can show that (4.36) holds because where L is the Lipschitz constant of ϕ, and hence Take now m in Lemma 4.4 as Then, m(v) is a seminorm on H and a norm on the constants because of (4.36). Hence, we can apply Lemma 4.4 to p ε with the choices done above and, thanks to (4.34-4.35), we get the bound p ε L 2 (Q) ≤ c . Finally, by comparison in (4.16) and using the estimates (4.28) and (4.32-4.37), we deduce that also ∂ H J H,ε (β) is bounded in L 2 (0, T ; V ′ ). Then, testing (4.16) with Bβ ε and then by using again (4.28) and (4.32-4.37) and the monotonicity properties of α ε , we get also β ε j L ∞ (0,T ;W 2,2 (Ω)) ≤ c (j = 1, 2, 3). (4.38) Now it remains only to pass to the limit in (4.13-4.16) as ε ց 0. This will be the aim of the next subsection.

Passage to the limit and uniqueness
As we have just mentioned, we want to conclude the proof of Theorem 3.2 passing to the limit in the well-posed (cf. Subsection 4.1) Problem (P ε ) as ε ց 0 using the previous uniform (in ε) estimates on its solution (cf. Subsection 4.2) and exploiting some compactness-monotonicity argument. Let us list before the weak or weak-star convergence coming directly from the previous estimates and well-known weak-compactness results. Note that the following convergences hold only up to a subsequence of ε ց 0 (let us say ε k ց 0). We denote it again with ε only for simplicity of notation. From the estimates (4.28-4.38) and the property (4.1) of ∂ H J H,ε , we deduce that and weakly star in Note that (4.46-4.47) imply immediately the convergence γ ε (w ε ) → ϑ = γ(w) and τ (γ ε (w ε )) → τ (ϑ) a.e. in Q .
Moreover, the two convergences (4.45) and (4.48) along with the property (4.6) and [2, Thm. 3.66, p. 373] give immediately the identification of the maximal monotone graph 3 and a.e. in [0, T ] with ξ = (ξ 1 , ξ 2 , ξ 3 ) and ξ j (j = 1, 2, 3) are the weak limits defined in (4.45). All these convergences with the identifications made above make us able to pass to the limit (as ε ց 0 or at least for a subsequence of it) in Problem (P ε ) finding a solution to Problem (P) and concluding in this way the proof of Theorem 3.2. Note that the convergences hold for all subsequences ε k of ε tending to 0 because of uniqueness of solutions. Indeed we may prove it in this way.

Proof of Theorem 3.3
In this section we give the proof of Theorem 3.3. We will use here the same symbol c for some positive constants (depending only on the data of the problem), which may also be different from line to line.
Then, let us take two sets of data (u i 0 , w i 0 , β i 0 ), (R i , F i ) (i = 1, 2) of Problem (P) and let (u i , w i , β i 1 , β i 2 , β i 3 , p i ) (i = 1, 2) be two solutions of Problem (P) corresponding to these data.
Then, write two times equations (3.28-3.31) with (u i , w i , β i 1 , β i 2 , β i 3 , p i ), make the difference between the two equations (3.30), and test the result with 2(w 1 − w 2 ). Make the difference between the two equations (3.28), test the result with 2(u 1 − u 2 ) t . Make the difference between the two equations (3.31), written for i = 1 and i = 2, and test the resulting vectorial equation by the vector (β 1 Finally, summing up the three resulting equations, integrating over (0, t), with t ∈ [0, T ], and exploiting the Lipschitz continuity of α (cf. assumption (3.36)), we get the following inequality Let us notice that we have estimated the terms containing the nonlinearity τ in (3.28) and (3.31) on the right hand side using (3.29) and the fact that γ , defined in (3.23), is a locally Lipschitz continuous function, ϑ i = γ(w i ) (i = 1, 2), and w i are bounded in L ∞ (Q) (cf. (3.25)), as follows Moreover, by adding to both sides in the inequality (5.1) we get the following inequality Applying now a standard version of Gronwall's lemma (cf. [16,Lemme A.4,p. 156]), we get the desired continuous dependence estimate (3.37). This concludes the proof of Theorem 3.3.