On the role of compressibility in poroviscoelastic models

1.   Introduction

Poroviscoelastic models of multi-component mixtures are often utilized in biological applications to describe the flow of fluids within the pores of a deformable solid skeleton, see e.g., ^[1,2,3,4]. Skeleton viscoelasticity is often due to the complex structures including extracellular matrix, collagen and elastin that are present in biological tissues. Specifically, our work focuses on the bio-fluid-mechanical response of poroviscoelastic media to non-smooth data, since this aspect is crucial in understanding the mechanisms leading to tissue damage in the optic nerve head, and consequent vision loss, associated with glaucoma ^[5,6,7,8].

In the absence of viscoelasticity, we have recently shown that time irregularities in the volumetric and/or boundary sources of linear momentum lead to a blow-up in the solution of poroelastic models ^[1,4]. Interestingly, the blow-up can be prevented by including structural viscoelasticity ^[4]. From the application viewpoint, examples of time-irregularities in the data are discontinuities in intraocular pressure, which acts as a boundary traction for the optic nerve head tissue, or discontinuities in the gravitational force, which acts as a volumetric source of linear momentum. The intraocular pressure exhibits rapid changes every time we rub our eyes or we change posture ^[9], whereas rapid changes in the gravitational acceleration are experienced by jet pilots and astronauts during flights ^[10,11]. Since tissue viscoelasticity has been shown to decrease with age and/or disease conditions, the solution blow-up identified by our theoretical work led to hypothesize that rapid changes in intraocular pressure and gravitational acceleration, even if within physiological ranges, could damage the optic nerve head tissue if its viscoelasticity was pathologically reduced.

It is important to emphasize that our previous work was built on the assumption that the poroviscoelastic medium under consideration was made of incompressible components. The incompressibility assumption is quite common in biological applications, since tissues are mostly made of water. However, compressibility is always present in real tissues, and this leads to wonder whether and to what extent compressibility of the mixture components would affect the tissue response to non-smooth data. The present work aims at addressing this question.

Proceeding as in ^[4], we assume: (i) a one-dimensional (1D) geometry; (ii) negligible inertial terms in the linear momentum balance equation; (iii) a Kelvin-Voigt model for the effective stress tensor and (iv) a constant hydraulic permeability of the porous medium. Then, we express the fluid pressure and the solid stress as functions of the sole solid phase displacement, and we obtain an initial-boundary value problem (IBVP) of pseudo-parabolic type for the solid phase displacement. For this equivalent formulation we are able to construct an analytical solution and prove the well-posedness of weak solutions. Moreover, we recover analytical formulas for fluid pressure and discharge velocity, and discuss their regularity in terms of the regularity of the data. Finally, we analyze the behavior of all the solutions for various continuous and discontinuous boundary loads, which are of particular interest in understanding how changes in intraocular pressure would impact the bio-fluid-mechanics of ocular tissues.

The main conclusion of our analysis is that the compressibility of the components of a poroviscoelastic mixture does not give rise to unbounded responses to non-smooth traction data. Interestingly, compressibility provides an additional mechanism that limits the maximum of the discharge velocity when the imposed boundary traction is irregular in time. This mechanism originates from the capability of the system to store potential energy as its components are elastically compressed, thereby delaying the transmission of irregularities in the linear momentum from the solid to the fluid. As a result, the fluid has the time to accommodate for sudden changes, resulting in bounded velocities.

Our results fit well with other poroviscoelastic studies motivated by geomechanical applications, where viscoelasticity was found to play a crucial role in the response of the medium to impulsive loads. Specifically, the studies focused on evaluating the consequences that different choices for the viscoelastic models would have on the medium response to external loads. For example, in ^[12], Schanz and Cheng considered a Kelvin-Voigt viscoelastic model and investigated the consequences of adopting it for the bulk compression modulus, the shear modulus and the compression modulus of the solid material. In ^[2,3], Huang et al. utilized the quasi-linear viscoelastic theory to study the response of articular cartilage under compression and tension experiments. However, these studies did not consider how the mechanical responses of the poroviscoelastic mixture would change depending on the level of compressibility of the mixture components, which is the main focus of the present work.

The outline of this article is as follows. In Section 2 we describe the compressible poro-visco-elastic model under consideration and discuss the finite compressibility of the components of the mixture. Section 3 introduces the energy identity associated with the system. In Section 4 we present an equivalent form of the fluid-solid mixture system, written solely in terms of the solid phase displacement. To simplify the theoretical analysis, in Section 5 we reduce the poroviscoelastic model into dimensionless form. Section 6 studies the well-posedness and regularity of solution for the IBVP introduced in Section 4, and provides analytical formulas for the elastic displacement, fluid pressure and discharge velocity. In Section 7 we compute and analyze the behavior of the solid displacements, the fluid pressures and the discharge velocities associated with different continuous and discontinuous boundary sources. We also display the energies, as well as the dissipation and source terms associated with the various cases. We conclude the article with Section 9, where we discuss our results and draw our final conclusions.

2.   The compressible poroviscoelastic model

In this paper, we focus on a viscoelastic, compressible Biot model in one spatial dimension. Let $x$ and $t$ denote the spatial and temporal coordinates, respectively. In the case where inertial terms are negligible, displacements are infinitesimal and external sources of linear momentum and mass are absent, the balance equations to be solved in the spatial interval $\left(0, L\right)$ and in the temporal observational interval $(0, T]$ can be written as:

where $\sigma$ is the total stress, $\zeta$ is the fluid content and $v$ is the discharge velocity. Equations (2.1) are complemented with the following constitutive laws:

Equation (2.2a) expresses the fact that the total stress is the sum of a contribution due to the interstitial fluid pressure (or pore pressure) $p$ and one due to the effective stress $\sigma _{0}$, which is assumed to be characterized by viscoelastic behavior of Kelvin-Voigt type (see Eq (2.2b) where $u$ denotes the solid phase displacement). Equation (2.2c) expresses the fact that the fluid content is altered by changes in fluid pressure and solid deformation. Equation (2.2d) is Darcy's law relating the discharge velocity $v$ with the pressure gradient by means of the hydraulic permeability $\mathcal{K}$. The Biot coefficient $\alpha$, the storage coefficient $c_{0}$, the material parameters $\theta$ and $\eta$ and the permeability $\mathcal{K}$ are assumed to be given positive constants.

In the classical Biot's theory, the material parameter $\theta$ can be expressed as $\theta = K - (4/3)\, G$, where $K$ is the drained bulk compression modulus and $G$ is the shear modulus. In the following, we will simply refer to $\theta$ and $\eta$ as the elastic and viscoelastic parameters, respectively ^[12]. In addition, we will use the notation $\mathcal{K} = \mathcal{K}_{0}$ to emphasize the fact that the permeability is assumed to be a given constant.

The problem must be equipped with appropriate boundary conditions. In the following, we will assume that the boundary located at $x = 0$ is fixed and impermeable, namely:

and that the boundary located at $x = L$ experiences an external stress $P$ (a compressive stress if $P$ is positive, a tensile stress if $P$ is negative) that is supported entirely by the solid component of the mixture (condition of exposed pores ^[13]), namely:

Finally, we complete the formulation of the problem by prescribing the following initial conditions:

We notice that (2.4a) and (2.4b) also imply the following initial conditions for the dilation of the solid material and the discharge velocity, respectively:

Remark 1. The case of incompressible mixture components can be obtained by setting $c_0 = 0$ and $\alpha = 1$ in the model described above, as detailed in ^[14,15]. The study of analytical solutions for the incompressible model was addressed in ^[1,4].

3.   Energy identity

The mathematical system described in Section 2 satisfies an energy identity that helps provide a physical interpretation of the solutions, assuming they exist. To this end, let us multiply (2.1a) by $\partial u/\partial t\in L^{2}(0, L)$; integrating over $(0, L)$ and using the constitutive laws (2.2a) and (2.2b) and the boundary conditions (2.3a), (2.3c) and (2.3d), we obtain

Multiplying (2.1b) by $p\in L^{2}(0, L)$, integrating over $(0, L)$ and using the constitutive laws (2.2c) and (2.2d) and the boundary conditions (2.3b) and (2.3c), we obtain

Subtracting (3.6) from (3.7) we get the following evolution equation for the total energy stored in the poroviscoelastic system

where the energy functional $\mathcal{E}_{tot} = \mathcal{E}_{tot}(t)$, the dissipation functional $\mathcal{D} = \mathcal{D}(t)$ and the force term $\mathcal{F} = \mathcal{F}(t)$ are defined as:

Remark 2. The physical units of $\mathcal E_{tot}$ are Joules per unit area, namely J m$^{-2}$. This is due to the fact that we are considering a one-dimensional problem in space and, as a consequence, all the problem variables are assumed to be constant on every plane perpendicular to the chosen direction $x$. Mathematically, $\mathcal E_{tot}$ is obtained by integrating in $x$ between 0 and L the energy density $\varepsilon_{tot}$ defined as

The units of $\varepsilon_{tot}$ are J m$^{-3}$. Analogously, the units of $\mathcal D$ and $\mathcal F$ are J m$^{-2}$ s$^{-1}$.

Since $\mathcal{E}_{tot}\geq 0$ and $\mathcal{D}\geq 0$, in absence of forcing terms (i.e., $P = 0$) the energy decreases in time. It is important to emphasize that the terms proportional to the storage coefficient $c_{0}$ and the elastic parameter $\theta$ contribute to the total energy in the form of potential energy, so that we can write

Conversely, the terms proportional to $\mathcal{K}_{0}$ and $\eta$ contribute to dissipate energy via viscous effects within the fluid and solid components. We notice that $\mathcal{F}$ does not have a definite sign since it depends on the boundary terms. The energy identity (3.8) will be very useful in interpreting the results presented in Section 7.

4.   The 1D poroviscoelastic model in displacement form

Now we express problem (2.1)–(2.4) solely in terms of the solid displacement. Combining (2.1a) and (2.3d), we obtain that the total stress is given by

and, moreover, the fluid pressure $p$ and the discharge velocity $v$ can be written in terms of the solid displacement $u$ as:

where $\sigma _{0} = \sigma _{0}\left(u\left(x, t\right) \right)$ is given by (2.2b). Let us now derive the problem satisfied by $u$. Integrating Eq (2.1b), where $\zeta$ is given by (2.2c), over the space interval $\left[0, x\right]$ and taking the boundary conditions (2.3) into account, yields

Now we substitute $p$ and $v$ by expressions (4.10) and note that

by virtue of (2.3a). Then we obtain

The boundary conditions are (2.3a) and $\left. \sigma _{0}\right\vert _{x = L} = -P\left(t\right)$ (coming from (2.3c) and (4.10a)). The initial conditions are (2.4a) and, on using (2.4b), (4.10a), (2.2b) and (2.5a),

or equivalently

$u_{1}$ being an arbitrary constant. Note that $u_{1} = 0$ if and only if the compatibility condition with (2.3a)

is satisfied. Summing up:

Find $u = u\left(x, t\right)$ such that:

Remark 3. The solution $u\left(x, t\right)$ of the problem (4.12) is the sum $u_{1}\left(x, t\right) +u_{0}\left(x, t\right)$ of the solution $u_{1}\left(x, t\right)$ of (4.12) where $P = 0$, and of the solution $u_{0}\left(x, t\right)$ of (4.12) where $u_{1} = 0$. In the first case, a perturbation is generated at time zero, with a certain initial velocity of propagation ($u_{1}\neq 0$): then $u_{1}\left(x, t\right)$ measures how the solid displacement changes through the medium when no disturbance is created at the boundary ($P = 0$). On the other hand, if there is no ''initial impulse" (i.e., $u_{1} = 0$) and a disturbance is generated at the boundary ($P\neq 0$), then the corresponding change of the solid displacement is represented by $u_{0}\left(x, t\right)$. Therefore, the compatibility condition (4.11) simply means to consider the response of the system to the sole external stress $P$.

Remark 4. We assume throughout the article that $c_{0} > 0$ and $\eta > 0$, i.e. the system is characterized by finite compressibility and structural viscoelasticity.

Remark 5. The IBVP (4.12) has a markedly different character compared to the linear poroviscoelastic system studied in ^[4] because of the second-order time derivative on the left-hand side of (4.12a).

Remark 6. If we define the following quantities:

then the partial differential equation (4.12a) that describes the dynamics of the solid phase displacement can be written as

Such an equation can be formally interpreted as a linear momentum balance equation for a single phase solid material whose dynamics is equivalent to that of the fluid-solid mixture under the assumptions of Section 2. In particular, we see that the finite compressibility of the mixture components, corresponding to $c_{0} > 0$, provides:

(i) an inertia-like term for the equivalent solid, even though inertial terms were neglected for the original solid phase within the mixture;

(ii) an additional term in the stress tensor introducing nonlocal effects in space;

(iii) a volumetric forcing term that results from the load applied as a boundary condition.

5.   The linear 1D model in dimensionless form

In order to simplify the theoretical analysis, in this section we reduce the 1D model (4.12) into dimensionless form. With this aim, for any variable $Y$, we define the corresponding non-dimensional variable by

where $[Y]$ is a suitably chosen scaling factor that has the same units as $Y$. The selection of the scaling factor is not unique and, in general, not trivial. In this article we generalize the choice made in ^[4], by introducing the following scaling factors:

We also define the non-dimensional quantity

Note that $0 < \gamma < 1$. We recover the same definitions of the scaling factors as in ^[4] by setting $\alpha = 1$ and $c_{0} = 0$ in (5.1).

Remark 7. For notational simplicity we will drop the '$\widehat{\cdot }$' notation for the dimensionless variables and instead use the same symbols adopted for the dimensional variables.

The linear 1D model for a poroviscoelastic mixture in dimensionless form reads:

Once $u\left(x, t\right)$ is known, the functions $\sigma$, $p$ and $v$ can be computed as follows:

Lastly, the dimensionless energy equation is written again in the form (3.8) where

6.   Well-posedness and Regularity of Solution

We make the following assumption on the boundary traction in (5.3).

Assumption 8. $P\left(t\right)$ is a piecewise smooth function on $\left[0, T\right]$.

We recall that a function $P\left(t\right)$ is piecewise smooth on $\left[0, T\right]$ if both $P$ and its derivative $P^{\prime }$ are continuous on $[0, T]$, except possibly at a finite number of points in $\left(0, T\right)$, where they have simple jump discontinuities.

6.1.   Auxiliary Problem

Define

Using Assumption 8 we have that $U\left(t\right)$ is absolutely continuous on $\left[0, T\right]$ and

Now we introduce the change of variable

and note that $w$ solves the following auxiliary IBVP with homogeneous boundary data:

where the volumetric source and initial datum are given by:

We shall prove the existence and uniqueness of the solution $w\left(x, t\right)$ for a very general class of data $f\left(x, t\right)$ and $\varphi \left(x\right)$. For sake of exposition we write $H = L^{2}\left(0, 1\right)$, and define the real Hilbert space

endowed with the equivalent norm $\left\Vert v\right\Vert _{V} = \left\Vert \partial v/\partial x\right\Vert _{L^{2}\left(0, 1\right) }$, due to Poincaré's inequality.

Remark 9. Sobolev's Embedding Theorem gives $W^{1, 2}\left(0, 1\right) \subset \mathcal{ C}^{0}\left[0, 1\right]$ so that $v\left(0\right) = 0$ holds in a strong sense for every $v\in V$.

Now we make the following assumptions on the functions $f\left(x, t\right)$ and $\varphi \left(x\right)$:

and write $w\left(t\right) = w\left(\cdot, t\right)$, $w^{\prime }\left(t\right) = \partial w\left(\cdot, t\right) /\partial t$, etc. We then define weak solutions of (6.9) as follows.

Definition 10. [Weak solution of problem (6.9)] A function $w:\left[0, T\right] \rightarrow V$ is a weak solution of the auxiliary problem (6.9) if:

D1 $w\in W^{1, 2}\left(0, T; V\right)$ and $w^{\prime \prime }\in L^{2}\left(0, T; V^{\prime }\right)$;

D2 for every $v\in V$ and for $t$ pointwise a.e. in $\left(0, T\right)$

or, equivalently,

D3 $w\left(0\right) = 0$ and $w^{\prime }\left(0\right) = \varphi.$

Remark 11. Condition [D1] implies that $w\in C^{0}\left(\left[0, T\right]; V\right)$ and $w^{\prime }\in C^{0}\left(\left[0, T\right]; H\right)$, and thus condition [D3] is well defined. The Dirichlet boundary condition in (6.9) at $x = 0$ is included in the regularity requirement that $w\left(t\right) \in V$, whereas the boundary condition at $x = 1$ is satisfied in a weak sense.

6.2.   A-priori Estimates

Lemma 12. Let $w$ be a weak solution of (6.9).Then there are constants $C_{i}$'s, depending only on $\gamma$, $\eta$ and $T$, such that the following estimates hold for $t$ pointwise in $\left[0, T \right]$:

Proof. Using $\left(\gamma \eta w^{\prime }+w\right) \in L^{2}\left(0, T; V\right)$ as multiplier in (6.14), we get

Integrating in time over $(0, t)$ and using the given initial conditions, we obtain

Using Cauchy-Schwarz and Young's Inequalities for the last term on the right-hand side, we obtain the following estimate:

Now, dropping the second and third terms on the left-hand side of (6.16) and using Gronwall's Inequality, estimate (6.15a) follows. Similarly, by dropping the first and third terms on the left-hand side and using (6.15a), we get (6.15b).

Lastly, we use triangle inequality and the embedding $V\hookrightarrow H$ to write

Then, estimate (6.15c) follows from (6.15a) and (6.15b).

Lemma 13. Let $w$ be a weak solution of (6.9).Then there is a constant $C$, depending only on $\gamma$, $\eta$ and $T$, such that

Proof. From (6.15b) it immediately follows

In addition, Eqs (6.16) and (6.15a) give

so that

for a suitable $C = C\left(\gamma, \eta, T\right)$.

Lastly, using the definition of a weak solution, we have that for $v\in V$

This implies that

from which, using estimates (6.15b) and (6.15a), we obtain

for a suitable $C = C\left(\gamma, \eta, T\right)$.

The following corollary is an immediate consequence of Lemma 13.

Corollary 14. (Uniqueness and continuous dependence on data) The weak solution to problem (6.9) is unique and depends continuously on the data.

6.3.   Existence of Solution

The IBVP (6.9) can be solved formally using separation of variables. If we look for solutions of the form $w(x, t) = T(t)X(x)$, then the associated regular Sturm-Liouville Problem is

with eigenvalues and eigenfunctions given by:

Remark 15. The sequence of functions $\left\{ \sqrt{2}X_{n}\left(x\right) \right\}$ forms a Hilbert space basis for $H$, whereas the sequence of functions $\left\{ \sqrt{\frac{2}{\lambda _{n}}}X_{n}\left(x\right) \right\}$ forms a Hilbert space basis for $V$ (see ^[17]).

The solution $w$ of (6.9) has the expansion

where $T_{n}(t)$ can be recovered using the data. Similarly, we use the basis $\left\{ X_{n}\left(x\right) \right\}$ to represent $\varphi$ and $f$ as follows:

where $\varphi _{n}$ and $f_{n}(t)$ are the Fourier coefficients of $\varphi$ and $f(\cdot, t)$ with respect to $X_{n}(x)$, respectively. Parseval's identity (consequence of Remark 15) provides the following relations:

Note that $T_{n}(t)$ satisfies the following ordinary differential equation (ODE) for all $n\geq 0$:

and initial conditions

The characteristic equation for the homogeneous counterpart of (6.27) is given by

Since the discriminant of the equation is

then for each $n\geq 0$ the characteristic equation has two real, negative, distinct roots $r_{1n} = -\Lambda _{1n}$ and $r_{2n} = -\Lambda _{2n}$, where:

Remark 16. $\Lambda _{1n}$ and $\Lambda _{2n}$ satisfy the following relations:

Therefore the solution for the homogeneous ODE (6.27) is given by $T_{n}^{0}(t) = a_{n}e^{-\Lambda _{1n}t}+b_{n}e^{-\Lambda _{2n}t}$. The particular solution is obtained from variation of parameters formula. The Wronskian is given by $W(t) = (\Lambda _{2n}-\Lambda _{1n})e^{-\Lambda _{1n}t}e^{-\Lambda _{2n}t}$, so that the particular solution has the following form

We introduce the following notation

and then write the solution to (6.27) as

where we used the following formula for convolution

Now we use (6.32) back into (6.22) and impose the initial conditions. We have:

which yields $a_{n} = -b_{n} = \varphi _{n}/\left(\Lambda _{2n}-\Lambda _{1n}\right)$. In conclusion, we obtain:

We note that the terms $G_{n}(t)$ given in (6.31) satisfy the following estimates.

Lemma 17. There exists a constant $C > 0$ such that for all $n\geq 0$ and for all $t\in \left[0, T\right]$ we have:

Proof. Since $0 < \Lambda _{1n} < \Lambda _{2n}$, then $0 < \exp \left(-\Lambda _{2n}t\right) \leq \exp \left(-\Lambda _{1n}t\right) \leq 1$ and thus for all $t\in \left[0, T\right]$ we have

and we get (6.34a) since the sequence $\frac{\lambda _{n}}{ \Lambda _{2n}-\Lambda _{1n}} = \frac{\gamma \eta \lambda _{n}}{\sqrt{(\eta \lambda _{n}+1)^{2}-4\gamma \eta \lambda _{n}}}\rightarrow \gamma$ as $n\rightarrow \infty$. Moreover, we have that for all $t\in \left[0, T\right]$

and (6.34b) follows since $\frac{\Lambda _{2n}}{\Lambda _{2n}-\Lambda _{1n}}\rightarrow 1$ as $n\rightarrow \infty$.

Now we can state and prove our well-posedness result.

Theorem 18. (Well-posedness of problem (6.9)) For every $f\in L^{2}(0, T; H)$ and $\varphi \in H$ there is a unique weak solution of (6.9), in the sense of Definition 10. Moreover, the solution depends continuously on the data.

Proof. Uniqueness and continuous dependence of weak solution have already been proved, see Corollary 14, so that it remains to prove existence, i.e. that $w(x, t)$ given in (6.33b) satisfies conditions (D1)–(D3) of Definition 10.

(D1) (A) First we show that $w\in L^{2}\left(0, T; V\right)$. For all $t\in \lbrack 0, T]$, we have

Using Lemma 17, we obtain^*

^*In what follows, for notational convenience, $C$ will denote possibly different constants depending only on $\gamma, \eta, T$.

Therefore $w\in L^{\infty }\left(0, T; V\right) \subset L^{2}\left(0, T; V\right)$.

(B) The weak derivative of $w$ with respect to time is given by

and we obtain the following estimate

Again by Lemma 17, we obtain

Thus $w^{\prime }\in L^{\infty }\left(0, T; H\right) \subset L^{2}\left(0, T; H\right)$. In order to show that $w^{\prime }\in L^{2}\left(0, T; V\right)$, we first use the Monotone Convergence Theorem to write

We use multiplier $T_{n}^{\prime }$ in (6.27) and the initial conditions (6.28) and obtain the following estimate:

Now we use (6.38) back into (6.37), and take advantage of estimate (6.36) to obtain

and this gives the desired conclusion that $w^{\prime }\in L^{2}\left(0, T; V\right)$.

(C) Left to show that $w^{\prime \prime }\in L^{2}\left(0, T; V^{\prime }\right)$. For any $v\left(x\right) = \sum_{n = 0}^{\infty }c_{n}X_{n}\left(x\right) \in V$ and for $t\in \lbrack 0, T]$, we use (6.27) to define the linear functional $w^{\prime \prime }\left(t\right)$ as follows

Hence, by virtue of the embedding $V\hookrightarrow H$, we obtain that

This shows that $w^{\prime \prime }(t)\in V^{\prime }$ and

or, equivalently,

In conclusion, using the estimates in part (A) and part (B), we obtain that $w^{\prime \prime }\in L^{2}\left(0, T; V^{\prime }\right)$.

(D2) Now we show that $w$ satisfies condition (D2). Since the $\left\{ \sqrt{\frac{2}{\lambda _{n}}}X_{n}\right\}$ is a basis in $V$, it suffices to consider the test function $v = X_{n}$. For $t$ pointwise almost everywhere in $\left(0, T\right)$, we use (6.27) and obtain

(D3) As stated in Remark 11, condition (D1) implies that $w\in C^{0}\left(\left[0, T\right]; V\right)$ and $w^{\prime }\in C^{0}\left(\left[0, T\right]; H\right)$. Moreover, solution $w(x, t)$ satisfies the initial conditions (D3) of Definition 10 by virtue of (6.28). This concludes the proof of our well-posedness theorem.

6.4.   Regularity of the Solution

We have established the existence and uniqueness of the weak solution $w$ to the auxiliary problem (6.9). Now we examine its regularity.

Proposition 19. (Regularity of $w$ and $w_{xx}$) Under the hypotheses in Theorem 18, the following is true:

Proof. (A) The result follows from Lemma 17, the fact that $|X_{n}(x)|\leq 1$ and Weierstrass Test. Indeed, for all $x\in \lbrack 0, 1]$ and $t\in \lbrack 0, T]$, and every $n\geq 0$, we have

Now, due to Lemma 17, we have

and

Applying Weierstrass Test, we obtain that the series $\sum_{n = 0}^{\infty }\ T_{n}(t)X_{n}\left(x\right)$ converges absolutely and uniformly in $Q = [0, 1]\times \lbrack 0, T]$. Moreover, we note that $G_{n}\left(t\right) \varphi _{n}X_{n}(x)+\frac{1}{\gamma \eta }(G_{n}\ast f_{n})\left(t\right) X_{n}(x)\in C^{0}(Q)$, for every $n\geq 0$, and thus $w(x, t)\in C^{0}(Q)$.

(B) The second order weak derivative in space of $w$ is given by

Then we use estimate (6.34a) in Lemma 17 and obtain

so that $w_{xx}\in L^{\infty }\left(0, T; H\right)$.

If the data are more regular with respect to $x$ the weak solution $w(x, t)$ enjoys stronger regularity properties.

Proposition 20. (Regularity of $w_{t}$, $w_{xx}$ and $w_{txx}$) In addition to the hypotheses in Theorem 18, we make the following ones:

Then the following is true:

Proof. (A) For all $x\in \lbrack 0, 1]$ and $t\in \lbrack 0, T]$, and every $n\geq 0$, we have

Now, due to Lemma 17, we have

and

Applying Weierstrass Test, we obtain that the series $\sum_{n = 0}^{\infty }\ T_{n}^{\prime }(t)X_{n}\left(x\right)$ converges absolutely and uniformly in $Q = [0, 1]\times \lbrack 0, T]$ hence $w_{t}(x, t)\in C^{0}(Q)$.

(B) Like in the proof of (6.40a), we have

so that $w_{xx}\in L^{\infty }\left(0, T; V\right)$.

(C) Since the second order weak derivative in space of $w^{\prime }$ is given by

then

In order to estimate the right hand side, we use multiplier $\lambda _{n}T_{n}^{\prime }$ on the ODE (6.27) that $T_{n}$ solves to obtain

Now we integrate from 0 to $t$ and use the initial conditions (6.28), to get

Thus, using (6.43) back into (6.42), we obtain

and the assertion follows from estimate (6.17).

6.5.   Analytical formulas for solid displacement and discharge velocity

Now we are in a position to return to our original linear 1D problem (5.3). The solid displacement solution $u\left(x, t\right)$ of (5.3) is the sum

where we recall that

and $w\left(x, t\right)$ is the unique weak solution of the auxiliary problem (6.9) with special data (6.10). From the Fourier expansions

as well as the Fourier series (6.23) of $\varphi$ and $f$ with respect to the basis $\left\{ X_{n}\right\}$, their Fourier coefficients are given by:

To summarize, we can write the solution $u$ as a sum

where:

Remark 21. The solid displacement solution $u\left(x, t\right) \in C^{0}\left(Q\right)$ since the physical data (6.10) satisfy $f\in$ $L^{2}\left(0, T; V\right)$ and $\varphi \in H$. Notice, however, that our special $\varphi$ belongs to $V$ if and only if $u_{1} = 0$, i.e. $\varphi = 0$ and $u_{1}\left(x, t\right) \equiv 0$.

By virtue of (6.46), the discharge velocity (5.4c) is given by

where:

Remark 22. It should be stressed that the Fourier coefficients of $v_{1}\left(x, t\right)$ are $O\left(n\right)$, hence the component of the discharge velocity due to the presence of an initial impulse $u_{1}$ always exibits a blow-up whatever boundary traction $P\left(t\right)$ is considered.

7.   Continuous vs. discontinuous boundary sources

In this section, we analyze the behavior of the solutions obtained for model (5.3) in the case where the compatibility conditions are satisfied (i.e., $u_1 = 0$), and the boundary traction $P(t)$ is characterized by continuous or discontinuous waveforms. To graphically represent the model solutions we proceed as follows: (1) we define the numerical values of model parameters (in dimensional form) consistently with the experimental data illustrated in ^[16]; (2) we perform the non-dimensionalization of the model equations according to the procedure described in Section 5; (3) we compute the model solutions in non-dimensional form; (4) we multiply the non-dimensional input data and model solutions by the scaling factors introduced in (5.1) and plot the obtained results for subsequent analysis. For each considered waveform of $P(t)$, we provide the non-dimensional expressions of the solutions $u$, $p$ and $v$. These latter expressions allow us to compute the non-dimensional energies $\mathcal E_\theta$, $\mathcal E_{c_0}$ and $\mathcal E_{tot}$, the non-dimensional dissipation functional $\mathcal{D}$ and the non-dimensional force term $\mathcal F$ using the expressions (5.5). The resulting non-dimensional energies are then multiplied by the scaling factor $[\mathcal{E}_{tot}]$ introduced in (5.1j) and the same is done for the non-dimensional dissipation and force terms that are multiplied by the scaling factor $[\mathcal{D}]$ introduced in (5.1i).

Table 1 reports the numerical values of model parameters (in dimensional form) that are used in the next sections. The values for $L$, $T$, $\theta$, $K_0$ and $P_{\text{ref}}$ have been chosen consistently with the experimental data illustrated in ^[16]. To observe the asymptotic behavior of the modeled system, $T$ has been taken equal to $10^5 \, \mathrm{{s}}$ in the example illustrated in Section 7.1. In this example, the corresponding value of $\eta$ is $4.85 \cdot 10^9 \, \mathrm{{N \, s \, m^{-2}}}$. The value of the Biot coefficient $\alpha$ has been set equal to 0.9, whereas in ^[16] it was equal to 1 since both mixture components were assumed to be incompressible. The values of the viscoelastic parameter $\eta$ and of the compressibility parameter $c_0$ (which were not present in the model studied in ^[16]) have been set equal to $\theta \tau_e$ and $1/P_{\text{ref}}$, respectively, where $\tau_e$ is an elastic time constant that has been set equal to $T/20$.

The numerical values of the parameters (in dimensionless form) that are used in the computations discussed in the next sections are reported in Table 2. These values have been obtained by applying the scaling procedure described in Section 5 to the values in Table 1. With a slight abuse of notation, the symbols used to denote the dimensionless parameters are the same ones that we used for the corresponding parameters expressed with their physical units.

For convenience, we recall here the formulas that we use to compute the solid displacement, the fluid pressure and the discharge velocity in dimensionless form:

where

Recalling formula (6.31) and using the fact that $(G_n \ast U^{\prime})^{\prime } = G_n^{\prime }\ast U^{\prime }$, we obtain the following simplification for the coefficients $B_{n}\left(t\right)$

7.1.   The case of step pulse at $t=t^*$

Let $t^* \in [0, T)$. We consider the following boundary source

Figure 1 gives a graphical representation of $P(t)$ for $t^* \in (0, T)$.

Replacing (7.2) into (6.6) we obtain

We can now compute the convolution needed in (7.1a)

We observe that for $t \geq t^*$,

and therefore the coefficients (7.1d) present in the expansions of pressure and discharge velocity become

Then displacement, pressure and discharge velocity are all zero for $t < t^*$, and have the following representations for $t \geq t^*$:

Note that all three series in the expressions (7.5) converge absolutely and uniformly on $[0, 1] \times [t^*, T]$, and therefore $p, v \in C^0([0, 1] \times [t^*, T])$. Moreover, we observe that $p(x, t) \to 0$ and $v(x, t) \to 0$, as $t \to t^*$.

Remark 23. Note that if $t^* = 0$, then this is the case of a continuous constant boundary source $P(t) = 1$, for all $t\geq 0$. Then the solid displacement, the fluid pressure and the discharge velocity have the following representations for $t \geq 0$:

All three series in (7.6) converge absolutely and uniformly on $[0, 1] \times [0, T]$, and therefore $p, v \in C^0([0, 1] \times [0, T])$.

7.1.1.   The case $t^*=0$

Figure 2 illustrates displacement, fluid pressure and discharge velocity (left, middle and right panel, respectively) as a function of $t$. Displacement and velocity are evaluated at $x = 1$ (right boundary) whereas the pressure is evaluated at $x = 0$ (left boundary). Within the observational time interval, the displacement decreases monotonically till it reaches an asymptotic value of approximately $50\; \mathrm{{\mu m}}$. Conversely, pressure and discharge velocity exhibit a nonmonotonic behaviour, characterized by a rapid increase, attaining a maximum value of approximately $2.3\; \mathrm{{kPa}}$ and $0.002\; \mathrm{{\mu m s^{-1}}}$, respectively, followed by a monotonic decrease that approaches zero asymptotically The behavior of the simulated solutions are consistent with those reported in ^[12,16].

Figure 3 illustrates displacement, fluid pressure and discharge velocity (left, middle and right panel, respectively) as a function of $x$ and $t$. We see that, after an initial transient due to structural viscoelasticity, the solid displacement tends to a linear behavior along the domain length, being maximum (in absolute value) at $x = L$. Mathematically, this behavior is due to the fact that, for long times, the viscoelastic terms in the series in (7.6a) become negligible with respect to the linear elastic component. Both pressure and discharge velocity distributions decrease as time $\rightarrow +\infty$, in agreement with the fact that definition (6.31) implies that $\lim\limits_{t\rightarrow +\infty} G_n(t) = 0$ for all $n\geq 0$.

Figure 4 illustrates the energies per unit area $\mathcal{E}_\theta$, $\mathcal{E}_{c_0}$ and $\mathcal{E}_{tot}$ (left, middle and right panels, respectively) as a function of $t$. In agreement with the previous analysis for the displacement and pressure profiles, we see that $\mathcal{E}_\theta$ tends to a constant value as time increases because the deformation of the mixture tends to become uniform in all the domain. At the same time, $\mathcal{E}_{c_0}$ decreases in time following the decrease of the pressure, being significant only during the initial transient. As a result, the total potential energy $\mathcal{E}_{tot}$ of the mixture tends to coincide with $\mathcal{E}_\theta$, as demonstrated by the right panel of Figure 4.

Figure 5 illustrates the dissipation $\mathcal{D}$ and force term $\mathcal{F}$ (left and right panel, respectively) as a function of $t$. The temporal profiles of both functions rapidly decay to zero as time increases. This is consistent with the fact that the domain of the mixture tends to deform in a uniform manner as time gets large so that its time variation becomes rapidly negligible.

7.1.2.   The case $t^* > 0$

We assume that the external traction applied at the right boundary $x = L$ has a jump discontinuity at $t^* = 0.25 T$.

Figure 6 illustrates the computed displacement, fluid pressure and discharge velocity (left, middle and right panel, respectively) as a function of $t$. Displacement and velocity are evaluated at $x = L$ (right boundary) whereas the pressure is evaluated at $x = 0$ (left boundary). We notice that the three graphs are the translation of the corresponding graphs in Figure 2. In particular, we see that $u$, $p$ and $v$ are continuous at $t = t^*$, where their value is equal to zero.

Figure 7 illustrates the displacement, fluid pressure and discharge velocity (left, middle and right panel, respectively) as a function of $x$ and $t$.

Figure 8 illustrates the energies per unit area $\mathcal{E}_\theta$, $\mathcal{E}_{c_0}$ and $\mathcal{E}_{tot}$ (left, middle and right panels, respectively) as a function of $t$.

Figure 9 illustrates the dissipation $\mathcal{D}$ and force term $\mathcal{F}$ (left and right panel, respectively) as a function of $t$. We notice that both $\mathcal D$ and $\mathcal F$ are discontinuous at $t = t^*$ where they experience a finite jump.

7.2.   The case of unbounded ramp pulse

Let $P\left(t\right)$ be the dimensionless unbounded ramp pulse of unit-slope starting at $t = 0$, represented in Figure 10

In this case, we have

so that:

Summing the two above expressions we obtain

Note that the above expression can be rewritten as

The time derivative of expression (7.10) is given by: for $t \geq 0$,

Note that

and the coefficients $B_n(t)$ become

Then the displacement, the pressure and the discharge velocity have the following representations for $t \geq 0$:

Figure 11 illustrates the displacement, fluid pressure and discharge velocity (left, middle and right panel, respectively) as a function of $t$. Displacement and velocity are evaluated at $x = L$ (right boundary) whereas the pressure is evaluated at $x = 0$ (left boundary). Unlike the case presented in Section 7.1, here the external pressure load continues to increase linearly with time, thereby inducing a continuous increase in displacement, pressure and velocity as time goes by. Over the observational time interval, the magnitude of the external load is smaller than what considered in Section 7.1 and this leads to a smaller (absolute) value of the maximum displacement, which is about $10\; \mathrm{{\mu m}}$ in this case.

Figure 12 illustrates the displacement, fluid pressure and discharge velocity (left, middle and right panel, respectively) as a function of $x$ and $t$. The displacement profile exhibits an approximately bilinear variation with respect to temporal and spatial coordinates, whereas pressure and velocity display a nonlinear behavior in the space-time domain. All dependent variables tend to increase in magnitude as a function of time at any spatial position of the mixture, the discharge velocity being closer to reach a stationary condition than the pressure.

Figure 13 illustrates the energies per unit area $\mathcal{E}_\theta$, $\mathcal{E}_{c_0}$ and $\mathcal{E}_{tot}$ (left, middle and right panels, respectively) as a function of $t$. Interestingly, $\mathcal{E}_\theta$ exhibits a nonlinear increase with respect to time in accordance with the fact that deformation is not constant in space. $\mathcal{E}_{c_0}$ follows a similar pattern because of the nonlinear trend of the pressure, but has a much smaller value, so that the $\mathcal{E}_{tot}$ almost coincides with $\mathcal{E}_\theta$ for all times.

Figure 14 illustrates the dissipation $\mathcal{D}$ and forcing term $\mathcal{F}$ (left and right panel, respectively) as a function of $t$. Dissipation increase with time is characterized by two markedly different slopes. In a first time interval, approximately equal to $T/5$, dissipation increases rapidly and is mainly determined by the fluid component of the mixture. In the remaining part of the observational time interval, dissipation increases less rapidly and is mainly determined by the structural viscoelasticity of the mixture. The forcing term increases linearly with time in accordance with the trend of the temporal variation of solid displacement at the right boundary of the domain, as shown in Figure 12.

7.3.   The case of bounded ramp pulse

Let $P^{\left(\varepsilon \right) }\left(t\right)$ be the dimensionless ramp pulse of unit amplitude and finite rise time ($= \varepsilon > 0$) starting at $t = 0$

represented graphically in Figure 15.

Using the linear superposition principle, the solid displacement, fluid pressure, and discharge velocity are given by:

where the expressions of $u_{0}\left(x, t\right)$, $p_{0}\left(x, t\right)$ and $v_{0}\left(x, t\right)$ are given by (7.13).

Notice that $p_{0}^{\varepsilon }\left(x, t\right), v_{0}^{\varepsilon }\left(x, t\right) \in C^{0}\left(Q\right)$ since $p_{0}\left(x, t\right), v_{0}\left(x, t\right) \in C^{0}\left(Q\right)$ and $p_{0}\left(x, 0\right) = v_{0}\left(x, 0\right) = 0$.

Figure 16 illustrates the displacement, fluid pressure and discharge velocity (left, middle and right panel, respectively) as a function of $t$. Displacement and velocity are evaluated at $x = L$ (right boundary) whereas the pressure is evaluated at $x = 0$ (left boundary). We see that the displacement increases in magnitude almost linearly during the increase in time of the externally applied pressure. Then, it rapidly tends to stationary conditions once the external applied pressure becomes constant. The asymptotic value of the displacement is the same as in the case of the step pulse illustrated in Section 7.1. The profile of fluid pressure increases in time and the externally applied pressure increases; once the solid deformation attains stationary conditions, the fluid pressure decreases. A similar trend is shown by the discharge velocity. The maximum value of pressure is the same as in the case of a step pulse external pressure whereas the maximum value of the discharge velocity is slightly smaller and coincides with the value of the velocity at $t = T/2$ that is obtained in the case of an unbounded ramp external pressure (cf. Figure 11, right panel).

Figure 17 illustrates the displacement, fluid pressure and discharge velocity (left, middle and right panel, respectively) as a function of $x$ and $t$. The spatial variation of the displacement becomes linear after the external pressure ceases to increase, whereas, in the same time interval, both fluid pressure and discharge velocity exhibit a spatial decrease.

Figure 18 illustrates the energies per unit area $\mathcal{E}_\theta$, $\mathcal{E}_{c_0}$ and $\mathcal{E}_{tot}$ (left, middle and right panels, respectively) as a function of $t$. Both $\mathcal{E}_\theta$ and $\mathcal{E}_{c_0}$ increase in time until the external pressure increases. Then, $\mathcal{E}_\theta$ becomes constant since deformation is constant in time, whereas $\mathcal{E}_{c_0}$ decreases because of the decrease of fluid pressure. As in all previously considered examples, the contribution of $\mathcal{E}_{c_0}$ is much smaller than that of $\mathcal{E}_{\theta}$ so that the $\mathcal{E}_{tot}$ almost coincides with $\mathcal{E}_\theta$.

Figure 19 illustrates the dissipation $\mathcal{D}$ and forcing term $\mathcal{F}$ (left and right panel, respectively) as a function of $t$. Both terms exhibit an increase with respect to time during the increase of the externally applied pressure. Then, they both experience a sudden decay once the externally applied pressure becomes constant. Dissipation tends to an asymptotic value that is much smaller than the value attained at $t = T/2$ whereas the force term tends to zero since the boundary displacement is almost constant in time after $t = T/2$.

8.   Dependence of the solution on compressibility

In this section we investigate the dependence of the solution of model (4.12) on the compressibility parameter $c_0$. In terms of the dimensionless equation system (5.3), this amounts to analyzing the solutions as a function of the quantity $\gamma$ defined in (5.2). We denote by $c_{0, ref} = 1.67 \cdot 10^{-5} \mathrm{{m^2 N^{-1}}}$ the reference value of the compressibility parameter. This value has been used in all the computations illustrated in Section 7. Then, we let $c_0$ assume the following values:

and we compute the solid displacement and discharge velocity at $x = L$, the fluid pressure at $x = 0$ and the energies $\mathcal E_{\theta}$, $\mathcal E_{c_0}$ and $\mathcal E_{tot}$, as functions of time in the interval $[0, \, T]$, $T = 10000 \mathrm{{s}}$ being the width of the observational window considered in Section 7. All the other model parameters have been set equal to the values adopted in Section 7. The applied boundary traction is a step pulse of amplitude $P_{ref}$ at $t^* = 0$.

Figure 20 shows the profiles of solid displacement (left panel), fluid pressure (middle panel) and discharge velocity (right panel) as functions of time and of the compressibility parameter $c_0$. We notice that compressibility has a significant impact on the quantitative values attained by the solution variables. In particular we see that increasing $c_0$ gives rise to:

● an increase of the magnitude of the solid displacement;

● a decrease of pressure and velocity;

● a right shift of the peak of pressure and velocity.

These behaviors are indicative of the fact that the compressibility of the mixture components allows the body to deform more under the same pressure load, thereby reducing the internal level of fluid pressure and limiting the impact on the fluid velocity. We also notice that decreasing $c_0$ has a much more significant impact than increasing $c_0$ on solution range variation. From the theoretical viewpoint, this is due to the fact that $\gamma \rightarrow 1^-$ for large values of $c_0$ and, as a consequence, the ratio $(1-\gamma)/(\gamma \eta)$, that characterizes the mathematical form of the solution expressions (7.5), tends to 0. From the physical viewpoint, decreasing the value of $c_0$ implies that the body cannot significantly deform upon applying a pressure load, thereby inducing higher levels of fluid pressure whose gradients lead to larger fluid velocities.

Figure 21 illustrates the discharge velocity computed by the model studied in this work and the discharge velocity computed by the model studied in ^[4]. The principal difference between the two models is that no compressibility was included in ^[4]. We see that:

● the velocity predicted by the model of ^[4] is a (very sharp) upper bound for all the velocities predicted by the model studied in the present work, when $c_0 \rightarrow 0^+$;

● for increasing values of $c_0$, the velocities predicted by the present model differ substantially from the upper bound velocity yielded by the model in ^[4]. In particular, we see that increasing $c_0$ gives rise to a decrease and to a right shift of the peak in the velocity profile.

Figure 22 shows the profiles of $\mathcal E_{\theta}$, $\mathcal E_{c_0}$ and $\mathcal E_{tot}$ as a function of $t$ and of $c_0$. We see that increasing $c_0$ has the effect of increasing substantially and monotonically the magnitude of $\mathcal E_\theta$ because the deformation profile gets larger as the components become more compressible, in accordance with Figure 20, left panel. On the contrary, $\mathcal E_{c_0}$ exhibits a nonmonotonic dependence on $c_0$. Specifically, for $c_0 < c_{0, ref}$ we see that increasing $c_0$ gives rise to an increase of $\mathcal E_{c_0}$ with a right shift of its peak; then, for $c_0 \geq c_{0, ref}$, we see that increasing $c_0$ leads to a significant monotonic decrease of $\mathcal E_{c_0}$. For the theoretical viewpoint, this is due to the fact that the ratio $\gamma/(1-\gamma) \to + \infty$ as $c_0 \rightarrow +\infty$ and the pressure profile tends to decrease in accordance with Figure 20, middle panel. The physical meaning of these results can be appreciated by observing that the predicted total energy is mainly determined by $\mathcal E_{\theta}$ for large values of $c_0$, since larger deformations can occur for more compressible components, whereas the contribution to the total energy given by $\mathcal E_{c_0}$ becomes more significant for smaller values of $c_0$, since larger pressures develop for less compressible components.

9.   Conclusions

The analysis presented in this work shows that the solutions of a poroviscoelastic model with compressible components remain bounded even in the case when the imposed boundary traction is irregular in time. In particular, given a certain functional form for the boundary traction and a certain level of structural viscoelasticity, the discharge velocity attains a maximum value that is lower when compressibility is higher. By investigating the dynamics of the energy functionals characterizing the system, we showed that this limiting effect is due to the capability of the system to store potential energy as its components are elastically compressed, thereby delaying the transmission of irregularities in the linear momentum from the solid to the fluid. As a result, the fluid has the time to accommodate for sudden changes, resulting in bounded velocities. This mechanism is very different from that provided by structural viscoelasticity, whose limiting effect on the discharge velocity is due to increased viscous dissipation, as shown in ^[4].

Ultimately, this work elucidates the specific role that compressibility plays in the control of fluid flow through complex deformable porous structure, which finds numerous applications in science and engineering. The work presented here offers many future directions of research. For example, it would be very interesting to investigate how the findings concerning the role of compressibility would translate to a more realistic three-dimensional setting. Furthermore, different viscoelastic models could be considered and the specific roles of fluid and solid viscosities could be investigated and compared.

Acknowledgments

Dr. Bociu has been partially supported by NSF CAREER DMS-1555062. Dr. Guidoboni has been partially supported by the award NSF DMS-1853222. Dr. Sacco has been partially supported by Micron Semiconductor Italia S.r.l., SOW nr. 4505462139.

Conflict of interest

The authors declare there is no conflict of interest.