On unsteady flows of pore pressure-activated granular materials

We investigate mathematical properties of the system of nonlinear partial differential equations that describe, under certain simplifying assumptions, evolutionary processes in water-saturated granular materials. The unconsolidated solid matrix behaves as an ideal plastic material before the activation takes place and then it starts to flow as a Newtonian or a generalized Newtonian fluid. The plastic yield stress is non-constant and depends on the difference between the given lithostatic pressure and the pressure of the fluid in a pore space. We study unsteady three-dimensional flows in an impermeable container, subject to stick-slip boundary conditions. Under realistic assumptions on the data, we establish long-time and large-data existence theory.


Introduction
The purpose of this study is to investigate mathematical properties of a system of nonlinear partial differential equations (PDEs) developed in [2] to describe processes, such as static liquefaction or enhanced oil recovery, in water-saturated (geological) materials. Such materials can be viewed as two component mixtures consisting of a granular solid matrix and a fluid filling the interstitial pore space. More specifically, we investigate the following system of PDEs: where S and Dv satisfy with τ (p f ) := q * (p s − p f ) + . (1.1e) The system (1.1) coincides with equations (2.24)-(2.26) stated in [2] (see also [7]) provided that we set δ * = 0 in (1.1e) and we identify the symbols v and p f with v s and p t f used in [2]. In [2], equations (1.1) are summarized at the end of Section 2 as the outcome of derivation starting from the general principles of the theory of interacting continua, also using several well-motivated simplifying assumptions. In (1.1), v represents the velocity of the granular solid matrix, v f is the velocity of the interstitial fluid, p stands for the total pressure of the whole mixture and p f is the pressure of the fluid in a pore space.
The vector and scalar fields v, v f , p and p f represent the unknowns, the other quantities are given material functions/parameters. More precisely, φ = φ(p − p f ) is the porosity given as a function of the "effective" pressure p − p f , ̺ m s and ̺ m f are the constant material densities of the solid and the fluid, b represents given external forces, α is the drag coefficient, ν * the viscosity of the fluid, K is a constant coefficient related to permeability, δ * is the non-zero activation parameter and p s is the given lithostatic pressure. Note that v f appears only in (1.1d) and can be always obtained a posteriori once v, p f and p are obtained from (1.1a)-(1.1c). Consequently, in what follows, we consider system (1.1) without equation (1.1d). It is worth observing that the constitutive relation (1.1e) can be rewritten in a more compact way as an implicit constitutive relation (see Figure 1): We will exploit formulation (1.2) in our analysis. A systematic study of implicit constitutive equations go back to the original works [9] and [10]. We study the system of PDEs (1.1a)-(1.1c) and (1.2) in time-space cylinder (0, T ) × Ω, where T > 0 and Ω ⊂ R 3 is a bounded open connected set with Lipschitz boundary ∂Ω. We complete the system by considering the following boundary and initial conditions: v · n = 0 and ∇p f · n = 0 on (0, T ) × ∂Ω, (1.3a) Here, we used the following notation: n : ∂Ω → R 3 stands for the unit outer normal vector, Figure 1. Representation of the material response described by (1.2) while for any vector z defined on ∂Ω, z τ := z − (z · n)n denotes the tangential component of z, in particular, s := −(Sn) τ , and γ * , β * , s * are non-negative constants. Condition (1.3b) describes the shifted stickslip (or threshold slip) and it is analogous to that for the stress tensor in the bulk (see (1.2)). It includes as special cases, the stick-slip by taking β * = 0 while s * , γ * > 0, Navier's slip s = γ * v τ by taking s * , β * = 0 while γ * > 0, and perfect slip s = 0 by setting s * = γ * = 0. Note that the no-slip condition is obtained by letting either s * → +∞ or by setting β * = 0 and letting γ * → +∞.
The main purpose of this study is to establish long-time and largedata theory to the initial-and boundary-value problem described by (1.1a)-(1.1c), (1.2), (1.3), see Theorem 2.1 below. The novelties consist not only in incorporating a more general model with δ * ≥ 0, but more importantly in providing a different proof for more general class of data (particularly for b that is merely L 2 -integrable). More precisely, we can avoid using L ∞ -estimates for p f needed in [2]. Consequently, the main tool for taking the limit in the constitutive equations cannot be applied in the form given in [2,Proposition 5.3], but has to be modified in an essential way due to a lower integrability of p f , but also a more complicated material response. The novel key tool regarding the attainment of the constitutive equations by the limiting objects is proved separately in Proposition 3.1.
The key assumption of this proposition, namely (3.5) and (3.9), call for taking v n − v as a test function in the weak formulation of balance of linear momentum. However, v n − v is not admissible test function in the setting considered here. This difficulty can be overcome by using the L ∞ -truncation method which requires to introduce an integrable pressure, as the truncations (v n −v) ∞ are not divergenceless. Following the approach originally developed in [6] (see also [5]), we overcome such difficulty by considering slipping boundary conditions (1.3a)-(1.3b). As pointed out in [3], the analysis for unsteady flows changes remarkably when the no-slip condition is considered.
We use the L ∞ -truncation method in proving Theorem 2.1 below. While the truncations (v n − v) ∞ are difficult to make solenoidal, the authors of [4] succeeded to make the Lipschitz approximations (v n − v) 1,∞ divergenceless and they thus developed a solenoidal version of the Lipschitz truncation method. This tool allows one to avoid the presence of the pressure in the setting, therefore one may include more general responses as well as boundary conditions. As a matter of fact, we present new results available for systems describing materials that behave after activation |Dv| > δ * , as a power-law fluid, i.e. the constitutive equation (1.2) is replaced by (see Figure 2) The available results are presented in Theorem 2.2. We are not providing the proof of these results as they can be deduced from the approach used when proving Theorem 2.1 and from the methods used recently, for example, in [3]. Note that the latter results are restricted to models (1.4) with q > 6 5 (in three dimensions). Recently, another concept of dissipative solution was introduced in [1] and, its long-time and large-data existence is proved independently of what is the value of q (in particular also for q ∈ [1, 6/5]). In fact in the theory developed in [1] the stress tensor can be merely subdifferential of a convex potential depending on Dv, whose growth is at least linear. There are other approaches to analyze the mathematical properties of Bingham fluids (see e.g. [8] and [11]), but they are usually based on regularity techniques requiring smoother data.

Preliminaries and main results
For the sake of simplicity in the right-hand side of (1.1c), which has the form g := ∂ t p s − div b, we omit the effect of ∂ t p s as it plays the role of a given external force and it can be easily incorporated into the analysis. We also set without loss of any generality ̺ m s = ̺ m f = K = 2ν * = γ * = q * = 1, while we assume δ * , s * , β * ≥ 0. Finally, to shorten the notation we set Q := (0, T ) × Ω and Σ := (0, T ) × ∂Ω, where we fix Ω to be a bounded open set in R 3 with either Lipschitz or C 1,1 boundary ∂Ω; such sets are denoted either Ω ∈ C 0,1 or Ω ∈ C 1,1 .
Before stating the main results, let us summarize the notation. The symbol Dϕ stands for the symmetric part of the gradient of a vector-valued function ϕ, i.e.
The symbols (L q (Ω), · q ) and (W 1,q (Ω), · 1,q ) with q ∈ [1, ∞], stand respectively for the Lebesgue spaces, the Sobolev spaces with their own norms. If X is a Banach space of scalar functions, then X 3 and X 3×3 denote the space of vector-valued functions having three components and the space of tensor-valued functions respectively, with each component belonging to X. For a Banach space X, L q (0, T ; X) denotes a corresponding Bochner space. We make use of the following function spaces * are the dual spaces to W 1,q n and W 1,q n,div respectively. In particular, assuming Ω ∈ C 1,1 the following Helmholtz decomposition holds W 1,2 n = W 1,2 n,div ⊕ {∇ϕ; ϕ ∈ W 2,2 (Ω), ∇ϕ · n = 0 on ∂Ω}. Note that such decomposition is not valid for W 1,2 0 (Ω) 3 .
We are ready to enunciate the first result, which is the existence of weak solutions to system (1.1a)-(1.1c), (1.2), (1.3), proved in Section 5.
satisfying the following weak formulations: and the following constitutive equations:

Theorem 2.2.
Let Ω ∈ C 0,1 , T > 0, and q > 6 5 . Set m := max{2, q ′ } and r := max q, 5q 5q−6 . For satisfying the initial conditions (2.5), the following weak formulations: and the following constitutive equations: This result is stated without the proof here. The proof can be however achieved in the spirit of Theorem 2.1 by employing a solenoidal version of the Lipschitz-truncation method developed in [4], and by using the approximation scheme presented in [3, Theorem 3.3].

Attainment of the constitutive equations
In this section, we establish a new scheme how to take the limit in the constitutive equations needed when proving Theorem 2.1.  In addition, assume that {V n } +∞ n=1 is a sequence such that Proof. First, note that by virtue of (3.4) and the Lipschitz-continuity of τ , it follows that Now, as for any A ∈ L 2 (U ), A = O, it holds (3.12) integrating this inequality over U , subtracting and adding Z n and using (3.1), we get Taking limsup as n → ∞ and employing the facts that the first integral converges to zero due to (3.11), the term D n |D n |+ 1 n is uniformly bounded in L ∞ (U ) and the sequence D n − A is bounded in (3.14) Referring then to the convergences (3.2) and (3.3) and using also (3.5), we conclude that Now, for any δ > 0, ε ∈ (0, δ) and for arbitrary matrices C and B 1 bounded in L 2 (U ) 3×3 and satisfying |C| ≤ 1 and B 1 = O, consider Note that such A's are non-zero in U . Inserting them into (3.15) we obtain Letting first ε → 0 in (3.16), we observe that for any B 1 = O. Consider, for any a > 0 and ω ⊂ U , the matrix B 1 of the form It then follows from (3.16) that with C positive constant, which implies letting a → 0 Since ω is arbitrary, we conclude that |Z| ≤ τ (p f ) on the set {|D| = 0}. Finally, letting δ → 0, we get, for arbitrary C, This implies The latter and (3.18) are equivalent to (3.6).
It remains to prove (3.10), which however follows from standard Minty's argument. Indeed, by the monotonicity, we have lim sup for any A ∈ L 2 (U ) 3×3 . By virtue of (3.9) and of convergences (3.8) and (3.3) we get Choosing A := D ± εC, with arbitrary C ∈ L 2 (U ) 3×3 and ε > 0, and after the limit as ε → 0 we obtain for any C, which implies (3.10).
Note that here we provided a proof of (3.6), which is simplified and shorter than the one given in [2, Proposition 5.3].

Approximations
In this section, we prepare all the needed tools in order to prove Theorem 2.1. For any n ∈ N, we introduce the following approximating system div v = 0 in Q, where G n : R → R is a smooth function such that G n (u) = 1 if |u| ≤ n, G n (u) = 0 if |u| ≥ 2n and |G ′ n | ≤ 2 n . Next, we consider the following regularization of the constitutive equations (both in the bulk and on the boundary) see Lemma B.1 in [3]. Therefore, due to the presence of the truncation in the convective term and the introduced approximations in the constitutive equations, the existence of weak solutions to system (4.1)-(4.5) can be proved through standard techniques of monotone operators, following also the spirit of the proof in [2, Proposition 5.1]. We enunciate the relevant result below and for the reader's convenience the proof can be found in Appendix.

Proof of Theorem 2.1
The proof is split in the following steps.
Step 1. Approximations. From Proposition 4.1 and following the reconstruction of the pressure in [2, Theorem 4.1, Step 2 of the proof], we get for each n ∈ N the existence of (v n , p n f , p n , S n , s n ), with p n ∈ L 2 (Q), satisfying with S n , s n fulfilling (4.13), (4.14) respectively, and satisfying also (4.15).
Step 2. Uniform estimates with respect to n and limit as n → +∞. Setting w := v n in (5.1) and z := p n f in ( v n → v strongly in L q (Q) 3 for all q ∈ 1, 10 3 , it follows from (5.12) that Finally, with the obtained convergences it is standard to prove that : ∇z for all z ∈ L 4 (0, T ; W 1,2 (Ω)). (5.16) Step 3. Attainment of the constitutive equations on the boundary. Using that s n ⇀ s weakly in L Thus, a suitable adjustment of Proposition 3.1 implies that (2.4) is fulfilled.
Step 4. Attainment of the constitutive equations in the bulk. In order to employ Proposition 3.1 we need to prove the limsup property (3.5), but as the solution itself can not be used as test function in (5.15), we follow the strategy as in [2] and perform the L ∞ -truncation method. To this aim, we introduce where λ n ∈ [A, B] with 0 < A < B < ∞ will be suitably chosen numbers independent of n, but depending on parameter N tending to +∞, see details below. For the reader's convenience, we recall all the properties of w n below, w n → 0 strongly in L s (Q) for every s ∈ [1, +∞), (5.17) w n → 0 strongly in L 2 (Σ), (5.18) w n ⇀ 0 weakly in L 2 (0, T ; W 1,2 n ), if |v n − v| > λ n . Now, let us define where I n := |p n 1 | 2 + |Z n | 2 + |Z| 2 + |V n | 2 + |V| 2 + |∇v n | 2 + |∇v| 2 . Note that it holds (see formula (6. Finally, all assumptions of Convergence Lemma 4.1 are fulfilled, thus (1.2) holds a.e. in U . Since |Q \ U | ≤ ε we can let ε → 0 and obtain that (1.2) holds a.e. in Q. Theorem 2.1 is proved.
Step 2. Uniform estimates. Multiplying (6.2) by c m r (t) and (6.5) by d m r (t) and taking the sum over r = 1, . . . , m, we obtain Adding {|Dv m |≤δ * } |Dv m | 2 to both sides of (6.6) and then by Young's inequality, we get Step 3. Limit. By virtue of uniform estimates established above there exist subsequences of {v m }, {p m f }, {S m } and {s m }, converging respectively weakly (or *-weakly) to v, p f , S and s in the corresponding function spaces. Furthermore, Aubin-Lions compactness Lemma and its variant including the trace theorem imply the following strong convergences: v m → v a.e. in Q and strongly in L q (Q) 3 for any q ∈ 1, 10 3 , (6.18) p m f → p f a.e. in Q and strongly in L q (Q) 3 for any q ∈ 1, 10 3 , (6.19) v m τ → v τ a.e. in Σ and strongly in L q (Σ) 3 for any q ∈ 1,