ON THE ROLE OF COMPRESSIBILITY IN POROVISCOELASTIC MODELS

. In this article we conduct an analytical study of a poroviscoelas-tic mixture model stemming from the classical Biot’s consolidation model for poroelastic media, comprising a ﬂuid component and a solid component, coupled with a viscoelastic stress-strain relationship for the total stress tensor. The poroviscoelastic mixture is studied in the one-dimensional case, corresponding to the experimental conditions of conﬁned compression. Upon assuming ( i ) negligible inertial eﬀects in the balance of linear momentum for the mixture, ( ii ) a Kelvin-Voigt model for the eﬀective stress tensor and ( iii ) a constant hydraulic permeability, we obtain an initial value/boundary value problem of pseudo-parabolic type for the spatial displacement of the solid component of the mixture. The dimensionless form of the diﬀerential equation is characterized by the presence of two positive parameters γ and η , representing the contributions of compressibility and structural viscoelasticity, respectively. Explicit solutions are obtained for diﬀerent functional forms characterizing the boundary traction. The main result 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 allows the system to store potential energy as its components are elastically compressed, thereby providing an additional mechanism that limits the maximum of the discharge velocity when the imposed boundary traction is irregular in time.

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.[2,8,9,17].Skeleton viscoelasticity is often that are present in biological tissues.Specifically, our work focuses on the bio-fluidmechanical 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 [3,5,7,12].
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 blowup in the solution of poroelastic models [2,17].Interestingly, the blow-up can be prevented by including structural viscoelasticity [17].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 [10], whereas rapid changes in the gravitational acceleration are experienced by jet pilots and astronauts during flights [6,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 [17], 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 nonsmooth 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 [13], 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 [8,9], 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 (0, L) and in the temporal observational interval (0, T ] can be written as: where σ is the total stress, ζ is the fluid content and v is the discharge velocity.Equations (1) are complemented with the following constitutive laws: Eq. (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 σ 0 , which is assumed to be characterized by viscoelastic behavior of Kelvin-Voigt type (see Eq. (2b) where u denotes the solid phase displacement).Eq. (2c) expresses the fact that the fluid content is altered by changes in fluid pressure and solid deformation.Eq. ( 2d) is Darcy's law relating the discharge velocity v with the pressure gradient by means of the hydraulic permeability K.The Biot coefficient α, the storage coefficient c 0 , the material parameters θ and η and the permeability K are assumed to be given positive constants.
In the classical Biot's theory, the material parameter θ can be expressed as θ = 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 θ and η as the elastic and viscoelastic parameters, respectively [13].In addition, we will use the notation K = 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 [14]), namely: Finally, we complete the formulation of the problem by prescribing the following initial conditions: We notice that (4a) and (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 α = 1 in the model described above, as detailed in [1,4].The study of analytical solutions for the incompressible model was addressed in [2,17].
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 (1a) by ∂u/∂t ∈ L 2 (0, L); integrating over (0, L) and using the constitutive laws (2a) and (2b) and the boundary conditions (3a), (3c) and (3d),we obtain Multiplying (1b) by p ∈ L 2 (0, L), integrating over (0, L) and using the constitutive laws (2c) and (2d) and the boundary conditions (3b) and (3c), we obtain Subtracting ( 6) from (7) we get the following evolution equation for the total energy 5 stored in the poroviscoelastic system where the energy functional E tot = E tot (t), the dissipation functional D = D(t) and the force term F = F(t) are defined as: Remark 2. The physical units of 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, E tot is obtained by integrating in x between 0 and L the energy density ε tot defined as The units of ε tot are J m −3 .Analogously, the units of D and F are J m −2 s −1 .
Since E tot ≥ 0 and D ≥ 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 θ contribute to the total energy in the form of potential energy, so that we can write Conversely, the terms proportional to K 0 and η contribute to dissipate energy via viscous effects within the fluid and solid components.We notice that F does not have a definite sign since it depends on the boundary terms.The energy identity (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 (1)-( 4) solely in terms of the solid displacement.Combining (1a) and (3d), we obtain that the total stress is given by σ(x, t) = −P (t), for all x ∈ (0, L] and for all t ∈ (0, T ] , and, moreover, the fluid pressure p and the discharge velocity v can be written in terms of the solid displacement u as: where σ 0 = σ 0 (u (x, t)) is given by (2b).Let us now derive the problem satisfied by u.Integrating Eq. (1b), where ζ is given by (2c), over the space interval [0, x] and taking the boundary conditions (3) into account, yields Now we substitute p and v by expressions (10) and note that x 0 σ 0 dy = θu + η ∂u ∂t by virtue of (3a).Then we obtain The boundary conditions are (3a) and σ 0 | x=L = −P (t) (coming from (3c) and (10a)).The initial conditions are (4a) and, on using (4b), (10a), (2b) and (5a), or equivalently u 1 being an arbitrary constant.Note that u 1 = 0 if and only if the compatibility condition with (3a) ∂u ∂t (0, 0) = 0 (11) is satisfied.Summing up: Find u = u (x, t) such that: Remark 3. The solution u (x, t) of the problem (12) is the sum u 1 (x, t) + u 0 (x, t) of the solution u 1 (x, t) of ( 12) where P = 0, and of the solution u 0 (x, t) of (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 = 0): then u 1 (x, t) 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 = 0), then the corresponding change of the solid displacement is represented by u 0 (x, t).Therefore, the compatibility condition (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 η > 0, i.e. the system is characterized by finite compressibility and structural viscoelasticity.
Remark 5.The IBVP (12) has a markedly different character compared to the linear poroviscoelastic system studied in [17] because of the second-order time derivative on the left-hand side of (12a).
Remark 6.If we define the following quantities: then the partial differential equation (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 fluidsolid 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 ( 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 [17], by introducing the following scaling factors: We also define the non-dimensional quantity Note that 0 < γ < 1.We recover the same definitions of the scaling factors as in [17] by setting α = 1 and c 0 = 0 in (13).
Remark 7.For notational simplicity we will drop the ' • ' 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 (x, t) is known, the functions σ, p and v can be computed as follows: Lastly, the dimensionless energy equation is written again in the form (8) where 6. Well-posedness and Regularity of Solution.We make the following assumption on the boundary traction in (15).
We recall that a function P (t) is piecewise smooth on [0, T ] if both P and its derivative P are continuous on [0, T ], except possibly at a finite number of points in (0, T ), where they have simple jump discontinuities.
6.1.Auxiliary Problem.Define Using Assumption 8 we have that U (t) is absolutely continuous on [0, T ] and Now we introduce the change of variable and note that w solves the following auxiliary IBVP with homogeneous boundary data: w(x, 0) = 0 in (0, 1), (21d) where the volumetric source and initial datum are given by: We shall prove the existence and uniqueness of the solution w (x, t) for a very general class of data f (x, t) and ϕ (x).For sake of exposition we write H = L 2 (0, 1), and define the real Hilbert space endowed with the equivalent norm v V = ∂v/∂x L 2 (0,1) , due to Poincaré's inequality.
Now we make the following assumptions on the functions f (x, t) and ϕ (x): and write w (t) = w (•, t), w (t) = ∂w (•, t) /∂t, etc.We then define weak solutions of (21) as follows.
or, equivalently, x = 0 is included in the regularity requirement that w (t) ∈ V , whereas the boundary condition at x = 1 is satisfied in a weak sense.

A-priori Estimates.
Lemma 12. Let w be a weak solution of (21).Then there are constants C i 's, depending only on γ, η and T , such that the following estimates hold for t pointwise in [0, T ]: Proof.Using (γηw + w) ∈ L 2 (0, T ; V ) as multiplier in (26), 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 (28) and using Gronwall's Inequality, estimate (27a) follows.Similarly, by dropping the first and third terms on the left-hand side and using (27a), we get (27b).Lastly, we use triangle inequality and the embedding V → H to write Then, estimate (27c) follows from (27a) and (27b).
Lemma 13.Let w be a weak solution of (21).Then there is a constant C, de-10 pending only on γ, η and T , such that Proof.From (27b) it immediately follows In addition, Eqs. ( 28) and (27a) give for a suitable C = C (γ, η, T ).Lastly, using the definition of a weak solution, we have that for This implies that from which, using estimates (27b) and (27a), we obtain The following corollary is an immediate consequence of Lemma 13.
Corollary 14. (Uniqueness and continuous dependence on data) The weak solution to problem (21) is unique and depends continuously on the data.

Existence of Solution.
The IBVP (21) 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 √ 2X n (x) forms a Hilbert space basis for H, whereas the sequence of functions 2 λn X n (x) forms a Hilbert space basis for V (see [16]).
The solution w of (21) has the expansion where T n (t) can be recovered using the data.Similarly, we use the basis {X n (x)} to represent ϕ and f as follows: where ϕ n and f n (t) are the Fourier coefficients of ϕ and f (•, 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 ≥ 0: and initial conditions The characteristic equation for the homogeneous counterpart of (39) is given by Since the discriminant of the equation is then for each n ≥ 0 the characteristic equation has two real, negative, distinct roots r 1n = −Λ 1n and r 2n = −Λ 2n , where Remark 16.Λ 1n and Λ 2n satisfy the following relations: Therefore the solution for the homogeneous ODE (39) is given by T b n e −Λ2nt .The particular solution is obtained from variation of parameters formula.The Wronskian is given by W (t) = (Λ 2n − Λ 1n )e −Λ1nt e −Λ2nt , so that the particular solution has the following form We introduce the following notation 14LORENA BOCIU * , GIOVANNA GUIDOBONI, RICCARDO SACCO, AND MAURIZIO VERRI and then write the solution to (39) as where we used the following formula for convolution Now we use (44) back into (34) and impose the initial conditions.We have w(x, 0) = a n + b n = 0, and In conclusion, we obtain We note that the terms G n (t) given in (43) satisfy the following estimates.
Now we can state and prove our well-posedness result.
Theorem 18. (Well-posedness of problem (21)) For every f ∈ L 2 (0, T ; H) and ϕ ∈ H there is a unique weak solution of (21), 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 (45b) satisfies conditions (D1)-(D3) of Definition 10.
(D1) (A) First we show that w ∈ L 2 (0, T ; V ).For all t ∈ [0, T ], we have Using Lemma 17, we obtain1 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 ∈ L ∞ (0, T ; H) ⊂ L 2 (0, T ; H).In order to show that w ∈ L 2 (0, T ; V ), we first use the Monotone Convergence Theorem to write We use multiplier T n in (39) and the initial conditions (40) and obtain the following estimate: Now we use (50) back into (49), and take advantage of estimate (48) to obtain and this gives the desired conclusion that w ∈ L 2 (0, T ; V ).
16LORENA BOCIU * , GIOVANNA GUIDOBONI, RICCARDO SACCO, AND MAURIZIO VERRI (C) Left to show that w ∈ L 2 (0, T ; V ).For any v (x) = ∞ n=0 c n X n (x) ∈ V and for t ∈ [0, T ], we use (39) to define the linear functional w (t) as follows Hence, by virtue of the embedding V → H, we obtain that This shows that w (t) ∈ V and or, equivalently, In conclusion, using the estimates in part (A) and part (B), we obtain that w ∈ L 2 (0, T ; V ).
(D2) Now we show that w satisfies condition (D2).Since the 2 λn X n is a basis in V , it suffices to consider the test function v = X n .For t pointwise almost everywhere in (0, T ), we use (39) and obtain where Proof.(A) The result follows from Lemma 17, the fact that |X n (x)| ≤ 1 and Weierstrass Test.Indeed, for all x ∈ [0, 1] and t ∈ [0, T ], and every n ≥ 0, we have Now, due to Lemma 17, we have and Applying Weierstrass Test, we obtain that the series for every n ≥ 0, and thus w(x, t) ∈ C 0 (Q).

(B)
The second order weak derivative in space of w is given by Then we use estimate (46a) in Lemma 17 and obtain 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: 18LORENA BOCIU * , GIOVANNA GUIDOBONI, RICCARDO SACCO, AND MAURIZIO VERRI Proof.(A) For all x ∈ [0, 1] and t ∈ [0, T ], and every n ≥ 0, we have Now, due to Lemma 17, we have and Applying Weierstrass Test, we obtain that the series (B) Like in the proof of (52a), we have (C) Since the second order weak derivative in space of w is given by In order to estimate the right hand side, we use multiplier λ n T n on the ODE (39) that T n solves to obtain Now we integrate from 0 to t and use the initial conditions (40): Thus, using (55) back into (54), we obtain and the assertion follows from estimate (29).
6.5.Analytical formulas for solid displacement and discharge velocity.Now we are in a position to return to our original linear 1D problem (15).The solid 5 displacement solution u (x, t) of ( 15) is the sum where we recall that and w (x, t) is the unique weak solution of the auxiliary problem (21) with special data (22).From the Fourier expansions as well as the Fourier series (35) of ϕ and f with respect to the basis {X n }, 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 (x, t) ∈ C 0 (Q) since the physical data (22) satisfy f ∈ L 2 (0, T ; V ) and ϕ ∈ H. Notice, however, that our special ϕ belongs to V if and only if u 1 = 0, i.e. ϕ = 0 and u 1 (x, t) ≡ 0.
20LORENA BOCIU * , GIOVANNA GUIDOBONI, RICCARDO SACCO, AND MAURIZIO VERRI By virtue of (58), the discharge velocity (16c) is given by where Remark 22.It should be stressed that the Fourier coefficients of v 1 (x, t) are O (n), 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 (t) is considered.
7. Continuous vs. discontinuous boundary sources.In this section, we analyze the behavior of the solutions obtained for model (15) 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 [15]; (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 ( 13) 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 E θ , E c0 and E tot , the non-dimensional dissipation functional D and the non-dimensional force term F using the expressions (17).The resulting non-dimensional energies are then multiplied by the scaling factor [E tot ] introduced in (13j) and the same is done for the non-dimensional dissipation and force terms that are multiplied by the scaling factor [D] introduced in (13i).
symbol value units L 0.81 Table 1.Numerical values of model parameters utilized in the numerical simulations.
Table 1 reports the numerical values of model parameters (in dimensional form) that are used in the next sections.The values for L, T , θ, K 0 and P ref have been chosen consistently with the experimental data illustrated in [15].To observe the asymptotic behavior of the modeled system, T has been taken equal to 10 5 s in the example illustrated in Section 7.1.In this example, the corresponding value of η is 4.85 • 10 9 N s m −2 .The value of the Biot coefficient α has been set equal to 0.9, whereas in [15] it was equal to 1 since both mixture components were assumed to be incompressible.The values of the viscoelastic parameter η and of the compressibility parameter c 0 (which were not present in the model studied in [15]) have been set equal to θτ e and 1/P ref , respectively, where τ 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.
Table 2. Numerical values of model parameters in dimensionless form.
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 (43) and using the fact that (G n * U ) = G n * U , we obtain the following simplification for the coefficients B n (t) 7.1.The case of step pulse at t = t * .Let t * ∈ [0, T ).We consider the following boundary source: 22LORENA BOCIU * , GIOVANNA GUIDOBONI, RICCARDO SACCO, AND MAURIZIO VERRI Figure 1 gives a graphical representation of P (t) for t * ∈ (0, T ).Replacing ( 63) into ( 18) we obtain We can now compute the convolution needed in (62a): We observe that for t ≥ t * , and therefore the coefficients (62d) 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 ≥ t * : Note that all three series in the expressions (66) converge absolutely and uniformly on [0, 1] × [t * , T ], and therefore p, v ∈ C 0 ([0, 1] × [t * , T ]).Moreover, we observe that Remark 23.Note that if t * = 0, then this is the case of a continuous constant boundary source P (t) = 1, for all t ≥ 0. Then the solid displacement, the fluid pressure and the discharge velocity have the following representations for t ≥ 0:  Conversely, pressure and discharge velocity exhibit a nonmonotonic behaviour, characterized by a rapid increase, attaining a maximum value of approximately 2.3kPa and 0.002µms −1 , respectively, followed by a monotonic decrease that approaches zero asymptotically The behavior of the simulated solutions are consistent with those reported in [15,13].
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 (67a) become negligible with respect to the linear elastic component.Both pressure and discharge velocity distributions decrease as time → +∞, in agreement with the fact that definition (43) implies that lim t→+∞ G n (t) = 0 for all n ≥ 0.
Figure 4 illustrates the energies per unit area E θ , E c0 and 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 E θ 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, E c0 decreases in time following the decrease of the pressure, being significant only during the initial transient.As a result, the total potential energy E tot of the mixture tends to coincide with E θ , as demonstrated by the right panel of Figure 4.   Figure 5 illustrates the dissipation D and force term 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.25T .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.
In this case, we have Summing the two above expressions we obtain Note that the above expression can be rewritten as The time derivative of expression (71) is given by: for t ≥ 0, in displacement, pressure and velocity as time goes by.Over the observational time 32LORENA BOCIU * , GIOVANNA GUIDOBONI, RICCARDO SACCO, AND MAURIZIO VERRI 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µ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 E θ , E c0 and E tot (left, middle and right panels, respectively) as a function of t.Interestingly, E θ exhibits a nonlinear increase with respect to time in accordance with the fact that deformation is not constant in space.E c0 follows a similar pattern because of the nonlinear trend of the pressure, but has a much smaller value, so that the E tot almost coincides with E θ for all times.
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 (x, t), p 0 (x, t) and v 0 (x, t) are given by (74).
36LORENA BOCIU * , GIOVANNA GUIDOBONI, RICCARDO SACCO, AND MAURIZIO VERRI Notice that p ε 0 (x, t) , v ε 0 (x, t) ∈ C 0 (Q) since p 0 (x, t) , v 0 (x, t) ∈ C 0 (Q) and p 0 (x, 0) = v 0 (x, 0) = 0.  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).As in all previously considered examples, the contribution of E c0 is much smaller than that of E θ so that the E tot almost coincides with E θ .
Figure 19 illustrates the dissipation D and forcing term 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 ( 12) on the compressibility parameter c 0 .In terms of the dimensionless equation system (15), this amounts to analyzing the solutions as a function of the quantity γ defined in (14).We denote by c 0,ref = 1.67 • 10 −5 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:   40LORENA BOCIU * , GIOVANNA GUIDOBONI, RICCARDO SACCO, AND MAURIZIO VERRI 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 γ → 1 − for large values of c 0 and, as a consequence, the ratio (1 − γ)/(γη), that characterizes the mathematical form of the solution expressions (66), 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 [17].The principal difference between the two models is that no compressibility was included in [17].We see that: • the velocity predicted by the model of [17] is a (very sharp) upper bound for all the velocities predicted by the model studied in the present work, when c 0 → 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 [17].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 E θ , E c0 and 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 E θ because the deformation profile gets larger as the components become more compressible, in accordance with Figure 20, left panel.On the contrary, E c0 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 E c0 with a right shift of its peak; then, for c 0 ≥ c 0,ref , we see that increasing c 0 leads to a significant monotonic decrease of E c0 .For the theoretical viewpoint, this is due to the fact that the ratio γ/(1 − γ) → +∞ as c 0 → +∞ 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 E θ for large values of c 0 , since larger deformations can occur for more compressible components, whereas the contribution to the total energy given by E c0 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 42LORENA BOCIU * , GIOVANNA GUIDOBONI, RICCARDO SACCO, AND MAURIZIO VERRI 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 [17].
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.

Definition 10 .
[Weak solution of problem (21)] A function w : [0, T ] → V is a weak solution of the auxiliary problem (21) if:

6 . 4 .
D3)-(D4) As stated in Remark 11, condition (D1) implies that w ∈ C 0 ([0, T ] ; V ) and w ∈ C 0 ([0, T ] ; H).Moreover, solution w(x, t) satisfies the initial conditions (D3)-(D4) by virtue of (40).This concludes the proof of our well-posedness theorem.Regularity of the Solution.We have established the existence and uniqueness of the weak solution w to the auxiliary problem (21).Now we examine its regularity.Proposition 19.(Regularity of w and w xx ) Under the hypotheses in Theorem 18, the following is true:

Figure 1 .
Figure 1.The step pulse at t = t * .

Figure 2 . 7 . 1 . 1 .
Figure 2. Left panel: solid displacement u at x = L as a function of t.Middle panel: fluid pressure p at x = 0 as a function of t.Right panel: discharge velocity v at x = L as a function of t.The applied boundary traction is a step pulse of amplitude P ref at t * = 0.7.1.1.The case t * = 0. Figure2illustrates 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µm.Conversely, pressure and discharge velocity exhibit a nonmonotonic behaviour, characterized by a rapid increase, attaining a maximum value of approximately 2.3kPa and 0.002µms −1 , respectively, followed by a monotonic decrease that approaches zero asymptotically The behavior of the simulated solutions are consistent with those reported in[15,13].

Figure 3 .
Figure 3. Top panel: solid displacement u as a function of x and t.Middle panel: fluid pressure p as a function of x and t.Bottom panel: discharge velocity v as a function of x and t.The applied boundary traction is a step pulse of amplitude P ref at t * = 0.

Figure 4 .
Figure 4. Simulated profiles of E θ (left), E c0 (middle) and E tot as a function of t.The applied boundary traction is a step pulse of amplitude P ref at t * = 0.

Figure 5 .
Figure 5. Left panel: dissipation as a function of t.Right panel: force term as a function of t.The applied boundary traction is a step pulse of amplitude P ref at t * = 0.

Figure 6 .
Figure 6.Left panel: solid displacement u at x = L as a function of t.Middle panel: fluid pressure p at x = 0 as a function of t.Right panel: discharge velocity v at x = L as a function of t.The applied boundary traction is a step pulse of amplitude P ref at t * = 0.25T .

Figure 8
Figure8illustrates the energies per unit area E θ , E c0 and E tot (left, middle and right panels, respectively) as a function of t.

Figure 7 .
Figure 7. Top panel: solid displacement u as a function of x and t.Middle panel: fluid pressure p as a function of x and t.Bottom panel: discharge velocity v as a function of x and t.The applied boundary traction is a step pulse of amplitude P ref at t * = 0.25T .

Figure 8 .
Figure 8. Simulated profiles of E θ (left), E c0 (middle) and E tot as a function of t.The applied boundary traction is a step pulse of amplitude P ref at t * = 0.25T .

Figure 9 .
Figure 9. Left panel: dissipation as a function of t.Right panel: force term as a function of t.The applied boundary traction is a step pulse of amplitude P ref at t * = 0.25T .

Figure 10 .
Figure 10.The unbounded ramp pulse.7.2.The case of unbounded ramp pulse.Let P (t) be the dimensionless unbounded ramp pulse of unit-slope starting at t = 0, represented in Figure10

Figure 12 .
Figure 12.Top panel: solid displacement u as a function of x and t.Middle panel: fluid pressure p as a function of x and t.Bottom panel: discharge velocity v as a function of x and t.The applied boundary traction is an unbounded ramp pulse.

Figure 13 .
Figure 13.Simulated profiles of E θ (left), E c0 (middle) and E tot as a function of t.The applied boundary traction is an unbounded ramp pulse.

Figure 14 .
Figure 14.Left panel: dissipation as a function of t.Right panel: force term as a function of t.The applied boundary traction is an unbounded ramp pulse.

Figure 14
Figure 14 illustrates the dissipation D and forcing term F (left and right panel, respectively) as a function of t.Dissipation increase with time is characterized by

Figure 16 .
Figure 16.Left panel: solid displacement u at x = L as a function of t.Middle panel: fluid pressure p at x = 0 as a function of t.Right panel: discharge velocity v at x = L as a function of t.The applied boundary traction is a bounded ramp pulse with ε = 0.5T .

Figure 16
Figure16illustrates 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

Figure 17 .
Figure 17.Top panel: solid displacement u as a function of x and t.Middle panel: fluid pressure p as a function of x and t.Bottom panel: discharge velocity v as a function of x and t.The applied boundary traction is a bounded ramp pulse with ε = 0.5T .

Figure 17
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 19 .
Figure 19.Left panel: dissipation as a function of t.Right panel: force term as a function of t.The applied boundary traction is a bounded ramp pulse with ε = 0.5T ,

Figure 20 .
Figure 20.Left panel: solid displacement u at x = L as a function of t.Middle panel: fluid pressure p at x = 0 as a function of t.Right panel: discharge velocity v at x = L as a function of t.The applied boundary traction is a step pulse of amplitude P ref at t * = 0.The compressibility parameter c 0 varies in the range [10 −3 , 10 3 ]c 0,ref , where the reference value c 0,ref is set equal to 1.67 • 10 −5 m 2 N −1 as in all the computations illustrated in Section 7.

Figure 21 .
Figure 21.Discharge velocity v at x = L as a function of t.The applied boundary traction is a step pulse of amplitude P ref at t * = 0.The compressibility parameter c 0 varies in the range [10 −3 , 10 3 ]c 0,ref , where the reference value c 0,ref is set equal to 1.67 • 10 −5 m 2 N −1 as in all the computations illustrated in Section 7. The black solid line represents the discharge velocity computed by the model studied in [17].

Figure 22 .
Figure 22.Simulated profiles of E θ (left), E c0 (middle) and E tot as a function of t.The applied boundary traction is a step pulse of amplitude P ref at t * = 0.The compressibility parameter c 0 varies in the range [10 −3 , 10 3 ]c 0,ref , where the reference value c 0,ref is set equal to 1.67 • 10 −5 m 2 N −1 as in all the computations illustrated in Section 7.