Global weak solutions for a coupled chemotaxis non-Newtonian fluid

This paper focuses on the mathematical analysis of a self-suggested model arising from biology, consisting of dynamics of oxygen diffusion and consumption, chemotaxis process and viscous incompressible non-Newtonian fluid in a bounded domain \begin{document}$Ω \subset \mathbb{R}^d$\end{document} , with \begin{document}$d = 2, 3.$\end{document} The viscosity of the studied fluid is supposed to be non constant and depends on the shear-rate \begin{document}$|{\bf{D}}\boldsymbol{v}|$\end{document} as well as the cell density \begin{document}$m$\end{document} and the oxygen concentration \begin{document}$c$\end{document} . Nonlinearities are also considered in the diffusion terms for the convection-diffusion equations corresponding to \begin{document}$m$\end{document} and \begin{document}$c$\end{document} . Under the choice of suitable structures and convenient assumptions for the nonlinear fluxes, we prove global existence of weak solutions, in the case of a smooth bounded domain subject to Navier's slip conditions at the boundary and for large range of initial data.


1.
Introduction. Chemotaxis is a biological process which consists of random movement of biological individuals and living cells (e.g. bacteria) at macroscopic scale as response to gradients of a chemical substance. This behavior can either be towards the chemical stimulus or away from. It was pointed out in [14] that this phenomenon is widely present in many prototypical biological situations, such as auto-organization of cells during embryonic development [28], migration of capillary sprouts in response to a chemo-attractant field set up by a tumour-released angiogenic factor [8], and many other situations.
In the 1970s, Keller and Segel [19,20] investigated the chemotactic aggregation of Dictyostelium discoideum bacteria. Based on previous work of Patlak [30], they suggested a mathematical model of such a chemotaxis mechanism consisting of four coupled reaction-advection-diffusion equations. Under quasi-steady-state assumptions, the system was later simplified by Hortsmann [18] to ∂ t n − ∆n = −div (n ∇c), where n and c are scalar functions denoting the cell density and the concentration of the chemical signal, respectively.
System (1) consists of two parabolic reaction-diffusion equations which describe the evolution of n and c within a certain chemotaxis process.
The Keller-Segel equations (1) have received a lot of interest by mathematicians. Diverse theoretical studies, such as global existence, finite time blow-up, pattern formation have been conducted. A survey of important results related to (1) and some of its variants can be found in [1,18,17]. Later, experimental studies [11,36] investigating the situation of swimming bacteria "Bacillus Subtilis" when suspended into sessile drops of water revealed that an aggregation of cells like plumes may occurs, and a possible spontaneous emergence of large-scale fluid motion and convection patterns.
On the basis of these conclusions, Truval et al. [36] proposed a coupling of Keller-Segel system with the Navier-Stokes equations modelling the dynamics of populations of aerobic bacteria within the flow of an incompressible Newtonian fluid. The arising governing system of equations resumes in In addition to the unknowns described in (1), the vector field u stands for the fluid velocity evolution, with associated pressure p. The third equation of (2) has the ability to describe the flow governed by either the incompressible Navier-Stokes system ( if κ = 0), or the Stokes system for κ = 0. The scalar function Φ stands for the gravitational potential forcing generated by aggregation of cells onto the fluid. The chemotactic sensitivity χ(c) is not yet constant. This is due to the fact that, in some chemotaxis situations, migrations of cells are not necessarily in parallel directions to the gradient of the signal. On the other hand, the source term f (c) comes from the oxygen consumption rate of the chemical by the cells. Both χ(c) and f (c) are non negative functions.
Observe that in (2), both the population density n and the chemical concentration c diffuse through the fluid and also are transported by it. Note that the oxygen is consumed proportional to a cut-off function f and the density of cells n which in turn is moving in the direction of the chemical gradient.
The chemotaxis Navier-Stokes system governed by (2) has recently drawn the attention of many mathematical studies. The first works were leaded by Lorz in [23] where he constructed local weak solution using Schauder's fixed point theory for the case of a bounded domain with constant chemotactic sensitivity and f fulfilling f (0) = 0 and also some monotonicity condition. Unique global classical solution result for the two dimensional case under large initial data was obtained by Winkler [38]. He also proved later, in [39], that the obtained solution (n, c, u) converges (as time goes to infinity) to a constant state (n ∞ , 0, 0). Recently, Braukhoff [4] established existence of a unique global classical solution when d=2, as well as existence of a global weak solution for the 3d case, when considering an additionally logistic growth of the bacteria population.
Considering slow fluid motion by vanishing the convective term (κ = 0) returns out (2) to the chemotaxis-Stokes system. In that situation, Duan et al. [12] proved the existence of global weak solutions for the Cauchy problem in R 2 . In the case of the three dimensional space, similar results were obtained by Winkler [38] for a bounded convex domain and Peng and Xiang [31] in an unbounded domain with boundary, assuming large initial cell density and velocity. We also mention a very recent work of Winkler [40] who obtained, for a nonlinear chemotaxis-Stokes system, global bounded weak solutions that approach a spatially homogeneous steady state in the large time limit.
In [17], Hillen and Painter reviewed within previous works [20,29] (and others) the possibility and even more the convenience of incorporating nonlinear diffusion terms into the cell dynamics depending on the signal concentration and/or the cell density. Nonlinearities have also been introduced in the diffusion term of bacteria's evolution equation (2) 1 by replacing ∆n by ∆n m . We review some mathematical results dealing with.
In the case of the chemotaxis-Stokes flow, Di Francesco et al. [10] established global-in-time existence of bounded solutions in a bounded domain of R 2 when the range of the exponent m is ( 3 2 , 2]. Tao and Winkler [34] extended the former result for all m > 1. The Cauchy problem in R 2 was investigated in [22] where weak solution are obtained for m = 4 3 . For the three dimensional case, Tao and Winkler [35] generalized the last stated result for smooth bounded convex domain for all m > 8 7 . A result which has been extended in [13] by Duan and Xiang to cover the range [1, 8 7 ]. When considering the chemotaxis Navier-Stokes system (κ = 1), results are fewer than in the former case. Nevertheless, Duan and Xiang [13] obtained similar results as described above. Vorotnikov [37] showed under restrictive conditions on the potential function Φ, existence of weak solution for the initial-boundary value problem for m > 4 3 . It is worth noticing also the recent work of Zhang and Li [42] where they proved global weak solution when m ≥ 2 3 for the general form of (2) i. e. κ ∈ R and arbitrary. As it was indicated, system of type (2) was proposed for modeling populations of aerobic bacteria when suspended into sessile drops of water [36] or more generally in a Newtonian fluid governed by the Navier-Stokes model. It is one of the most widely used models when describing flows of incompressible Newtonian fluids.
Recall that a fluid is called Newtonian fluid if its Cauchy stress is linearly related to the symmetric part of the velocity gradient via some constant µ called the viscosity of the fluid through the constitutive equation Here, S stands for the constitutively determined part of the Cauchy stress T, I is the identity matrix and p is the pressure of the fluid whereas the matrix Dv is the symmetric part of the velocity gradient defined through Dv = 1 2 (∇v + (∇v) t ). But in nature, when modeling some flows in many areas of engineering sciences such as polymer mechanics, glaciology, blood and food rheology... whose behaviors can not be captured by the Newtonian's law of viscosity (3), i. e. the viscosity µ is not still constant and may depend on the shear-rate |Dv| (i. e. the Euclidean norm of Dv) and even more on some parameters affecting the fluid motion like pressure, density, concentration, temperature, etc... These fluids are then referred to as non-Newtonian fluids. For example fluids with shear rate dependent viscosity are capable of capturing shear-thinning or shear-thickening behaviors are modeled by where µ(|Dv|) ∼ µ 0 (1 + |Dv| 2 ) is the generalized viscosity of the power-law type, dependent of the shear-rate. The interested reader is invited to review the works of Málek, Nečas, Ruzička, Bulíček, Rajagopal and many others.
During the three last decades, mathematical studies of non-Newtonian fluids have became very popular and challenging. Mainly, medical, biological and engineering fields manifested a great interest to the theoretical and applied mathematical analysis on biological fluids (such as blood, polymeric solutions), and fluids of food industry and agriculture, etc. Here, we will be content of reviewing only some of the most famous theoretical results dealing with generalized Navier-Stokes equations where the viscosity is depending on the shear-rate. The credit of the pioneering result goes back to Ladyzhenskaya [21] where the author established existence of weak solution for r ≥ 3d d+2 (r is the power-law index) using the theory of monotone operators and compactness arguments. Since the end of 90's, exhaustive studies on existence, regularity of weak and strong solutions were conducted by Bulíček, Málek, Rajagopal and coworkers ... for both steady and evolutionary case. Sometimes other parameters like pressure, temperature, density, and concentration were incorporated to the dependence of the viscosity function. A variety of papers dealing with these subjects is listed here [6,7,15,24,25,26].
In the second half of this first section we state the governing equations of our studied problem.
1.1. Formulation of the problem. In this work, we shall consider the dynamics of bacteria or micro-organisms living in an incompressible flow of a non-Newtonian fluid. We assume that the flow in question takes place in a bounded closed environment which does not let the opportunity to the bacteria to come or leave the system. Therefore we can suppose that the total mass of bacteria in the system does not change in time. So the container will be identified with a bounded open connected set Ω ⊂ R d (d=2,3) with sufficiently smooth boundary ∂Ω, ∂Ω ∈ C 1,1 .
Let T > 0 be an arbitrary fixed real number. The system of equations consists of a generalized homogeneous incompressible Navier-Stokes equations coupled with two convection-diffusion equations for the oxygen concentration and the cell population, taking place in Ω and a long the time interval [0, T ].
The considered flow is then governed by the following system of partial differential equations where Q := [0, T ] × Ω is the space-time cylinder and which will be associated to the lateral surface Σ := [0, T ] × ∂Ω. The unknowns of the problem are : is the velocity field representing the fluid motion, p : Q −→ R is the associated pressure of the fluid, m : Q −→ R + is the cell density, c : Q −→ R + is the oxygen concentration.
Unlike mostly of the previous works dealing with chemotaxis Navier-Stokes equations where the cell density is notified by the small letter " n ", we have chosen to replace it by the small letter " m " to avoid possible misunderstandings and confusion with the unit outward normal vector n, and the inconvenience with the superscript N of the constructed sequences, which will be widely used throughout the paper.
In (5) 3 , the extra stress tensor is given by S(m, c, Dv) := µ(m, c, |Dv| 2 )Dv where µ is the viscosity function of the fluid and manifests dependence on the cell density, the concentration and the shear-rate as it will be discussed later. We prefer to use S instead of S(m, c, Dv) as an abbreviate notation of the stress tensor. The cell density and concentration fluxes are given by ϑ(m)∇m and k(c, |Dv|)∇c, respectively. The former flux vector is only depending on the cell population whereas the later depends on both c and |Dv|.
Our problem is supplemented with the following boundary conditions. We shall suppose that the velocity field satisfies a Navier's slip boundary conditions v.n = 0 and α v tan + (S n) tan = 0 on Σ, where n is the unit outward normal vector to ∂Ω, α ≥ 0 is the domain-wall friction coefficient and called also slip parameter, the subscript tan refers to the tangential component at the boundary given by (u) tan = u − (u . n)n for any vector field u.
Regarding the fluid-solid boundary conditions (6), we have assumed in the first half that the boundary is impermeable and in the second part that the fluid has the ability of slipping at the solid boundary.
Zero-flux on the boundary will be assumed for the cell density and the oxygen concentration : k(c, |Dv|)∇c . n = 0 on Σ.
We shall close the PDE's system (5) by specifying the initial data as follows In addition we shall assume that the initial cell density and oxygen concentration are uniformly positive, i.e. there exist constants δ m , δ c ∈ (0, ∞) such that Further information on the given m 0 , c 0 and v 0 will be stated later.
The considered system (5) attempts to describe the behavior of this complex mixture of biological dynamics and the flow of a non-Newtonian fluid. We are unaware of any experimental or theoretical biological fluid dynamics studies which play the role of support and consolidation to the consideration we have referred to. Also, our problem as it was assumed by incorporating dependance in a nonlinear way of the viscosity function, the concentration and the cell population fluxes on the cell density, the oxygen concentration and the shear rate (even if with less dependance factors) seems to be never investigated before. Thereby, our obtained result looks to be the first of its kind.
Generalized NSE with Navier's slip boundary condition have already been the subject of multiple studies of Bulíček, Málek et al. [6,7], where the considered viscosity is mainly shear-rate depending, and even sometimes concentration or pressure or temperature depending also. The common feature of these works is the achievement of existence of weak solutions.
The main result of this paper is to prove global existence of weak solutions for the model (5) under convenient assumptions on the nonlinear diffusive terms.
Let us mention that when multiplying formally the dissipative term of the balance of linear momentum equation with an admissible test vector having zero normal component we get with the help of Green's formula The integrals evolved at the boundary could be simplified by dividing the vector field S n onto its normal and tangential components and with the help of (6) as follows Similarly, for the nonlinear diffusion functions ϑ and k, and according to the boundary conditions (7) and (8) we see that for an admissible test function ψ as well as As far as the structure of this paper goes, next we are about to introduce tools which will be of use, such as spaces and objects notations, technical lemmas and the necessary assumptions on the viscous stress tensor S and the diffusion functions k and ϑ, which are of important request for the analysis of the problem. We will conclude this section by stating the main result of this work. Section 3 is devoted to the resolution of an η-approximative problem which returns to (5) as far as η −→ 0 + where we establish existence of approximate weak solutions. In the last section, we prove that the obtained solutions converge to some weak solution of the original problem as η vanishes.
2. Tools and principle result.

Notations and function spaces.
In this paragraph we fix notations and provide useful definitions of function spaces. Given a Banach space X of scalar functions then X d and X d×d are the space of vector-valued and tensor-valued functions with d and d 2 components whose each of them belongs to X, respectively. The space X * stands for the dual space of X and for the duality pairing we use the notation < . , . > X,X * . But sometimes, for simplicity, the latter symbol will be skipped if we are sheltered of any possible confusion. Moreover, if X is reflexive then X weak denotes the space equipped with the weak topology.
The punctuation marks . and : stand for the scalar product of vector fields and tensors (which will take the bold character), respectively. While for the product of scalar functions no particular symbol will be made for that. We denote (. , .) and (. , .) ∂Ω the inner product in L 2 (Ω) and L 2 (∂Ω), respectively.
For q ∈ [1, ∞] and k ∈ N, we use the standard notation for the usual Lebesgue and Sobolev spaces (L q (Ω), . q ) and (W k,q (Ω), . k,q ), respectively. The space of continuously differentiable functions on Ω is denoted C ∞ (Ω). Recall that C ∞ (Ω) is dense in both L q (Ω) and W k,q (Ω).
We also set the following spaces : Recall the conjugate exponent q := q q−1 relatively to the notion of dual Sobolev spaces, so we sometimes use the following notations : (Ω) and W k,q n,div (Ω) * := W −k,q n,div (Ω). We also introduce relevant spaces of Bochner-type, which will be of further use.
Auxiliary lemmas. In this part, we recall some useful inequalities and technical results that will be of interest in the sequel. Let Ω be a bounded domain of class C 0,1 and q ∈ (1, ∞). Then there exists a positive constant K, depending only on Ω and q, such that for all u ∈ W 1,q (Ω) d ∩ L 2 (∂Ω) d the following inequality holds For the proof see for example [6]. Another classical inequality of Poincaré needs to be recalled regarding the context of the work.

Lemma 2.2 (Poincaré's Inequality).
Let Ω ⊂ R d be a bounded domain and q ∈ (1, 2). Then there exists a positive constant C P,q , depending only on Ω, d and q, such that for all u ∈ W 1,q (Ω) d the following inequality holds where the Poincaré's constant is given by We refer the reader to the second chapter in [16].
A key point of obtaining direct compactness results is the so-called Aubin-Lions lemma.
Lemma 2.3 (Aubin-Lions). Let 1 < α < ∞, 1 ≤ β ≤ ∞ and X 0 , X 1 , X 2 be Banach reflexive separable spaces such that See [32] and [33] for the proof. Next, we state the Helmholtz decomposition for vector fields and some related estimates for the Neumann problem of the Laplace operator. Let with Ω ∈ C 1,1 , then the Helmholtz decomposition is given by where h u is solution of the following homogeneous Neumann auxiliary problem Clearly, from (17) and (18), we have div u div = 0. Important results from the standard L q -theory of the solvability of problem (18) are stated below: According to the Helmholtz decomposition, it yields the direct sum of the following Lemma 2.4.
Let Ω ∈ C 0,1 and q 1 ≥ 1 and r, q 2 ∈ (1, ∞). Let S be the set defined by The result is due to Bulícek et all. in [6], see lemma 1.12 and corollary 1.13 that follows.
3. Assumptions on non linearities and main result. In this paragraph, we formulate the assumptions involved in our study regarding the non linear functions ϑ, k, S, χ and f . Later we state our main theorem. A 1. Assumption on ϑ. We suppose that ϑ : R + → R + is continuous and verifying for some constants α 1 , α 2 ∈ (0, ∞) Secondly, we state our assumptions on the diffusion function k of the c-equation. A 2. Assumption on k. We shall assume that k : R + × R + −→ R + is a continuous function and satisfying where β 1 , β 2 are fixed positive constants. Lastly, conditions on the stress tensor S of the stress tensor of the momentum equation are given below. A 3. Assumptions on S. As it is mentioned above, we have considered that the viscosity of the fluid is a nonlinear function of m, c, and |Dv|.
Moreover, we shall assume that S : for some non negative constants γ 1 , γ 2 . See [7,6,25,24] for an extensive discussion about similar considered structure of the shear stress.
Note that (24) expresses the coercivity of the stress tensor S while (25) expresses polynomial growth's behavior. The last inequality (26) indicates that it is strictly monotone.
Moreover, suppose that the potential Φ satisfies We say that a triplet (m, c, v) is a weak solution to the problem (5)-(9) if m ∈ L ∞ (0, T ; L 2 (Ω)) ∩ L 2 (0, T ; W 1,2 (Ω)), and satisfying for all ζ ∈ L 2 (0, T ; W 1,4 0 (Ω)), ξ ∈ L ∞ (0, T ; W 1,∞ 0 (Ω)) and ϕ ∈ X r, To the best of our knowledge, our work represents the first mathematical study for a non linear chemotaxis NSE (5) where material coefficients depend on the cell density, the oxygen concentration and the shear-rate. Besides, this is accompanied by a long time existence of weak solutions result which holds for large data and smooth bounded domains fulfilling Navier's slip boundary conditions and suitable nonlinear functions.
We would like to mention here that the lower bound for the parameter r ( 2(d+1) d+2 ) comes from the lemma 2.4 where compactness of the velocity field at the boundary in L 2 (0, T ; L 2 (∂Ω) d ) is verified only over this value.
4. Approximative problem. We will adopt the Faedo-Galerkin method. In order to guarantee the convergence from the Galerkin approximations to the original problem we need to test the balance of linear momentum (5) 3 by the velocity field v. But we remark that there is a default of integrability in the convective term for r < 2, there is no need for such approximation for r ≥ 2. So to overcome this difficulty we strongly think to regularize v in the convective term using a suitable divergenceless mollifier v η in order to achieve integrability of the concerned term. Manifestly, this mollification will also be used in the two other convective terms, those of the equation of the cell density and of the chemical concentration.
Hence, for this reason we introduce the standard regularization by means of mollifiers.
Let η > 0 be a (small) real fixed number. Recall a standard regularization kernel through a function ρ satisfying : ρ ∈ C ∞ (R d ) non-negative, radially symmetric and We also consider the function ϕ η : Ω → {0, 1} fulfilling Given v ∈ W 1,r n (Ω), we define v η := ((ϕ η v) * ρ η ) div where u * ρ η is the standard convolution product of an integrable function u with the kernel ρ η whose the support is located in a ball of radii η. The subscript div comes from the Helmholtz decomposition.
Observe that for an admissible v ∈ L r (0, T, L r n,div (Ω) d ) then we have v η −→ v in L r (0, T, L r n,div (Ω) d ) as η → 0 + . We would like to briefly draw the attention of the reader to some simple computations in relationship with the constructed mollification v η .
Using the Green's formula, we see that

ON A COUPLED CHEMOTAXIS NON-NEWTONIAN FLUID 917
Similarly, and by the same way we get Now, we are in a position to define our approximative problem : subject to the same boundary conditions (6)-(8) and fulfilling the same initial conditions (9).
We look for Galerkin approximations (m N , c N , v N ) being of the form solution of the following system of differential equations We also define the extensions Let θ 1 N be a regularization kernel of radii 1 N . We define the initial data m N 0 and c N 0 for the above system by regularizingm 0 (x) andc 0 (x) as follows Note that with such a definition we have On the basis of (10) and the definitions (53) we have In addition, we require that v N is subject to the initial condition Existence of solution to the system of equations (50), (53) and (57) can be handled for at least a short time interval [0, t] (for t < T ) with the help of the classical Carathéodory theory (see Chapter 30 of [41]). The uniform estimates that we will derive in the next subsection enable us to extend the solution onto the whole time interval [0, T ].

4.1.
First a priori estimates. Next, we establish a priori estimates for energy inequalities. We shall begin with the convection diffusion equation of the oxygen concentration.
In this temporary context, for t ∈ (0, T ], we denote Q t := [0, t] × Ω and Σ t := [0, t] × ∂Ω. Multiplying the j-th equation in (51) by b N j , using (23) and (27) 1 and taking the sum over j = 1, ..., N we get Here, the convective term vanish with the same manner as in (39). Secondly, multiplying the j-th equation in (50) by a N j , summing over the indices j, we obtain with the help of assumptions (22) and (27) Adding (58) to (59) and using the Hölder's inequality and the Poincaré's inequality (15) we obtain where C P,2 is the Poincaré's constant given in (16). On the basis of the conditions (28), one can after integrate over the time interval (0, t) and absorbing get the following estimate Consequently {m N } is uniformly bounded in L ∞ (0, t; L 2 (Ω)) ∩ L 2 (0, t; W 1,2 (Ω)).
The same uniform bound can be handled for the sequence {m N }.

HAFEDH BOUSBIH
Finally, in order to handle estimates for the velocity field approximations, we multiply each k-th equation in (52) by d N k and take the sum of them. We obtain 1 2 (65) Note that we have used similar results to those established in (11) and (38), in addition to the r-coercivity property (24) of the stress tensor S.
Thus, with the help of (24), Hölder's, Korn's and Young's inequalities and Sobolev embedding we infer that , choosing ε small enough, absorbing the second term of the right hand side and integrating over time interval (0, t) we obtain thanks to (61) and the assumption (43) the Korn's inequality which yields {v N } is uniformly bounded in L ∞ (0, t; L 2 (Ω) d ) ∩ L r (0, t; W 1,r n,div (Ω) d ), (67) {v N } is uniformly bounded in L 2 (0, t; L 2 (∂Ω) d ). (68) Note that when d = 3, then using interpolation inequality (see [2]) with p 1 = 6r−6q+2rq q(5r−6) and p 2 = 3r(q−2) q(5r−6) for q ∈ In the case of d = 2, if r ∈ [1, 2) then 2r 2−r ≥ 2, so we have (see also [2]) Then with the help of (67) we conclude that {v N } is uniformly bounded in L 2r (0, t; L 2r (Ω) 2 ). Consequently, and the same estimates hold for v N η ⊗v N thanks to the regularisation effect. Remark that the last two bounds can be unified into Using the polynomial growth properties (23) and (25) for k and S, we can see that Collecting all the estimates established above, one can deduce that there exist (m η , c η , v η ) and a not relabeled subsequences (m N , c N , v N ) such that for some S η ∈ L r (0, t; L r (Ω) d×d ) and K η ∈ L 1 (0, t; L 1 (Ω) d ).
At this stage we need to control the time-derivative of m N , c N and v N in order to ensure compactness by means the famous Aubin-Lions lemma.
On the basis of the estimates established above we can easily handle (uniformly with respect to η) that Let us mention that the large choice of the space W 3,2 (Ω) d comes from the continuous injection W 3,2 (Ω) d → W 1,∞ (Ω) d . Since we need to bound the gradient of the test function in L ∞ in space. Therefore, this implies Hence, on the basis of the obtained estimates and applying the Aubin-Lions lemma (2.3) we get c n → c η strongly in L p (0, T ; L p (Ω)), and Thanks to (63) and (95) we can prove that On the other hand using lemma 2.4 we can reach Thus due to the continuity of the functions ϑ, χ and f , the convergence results (80), (91-96) and (100) we can ensure the establishment of the weak formulations (34), (35) and (36) except for the terms involving the stress tensor of the momentum equation and the concentration-flux. This is the subject of the next paragraph.
First, let us identify the limit in (52). Note that on the basis of (87), (93), (94), (95), (96) and (100) it follows that for all ϕ ∈ X r, (d+2)r 2d div . Observe that the last integral in (101) is well defined since ϕ is free divergence and a control of that term can be handled by the same way as in lat term of (66).
The vector field v η being an element of X r, (d+2)r 2d div then it is absolutely legitime that we can replace ϕ by v η in (101). Doing so, we get Now, reviewing (65) and taking the limit as N −→ ∞ we deduce with the help of weak lower semicontinuity of norms that lim sup On the other hand, using (25), the continuity of S regarding m and c and the strong convergences (94) Recall that S N := S(m N , c N , Dv N ). Using (26) and the continuity of S regarding the first and the second variable, we have Observe now that choosing A = Dv η ± λDu for arbitrary u ∈ W 1,r n (Ω) d and λ > 0, dividing the result by λ and after that taking the limit as λ −→ 0 we obtain which means that Consequently, on the basis of the strong monotonicity property (26) of S Dv N −→ Dv η a. e. in Q T .
Now, it remains to show that the weak limit K η obtained in (86) corresponds really to k(c η , |Dv η |)∇c η . That is Let us mention that the flux-vector can be rewritten as We shall prove that strongly in L 2 (0, T ; L 2 (Ω)).
Observe that by (23) and (63) we have Indeed having already point-wise convergence of the sequences {(m N , c N , Dv N )} N , the Vitali theorem enables us to check up (109). Consequently, with the help of (99) we deduce (108) which finishes the proof of Lemma 4.2.
5. Proof of Theorem 3.2. Now, we launch the proof of the main theorem of this paper. In the previous section, we have already collected important information for the constructed approximate solution (m η , c η , v η ) of the problem (41). Our goal here is to show that the obtained solution converges (as η −→ 0 + ) to a triplet (m, c, v) which should be solution of the original problem (5)- (9).
On the basis of the uniform estimates established in the last section and standard tools such as weak lower semi-continuity of norms and the Fatou's lemma we have uniformly with respect to η. Consequently, with the help of the Aubin-Lions lemma we can deduce that up to (not relabeled) subsequences S η S weakly in L r (0, T ; L r (Ω) d×d ), where the values of p are given by (97). The continuity of the nonlinear functions ϑ, χ, and f , together with (120) and (121) yield the required result.
It remains to show that the S := S(m, c, Dv) and K := k(c, |Dv|)∇c. The key point to fulfill that purpose is to establish the point-wise convergence of Dv η . This is the subject of the next paragraph.
Certainly g η is positive for all η. From (113) we deduce, uniformly to η and up to a certain fixed real number K > 1, that A crucial property useful for the proof our goal was proved in [6]. Here we give the statement.
Moreover, we can handle From the Helmholtz decomposition (17) the vector field u j could be written as u j = u j div + ∇h u j . Moreover the estimate (19) 1 yields that T 0 h u j r 2,r dt ≤ Cε.
By virtue of (132) and (19) we infer that u j div → 0 strongly in L p (0, T ; L p (Ω) d ) for all p < ∞, which with the help of (130)-(132) implies u j 0 weakly in X r, Setting S j = S(m j , c j , Dv j ), reconsidering the monotonicity property (26), it follows that Now, observe that, letting j grows great enough ( in such a way that one can ensure that |v j −v| l < 1), the second term in the right hand side of (137) may be decomposed as follow First, setting ϕ = u j div in the weak formulation (36) of v j , we get Let us investigate the T i 's.
The limit of the first term of the right hand side of (140) vanishes, see the footnote at the page 79 of [6].
For T 2 , we have ∇v j L r (L r ) .
The last two norms are uniformly bounded, so by (135) we see that lim j→∞ T 2 = 0.
Similarly, thanks to (117) and (131) we have Lastly, we can see that T 4 have the same behavior as the previous terms

Consequently, we infer
and thus the proof of Theorem 3.2 is complete. We finish by mentioning that using the coercivity assumptions of the flux vectors of the cell density and oxygen concentration (22) 1 and (23) 1 and by means of the maximum principle we have m(t, x) ≥ δ m and c(t, x) ≥ δ c a.e. (t, x) ∈ Q T . (144)