On the three dimensional Kelvin-Voigt fluids: global solvability, exponential stability and exact controllability of Galerkin approximations

In this work, we consider the three-dimensional viscoelastic fluid flow equations, arising from the motion of Kelvin-Voigt fluids in bounded and unbounded domains. We investigate the global solvability results, asymptotic behavior and also address some control problems of such viscoelastic fluid flow equations with "fading memory" and "memory of length \begin{document}$ \tau $\end{document} ". A local monotonicity property of the linear and nonlinear operators and a localized version of the Minty-Browder technique are used to obtain global solvability results. Since we are not using compactness arguments in the proofs, the global solvability results are also valid in unbounded domains like Poincare domains. We also remark that using an \begin{document}$ m $\end{document} -accretive quantization of the linear and nonlinear operators, one can establish the existence and uniqueness of strong solutions for the Navier-Stokes-Voigt equations and avoid the tedious Galerkin approximation scheme. We examine the asymptotic behavior of the stationary solutions and also establish the exponential stability results. Finally, under suitable assumptions on the Galerkin basis, we consider the controlled Galerkin approximated 3D Kelvin-Voigt fluid flow equations. Using the Hilbert uniqueness method combined with Schauder's fixed point theorem, the exact controllability of the finite dimensional Galerkin approximated system is established.

1. Introduction. The mathematical theory of the three-dimensional viscoelastic fluid flow equations, arising from the Kelvin-Voigt model for the non-Newtonian fluid flows and related models has been examined by various mathematicians over the past several decades starting from the works of Oskolkov (see for example [32,33,34,37,35,38,39,8,36,42,14,11,16], etc). The work under our consideration is the initial-boundary value problem for the three dimensional Navier-Stokes-Voigt equations or Oskolkov equations with "fading memory" and "memory of length τ " (see (9) below). This is used to model the viscoelastic fluid flow motions. One can also view this model as a regularization of the equations of motion arising in Oldroyd fluids (see for instance [39,26], etc). The interested readers are referred to see [42] for extensive mathematical studies on the model. 1.1. Linear viscoelastic fluids. When viscoelastic materials (for example like, polymers, metals at very high temperatures, etc) are exposed to strain (or deformation), stresses are generated in the material. The generated stresses are gradually released and the material experiences a process along time called "relaxation"; a kind of recovery to the initial conditions. If the system is constantly strained, the stresses decrease with time. The retardation function describes, how the strain is retarded and develops in time. Both relaxation and retardation functions characterize the behavior of linear viscoelastic bodies.
A linear viscoelastic fluid with a finite discretely distributed relaxation times {λ l } and retardation times {κ m }, we mean a fluid whose defining (or rheological) equation, connecting the deviator of the stress tensor σ and the strain tensor D, has the form (cf. [33,34,37,5,24]): where the numbers L and M are connected by one of the following three conditions: In the sequel, we shall assume that the relaxation times {λ l } satisfy the following conditions: • the roots {α l } of the polynomial Q(p) = 1 + where for the L , κ 0 = 1, λ −1 = 0 and µ 0 = (κ L − µ 1 λ L−1 )λ −1 L . Without loss of generality, we assume that σ(x, t) ≡ 0, for t < 0 and we apply the Laplace transform L with respect to the variable t in (1) to obtain where From (3) and (4) • for Kelvin-Voigt fluids of order L = 1, 2, . . .
Substituting (5)- (7) into the equations of the motion of a continuous incompressible medium in Cauchy form, that is, in the following equation: ∂u ∂t + (u · ∇)u + ∇p = ∇ · σ + f (x, t), ∇ · u = 0, we obtain the integro-differential equations of the motion of Maxwell fluids of order L, Oldroyd fluids of order L, and Kelvin-Voigt fluids of order L = 1, 2, . . . , respectively as: In this work, we consider the Kelvin-Voigt fluids of order one with "fading memory" as where

MANIL T. MOHAN
The Kelvin-Voigt fluid flow motion of order one with "memory of length τ " can be written as (see [36]) for 0 < τ ≤ t. Let us introduce w(t) := t t−τ e −δ(t−s) u(s)ds, so that the system (10) can be rewritten as The system (9) or (10) with γ = 0 is known as Navier-Stokes-Voigt equations or Oskolkov equations. The global solvability results of the system (9) in bounded domains are available in the literature, cf. [38,39,36,42], etc. The method we followed in this work to establish the existence and uniqueness of weak solutions to the system (9) can also be applied to unbounded domains like Poincaré domains, as we are not using any compactness arguments in the proofs. Since, the global solvability results are true for all γ ≥ 0, the results obtained in this paper are true for Navier-Stokes-Voigt equations also (cf. [32,15,2], etc for solvability results regarding Navier-Stokes-Voigt equations and references therein). We study the asymptotic behavior of the steady state solutions and also establish the exponential stability results. Under suitable assumptions on the Galerkin basis, we consider the controlled Galerkin approximated 3D Kelvin-Voigt fluid flow equations and establish the exact controllability results, using the Hilbert uniqueness method combined with Schauder's fixed point theorem (see [18,19,20,1,22], etc for different models). Moreover, one can extend these results to the Kelvin-Voigt fluid flow equations with "memory of length τ ", under some additional assumptions.

1.2.
Organization of the paper. The rest of the paper is organized as follows. In the next section, we provide an abstract formulation of the problem (9) and describe the necessary function spaces needed to obtain the global solvability results. We also explain some important properties of the memory kernel appearing in the system (9) (Lemmas 2.6 and 2.9). In section 3, we establish the global existence and uniqueness of a weak solution to the system (9), by first considering a Galerkin approximated system (Theorem 3.5). The local monotonicity property of the linear and nonlinear operators (Theorems 3.2 and 3.3), and a localized version of the Minty-Browder technique are exploited in the proofs. As we are not using any compactness arguments in the proofs, the global solvability results remain valid in unbounded domains (like Poincaré domains) also (Remark 7). Under suitable regularity on the initial data and external forcing, we also obtain further regularity of the weak solutions (Remark 8). An important remark has been given in this section that using an m-accretive quantization of the linear and nonlinear operators, the global solvability results for the Navier-Stokes-Voigt equations can be established, avoiding the tedious Galerkin approximation scheme (Remark 9). We consider the stationary solutions in section 4 and establish the exponential stability results of the steady state solution (Theorem 4.3). In the final section, under suitable assumptions on the Galerkin basis, we consider the controlled Galerkin approximated 3D Kelvin-Voigt fluid flow equations. Using the Hilbert uniqueness method combined with Schauder's fixed point theorem, we establish the exact controllability of the finite dimensional Galerkin approximated system (Theorem 5.2).
2. Mathematical formulation. In this section, we describe the necessary function spaces needed to obtain the global solvability results of the system (9) in three dimensional bounded domains. We consider the following three dimensional Kelvin-Voigt fluid flow equations of order one: with the conditions where Ω is a bounded open domain in R 3 with smooth boundary ∂Ω, u(t, x) ∈ R 3 denotes the velocity field at time t and position x, p(·, ·) denotes the pressure field, f (·, ·) is an external forcing, and The final condition in (12) is introduced for uniqueness of the pressure p.
2.1. Functional setting. Let us denote is the space of all infinitely differentiable functions with compact support in Ω. Let H and V be the completion of V in L 2 (Ω; R 3 ) and H 1 (Ω; R 3 ) norms respectively. Then under some smoothness assumptions on the boundary, we characterize H and V as where n is the outward normal to ∂Ω and The inner product in the Hilbert space H is denoted by (·, ·) and the induced duality, for instance between the spaces V and its dual V by ·, · . Since V is densely and continuously embedded into H and H can be identified with its dual H , we have the following Gelfand triple: If Ω is a bounded domain, then the embedding of V ⊂ H is compact. In the sequel, we also use the notation H 2 (Ω) := H 2 (Ω; R 3 ) for the second order Sobolev spaces.
2.2. Linear operator. Let P H : L 2 (Ω) → H be the Helmholtz-Hodge orthogonal projection. Then, we define The operator A is a non-negative self-adjoint operator in H with V = D(A 1/2 ) and Note that for a bounded domain Ω, the operator A is invertible and its inverse A −1 is bounded, self-adjoint and compact in H. Hence the spectrum of A consists of an infinite sequence 0 < λ 1 ≤ λ 2 ≤ . . . ≤ λ k ≤ . . . , with λ k → ∞ as k → ∞ of eigenvalues, and there exists an orthogonal basis {e k } ∞ k=1 of H consisting of eigenvectors of A: Ae k = λ k e k , for all k ∈ N. We know that u can be expressed as u, e j e j and Au = ∞ j=1 λ j u, e j e j . Thus, it is immediate that If u, v are such that the linear map b(u, v, ·) is continuous on V, the corresponding element of V is denoted by B(u, v). We also denote (with an abuse of notation) B(u) = B(u, u) = P H (u · ∇)u. An integration by parts yields Let us next give a particular case of the Gagliardo-Nirenberg inequality, which is used frequently in the sequel. Lemma 2.1 (Gagliardo-Nirenberg inequality, Theorem 2.1, [31]). Let Ω ⊂ R n and u ∈ W 1,p 0 (Ω; R n ), p ≥ 1. Then for any fixed number q, r ≥ 1, there exists a constant C > 0 depending only on n, p, q such that where the numbers p, q, r and θ satisfy the relation Let us take r = n = 3 and p = q = 2 in (15) to get θ = 1 2 and u L 3 ≤ C ∇u 1/2 Now if we consider r = 4, n = 3 and p = 2q = 2 in (15), then we find θ = 3 4 and u L 4 ≤ C ∇u 3/4 where the constant C = √ 2 (see Lemma 2, Chapter 1 [17]). If we take r = 6 n = 3 and p = q = 2 in (15), then we obtain θ = 1 and where the constant C = 48 1 6 (see Lemma 2, Chapter 1, [17]). Next, we state an important inequality due to Ladyzhenskaya (see Lemma 1 and 2, Chapter 1, [17]), which is valid in unbounded domains also: where C = 2, 4, for n = 2, 3, respectively.
We also need the following general version of the Gagliardo-Nirenberg interpolation inequality for higher order estimates. Let Ω ⊂ R n , u ∈ W m,p (Ω; R n ), p ≥ 1 and fix 1 ≤ p, q ≤ ∞ and a natural number m. Suppose also that a real number θ and a natural number j are such that and j m ≤ θ ≤ 1. Then for any u ∈ W m,p (Ω; R n ), we have where s > 0 is arbitrary and the constant C depends upon the domain Ω, m, n.
If 1 < p < ∞ and m − j − n p is a non-negative integer, then it is necessary to assume also that θ = 1. It should be noted that for u ∈ W 1,p 0 (Ω; R n ), Lemma 2.1 is a special case of the above inequality, since for j = 0, m = 1 and 1 s = θ p + 1−θ q in (21), and an application of the Poincaré inequality yields (15). Note also that (21) can also be written as for any u ∈ W m,p (Ω; R n ). If we take j = 1, n = 3, p = q = 2 and r = 3 in (22), then we obtain θ = 7 8 and ∇u L 3 ≤ C Au for all u ∈ H 2 (Ω) ∩ H 1 0 (Ω). In the trilinear form, an application of Hölder's inequality yields The trilinear map b : V × V × V → R has a unique extension to a bounded trilinear map from L 4 (Ω) × (L 4 (Ω) ∩ H) × V to R. It can also be seen that B maps L 4 (Ω) ∩ H (and so V) into V and for all v ∈ V, so that using the Poincaré inequality. Note that B also maps L 6 (Ω) ∩ H into V and so that once again, we find (24). Making use of (24), we also have so that B(·) is a locally Lipschitz operator from V to V .
Remark 1. One can show that the norms u V and (I+µ 1 A) 1/2 u H are equivalent. It can be easily seen that using the Poincaré inequality. From (26), it is also clear that Combining (26) and (27), we obtain the required result.

2.4.
Properties of the kernel. We now consider a general kernel β(·) and discuss some of its properties. Later, we specialize it to the kernel β(t) = γe −δt . Let us define  Then, we have the following result on the positivity of the kernel β(·): Then, β(t) defines a positive kernel.
Lemma 2.6 (Lemma 2.8, [27]). Let β ∈ L 1 (0, T ), f, g ∈ L 2 (0, T ), for some T > 0. Then, we have Proof. With a change of variable and change of integrals, it can be easily seen that which completes the proof.
For completeness, we give a proof for an elementary version of the above lemma. Let us first provide a result on the positivity of a class of Riemann-Stieltjes integrals (see [21] also). Remember that if a function g is of bounded variation, then there exist nondecreasing functions g + = max{g, 0} and g − = − min{g, 0} such that g = g + − g − . Therefore, finite lim x→y + g(x) and lim x→y − g(x) exist for any y ∈ [a, b].
is positive.
Proof. Since f is continuous and f : Since g(a) = 0, we may assume that g + and g − are nonnegative and satisfy g ± (a) = 0. If g + (y) > 0 and g − (y) = 0, then the integral (37) is at least which completes the proof.
Then, for any continuous function f : [0, T ] → [0, ∞) and for a.e. continuous function g : Proof. We know that a.e. t ∈ [0, T ] and the derivative is continuous a.e. Hence the function ((a * g)(t)) 2 is differentiable a.e. t ∈ [0, T ] and the derivative is integrable on [0, T ], and hence is an absolutely continuous function. Thus, it is of bounded variation also. Using Clearly the right hand side of the above equality is greater than or equal to 0 using the positivity of f , absolute continuity of ((a * g)(·)) 2 (of bounded variation) and Lemma 2.8 (note that (a * g)(0) = 0).

Remark 4.
Changing the order of integration, we find Abstract formulation. Our next aim is to provide an abstract formulation of the system (11)- (12). Let us take the Helmholtz-Hodge orthogonal projection P H in the equation (11) to write down an abstract form of the system (11)-(12) as: where f ∈ L 2 (0, T ; V ) (we used P H f = f for simplicity). The system (42) is equivalent to the following system: for t ∈ (0, T ). Let us now give the definition of weak solution for the system (42).
and the energy equality Let us now state one of the main results of this work, the existence and uniqueness of weak solutions to the system (42).
Theorem 2.11. Let u 0 ∈ V and f ∈ L 2 (0, T ; V ) be given. Then, there exists a unique weak solution u(·) to the system (42) in the sense of Definition 2.10.
Even though a proof of Theorem 2.11 using compactness methods is available in the literature (cf. [38,39,36,42], etc) for bounded domains, we prove Theorem 2.11 in the following section using a localized version of the Minty-Browder technique, which avoids compactness arguments.
3. Global existence and uniqueness of weak solution. In this section, we prove the global existence and uniqueness of weak solutions to the system (42) using local monotonicity property of the linear and nonlinear operators, and a localized version of the Minty-Browder technique. Several authors have used this method to prove global solvability results of various models (cf. [13,30,28,29,26], etc.) 3.1. Local monotonicity. In this subsection, we establish a local monotonicity property of the operator F(u) = µ 0 Au + B(u). Before proceeding further, let us provide the definition of monotone, hemicontinuous and demicontinuous operators. [3]). Let X be a Banach space and let X be its topological dual. An operator F : Clearly, demicontinuity implies hemicontinuity.
We prove two kinds of local monotonicity results, that is, we prove that the operator is locally monotone in a V-ball as well as in an L 4 -ball.
Theorem 3.2. For the operator F(u) = µ 0 Au + B(u), we have for all v ∈ B N , where B N is a ball of radius N , i.e., Proof. Using an integration by parts, it is immediate that We apply Hölder's, Gagliardo-Nirenberg and Young's inequalities to estimate Combining (45) and (46), we find Since v ∈ B N , we finally obtain (44).

Remark 5.
Using (17), one can also estimate B that is we can take the constant C appearing in (44) as 16.
Proof. We use an integration by parts, Hölder's, Gagliardo-Nirenberg and Young's
is demicontinuous and hence is hemicontinuous.
Proof. Let us take u n → u in L 4 (0, T ; V) as n → ∞. For v ∈ V, we consider where we used an integration by parts and Hölder's inequality. Since u n , u ∈ L 4 (0, T ; V) and u n → u in L 4 (0, T ; V), using Ladyzhenskaya and Poincaré inequalities, we get u n → u in L 4 (0, T ; L 4 (Ω)). It can be easily seen that where we used Hölder's inequality and (35). Similarly, we have . Moreover using Hölder's inequality, we also have . Hence the right hand side of the inequality (53) tends to zero as n → ∞, for all and hence the operator is demicontinuous. Since demicontinuity implies hemicontinuity, we also have the operator F(·) is hemicontinuous.
A calculation similar to above also shows that is demicontinuous and hence is hemicontinuous.
3.2. Energy estimates. Let {e 1 , . . . , e n , . . .} be a complete orthonormal system in H belonging to V and let H n be the n-dimensional subspace of H. Let P n denote the orthogonal projection of V to H n , that is, x, e i e i . Since every element x ∈ H induces a functional x * ∈ H by the formula x * , y = (x, y), y ∈ V, then P n H , the orthogonal projection of H onto H n is given by Hence in particular, P n is the orthogonal projection from H onto span{e 1 , . . . , e n }. We define B n (u n ) = P n B(u n ) and f n = P n f .
With the above setting, let us now consider the following system of ODEs: with u n 0 = P n u 0 , for all v ∈ H n . Equivalently, the system (54) can be written as Since B n (·) is locally Lipschitz (see (25)), the system (55) (or equivalently (54)) has a unique solution u n ∈ C([0, T ]; H n ), (existence due to Carathéodory's existence theorem and uniqueness is immediate from the local Lipschitz property (Picard's)).
Let us now establish the energy estimates satisfied by the system (54).
Proof. Let us take inner product with u n (·) to the system (54) to obtain 1 2 where we performed an integration by parts and used the fact that (B(u n ), u n ) = 0. Let us integrate the inequality from 0 to t and use (31) to find Using Cauchy-Schwarz and Young's inequalities, we estimate |(f n , u n )| as Thus, from (57), we infer that where we also used the fact that u n 0 V = P n u 0 V ≤ u 0 V and f n V ≤ f V . An application of Gronwall's inequality in (58) yields for all t ∈ [0, T ] and (56) follows. Note that the sequence u n (·) remains bounded in L ∞ (0, T ; V) and the estimate is independent of µ 0 also. That is, the energy estimate is valid for inviscid limit (i.e., µ 0 = 0) also.

Remark 6.
Since the estimate (59) is independent of µ 0 , the results obtained in this work can be extended to the models like three dimensional inviscid simplified Bardina turbulence models also (cf. [23,9], etc for more details).

Global existence and uniqueness.
Let us now show that the system (42) has a unique weak solution using Galerkin approximations and energy estimates obtained in the previous subsection. We make of use the local mononicity property of the linear and nonlinear operators, and a localized version of the Minty-Browder technique to obtain global solvability results. Similar methods are applied for establishing global solvability results for the 2D Navier-Stokes equations in Theorem 4.2, [13] and Oldroyd fluids of order one in Theorem 3.7, [26].
Proof. Part (1). Existence: We establish the existence of a weak solution to the system (42) by using a localized version of the Minty-Browder technique in the following steps.
Step (i): Finite-dimensional (Galerkin) approximation of the system (42): Let {e 1 , e 2 , . . .} be a complete orthonormal basis in H belonging to V (see section 3.2). Let us denote by H n , the n-dimensional subspace of H. We consider the following ODE in H n : for any v ∈ H n , where u n 0 is the orthogonal projection of u 0 in the span{e 1 , . . . , e n }. We define the operator F(u) := µ 0 Au+β * Au+B(u). Then the system (61) becomes for any v ∈ H n . By taking the test function as u n (·) in (62), we get the following energy equality: for any t ∈ (0, T ). One can also show that (see Theorem 3.2, [28]) and hence the following energy equality holds: e −2r(s) (F(u n (s)) − f n (s) +ṙ(s)u n (s), u n (s))ds for any t ∈ (0, T ]. The quantity r(t) given in (65) will be specified later.
Step (ii): Weak convergence of the sequences u n (·) and F(u n (·)): Remember that the dual of L 1 (0, T ; V ) is L ∞ (0, T ; V) (that is, L 1 (0, T ; V ) ∼ = L ∞ (0, T ; V)) and L 1 (0, T ; V ) is separable. Using the Banach-Alaoglu theorem (or Helly's theorem, note that we have energy estimates in L ∞ (0, T ; V), see Proposition 1), we can extract subsequences {u n k (·)} and {F(u n k (·))} such that (for notational convenience, we use {u n (·)} and {F(u n (·))}) The final convergence in (66) can be justified using (32) and (24) as and hence (I + µ 1 A)u n w − → (I + µ 1 A)u, in L 2 (0, T ; V ). It should also be noted that using (67). Thus, we have ∂ t u n w − → ∂ t u in L 2 (0, T ; V) and hence the energy estimates in Proposition 1 and (66) also guarantee that ∂ t (I + µ 1 A)u n w − → ∂ t (I + µ 1 A)u in L 2 (0, T ; V ). Note that f n → f in L 2 (0, T ; V ). On passing to limit n → ∞ in (62), the limit u(·) satisfies: One can easily show that the system (68) satisfies the following energy equality: for any t ∈ (0, T ]. A calculation similar to (65) yields the following estimate: for any t ∈ (0, T ]. Since u n (0) = P n u(0), it is immediate that the initial value u n (0) converges to u(0) strongly in V and H, i.e., we have is controlled using Ladyzhenskaya inequality (see Lemma 2.2). Using the local monotonicity result (see Corollary 1 and (51)), we have Let us apply the energy equality (65) in (73) to find Next we take liminf on both sides of (74) to obtain Using the lower semicontinuity property of the V and H norms, and the strong convergence of the initial data u n (0) (see for example (71)), one can easily deduce that the second term on the right hand side of the inequality (75) satisfies the following inequality: Let us now apply (76) in (75) to find where in the second step, we used the energy equality (70). It should be noted that the estimate (77) holds true for any v ∈ L ∞ (0, T ; V m ), m ∈ N, since the inequality given in (77) is independent of both m and n. One can also use a density argument that the inequality (77) remains true for any v ∈ L ∞ (0, T ; V). In fact, for any v ∈ L ∞ (0, T ; V), there exists a strongly convergent subsequence v m ∈ L ∞ (0, T ; V) that satisfies the inequality in (77). Let us now take v = u + λw, λ > 0, where w ∈ L ∞ (0, T ; V), and substitute for v in (77) to obtain T 0 e −2r(t) F(u(t) + λw(t)) − F 0 (t) +ṙ(t)λw(t), λw(t) dt ≥ 0.
Thus, the inequality (78) changes to Let us now divide the inequality in (79) by λ, use the hemicontinuity property of the operator A + β * A + B(·) (see Lemma 3.4), and then pass λ → 0 to obtain Note that the term Thus, from (80), we infer that for any w ∈ L ∞ (0, T ; V). Hence from (82), we finally deduce that F(u(t)) = F 0 (t).
Remark 7. 1. Since there are no compactness arguments involved in the proof of Theorem 3.5, it is valid in both bounded and unbounded domains like Poincaré domains (see [26] for more details).
2. It should also be noted that the Theorem 3.5 holds true for Kelvin-Voigt fluids of order L = 1, 2, . . . , also.
Remark 8 (Regularity). If we take, u 0 ∈ D(A) and f ∈ L 2 (0, T ; H), then one can show that u ∈ C([0, T ]; D(A)). Let us take inner product with Au(·) to the first equation in (42) to obtain Using Hölder's, Gagliardo-Nirenberg, and Young's inequalities, we estimate |(B(u), Au)| as (see (18) and (23)) Using Cauchy-Schwarz and Young's inequalities, we also get Let us combine (87)  Remark 9. If we take γ = 0, that is, if we consider Navier-Stokes-Voigt equations, then one can get the strong solution to the system (42) using the same techniques followed in [6] for the Navier-Stokes equations. Here we use the m-accretive quantization of the linear and nonlinear operators to get such a result (see [3,4] for more details on accretive operators). This method also avoids the tedious Galerkin approximation scheme, for example used in Theorem 3.5, Theorem 3.2, [2], etc. Since, the norms u V and (I + µ 1 A) 1/2 u H are equivalent, we can also consider and another the inner product on V as ((u, v)) V = ((I + µ 1 A) 1/2 u, (I + µ 1 A) 1/2 v), for all u, v ∈ V. We first define the quantized nonlinearity B N (·) : V → V as Let us now define the following quantized nonlinearity B N (·) : V → V also: Since v ∈ V is arbitrary, we get B N (u) V ≤ 2

then a calculation similar to (47) yields
Let us now consider the case, u V , u V > N . Without loss of generality, we may assume that v V ≥ u V . Then, using Hölder's, Ladyzhenskaya, Poincaré and Young's inequalities, we find A similar estimate can be obtained for the case u V ≤ N and v V > N . Thus, we have , and λ ≥ C N , and hence Γ N + λI is a monotone operator. An argument similar to Lemma 5.1, [6] yields the following result: There exists an α N > 0 such that Γ N + α N I is maximal monotone on V (or Γ N is quasi-m-accretive on V). Then using (91) For f ∈ W 1,1 (0, T ; D(A)) and u 0 ∈ D(A), the system (91) can also be written as: where B N is defined in (90), so that u ∈ C([0, T * ]; D(A)). Thus, it turns out that for large enough N , u N is the solution of the system (92) on a certain interval (0, T * ). Let us take inner product with u N (·) to the first equation in (92) to obtain But, we know that 4. Exponential stability. In this section, we examine some asymptotic behavior of the problem (42) and establish stability results for the steady state solution corresponding to the system (42). Similar results for the 2D Navier-Stokes equations and 2D Oldroyd models in H are established in [41,26], respectively, and for the 2D tidal dynamics system is proved in [25]. For the exponential stability results regarding the abstract Kelvin-Voigt equations or Oskolkov models, the interested readers are referred to see [32]. Here we assume that the external forcing f is independent of t and consider the following stationary system:
In particular, if u ∞ is a stationary solution of system (42), then u ∞ is called exponentially stable in V, provided that any weak solution to (42) converges to u ∞ at the same exponential rate α > 0. Theorem 4.3. Let f ∈ V and u ∞ be the unique weak solution of the system (94). If u(·) is any weak solution of the system (42) with u 0 ∈ V arbitrary, then u ∞ ∈ V is exponentially stable in V, and for µ 0 > 4 Proof. Since β(t) = γe −δt , δ > 0 and (96) holds true, we can chose a constant η > 0 such that Since u ∞ is the steady state solution, we know that Let us multiply the first equation in (98) by e ηt , take inner product with e ηt (u(t)− u ∞ ) and then integrate from 0 to t to obtain Remember that u ∞ is a stationary solution to (94), so that we have A calculation similar to (46) yields Using the definition of β(t) = γe −δt , we also have We denote the final two terms from the right hand side of the equality (103) as I 1 and I 2 . Thanks to Lemma 2.9, we have We apply Cauchy-Schwarz, Hölder's and Young's inequalities to estimate I 1 as which holds true for all δ > η. Using divergence free condition, Hölder's, Ladyzhenskaya and Poincaré inequalities, Using (106), we estimate |I 2 | as Let us substitute (104)-(107) in (103) to obtain If (96) holds true, then we have so that (95) follows and the stationary solution u ∞ is exponentially stable.
Next, we assume that f (·) depends on u(·) in such a way that f (·) is globally Lipschitz. That is, there exists a constant L > 0 such that and f (0) = 0. Then, one can establish the following theorems:

MANIL T. MOHAN
Theorem 4.4. If u(·) is any weak solution of the system (42) with u 0 ∈ V arbitrary, then we have Proof. Since µ 0 − L > 0, we can find an α such that Now we take −r(t) = αt in (70) with F 0 (·) = F(u(·)) to obtain t 0 e 2αs (β * ∇u(s), ∇u(s))ds Note that the final integral in the left hand side of the equality (112) is positive, using Lemma 2.9. Let us apply the Poincaré inequality, Lemma 2.9, Cauchy-Schwarz inequality and Young's inequality in (112) to find Hence under the condition given in (111), the solution u(t) of the system (42) converges to zero in V exponentially, since and hence zero solution is exponentially stable.
The above results motivated us to consider the problem of the existence of global and exponential attractors for the system (42) and it has been examined in the work [27].
Theorem 4.5. Let u and v be any two weak solutions of the system (42) with initial data u 0 and v 0 in V, respectively. Then, we have Proof. Let us define w := u − v and then w satisfies: Since µ 0 − 2L > 0, we can find an ζ such that t 0 e 2ζs (β * ∇w(s), ∇w(s))ds Once again a calculation similar to (106) gives Using the global Lipschitz property of f (·), we get Let us use the above estimates in (118) to find where we used Lemma 2.7. An application of Gronwall's inequality in (119) gives for ζ given in (117). We also have the following energy estimate: Thus, for µ 0 > 2L > L, we infer that From (120), we finally have and hence (115) is holds true.

5.
Exact controllability of Galerkin approximations. In this section, our main goal is to establish the exact controllability of Galerkin approximation of the controlled three dimensional Kelvin-Voigt fluids of order one. Let us first consider the following controlled system: where O is the control domain (the domain in which control is acting on, and is supposed to be as small as necessary), U is the distributed control function acting over the system and χ O denotes the characteristic function of O.
Let us take inner product with v ∈ V to the first equation in (122) to obtain the following variational formulation of (122): Let us take {e j } ∞ j=1 as basis of V such that The existence of this basis is guaranteed from the abstract results due to Lions  Since the system under our consideration is nonlinear, we employ the similar techniques adopted in [18,19,20,1,22] to obtain the exact controllability of Galerkin approximated system (125). The authors in [18] obtained the exact controllability of Galerkin approximated Navier-Stokes equations with distributed control and [19] obtained the same with boundary control. Exact controllability of Galerkin's approximations of micropolar fluids is established in [1] and of Oldroyd fluids is obtained in [22].
Next, we consider the following finite dimensional space: We consider Galerkin approximation of the system (42) with f = 0. Let us introduce the following Galerkin approximation of the variational formulation (see (123)): ((I − µ 1 ∆)u t , e) + µ 0 (∇u, ∇e) + ((β * ∇u), ∇e) + ((u · ∇)u, e) = (Uχ O , e), for all e ∈ E. We have already shown that the system (125) has a unique solution with u in C([0, T ]; E) (see section 3.2). We abused the notations in (125), for example u 0 in (125) needs to be interpreted as u 0 = n i=1 (u 0 , e i )e i . In the sequel, we use the notations · and (·, ·) for the norm and inner product, respectively on L 2 (Ω) and finite dimensional space E.

Exact controllability.
Let us first provide the definition of exact controllability of the system (125).
Definition 5.1. We say that Galerkin approximated system (125) is said to be exactly controllable at time T > 0, if for given u 0 , u T ∈ E, there exists a control U ∈ L 2 (0, T ; L 2 (O)) such that the solution u of (125) satisfies: In order to achieve this goal, we first consider the following cost functional: Let us now state and prove the main result of this section. Proof. In the second part of the Theorem, we need to show that the cost functional is bounded independent of the nonlinearity, so we consider the following system: with α ∈ R. In step 2 below, we prove some estimates for the system (128) independent of α. We establish the exact controllability results for the system (128), which easily gives the required result for the system (122) (by taking α = 1). Let us first introduce the following variational formulation of the system (128): for all e ∈ E. We provide rest of the proof of Theorem 5.2 in the following steps.
Step 1. Linear system. Let us take a function h such that h ∈ L 2 (0, T ; E), and consider the following linear system: Since the system (130) is linear, it is immediate that the system has a unique solution with u ∈ C([0, T ]; E). Note that the linearity of the system (130) also allow us to work with the null initial data. The result is still valid if one works with non-null initial data also, i.e., if we take u(0) = u 0 ∈ E. Let us now show that the system (130) is exactly controllable in any time T > 0, in the sense of Definition 5.1. In order to establish this, it suffices to prove that if g ∈ E satisfies (u(·, T ; U), g) = 0, for all U ∈ L 2 (0, T ; L 2 (O)), then g = 0.
(131) To achieve this goal, we first consider the following adjoint system: where g ∈ E and (β • ∆p)(t) = T t β(s − t)∆p(s)ds. Moreover, the variational formulation of the above system is given by for all e ∈ E. Due to the linearity of (132), it is immediate that the system has a unique solution with p ∈ C([0, T ]; E). Let us now take e = u in (133) to find where we used (41) and performed an integration by parts. Using (129) and the fact that p(T ) = (I − µ 1 ∆) −1 g in (135), we deduce that Now, if the assumption given in (131) holds true, then from (136), it can be easily seen that T 0 (U(t)χ O , p(t))dt = 0, for all U ∈ L 2 (0, T ; L 2 (O)).
Step 2. Estimates using duality argument. Using the results obtained in Step 1, we define the functional M : where U ad is the set of all admissible controls U ad = U ∈ L 2 (0, T ; L 2 (O)) : u is a solution of (129) satisfying (126) .
Our main aim is to establish that where K is a constant independent of both h and α. We establish the bound in (139) using a duality argument. Let us consider the continuous linear operator L : L 2 (0, T ; L 2 (O)) → E as L(U) = u(·, T ; U).
We also introduce the following functionals: With the above definitions, one can rewrite functional M given in (138) as follows: Making use of the Fenchel-Rockafellar duality theorem (see Theorem 31.1, [40]), we have where L * : E → L 2 (0, T ; L 2 (O)) is the adjoint of L. Using (136), it can be easily shown that We know that |p(x, t)| 2 dxdt and F * 2 (−p) = −(g, u T ), and then we have Since E is finite dimensional, all the norms on E are equivalent. Also, making use of the hypothesis on the basis of E, we have O |e(x)| 2 dx is a norm on E, so that for some positive constants c and C depending only on E, we obtain Using this in (140), we find Let us now take e = p(t) in (133) and integrate from t to T to obtain Next, we integrate from 0 to T to find (change the order of integration also) Since E and F are finite dimensional, we have for some constants c, C > 0, that depends only on E and F . Using a change of order of integration, Cauchy-Schwarz and Hölder's inequalities, and (32), we obtain Thus, from (143), we deduce that for some suitable constants c, C, depending only on E and F . From (141), we also infer that −M(h) ≥ inf g∈E cT 4 1 2 (1 + cµ 1 ) + CT µ 0 + γ δ g 2 − (g, u T ) .
Step 3. Nonlinear system. Our next aim is to show that the nonlinear system (128) is exactly controllable. Let h ∈ L 2 (0, T ; E) be given. For U ∈ L 2 (0, T ; L 2 (O)), we choose the unique element such that which defines a continuous mapping h → U from L 2 (0, T ; E) to L 2 (0, T ; L 2 (O)).
Let us now denote u(h) as the solution of (130) with the control U = U(h). Next, we take e = u(t) in (130) to obtain 1 2 d dt u(t) 2 + µ 1 ∇u(t) 2 + µ 0 ∇u(t) 2 = −((β * ∇u(t)), ∇u(t)) + (U(t)χ O , u(t)), since (B(u, u), u) = 0. Integrating from 0 to t, we find Our next goal is to establish that the map h → u(h) admits a fixed point in K. Remember that Schauder's fixed point theorem states that if K is a convex and closed subset of a Banach space X, then any continuous and compact map (i.e., bounded sets in K are mapped into relatively compact sets) f : K → K has a fixed point. In order to use Schauder's fixed point theorem, it is sufficient to establish that the range of u(h), when h spans over K is relatively compact in K. This easily follows from the following estimate: for all e ∈ E. Thus, it is immediate that Using Hölder's inequality, we estimate β * ∇u(t) as β * ∇u(t) ≤ .
Since E is finite dimensional, from (149), we also have which easily gives (147). Invoking Schauder's fixed point theorem, we infer that the map h → u(h) admits a fixed point in K. Thus, if h is a fixed point, then since the system (130) is exactly controllable in time T > 0, we easily obtain the exact controllability of the system (129). Note that the the system (129) is exactly controllable for any α ∈ R, and in particular for α = 1, which easily gives the exact controllability of the system (125). Furthermore, for any h, we have the uniform estimate (139), so that the cost functional given in (127) is bounded independent of the nonlinearity.
Remark 10. The results obtained in this paper can easily be extended to the Kelvin-Voigt fluids of order L = 1, 2, . . . , with "fading memory" (see the final equation in (8)) as well as "memory of length τ " also. The Kelvin-Voigt fluid flow equations of order L = 1, 2, . . . , with "memory of length τ " are described by the rheological equation (see [36]): σ(x, t) = µ 1 ∂D(x, t) ∂t + µ 0 D(x, t) + Thus, the integro-differential equations of the motion of Kelvin-Voigt fluids with "memory of length τ " of order L = 1, 2, . . . , can be written as: where we used (31) and (30). This implies that the conditions (96) and (111) need to be changed to respectively, where in the second condition, we take µ 0 > γ δ + L. A similar change has to be made in (117) also.