Analysis of the 3D non-linear Stokes problem coupled to transport-diffusion for shear-thinning heterogeneous microscale flows, applications to digital rock physics and mucociliary clearance

This study provides the analysis of the generalized 3D Stokes problem in a time dependent domain, modeling a solid in motion. The fluid viscosity is a non-linear function of the shear-rate and depends on a transported and diffused quantity. This is a natural model of flow at very low Reynolds numbers, typically at the microscale, involving a miscible, heterogeneous and shear-thinning incompressible fluid filling a complex geometry in motion. This one-way coupling is meaningful when the action produced by a solid in motion has a dominant effect on the fluid. Several mathematical aspects are developed. The penalized version of this problem is introduced, involving the penalization of the solid in a deformable motion but defined in a simple geometry (a periodic domain and/or between planes), which is of crucial interest for many numerical methods. All the equations of this partial differential system are analyzed separately, and then the coupled model is shown to be well-posed and to converge toward the solution of the initial problem. In order to illustrate the pertinence of such models, two meaningful micrometer scale real-life problems are presented: on the one hand, the dynamics of a polymer percolating the pores of a real rock and miscible in water; on the other hand, the dynamics of the strongly heterogeneous mucus bio-film, covering the human lungs surface, propelled by the vibrating ciliated cells. For both these examples the mathematical hypothesis are satisfied.

micro-devices or highly viscous flows (food industry). The complexity of the physics involved, especially in natural sciences such as environment and life, leads to a very wide class of models with non-linear coupling between meaningful quantities.
In particular, considering polymers modeling, there are multiple ways to take into account the complexity of the non-linear behavior of the material, that is to say its rheology. Among the most common models, one can mention visco-plastic models involving a yield stress (such as Bingham fluids [3]), visco-elastic involving relaxation of strain rate and/or stress tensor possibly by means of upper-convected time derivation (such as Oldroyd-B or Maxwell models [39]), potentially extended to non-linear strain-stress relationships (Giesekus model [27,46]). Rheology models also include non-linear viscosity models of shear-thinning or shear-thickening fluid, such as the Ostwald law (or power-law), the Carreau model [8], or Herschel-Bulkley [29] model involving a viscous-stress threshold (thus also featuring visco-plasticity).
The present study focuses on micro-scale modeling, involving the physics of highly viscous flows modeled by the stationary Stokes equation. The viscosity is assumed to follow the non-linear Carreau law (shear-thinning rheology), variable in space and time, modeling heterogeneities. This heterogeneous medium is in motion resulting from the velocity solution of the Stokes equation and also diffused. This results in a three-dimensional model coupling diffusion and transport, the Stokes equation, whose non-linear viscosity takes into account the rheology, with parameters varying with respect to space and time, due to the miscibility.
More precisely, a smooth immersed deformable body B(t) moving in the three-dimensional domain Ω ⊂ R 3 is considered interacting with a highly viscous, incompressible, miscible and shear-thinning fluid in a timedependent fluid domain Ω \ B(t). Such a prescribed domain motion is especially suitable for phenomena driven by the solid, when the fluid feedback can be neglected. This avoids at the same time the parametrization of the fluid/solid interaction, which is sometimes difficult to carry out or subject to large uncertainties (such as elasticity features). In order to set up models around these geometries, the following set of R 4 is introduced, for a final real time T > 0, defining the fluid domain: The shear-thinning feature induces a non-linear Stokes equation, whose solution provides the fluid velocity u and its pressure p. The heterogeneity means that the fluid parameters depend on a constituent concentration α and the miscible aspects are modeled by the diffusion of this constituent α. These modeling assumptions lead to the following system of partial differential equations: where D(u) = (∇u + t ∇u)/2 is the strain tensor and µ is the variable viscosity, with given initial conditions for α. The diffusion coefficient σ is a strictly positive real number. The boundary conditions are defined in the sets of conditions (2.4), (2.5) and are quite usual: Dirichlet, Neumann, a combination of both, or periodic conditions on u, and typically no-flux (homogeneous Neumann) on α on the body inner boundary ∂B(t) \ ∂Ω and on the domain boundary ∂Ω. For these conditions, the velocity is required to match a prescribed solid velocityū almost everywhere on the body: a one-way fluid structure interaction is assumed, which means that the counter force exerted by the fluid on the solid is neglected, which is valid for the targeted applications.
The classical example is the Ostwald law µ = Kγ q−2 /2 with the exponent q 2 and the consistency K > 0 (this law is also called power-law). It leads to a diffusion operator −div(|D(u)| q−2 D(u)) in equation (1.1), whose analysis is very close to the q-Laplacian theory. Its bounded version, exhibiting two physical Newtonian bounds and sometimes used as a regularization of the Ostwald law, is the following Carreau law: where 0 < µ ∞ µ 0 (α) and q(α) ∈]1, 2] for all α, where µ 0 , β and q are C ∞ functions of α, and where D is an arbitrary 3 × 3 matrix, in practice applied to the shear rate D(u). With respect to α, the viscosity at rest µ 0 is assumed to be increasing, q is assumed to be decreasing (N = q + 1 is called the fluid index), and the material time β is assumed to be monotonic, for anyγ. These conditions include constant values. In practice α is chosen in [0, 1] but is free to be in any closed interval of R + . It is noticeable that the need of an exponent q > 1 excludes the Herschel-Bulkley model of visco-plasticity, which exhibits a 1/γ singularity such as q = 1.
In order to go forward to practical aspects of real-world applications, this system of partial differential equations is approximated by means of the following penalization of B(t) to satisfy the prescribed divergencefree velocity u inside the body: where 1 B(t) denotes the characteristic function of B(t). We set ∂ n α ε (t, ·) = 0 on ∂B(t) \ ∂Ω for any time t, which ensures that there is no diffusion of fluid material inside the body. The related boundary conditions are detailed in the sets of conditions (2.4) and (2.5), as for the the non-penalized problem.
The penalized equations set (1.3) is definitely of interest in order to perform numerical simulations, as it makes possible the use of numerical methods which avoid to mesh the body boundary ∂B(t) and its displacements.
The present study provides an original analysis of the non-linear system (1.3). The convection-diffusion equation analysis is performed in a time-dependent domain where classical methods can not be applied, and a mapping to a stationary domain with time-dependent operators is provided. Indeed, Galerkin methods cannot be used in this framework (an evolutionary spectrum can lead to possible multiple eigenvalues by means of path crossing). Hence an innovative technique is presented involving the limit of operators with piece-wise constant coefficients. Moreover, the equations (1.3) involving the penalization of the time-dependent domain are shown to have a unique solution, converging toward the solution to the original system of partial differential equations (1.1). Furthermore, this non-linear coupling involves shear-thinning effects leading to a second non-linearity by means of a generalized Stokes operator. This is analyzed using p(x, t)-Laplacian analytical techniques (variables exponent Lebesgue spaces and Luxemburg norms typical from Orlicz spaces) which is a novelty for the analysis of the penalized formulation of fluids with rheology.
This article is consequently structured in four main parts, describing mathematically and numerically this coupled problem. The first part is a global overview which presents the model and the main mathematical results in Section 2.
The second part presents in Section 3 the analysis of each PDE involved in model (1.3), for time independent and time dependent domains (generalizing [14,18] to non-linear framework). This section introduces the transformation of the non-linear viscosity into a variable q-Laplacian problem, in the spirit of the work from Diening et al. [19][20][21] for electro-rheological fluid, where they introduced the theory of variable exponent Lebesgue spaces. Still about this problem, Sections 3.1 and 3.2 use p-Laplacian methods [33] to prove solutions for the non-coupled Stokes problem. In Section 3.3 we study the rather simple convection-diffusion equation written in the non simple case of a time-dependent domain Ω \B(t). For this study, in order to improve the readability of the proofs, various technical arguments are presented in Appendix A, including the Luxemburg norms and the required elements of non-linear analysis, and the regularity and L 2 estimates of the geometry lifting are provided in Appendix B.
The third part of this article focuses on the mathematical aspects of the coupling in Section 4, and shows the proofs of Theorems 2.4 and 2.6. They provide the existence of the solution of the penalized coupled problem (1.3), its convergence toward the solution of the initial problem (1.1) through the asymptotic analysis, and its solution's existence at the same time. To do so, an innovative technique of domain mapping is presented involving the limit of operators with piece-wise constant coefficients, since such a mapping leads to time-dependent differential operators and thus Galerkin methods cannot be used in this case (an evolutionary spectrum can lead to possible multiple eigenvalues by means of path crossing). To simplify the readability, some proofs are also provided in Appendix C. Such a problem has been studied by Boyer and Fabrie [5] for a coupling between Navier-Stokes and transport equations, and also by the present authors in [14] for the Newtonian case, that is to say for a viscosity independent of the shear-rate (i.e. a linear Stokes problem). The theoretical study of a fixed penalized domain was performed by Carbou and Fabrie [7] for the Navier-Stokes equation.
The fourth part is dedicated to numerical simulations of such problems in real world configurations, showing that all the hypotheses set up before are necessary and satisfied. The α-dependency of the parameters in (1.2) is subject to the mixing model. Firstly, an example from geosciences is provided in Section 5.1, focusing on the heterogeneous rheology in a non-moving complex geometry. Secondly, a life science example related to the mucociliary clearance of human lungs is performed in Section 5.2, and involves the shear-thinning mucus subject to mucin heterogeneity around epithelial ciliated cells in motion. The numerical parameters used in the simulations are available on Table 2.

Domains, equations and boundary conditions
Let T > 0, B(t) a domain of R 3 depending on time t ∈ [0, T ], and Ω a sufficiently smooth open set of R 3 . If one considers a problem with one periodic boundary condition at least, let Q be an open set of R 3 whose boundary is sufficiently smooth. An equivalence relationship ∼ (whose kernel is denoted G) is also defined such that the quotient space Ω ≡ Q/G, with a topology induced from Q, is sufficiently smooth.
Consequently, thanks to the quotient topology, the domain is smooth even if the faces involving periodic boundary conditions exhibit corners with the remaining part of the domain boundary. The Figure 1 illustrates the different configurations of acceptable and non acceptable domains with adequate boundary conditions.
One typical example of interest for the present study is a cube ]0, L[ 3 on which two periodic boundary conditions are set: in that case, Q = R 2 ×]0, L[ and G = LZ 2 × {0}, so that Ω ≡ Q/G is smooth.
Four sets of R 4 are now introduced: × Ω the whole domain in time and space coordinates, The problem (1.1) rewrites in these spaces: On the boundaries of B(t) and Ω, where the external normal field ν is defined, we introduce u ⊥ = (u · ν)ν and u = u − u ⊥ = −(u ∧ ν) ∧ ν, respectively the normal and the tangential parts of u ∈ R 3 (the hypothesis H1 introduced after will forbid the fact that ν may be undefined at ∂B(t) and ∂Ω intersection). In order to easily setup up the boundary conditions of problem (1.1), let Γ 1 , . . . , Γ I denote the I connected components of ∂Ω and their union 3) Figure 1. Acceptable domains Ω whose boundaries have to be sufficiently smooth, and how to manage corners in domain regularity. Pairs of symbols denote periodicity.
on which a set of θ i ∈ {0, 1}, constant on each connected component Γ i , is introduced in order to naturally switch between Neumann and Dirichlet boundary conditions on u . It can be noticed that due to the quotient topology, periodic boundary conditions do not appear among the Γ i and consequently are not part of the listing of boundary conditions. Nevertheless, the mention periodic otherwise is written in order to remind that periodic conditions are included in the quotient topology. This leads to the following set of initial and boundary conditions: . . , I}, u is periodic otherwise.  in Ω thanks to its definition as a quotient. The condition "periodic otherwise" is a reminder throughout the article of this definition. The viscosity µ and its regularity, with respect to α, is assumed to satisfy: , ∀D ∈ M 3,3 (R) and ∀α ∈ R, µ 0 ∈ C 1 (R) ∩ W 1,∞ (R) and ∀α ∈ R, µ 0 (α) µ ∞ > 0, β ∈ C 1 (R) ∩ W 1,∞ (R) and ∀α ∈ R, β(α) 0, q ∈ C 1 (R) ∩ W 1,∞ (R) and ∀α ∈ R, q(α) ∈ [q − , 2] where q − > 1.
According to the modeling, the penalized problem (1.3) can be rewritten: with its initial and boundary conditions following the equation set (2.4) and (2.5).
In the following the norm on the spaces L p (Ω; R n ) will be simply denoted by · L p (Ω) to ease the notations. Finally the following functional spaces are introduced: Definition 2.1. We let F(Ω; X) denote the set of functions defined on Ω with values in X and with the usual H 1 -norm, the periodicity of function and normal field ν being implicitly taken into account by definition of Ω. Thanks to the generalized Poincaré's inequality (see [4] Sect. III.6), there exists two constants depending only on Ω such that for all w ∈ V (Ω)

Main results
The following existence and uniqueness results require some regularity assumptions on the obstacle's time evolution (H1) and (H2) and on the velocityū inside the obstacle (H3). We first introduce a new definition: One can notice that ifB ∩ ∂Ω = ∅ and B has the uniform C k -regularity property then B has the C k -bi-regularity property. Moreover, in the mucociliary clearance application, there exists (l, l ) ∈ (R + ) 2 such that 0 < l < l < L and Γ × [0, l] ⊂ B ⊂ Γ × [0, l ], i.e. the bronchial wall is assumed to be the fixed domain Γ × [0, l] where cilia are attached. Cilia are also supposed to never reach the top of the domain Ω. There exists then domains B having the C k -bi-regularity property. Moreover the inner boundary Σ = ∂B \ ∂Ω of B is smooth and may be smoothly periodically extended. It splits Ω in two regular open sets.
Furthermore, we need existence and regularity of solutions for the convection-diffusion equation in the deformable domain, obtained by rewriting the equations in a stationary domain thanks to a diffeomorphism Ψ. The following assumptions on the regularity of Ψ are required for the continuity of field v and operator A in equations (3.7) and (3.8): Hypothesis 2 (H2). There exists a function Ψ ∈ C 1 ([0, T ]; C 2 (Ω; R 3 )) such that for all t ∈ [0, T ] -Ψ(t) preserves the C 2 -bi-regularity property.

Remarks:
-The domain B(t) has the C 2 -bi-regularity property. -In the following, Ψ(t)(ω) is also denoted Ψ(t, ω) or Ψ t (ω). -The preservation of the measure also implies that   A last technical assumption on the viscosity is required: Hypothesis 4 (H4). We assume that We can now define what is the meaning of a weak solution of the penalized problem (2.6) with boundary conditions (2.4) and (2.5), and then state in Theorem 2.4 that there is a unique weak solution of the penalized problem under the assumption (H1)-(H4).
Furthermore, many technical arguments are provided in Appendix A in order to build the proof of both theorems, detailed in Section 4.

The generalized Stokes problem in a non moving domain
We are now interested in the following non-coupled penalized Stokes problem where u satisfies the boundary conditions (2.4) and ε > 0. Since there is for instance no coupling we work with a viscosity µ that only depends on x and u. In this problem, the given data are (µ 0 , β) ∈ (F(Ω; R + )) 2 , q ∈ P(Ω) the set of all measurable functions on Ω with values in [1, +∞], see also Section A.1, B an open set in Ω,ū ∈ F(B; R 3 ) fulfilling divū = 0, f ∈ F(Ω; R 3 ), and the unknowns are D(u) = ∇u + t ∇u /2 with u ∈ F(Ω; R 3 ) and p ∈ F(Ω; R).
First a time-independent problem is considered. To obtain the following results (H1) is assumed and an additional assumption is needed onū (H5) which is a simpler time-independent version of (H3).
Hypothesis 5 (H5). In order to ensure well-posedness of equation (3.1),ū does not need to match the velocity of the B(t), which is accessible via the applicatin Ψ in hypothesis (H2). Nevertheless, the goal of the coupled model and assumption (2.5) is to avoid the flux of α into the body B(t), soū is required to match the body velocity field at the body boundary when the full systems (2.2) and (2.6) are considered. The remark holds for both the penalized problem (3.1) and the original non-penalized problem (3.5).

Existence and uniqueness
Definition 3.1. We associate to the problem (3.1) the following variational formulation for all v ∈ V (Ω), Then there exists a unique (u, p) ∈ V (Ω) × L 2 0 (Ω) solution of (3.1). Proof. This proof is inspired from various results from [4], Section II.3.2 and from [14] but is adapted to the non-linear case using Theorem A.14. We use the notation V = V (Ω). We let Then It follows that Since µ 0 µ ∞ , thanks to the Korn's and Poincaré's inequalities: Consequently is well defined, linear and continuous since so the function L is in V . Thanks to the non-linear Lax-Milgram Theorem A.14 there exists a weak solution The solution is then unique and the end of the proof is classical. Indeed, the first part of the proof implies Testing this function against a function φ ∈ H 1 0 (Ω) such that divφ = 0 we obtain: De Rham's theorem induces the existence of a function p ∈ L 2 0 (Ω) such that g = −∇p. We now verify that the last boundary conditions are fulfilled. The function g + ∇p rewrites From this expression we deduce that 2µ(x, u)D(u) − pI lies in and admits a normal trace in Since divφ = 0 and thanks to (3.2) expressed with v = φ and to the boundary conditions already fulfilled by u and φ we obtain Since Γ i is assumed to be flat (no curvature), so ν is constant and can be extended to a tubular neighborhood of Γ i , this writes: and we finally get, by means of tangential projection, that is the solution to the Stokes problem (3.1).

"Harmonic" prolongation ofū
To obtain precise estimates on (u, p) we need a prolongationũ ofū in the whole domain Ω. We study the following system: u satisfies the boundary conditions in (2.4).

Estimates on (u, p)
Under the assumptions of Theorem 3.2 and Proposition 3.3 we let (u, p) be the unique weak solution to (3.1) given by Theorem 3.2.
Since the mean of p is zero, Poincaré's lemma gives

The generalized Stokes problem in a time-dependent domain
Let T > 0. In the following the domain B depends on the time t ∈ [0, T ] and we need to track the dependency on time t of the constants (especially the ones in all Sobolev related theorems). Using classical trace and lifting theorems (on the domain B(t) or Ω \B(t)) there appear to be some constants depending on B(t) and Ω \B(t). Thanks to the Proposition A.6 we can estimate the time-behaviour of these constants.
Theorem 3.7. We assume (H1)-(H3). Let Then there exists a unique weak solution (u, p) to the problem We also have the following inequalities:
Proof. In order to deal with the variable domain we use the diffeomorphism Ψ defined in hypothesis (H2) to rewrite the equation in a fixed domain but with time-dependent differential operators. We let β(t, y) = α(t, Ψ(t, y)) for all t ∈ [0, T ] and y ∈ Ω \B. In the following we note Ψ t (y) instead of Ψ(t, y). β fulfills where ν is the outward unitary normal to Since we have assumed that Ψ preserves the measure we have | det J t (y)| = 1 for all (t, y) ∈ [0, T ]×Ω hence we get the expression of the operator A. This also allows us to work with the classical Lebesgue and Sobolev spaces instead of weighted ones (with a dependence on time). This assumption also leads to divJ −1 Instead of using a Galerkin's method where the time-dependence of the eigenfunctions of A will be problematic we discretize in time the equation to prove the existence of solutions. Let N ∈ N * , h = T N and t n = nh. We are looking for β n (y) approximation of β(t n , y), which is the solution of the following implicit numerical scheme: By taking the scalar product of the equation with β ∈ H 2 (Ω \B) fulfilling the same boundary conditions and integrating by parts, we get This formulation is well-posed in H 1 (Ω \B). Thanks to the regularity of Ψ , the operator L is linear and continuous on H 1 (Ω \B) and a is bilinear and continuous on Thanks to the assumptions on u and Ψ, (v n ) n∈{0,...,N } is bounded in ∞ ({0, . . . , N }; L 6 (Ω \B)), then for h independent of n and small enough, the form a is coercive. Thanks to Lax-Milgram theorem we then have the existence and uniqueness of a solution β n+1 ∈ H 1 (Ω \B) for all n ∈ {0. . .N }.
Thanks to all these convergences we deduce that β is a solution to (3.7) with all the expected boundary and initial conditions. Since we assume that β 0 ∈ H 2 (Ω \B) we have indeed more regularity on β. To prove this we write all the energy estimates we have on β, available in Appendix B, which finishes the proof.
Remark. Since we use Lemma A.2 the first inequality in Theorem 3.8 is not sharp. The second one is a very rough bound.

Study of the coupled penalized problem
To prove the existence of a solution to this coupled non-linear problem, we introduce the sequence (α n , u n , p n ) defined by (1) α 0 =α 0 the extension by 0 of α 0 outside Ω \B(0), (2) (u n , p n ) is the solution to × Ω, u n satisfies the set of boundary conditions (2.4).
In the following we prove that G is a contraction in a well-chosen space. One can notice that the extension of α n by 0 means α n+1 = 0 in O p .

Study of G
Let E = L ∞ (0, T ; V (Ω)). The operator G maps E to E thanks to Theorems 3.7 and 3.8. Let u and u two solutions to (3.6) with α and α respectively. We need the following proposition, whose proofs are available in Appendix C.  and  Proof. See Appendix C.4.
Proof. The proof is similar to the proof of Proposition 4.4, detailed in Section C.4.
Study of (u, p). Thanks to the weak convergence of u ε towards u in L 2 (0, T ; V (Ω)) and to the strong convergence of u ε towardsū in L 2 (O p ), we have u =ū in O p .
Since u ε u in L 2 (0, T ; V (Ω)) and A α u ε U in L 2 (0, T ; V (Ω)) and A α is pseudo-monotone, we have thanks to Proposition A.13 that A α u = U .

Applications to shear-thinning heterogeneous micro-scale flows
Two microfluidic real-world problems satisfying the coupled problem (1.2,1.3) are addressed in these sections. Firstly, in Section 5.1, the transport and diffusion of Xanthan polymer in a porous rock is studied at the scale of the pores. In this case the domain does not depend on time, but a strongly non-linear miscible fluid is in motion in a micro-channel exhibiting a complex geometry. Secondly, in Section 5.2, a mixture of mucins, the proteins involved in pulmonary mucus, are displaced by means of the vibratile motion of the ciliated epithelium cells covering the lungs. In that last case, the domain is moving and is the dominant motion effect. The full Stokes-transport coupling thus goes beyond the work done in [10] where only a prescribed α(x) was considered.
Dedicated computational algorithms have been developed to compute the studied coupled Newtonian problem [10][11][12], and for the uncoupled non-linear problem [9]. These algorithms involve a hybrid grid-particles framework in order to consider a suitable discretization of each phenomenon: Cartesian grids for the Stokes problem and diffusion, Lagrangian method for the transport.
The resolution of the Stokes problem is based on an iterative projection method. This ensures accurate computations of both the inviscid velocity and the non-linear effects as well as the penalized domain velocity and boundary conditions (an inherent problem of projection methods [28]). It has been shown in [11], and improved numerically in [12], that for a given function α the solution to the generalized Stokes problem is the limit of: where P is the projector on divergence-free fields P(v) = v − ∇ζ with ζ the solution to −∆ζ = −div v satisfying homogeneous Neumann or periodic boundary conditions, and where u * n is the sequence of functions defined as follows: with µ = µ (α, P(u * n )) satisfying the Carreau law (1.2). Adequate boundary conditions are setup for the limit function to satisfy the correct boundary conditions of the Stokes problem. Indeed, the non homogeneous Dirichlet condition based on the Richardson extrapolation leads to the boundary condition u = g. θ ∈ [0, 1] leads to the usual relaxation scheme, while θ = −1 gives the Richardson extrapolation, used in the present simulations. The sequence defined by equation (5.1) reads, once divided by µ µ ∞ > 0: This numerical strategy reduces the grid computations to a sequence of Poisson problems resolution computed with a FFT-based solver FishPack [45] for the numerical evaluation of the projector P, and an algebraic multigrid solver MudPack [1] for the non-separable Helmholtz equation (5.2). This ensures a quasi-linear computational cost with respect to the number of discretization points.
The idea of hybrid Particle-Grid methods is to use particles for the diffusion-transport and a grid-based method for the penalized Stokes problem. A particle is a set (α j , x j , v j ) with a position x j , a volume v j and a weight α j , all of them depending only on time. Particles are a mean to approximate a function in the following sense: for any regularizing kernel ρ. The transport equation can then be replaced by a set of differential equations for the particle weights and locations which, for an incompressible flow, reads: where the velocity field u is computed through the resolution of the penalized Stokes problem. This method in particular exhibits nice stability properties, in the sense that, unlike for traditional gridbased methods, the time-step is not constrained by the grid size (that is to say there is no transport CFL condition). For the high resolution simulations that we will consider, this significantly reduces the computational time. In order to transfer data from grid to particles and the way back, high order convolution with strongly compact supported kernels are used. For instance, we use the M 3 , M 5 and M 4 kernels (depending on whether the positivity is crucial or not) introduced in [36] and developed further in the context of particle methods in [16].
Since the convergence and consistency of this numerical method have been already investigated in [11,12], no further description is provided in the present article. The following section addresses the use of this method for two real-world problems treated for the first time, that is to say solving the full coupled PDEs ((1.2) and (1.3)).

Application to polymer dynamics in geoscience
We are interested in the transport of Xanthan mixed in water (called the solvent, a Newtonian fluid of constant viscosity µ ∞ ) by a viscous flow in a porous rock at its pore scale. This configuration is displayed on Figure 3 and is setup as follows.
First, the solid is defined by the rock core matrix surrounded by a cylindrical cell filling the outer part of the computational box (the solid domain B is the union of these two sets, as displayed on Fig. 3A). The rock core matrix geometry is a Bentheimer sandstone defined by a set of 257 3 voxels of 1 and 0, here obtained by MicroCT X-Ray tomography with a SkyScan 1172 (Bruker). This provides a straightforward characteristic function 1 B of the solid penalized region and leading to model (1.3), the penalized version of (1.1). Despite the fact that the domain here is not moving, this case provides a meaningful application case of (1.3) and illustrates how the hypothesis (H1)-(H4) are satisfied and what they mean in practice (especially (H4)).
Second, the fluid domain is the open complementary set of B (shown on Fig. 3B) and is filled with a concentration αC 0 of Xanthan (α is variable in space and time), where C 0 = 600 mg L −1 or equivalently C 0 /ρ = 0.0713%, where ρ = 0.842 kg L −1 is its density. The quantity α is transported and diffused in the fluid domain with a velocity u solution to the 3D generalized Stokes problem (1.1). When α reaches zero, there is no Xanthan at all at this location.
The Xanthan is a miscible polymer, used in Enhanced Oil Recovery (EOR). It exhibits a shear-thinning behavior, and is ommonly used as an emulsifier or thickener in human food industry, widely used in glutenfree food production. A usual way to produce Xanthan is by fermentation of corn sugar, wheat or soy with  . Xanthan viscosity at rest µ 0 (zero shear-rate) and inverse shear-rate cut-off β with respect to xanthan concentration (g L −1 ): experimental data (from [30]) and best exponential fits, with sigmoid transition from linear to exponential for the µ 0 . the bacteria Xanthomonas campestris. Its rheological properties are quite stable to temperature and acidity variations. These features are close to the Scleroglucan polymer, also used in EOR and reasonably clean since it can also be consumed by bacterial activity, but Xanthan exhibits no yield-stress and thus its micro-scale transport and diffusion is correctly modeled by our coupled Stokes-Tranport/Diffusion system (1.1) involving the Carreau law. The Carreau law (1.2) satisfied by this miscible Xanthan reads: whereγ = √ 2 |D(u)| and µ ∞ = 10 −3 is the solvent dynamic viscosity, here water. In this solvent, the polymer diffusion coefficient is σ = 4 10 −12 m 2 /s (from [44]) and its exponent q = 1.386 can be considered constant due to its very small variations [30], corresponding to a fluid index N − 1 = q − 2 = −0.614.
To our knowledge, for any miscible polymer, the index N and the viscosity at rest µ 0 (zero shear-rate) are slightly or moderately increasing, while the inverse cut-off shear-rate β can be constant, increasing [17] or decreasing [25] with respect to the concentration α. In our case, the Xanthan polymer exhibits an exponential growth in β [30], an almost constant exponent chosen as N = 0.613, an exponential growth of η 0 for large concentrations (αC 0 /ρ > 0.01%) and an affine transition to µ ∞ for small concentrations, modeled by a 1−e −x Figure 5. Xanthan shear-thinning feature: Viscosity with respect to the shear-rate. Two curves fit data from [48], with inverse shear-rate cutoff β obtained by the regressions displayed on Figure 4 and based on [30].
weighting, dominant for α close to 0: and β(α) = A β e αC0/ρB β , with a nominal time A β = 2.4 s. These curves, with respect to experimental values, are displayed on Figure 4. These five coefficients for salt-free and neutral pH solvent are displayed in the Table 1, and the best fitting curves are displayed on Figure 4. The final Carreau law obtained is displayed on Figure 5 and all the related parameters are mentioned in Table 1.
Considering these fluid features, on the one hand q is constant so q (α) = 0 and hypothesis (H4.1) holds, and on the other hand α remains in a compact interval I ⊂ R so µ 0 (α) is bounded, and one gets β /β = C 0 /ρB β so hypothesis (H4.2) is also satisfied. Furthermore, one can notice that the hypothesis (H4.1) is robustly satisfied, that is to say is still satisfied if q is only C 1 (I) and lightly deviating from its value such as m(α) = 1 − q(α)/2 remains positive. For any α ∈ I, let z = 2|D(u)| 2 =γ 2 be the square of the shear-rate and M > 0 bounding |µ 0 − µ ∞ |. Indeed, this gives directly that the expression of hypothesis (H4.1) is bounded by and thus is bounded for any D and α ∈ I. Furthermore, the initial condition of field α is based on a random field generated by FFTMA algorithm [38,41]. This method allows to build 3D periodic random scalar fields by the use of fast Fourier transform. A covariance field C is first computed on the domain Ω according to a Gaussian model. The associated random field of arithmetic average M , geometric average G and variance τ 2 is given bỹ where F (respectively F −1 ) denotes the Fourier transform in R 3 (respectively the inverse Fourier transform in R 3 ) and Z is a gaussian white noise. This randomα is finally cut smoothly on a few voxels from the domain boundary by means of a sigmoid function K so that the support of α =αK1 Ω\B does not intersect ∂Ω or B. Figure 6. Top picture: isosurfaces of Xanthan, colored by the velocity, at level α/α max = 30%. Initial value is displayed with transparency, while the value after 24 time steps is displayed with solid colors. Zooms exhibit locations with highest and lowest velocity regimes (respectively in red and blue). Bottom left picture: similar but colored by the shear rate. Bottom right picture: mean velocity of polymer at concentration higher than a level α p , quantifying the difference between water and diluted polymer. Figure 6 displays the motion of the polymer (initial concentration α is displayed in transparency, and its value after 24 iterations is displayed in solid). On this figure, coloring shows the shear rate on the one hand and the velocity on the other hand. One can see that some regions are transport dominant, while others at small velocity are diffusion dominant, due both to the shear-thinning feature and the domain shape. Thanks to that computation, we can conclude that the polymer is here 3-5 times slower than the water. Despite the fact that the domain is not time-dependent, our coupled model is crucial to capture these physical phenomena.

Application to the bio-mechanics of the human lung mucus
In the lung the tracheobronchial tree is protected from the outer world by the airway surface liquid (ASL). It forms a thin film (∼10−20 µm) on the bronchial walls and inhaled agents (dust, pathogens, pollution particles) are trapped therein to prevent lung contaminations. On the bronchial wall micro-metric cilia are constantly beating to propel this ASL film to the trachea, in order to evacuate inhaled agents. It is then swallowed in the stomach. This carpet of cilia is called the lung epithelium (or ciliated epithelium) and this ASL renewal mechanism is called the mucociliary clearance. When it fails the cough attempts to overcome the ASL accumulation to prevent the pathogens proliferation.
Due to the complex structure of the tracheobronchial tree and the micro-metric scale of the flow, a noninvasive experimental study of the mucociliary clearance is impossible. That is why it was intensively studied numerically in the last decades [43] -see [10] for a detailed and more exhaustive bibliography on recent results. Using experimental rheology, the shear-thinning behavior of the ASL was measured and identified, it allows real data inputs in the computational models [10,32,40]. These studies conclude that the flow can be modeled by a generalized Stokes problem with shear-thinning effects interacting with immersed obstacles: the beating cilia.
The ASL is a non-homogeneous fluid whose viscosity varies due to proteins: the mucins [6]. These proteins are released on the bronchial wall. They maturate and hydrate while they are transported in the ASL film, increasing the viscosity and changing the rheological properties. The viscosity of the ASL and the corresponding rheological parameters are assumed to be depending on a mucins ratio α, quantifying the mucins maturation: when α = 0 mucins have just been released, the fluid is newtonian and its viscosity is equal to water; when α = 1 mucins are fully polymerized and the fluid is non-Newtonian with shear-thinning effects. A continuous transition between these regimes is modeled with the Carreau rheology involving α − dependant rheological parameters, similarly to the previous digital rock physics application.
The parameter dependent Carreau law models the regions with a low density of mucins (the bottom part of the ASL), a high density of mucins (the upper part of the ASL), and the transition between both layers. The lower Newtonian part of the fluid baths the cilia environment and is called the periciliary fluid layer (denoted PCL thereafter) and the upper non-Newtonian part the mucus layer (denoted ML).
The function α is assumed to be the solution to a convection-diffusion equation and the Cilia-ASL interaction is handled using the penalization method, assuming a one-way fluid-structure interaction. This cilia motion can be divided in two parts: the effective stroke (when cilia are polymerizing to beat forward) and the recovery stroke (when cilia are moving backward). A cilium shape and velocity follow the damped wave equation described in [11] with the diameter shown on Table 2. The carpet of cilia is a collection of such cilia, with the metachronal delay used in [10].
Finally a no slip boundary condition is imposed on the bronchial wall and a free slip boundary conditions is imposed at the air mucus interface which is assumed to be flat [10]. This leads to the exactly same problem ((1.2) and (1.3)) analyzed in the previous sections. The rheology is here defined by the following Carreau law using the fluid index N (α) = q(α) − 1: where µ ∞ = 10 −3 is the PCL viscosity, equal to the water dynamics viscosity, and the material time β = 4×10 3 s is constant and independent of the mucin dilution. As established in [10], the fluid index is set to N (α) = αN ML + (1 − α) and the viscosity is set to µ 0 (α) = µ ∞ (K/µ ∞ ) α , where the consistency K and the mucus layer fluid index N ML are displayed on Table 2. Similarly to the digital rock physics configuration, for any α > 0 (and α 1), one gets q − 2 < 0 and q (α) = N ML − 1 so hypothesis (H4.1) is satisfied for α > 0 as the expression decreases with respect to |D|: the limit and bound from (5.6) still holds in the present case. If α = 0, µ 0 (0) = µ ∞ so the expression is zero and consequently hypothesis (H4.1) holds. The second condition (H4.2) is trivially satisfied since β ≡ 0.
Four computations of the ASL propelled by beating cilia are presented in this section to investigate the influence of both the transport and non-Newtonian effects. In order to compare these effects the following simulations are computed (see Tab. 2 for parameters meanings and values): Table 2. Reference parameters for the numerical simulations (gathered from [9,13,23,42,43] Notes. ‡ Only used in computations with mucins advection. † Only used in non-Newtonian computations. * Only used in Newtonian computations. -with both non-Newtonian and transport effects (analyzed in the previous sections), -with non-Newtonian effects but no transport, meaning α(x, t) = α 0 (x) ∀t 0, -with transport but without non-Newtonian effects, meaning µ(α, u) = µ(α) = µ ∞ (1 + κα(x, t)), -without transport and without non-Newtonian effects, meaning where the initial (or stationary, depending on the case considered) function α is defined by: In all these simulations the other parameters remain identical and are gathered in the Table 2. Each computation involves an array of 100 cilia beating asynchronously with a metachronal wave length of 100 µm. Several snapshots of the simulation are presented on Figure 7.
The mucociliary clearance efficiency is quantified by the mean velocity of the mucus layer (the upper part of the ASL). This mean velocity V (t) in the proximal direction (x) is computed at each time step, it is then averaged over six beating cycles (U ): On Figure 8 the evolution of the quantity V (t) is displayed with respect to t for each previous simulation. The associated quantity U is displayed in the legend.
For both simulations with transport (respectively without transport), one can observe that the behavior of V (t) is very similar since the curves are almost overlapping. Differences can be observed at the middle of the cilia strokes (effective and recovery) when extreme values of V (t) are reached (see zoomed parts of the Fig. 8).
The highest velocities are reached simultaneously and advected simulations give higher values, whereas lowest velocities are not reached at the same times: advected simulations are in advance with respect to non-advected. For Newtonian computations lowest velocities are similar while the Newtonian-advected simulation reach a lower one. A consequence of these observations is that U is similar for both non-Newtonian computations, it is higher than for Newtonian cases: a 11% reduction (respectively 32% reduction) is observed for simulations without advection (respectively with advection). The Newtonian simulation with advection presents the lowest U . Hence one can conclude that the non-Newtonian behavior of the mucus significantly increases the mucociliary clearance efficiency, so Newtonian computations underestimate the ASL transport. Moreover at Newtonian regimes the advection has an important role to play: mucociliary clearance is reduced by 21% when mucins transport is neglected.

Conclusion, perspectives and future work
In this study we have presented the analysis of a system of non-linear partial differential equation modeling the microscale dynamics of a miscible and heterogeneous shear-thinning fluid in its solvent, in a moving domain. The analysis has focused on the well-posedness and the convergence of its penalized version toward the original problem in the time-dependent domain.
The penalized problem is shown to be of interest in order to perform numerical simulation. Indeed, two kinds of heterogeneous flows involving miscible polymers at the micrometer scale have been presented, one in a non-moving domain relative to geoscience, and one in a moving domain relative to life science. It can be noticed that the present model is already valid for the heterogeneous mixing of two shear-thinning fluids moving and diffusing into one another, but there is no usual real-world situation that can easily exhibit this.
An interesting perspective is to introduce a relaxation time with second order stress, as described by the Giesekus model [27,46], which would describe a wide class of viscoelastic fluids by means of two parameters. Generalization to such models with time-relaxation (Maxwell, Oldroyd or Giesekus) could be of interest in order to give a mathematical framework of the heterogeneity of such fluids, but would involve a large amount of physical parameters which has little chance to be put easily in practice.
Nevertheless, two extensions of this study may be of direct interest. First by focusing on mass transfer with a diffusion σ dependent on the velocity u. This would be of practical interest for modeling the dispersion of the front between a fluid and a polymer at high Peclet number, as quantified for example in [17]. Second, the case of Herschel-Bulkley fluids involving yields stress, for which the viscosity exhibits a singularity scaling as 1/|D(u)|: the present results don't hold is this case and such study would require a different analysis.

Appendix A. Technical arguments
Lemma A.1. Let m ∈ N * . There exists a constant C depending only on Ω and m such that for all v ∈ H m (Ω; where ν is the outward unitary normal on ∂Ω. The result remains true if we consider a periodic domain.
Proof. These results are straight adaptations of classic ones: The first result comes from the regularity of the operator A t = I − divA(t, ·)∇ with the domain D(A t ) = u ∈ H 2 (Ω \B), (A(t, ·)∇u) · ν = 0 on Γ ∪ σ, u has periodic conditions .
The second estimate requires an adaptation of the proof of Lemma A.1 based on Hodge decomposition. We need to work with the weighted scalar product defined on L 2 (Ω \B; Proof. These two inequalities are deduced from the Sobolev embeddings of H 1/2 (Ω) and H 7/4 (Ω) respectively into L 3 (Ω) and L ∞ (Ω) (see [2], Sect. 4) and from the estimates of the H 1/2 (Ω) and H 7/4 (Ω) norms thanks to the interpolation space theory (see [35], Sect. I).
Let Ω be a regular open bounded set of R 3 and Ω t = Ψ(t, Ω).
Theorem A.5 (Aubin-Lions-Simon). Let B 0 ⊂ B 1 ⊂ B 2 three Banach spaces. We assume that the injection of B 1 into B 2 is bounded and that the injection of B 0 into B 1 is compact. Let p, r such that 1 p, r +∞. Let T > 0. We note Proof. See [4], Section II.5.5. Moreover, there exists a time-independent constant C > 0 such that for all t ∈ [0, T ] and for all g ∈ H 3/2 (∂B(t)) there exists a lifting v of g in H 2 (B(t)) such that v H 2 (B(t)) C g H 3/2 (∂B(t)) .
Let v ∈ H 2 (B(t)) and for all x ∈ B let w t (x) = v(Ψ(t, x)). There exists some constant K a,c > 0 depending only on (a, c) such that ∀t ∈ [0, T ] w t ∈ H 2 (B) and w t H 2 (B) K a,c v H 2 (B(t)) .
In the same way, by writing v(y) = w t (Ψ −1 (t, y)) for all y ∈ B t , there exists some constant K b,c > 0 such that Following [37], Section 3.8, the space H s (Ω), where Ω ⊂ R n and s = k + θ with k ∈ N and 0 < θ < 1, is defined by Using this definition for the space H 3/2 (∂B t ) we obtain in the same way the existence of two constants (C, We now apply the classical trace and lifting theorem to the function w t and obtain the desired results thanks to the previous inequalities. Remark. The lifting theorem can also be generalized in Ω \B(t) with time-independent constants to the case of periodic boundary conditions, v ⊥ = 0 and

A.1. Variable exponent Lebesgue spaces
To deal with the variable q(x)-laplacian type problem we introduce the following spaces (see [20], Sect. 3.1).
Back to the original equation on α we obtain the announced result.
Appendix C. Analysis tools for penalized problem C.1. Proof of Proposition 4.1 We take the difference between the two non-linear equations (3.6) and we get
Estimations in Theorem 3.7 conclude the proof.