An existence result for the fractional Kelvin-Voigt's model on time-dependent cracked domains

We prove an existence result for the fractional Kelvin-Voigt's model involving Caputo's derivative on time-dependent cracked domains. We first show the existence of a solution to a regularized version of this problem. Then, we use a compactness argument to derive that the fractional Kelvin-Voigt's model admits a solution which satisfies an energy-dissipation inequality. Finally, we prove that when the crack is not moving, the solution is unique.


Introduction
This paper deals with the mathematical analysis of the dynamics of elastic damping materials in the presence of external forces and time-dependent brittle fracture. In this framework, it is important to find the behavior of the deformation when the crack evolution is known. This is the first step towards the development of a complete model of dynamic crack growth in viscoelastic materials. From a mathematical point of view, this means solving the following dynamic system u(t) − div(σ(t)) = f (t) in Ω \ Γ t , t ∈ (0, T ). (1.1) In the equation above, Ω ⊂ R d represents the reference configuration of the material, the set Γ t ⊂ Ω models the crack at time t (which is prescribed), u(t) : Ω \ Γ t → R d is the displacement of the deformation, σ(t) the stress tensor, and f (t) is the forcing term.
In the classical theory of linear viscoelasticity, the constitutive stress-strain relation of the so called Kelvin-Voigt's model is given by σ(t) = Ceu(t) + Beu(t) in Ω \ Γ t , t ∈ (0, T ), (1.2) where C and B are two positive tensors acting on the space of symmetric matrices, and ev denotes the symmetric part of the gradient of a function v (which is defined as ev := 1 2 (∇v + ∇v T )). The local model associated to (1.2) has already been widely studied and we can find several existence results in the literature; we refer to [2,3,6,7,16,23] for existence and uniqueness results in the pure elastodynamics case (B = 0) and in the classic Kelvin-Voigt's one.
In recent years, materials whose constitutive equations can be described by non-local models are of increasing interest. In this context, by non-local we mean that the state of the stress at instant t depends not only on that instant, but also on the previous ones (long memory). For solid viscoelastic materials, some experiments are particularly in agreement with models using fractional derivative, see for example [10,11,22,24] and the reference therein.
In this paper, we focus on the fractional Kelvin-Voigt's model, i.e. we consider the following constitutive stress-strain relation σ(t) = Ceu(t) + BD α t eu(t) in Ω \ Γ t , t ∈ (0, T ), where D α t denotes a fractional derivative of order α ∈ (0, 1). In the literature we can find several definitions for the fractional derivative of a function g : (a, b) → R; here we focus on the most used ones which are Riemann-Liouville's derivative of order α at starting point a (r) (t − r) α dr. We recall that Γ denotes Euler's Gamma function; notice that in order to define Caputo's derivative the function g must be differentiable, while this is not necessary for Riemann-Liouville's derivative. Given g ∈ AC([a, b]), and t ∈ (a, b) we have the following relation between Riemann-Liouville's and Caputo's derivative (see, e.g., [12]): In particular, when g(a) = 0, these two notions coincide. For more properties regarding these two fractional derivatives, we refer for example to [4,15,19,20] and the references therein.
In this paper we use Caputo's derivative, which means we consider the dynamic system One of the quality of this definition for the fractional derivative is that the initial conditions can be imposed in the classical sense, see for example [15,19]. The choice of 0 as a starting point is due to the fact that we want to couple dynamic system (1.1) with the initial conditions at time t = 0. Dealing with (1.4) is very difficult, since in the definition of C 0 D α t eu(t) we need that eu is differentiable, which is a very strong request. Hence, we rephrase Caputo's derivative in a more suitable way. Thanks to (1.3) for g ∈ AC([0, T ]) we can write (t − r) α (g(r) − g(0)) dr. (1.5) This formulation of Caputo's derivative is well-posed in the distributional sense also when the function g is only integrable. We point out that formula (1.5) can be found in the recent literature on fractional derivatives, where it is used to define the notion of weak Caputo's derivative for less regular functions, see for example [9,14]. Thanks to formula (1.5), we can write system (1.4) in a weaker form (see Definition 2. (1.7) Notice that the scalar function ρ appearing in F is positive, decreasing, and convex on (0, ∞). Moreover, ρ ∈ L 1 (0, T ) for every T > 0, but it is not bounded on (0, T ). In particular, we can not compute the derivative in front of the convolution integral in (1.6). When there is no crack, existence results for this kind of system can be found for example in [1,5,13,18]. However, in the case of a dynamic fracture, the techniques used in the previous papers can not be applied and up to now there are no existence results in this setting.
To prove the existence of a solution to (1.6) we proceed into two steps. First we consider a regularized version of (1.6), where we replace the kernel F in (1.6) by a regular kernel G ∈ C 2 ([0, T ]). Then we prove the existence of a solution to the more regular system in Ω \ Γ t , t ∈ (0, T ), (1.8) and we show that this solution satisfies a uniform bound depending on the L 1 -norm of G. Finally, we consider a sequence of regular tensors G ǫ converging to F in L 1 and we take the solutions to (1.8) with G := G ǫ . By a compactness argument, we show that the sequence u ǫ converge to a function u * which solves (1.6). Moreover, we prove that this solution satisfies an energy-dissipation inequality. We conclude this paper by showing that, when the crack is not moving, the fractional Kelvin-Voigt's system (1.6) admits a unique solution. The paper is organized as follows: in Section 2 we fix the notation and the framework of our problem. Moreover, we give the notion of solution to the fractional Kelvin-Voigt's system involving Caputo's derivative (1.6) and we state our main existence result (see Theorem 2.4). Section 3 deals with the regularized system (1.8). First, by a time-discretization procedure in Theorem 3.13 we prove the existence of a solution to (1.8). Then, in Lemma 3.14 we derive the uniform energy estimate which depends on the L 1 -norm of G. In Section 4 we consider Kelvin-Voigt's system (1.6): we prove the existence of a generalized solution to system (1.6) and in Theorem 4.2 we show that such a solution satisfies an energy-dissipation inequality. Finally, in Section 5 we prove that, for a not moving crack, the solution to (1.6) is unique.

Notation and framework of the problem
The space of m × d matrices with real entries is denoted by R m×d ; in case m = d, the subspace of symmetric matrices is denoted by R d×d sym . Given a function u : R d → R m , we denote its Jacobian matrix by ∇u, whose components are (∇u) ij := ∂ j u i for i = 1, . . . , m and j = 1, . . . , d; when u : R d → R d , we use eu to denote the symmetric part of the gradient, namely eu := 1 2 (∇u + ∇u T ). Given a tensor field A : R d → R m×d , by div A we mean its divergence with respect to rows, namely (div We denote the d-dimensional Lebesgue measure by L d and the (d − 1)-dimensional Hausdorff measure by H d−1 ; given a bounded open set Ω with Lipschitz boundary, by ν we mean the outer unit normal vector to ∂Ω, which is defined H d−1 -a.e. on the boundary. The Lebesgue and Sobolev spaces on Ω are defined as usual; the boundary values of a Sobolev function are always intended in the sense of traces.
The norm of a generic Banach space X is denoted by · X ; when X is a Hilbert space, we use (·, ·) X to denote its scalar product. We denote by X ′ the dual of X and by ·, · X ′ the duality product between X ′ and X. Given two Banach spaces X 1 and X 2 , the space of linear and continuous maps from X 1 to X 2 is denoted by L (X 1 ; X 2 ); given A ∈ L (X 1 ; X 2 ) and u ∈ X 1 , we write Au ∈ X 2 to denote the image of u under A.
Moreover, given an open interval (a, b) ⊆ R and p ∈ [1, ∞], we denote by L p (a, b; X) the space of L p functions from (a, b) to X; we use W k,p (a, b; X) and H k (a, b; X) (for p = 2) to denote the Sobolev space of functions from (a, b) to X with k derivatives. Given u ∈ W 1,p (a, b; X), we denote byu ∈ L p (a, b; X) its derivative in the sense of distributions. When dealing with an element u ∈ W 1,p (a, b; X) we always assume u to be the continuous representative of its class; in particular, it makes sense to consider the pointwise value u(t) for every t ∈ [a, b]. We use C 0 w ([a, b]; X) to denote the set of weakly continuous functions from [a, b] to X, namely, the collection of maps u : Let T be a positive real number and let Ω ⊂ R d be a bounded open set with Lipschitz boundary. Let ∂ D Ω be a (possibly empty) Borel subset of ∂Ω and let ∂ N Ω be its complement. Throughout the paper we assume the following hypotheses on the geometry of the cracks: (H1) Γ ⊂ Ω is a closed set with L d (Γ) = 0 and H d−1 (Γ ∩ ∂Ω) = 0; (H2) for every x ∈ Γ there exists an open neighborhood U of x in R d such that (U ∩ Ω) \ Γ is the union of two disjoint open sets U + and U − with Lipschitz boundary; (H3) {Γ t } t∈[0,T ] is an increasing family in time of closed subsets of Γ, i.e. Γ s ⊂ Γ t for every 0 ≤ s ≤ t ≤ T . Thanks (H1)-(H3) the space L 2 (Ω \ Γ t ; R m ) coincides with L 2 (Ω; R m ) for every t ∈ [0, T ] and m ∈ N. In particular, we can extend a function u ∈ L 2 (Ω \ Γ t ; R m ) to a function in L 2 (Ω; R m ) by setting u = 0 on Γ t . To simplify our exposition, for every m ∈ N we define the spaces H := L 2 (Ω; R m ), H N := L 2 (∂ N Ω; R m ) and H D := L 2 (∂ D Ω; R m ); we always identify the dual of H by H itself, and L 2 ((0, T ) × Ω; R m ) by the space Notice that in the definition of U t we are considering only the distributional gradient of u in Ω\Γ t and not the one in Ω. By (H2) we can find a finite number of open sets U j ⊂ Ω \ Γ, j = 1, . . . m, with Lipschitz boundary, such that Ω \ Γ = ∪ m j=1 U j . By using second Korn's inequality in each U j (see, e.g., [17,Theorem 2.4]) and taking the sum over j we can find a constant C K , depending only on Ω and Γ, such that where eu is the symmetric part of ∇u. Therefore, we can use on the space U t the equivalent norm Furthermore, the trace of u ∈ H 1 (Ω \ Γ; R d ) is well defined on ∂Ω. Indeed, we may find a finite number of open sets with Lipschitz boundary V k ⊂ Ω \ Γ, k = 1, . . . l, such that ∂Ω \ (Γ ∩ ∂Ω) ⊂ ∪ l k=1 ∂V k . Since H d−1 (Γ ∩ ∂Ω) = 0, there exists a constant C, depending only on Ω and Γ, such that Hence, we can consider the set U D t := {u ∈ U t : u = 0 on ∂ D Ω} for every t ∈ [0, T ], which is a closed subspace of U t . Moreover, there exists a positive constant C tr such that u HN ≤ C tr u UT for every u ∈ U T . Now, we define the following sets of functions T ]}, in which we develop our theory. Moreover, we consider the Banach space represents the space of symmetric tensor fields, i.e. the collections of linear and continuous maps A : R d×d sym → R d×d sym satisfying Aξ · η = Aη · ξ for every ξ, η ∈ R d×d sym . We assume that the Dirichlet datum z, the Neumann datum N , the forcing term f , the initial displacement u 0 , and the initial velocity u 1 satisfy We consider a coercive tensor C ∈ B, which means that there exists γ > 0 such that C(x)ξ · ξ ≥ γ|ξ| 2 for every ξ ∈ R d and a.e. x ∈ Ω. (2.4) Moreover, let us take a time-dependent tensor F : (0, T + δ 0 ) → B, with δ 0 > 0, satisfying
Since in our existence result we first regularize the tensor F by means of translations (see Section 4) we need that F is defined also on the right of T . This is not a problem, because our standard example for F, which is (1.7), is defined on the whole (0, ∞).
In this paper we want to study the following problem (2.9) We give the following notion of solution to system (2.9): Definition 2.2 (Generalized solution). Assume (2.1)-(2.8). A function u ∈ C w is a generalized solution to system (2.9) if u(t) − z(t) ∈ U D t for every t ∈ [0, T ], u(0) = u 0 in U 0 ,u(0) = u 1 in H, and for every ϕ ∈ C 1 c the following equality holds Remark 2.3. The Neumann conditions appearing in (2.9) are only formal; they are used to pass from the strong formulation in (2.9) to the weak one (2.10).
The main existence result of this paper is the following theorem: Then there exists a generalized solution u ∈ C w to system (2.9).
The proof of this theorem requires several preliminary results. First, in the next section, we prove the existence of a generalized solution when the tensor F is replaced by a tensor G ∈ C 2 ([0, T ]; B). Then, we show that such a solution satisfies an energy estimate, which depends via G only by its L 1 -norm. In Section 4 we combine these two results to prove Theorem 2.4.

The regularized model
In this section we deal with a regularized version of the system (2.9), where the tensor F is replaced by a tensor G which is bounded at t = 0. More precisely, we consider the following system and we assume that G :
in H, and for every ϕ ∈ C 1 c the following equality holds Since the time-dependent tensor G is well defined in t = 0, we can give another notion of solution. In particular, the convolution integral is now differentiable, and we can write  (3.6). A function u ∈ C w is a weak solution to system (3.1) if u(t) − z(t) ∈ U D t for every t ∈ [0, T ], u(0) = u 0 in U 0 ,u(0) = u 1 in H, and for every ϕ ∈ C 1 c the following equality holds In this framework the two previous definitions are equivalent. Proof. We only need to prove that (3.8) is equivalent to (3.7). This is true if and only if the function u ∈ C w satisfies for every ϕ ∈ C 1 c the following equality Let us consider for t ∈ [0, T ] the function we deduce that p ∈ Lip([0, T ]). In particular, there existsṗ(t) for a.e. t ∈ (0, T ). Given t ∈ (0, T ) and h > 0 we can write Let us compute these three limits separately. We claim that for a.e. t ∈ (0, T ) we have Indeed, by the Lebesgue's differentiation theorem, for a.e. t ∈ (0, T ) we get Finally, for every t ∈ (0, T ) we get Therefore, by the identity and the previous computations we deduce (3.9).
In the particular case in which the tensor G appearing in (3.1) is the one associated to the Standard viscoelastic model, i.e.
with β > 0 and B ∈ B non-negative tensor, then the existence of weak solutions (and so generalized solutions) was proved in [21]. Here we adapt the techniques of [21] to a general tensor G satisfying (3.2)-(3.5).
3.1. Existence and energy-dissipation inequality. In this subsection we prove the existence of a generalized solution to system (3.1), by means of a time discretization scheme in the same spirit of [6]. Moreover, we show that such a solution satisfies the energy-dissipation inequality (3.40). We fix n ∈ N and we set τ n := T n , u 0 n := u 0 , u −1 n := u 0 − τ n u 1 , δz 0 n :=ż(0), δG 0 n := 0. Let us define U j n := U D jτn , z j n := z(jτ n ), G j n := G(jτ n ) for j = 0, . . . , n, Regarding the forcing term and the Neumann datum we pose For every j = 1, . . . , n let us consider the unique u j n ∈ U T with u j n − z j n ∈ U j n , which satisfies The existence and uniqueness of u j n is a consequence of Lax-Milgram's lemma. Notice that equation (3.10) is a sort of discrete version of (3.8), which we already know that is equivalent to (3.7).
We now use equation (3.10) to derive an energy estimate for the family {u j n } n j=1 , which is uniform with respect to n ∈ N. Proof. First, we notice Therefore, equation (3.10) can be written as . . , n}. By taking v := τ n (δu j n − δz j n ) ∈ U j n and summing over j = 1, . . . , i, we get the following identity By using the identity Moreover, we can write Finally, let us consider the term We can write By combining together (3.12)-(3.16), we obtain for i = 1, . . . , n the following discrete energy equality By our assumptions on G we deduce x ∈ Ω and every ξ ∈ R d and j = 0, . . . , n, x ∈ Ω and every ξ ∈ R d and j = 1, . . . , n, x ∈ Ω and every ξ ∈ R d and j = 2, . . . , n.
Hence, thanks to (3.17), for every i = 1, . . . , n we can write Let us estimate the right-hand side in (3.18) from above. We set Therefore, we have the following bounds Notice that the following discrete integrations by parts hold . . , n. In particular, since the right-hand side is independent of i, u 0 n = u 0 and δu 0 n = u 1 , there exists another constant C 2 = C 2 (z, N, f, u 0 , u 1 , C, G) for which we have K 2 n + E 2 n ≤ C 2 (1 + K n + E n ) for every n ∈ N. This implies the existence of a constant C = C(z, N, f, u 0 , u 1 , C, G) independent of n ∈ N such that δu j n H + eu j n H ≤ K n + E n ≤ C for every j = 1, . . . , n and n ∈ N, which gives (3.11).
A first consequence of Lemma 3.4 is the following uniform estimate on the family {δ 2 u j n } n j=1 . Proof. Thanks to equation (3.10) and to Lemma 3.4, for every j = 1, . . . , n and v ∈ U D 0 ⊆ U j n with v U0 ≤ 1 we have HN ) . By taking the supremum over v ∈ U D 0 with v U0 ≤ 1 we obtain HN ) . We multiply this inequality by τ n and we sum over j = 1, . . . , n to get (3.31).
We now want to pass to the limit into equation (3.10) to obtain a generalized solution to system (3.1). Let us recall the following result, whose proof can be found for example in [8].
Lemma 3.6. Let X, Y be two reflexive Banach spaces such that X ֒→ Y continuously. Then for t ∈ [(i − 1)τ n , iτ n ) and i = 1, . . . , n,ũ − n (T ) = δu n n , which approximate the first time derivative of the generalized solution. In a similar way, we define also f + n , N + n ,Ñ + n , z ± n ,z n ,z + n , G ± n ,G n ,G + n . Thanks to the uniform estimates of Lemma 3.4 we derive the following compactness result: and for every t Proof. Thanks to Lemma 3.4 and the estimate (3.31), the sequences Again, thanks to (3.11) and (3.31), for every t ∈ [0, T ] we get which imply (3.33). Finally, observe that for every t ∈ [0, T ] Therefore, u(t) ∈ U t for every t ∈ [0, T ] since U t is a closed subspace of U T . Hence, u ∈ C w . Let us check that the limit function u defined before satisfies the boundary and initial conditions.
Proof. We only need to prove that the function u ∈ C w satisfies (3.7). We fix n ∈ N and a function ϕ ∈ C 1 c . Let us consider ϕ j n := ϕ(jτ n ) for j = 0, . . . , n, δϕ j n := ϕ j n − ϕ j−1 n τ n for j = 1, . . . , n, and, as we did before for the family {u j n } n j=1 , we define the approximating sequences {ϕ + n } n and {φ + n } n . If we use τ n ϕ j n ∈ U j n as a test function in (3.10), after summing over j = 1, ..., n, we get where t n := † t τn £ τ n for t ∈ (0, T ) and ⌈x⌉ is the superior integer part of the number x. Thanks to (3.35) we deduce We use (3.32) and the following convergences Moreover, for every fixed t ∈ (0, T ) which together with (3.32) gives By (3.11) for every t ∈ (0, T ) we deduce Therefore, we can use the dominated convergence theorem to pass to the limit in the double integral of (3.36), and we obtain that u satisfies (3.7) for every function ϕ ∈ C 1 c . Now we want to deduce an energy-dissipation inequality for the generalized solution u ∈ C w of Lemma 3.7. Let us define for every t ∈ [0, T ] the total energy E(t) and the dissipation D(t) as Notice that E(t) is well defined for every time t ∈ [0, T ] since u ∈ C 0 w ([0, T ]; U T ) andu ∈ C 0 w ([0, T ]; H). Moreover, by the initial conditions we have  6). Then the generalized solution u ∈ C w to system (3.1) of Lemma 3.7 satisfies for every t ∈ [0, T ] the following energy-dissipation inequality where the total work is defined as . After setting t n := iτ n and using that δG 0 n = 0, we rewrite (3.17) as where r n := † r τn £ τ n for r ∈ (τ n , t n ), and the approximate total work W + n (t) is given by [(f + n (r),ũ + n (r) −z + n (r)) H + (N + n (r),ũ + n (r) −z + n (r)) HN + (u n (r),z + n (r)) H ] dr + tn 0 [(Ceu + n (r), ez + n (r)) H + (G − n (r)(eu + n (r) − eu 0 ), ez + n (r)) H ] dr + tn τn rn−τn 0 (G − n (r n − s)(eu + n (s) − eu + n (r)), ez + n (r)) H ds dr.

By (2.4), (3.3), and (3.33) we derive
Moreover, the estimate (3.11) imply Hence, we can argue as before to deduce In particular, we can use Fatou's lemma and the fact that t ≤ t n to obtain By arguing in a similar way, we can derive Let us consider the double integral in the left-hand side. We fix r ∈ (0, t) and by (3.5) for every s ∈ (0, r) we have Moreover, for a.e. s ∈ (0, r n − τ n ) by defining s n := † s τn Therefore, for a.e. s ∈ (0, r) we get (G(r − s)(eu(r) − eu(s)), eu(r) − eu(s)) H ≤ lim inf n→∞ (Ġ n (r n − s)(eu + n (r) − eu + n (s)), eu + n (r) − eu + n (s)) H , since s ∈ (0, r n − τ n ) for n large enough. If we apply again Fatou's lemma we conclude from which we derive (Ġ n (r n − s)(eu + n (r) − eu + n (s)), eu + n (r) − eu + n (s)) H ds.
Let us study the right-hand side of (3.42). Given that By arguing as before we deduce (Ṅ (r), u(r) − z(r)) HN dr, (3.52) thanks to Lemma 3.7 and to the following convergences:  Unfortunately, W bdry (t) is not well defined under our assumptions on u. In particular, the term involving the Dirichlet datum z is difficult to handle since the trace of the function Ceu(r)ν + d dr r 0 G(r − s)eu(s)ds ν on ∂ D Ω is not well defined. If we assume that u ∈ L 2 (0, T ; H 2 (Ω \ Γ; R d )) ∩ H 2 (0, T ; L 2 (Ω \ Γ; R d )) and that Γ is a smooth manifold, then the first term of W bdry (t) makes sense and satisfies  We conclude this subsection by showing that the generalized solution of Lemma 3.7 satisfies the initial conditions in a stronger sense than the ones stated in Definition 2.2.
It remains to study the last two terms, which are Therefore, since 0,T ;B) ). In particular, being the right-hand side independent of t ∈ [0, T ], we conclude This implies the existence of a constant M = M (C 0 , C 1 , C 2 ) for which (3.57) is satisfied. By equation (3.7) it is easy to see thatu ∈ H 1 (0, T ; (U D 0 ) ′ ) and thatü satisfies for a.e. t ∈ (0, T ) and for Hence, we derive . Therefore the bounds onü depends on G(0) B even when z ∈ W 2,1 (0, T ; U 0 ).
As explained in the previous remark, we can not deduce a uniform bound foru in H 1 (0, T ; (U D 0 ) ′ ) depending on G only via its L 1 -norm. On the other hand, the bound onu in H 1 (0, T ; (U D 0 ) ′ ) is useful if we want to prove the existence of a generalized solution u * to the fractional Kelvin-Voigt system (2.9), especially to show thatu * ∈ C 0 w ([0, T ]; H). To overcome this problem, we introduce another function that is related tou and for which is possible to derive a uniform bound. Let us consider the auxiliary function α : we have for every v ∈ U D 0 the following convergence The second convergence is true because Clearly For this function α is possible to find a uniform bound in H 1 (0, T ; (U D 0 ) ′ ) which depends on G L 1 (0,T ;B) . Proof. First, by Lemma 3.14 we have ) for every t ∈ [0, T ]. Moreover, by the definition of generalized solution, we deduce that for every ψ ∈ C 1 c (0, T ) and v ∈ U D 0 it holds This gives that there existsα ∈ L 2 (0, T ; (U D 0 ) ′ ) and , v) HN for every v ∈ U D 0 and for a.e. t ∈ (0, T ).

The fractional Kelvin-Voigt's model
In this section we prove the existence of a generalized solution to (2.9) for a tensor F which is not necessary bounded at t = 0, as it happens in (1.7). Here, we assume that our data z, N, f, u 0 , u 1 , C, and F satisfy the conditions (2.1)-(2.8). To prove the existence of a generalized solution to (2.9) under these assumptions, we first regularize F by a parameter ǫ > 0 and we consider system (3.1) associated to this regularization. Then, we take the solution u ǫ given by Theorem 3.13 and thanks to Lemma 3.14 and Corollary 3.16 we obtain a generalized solution to (2.9).
Let us regularize F by defining For every fixed ǫ ∈ (0, δ 0 ) we can consider the generalized solution u ǫ to system (3.1) with G replaced by G ǫ of Theorem 3.13. By Lemma 3.14 and Corollary 3.16 we deduce the following compactness result: For every ǫ ∈ (0, δ 0 ) let u ǫ be the generalized solution associated to system (3.1) with G replaced by G ǫ given by Theorem 3.13. Then there exists a function u * ∈ C w and a subsequence of ǫ, not relabeled, such that
We can now prove the main existence result of Theorem 2.4 for the fractional Kelvin-Voigt's system involving Caputo's derivative.
Proof of Theorem 2.4. It is enough to show that the function u * given by Lemma 4.1 is a generalized solution to (2.9). To this aim, it remains to prove that u * satisfies (2.10). For every ϕ ∈ C 1 c we know that the function u ǫ ∈ C w satisfy for every ǫ ∈ (0, δ 0 ) the following equality Let us pass to the limit as ǫ → 0 + . Clearly, by (4.1) we have It remains to study the behaviour as ǫ → 0 + of T 0 t 0 (G ǫ (t − r)(eu ǫ (r) − eu 0 ), eφ(t)) H dr dt.
We define for every ǫ ∈ (0, δ 0 ) the function v ǫ (t) : with M independent of ǫ being G ǫ L 1 (0,T ;B) ≤ F L 1 (0,T +δ0;B) . Notice that and thanks to (4.6) and to the fact that G ǫ → F in L 1 (0, T ; B) as ǫ → 0 + , we get On the other hand, since r → T r F(t − r)eφ(t) dt belongs to L ∞ (0, T ; H), we can write As a consequence, u * is a generalized solution to system (2.9).
By arguing in a similar way, we can derive For the term involvingF, we argue as we already did forḞ and by using two times Fatou's lemma we get It remains to study the right-hand side of (3.40) with the formulation of the total work as in (4.7). Thanks to Lemma 4.1 and the fact that G ǫ → F in L 1 (0, T ; B) we deduce