The role of structural viscoelasticity in deformable porous media with incompressible constituents: applications in biomechanics

The main goal of this work is to clarify and quantify, by means of mathematical analysis, the role of structural viscoelasticity in the biomechanical response of deformable porous media with incompressible constituents to sudden changes in external applied loads. Models of deformable porous media with incompressible constituents are often utilized to describe the behavior of biological tissues, such as cartilages, bones and engineered tissue scaffolds, where viscoelastic properties may change with age, disease or by design. Here, for the first time, we show that the fluid velocity within the medium could increase tremendously, even up to infinity, should the external applied load experience sudden changes in time and the structural viscoelasticity be too small. In particular, we consider a one-dimensional poro-visco-elastic model for which we derive explicit solutions in the cases where the external applied load is characterized by a step pulse or a trapezoidal pulse in time. By means of dimensional analysis, we identify some dimensionless parameters that can aid the design of structural properties and/or experimental conditions as to ensure that the fluid velocity within the medium remains bounded below a certain given threshold, thereby preventing potential tissue damage. The application to confined compression tests for biological tissues is discussed in detail. Interestingly, the loss of viscoelastic tissue properties has been associated with various disease conditions, such as atherosclerosis, Alzheimer's disease and glaucoma. Thus, the findings of this work may be relevant to many applications in biology and medicine.


Introduction
Fluid flow through deformable porous media is relevant for many applications in biology, medicine and bioengineering. Some important examples include blood flow through tissues in the human body [11,9] and fluid flow inside cartilages, bones and engineered tissue scaffolds [13,21,37,32,10]. The mechanics of biological tissues typically exhibits both elastic and viscoelastic behaviors resulting from the combined action of various components, including elastin, collagen and extracellular matrix [23,24,27,31]. Thus, from the mathematical viewpoint, the study of fluid flows through deformable porous biological structures requires the coupling of poro-elasticity with structural viscoelasticity, leading to poro-visco-elastic models.
The theoretical study of fluid flow through deformable porous media has attracted a lot of attention since the beginning of the last century, initially motivated by applications in geophysics and petroleum engineering. The development of the field started with the work of Terzaghi in 1925 [38], which focused on finding an analytic solution for a one-dimensional (1D) model. However, it was Biot's work in 1941 [6] that set up the framework and ignited the mathematical development for fluid flow through poro-elastic media. To date, several books and articles have been devoted to the mathematical analysis and numerical investigation of poro-elastic models, such as [37,14,12,15,16,40,26,35,36,39,28,29,30,1], with applications ranging from engineering and geophysics to medicine and biology. Recently, our team has developed a theoretical and numerical framework to study both poro-elastic and poro-visco-elastic models, as motivated by biological applications [7]. The study showed that structural viscoelasticity plays a crucial role in determining the regularity requirements for volumetric and boundary forcing terms, as well as for the corresponding solutions. Moreover, in [3] it has been shown that the solution of the fluid-solid mixture (elastic displacement, fluid pressure, and Darcy velocity) is more sensitive to the boundary traction in the elastic case than in the visco-elastic scenario. These theoretical findings are also supported by experimental and clinical evidences showing that changes in tissue viscoelasticity are associated with various pathological conditions, including atherosclerosis [25], osteoporosis [2], renal disease [20] and glaucoma [17].
Interestingly, the study in [7] provided numerical clues that sudden changes in body forces and/or stress boundary conditions may lead to uncontrolled fluid-dynamical responses within the medium in the absence of structural viscoelasticity. This finding led us to formulate a novel hypothesis concerning the causes of damage in biological tissues, namely that abrupt time variations in stress conditions combined with lack of structural viscoelasticity could lead to microstructural damage due to local fluid-dynamical alterations, as illustrated in Fig. 1.
The importance of the coupling between structural mechanics and fluid dynamics in the damage of deformable porous media has been investigated by several authors [24,22,34]. In the present study, we focus on a particular aspect of this coupling and we aim at characterizing and quantifying the influence of structural viscoelasticity on the biomechanical and fluid-dynamical responses to sudden changes in stress conditions. Biomechanical applications are characterized by the fact that tissues have a mass density that is similar to that of water. For this reason, we consider in these pages the case of deformable porous media constituted by incompressible solid and fluid components. To mathematically define this concept, we introduce the following relation between fluid pressure p, variation of fluid content ζ and volumetric dilation ∇ · u, u denoting the solid displacement (1a) ζ = c 0 p + α∇ · u, Figure 1. Illustration of a novel hypothesis on potential causes of microstructure damage in deformable porous media. In the case that (i) the medium is subjected to a time-discontinuous mechanical load and (ii) structural viscoelasticity is reduced or absent, then the fluid velocity within the porous medium will experience a blow-up, possibly leading to microstructural damage.
where c 0 is the storage coefficient and α is the Biot coefficient (see [16]). In the case of incompressible solid and fluid components, we have c 0 = 0 and α = 1 (see [16]) so that (1a) becomes Notice that unlike in standard elasticity theory, incompressibility of each component of a deformable porous medium does not mean that both solid displacement and fluid velocity are divergence-free, rather, that the volumetric deformation of the solid constituent corresponds to the variation of fluid volume per unit volume of porous material, with the convention that ∇ · u is positive for extension while ζ is positive for a "gain" of fluid by the porous solid. Building upon our theoretical and numerical results presented in [7], we devise a 1D problem for which we exhibit an explicit solution where the discharge velocity goes to infinity if the stress boundary condition is not sufficiently smooth in time and the solid component is not viscoelastic. Interestingly, this blow-up in the velocity occurs even in the simple case where the permeability of the medium is assumed to be constant, in comparison to the general case of nonlinear permeability depending on dilation [39,7] or pressure [36]. In addition, we perform a dimensional analysis that allows us to identify the parameters influencing the solution blow-up, thereby opening the path to sensitivity analysis on the system, and providing practical directions on how to control the biomechanical and fluid-dynamical response of the fluid-solid mixture, prevent microstructural damage and, in perspective, aid the experimental design of bioengineered tissues [18].
The article is organized as follows. Poro-visco-elastic models are described in Section 2, along with a summary of the related theoretical results. Section 3 focuses on a special 1D case, for which an explicit solution is derived and its well-posedness is studied in the presence or absence of viscoelasticity. The dimensional analysis of the 1D problem is carried out in Section 4. Solution properties are explored in detail for the particular case of boundary traction with a discontinuity in time, see Section 5, and with a trapezoidal time profile, see Section 6. The application of this analysis to the case of confined compression of biological tissues is discussed in Section 7. Conclusions and perspectives are summarized in Section 8.
2. Poro-Visco-Elastic Models: description and theoretical results Following the same notation as in [7], let Ω ⊂ R d , with d = 1, 2 or 3, be the region of space occupied by a deformable porous medium and let (0, T ) be the observational time interval. In the assumptions of small deformations, full saturation of the mixture, incompressibility of the mixture constituents and negligible inertia, we can write the balance of linear momentum for the mixture and the balance of mass for the fluid component as respectively, where σ is the total stress tensor, ζ is the fluid content, v is the discharge velocity and F and S are volumetric sources of linear momentum and fluid mass, respectively. The system must be completed by constitutive equations for σ, ζ and v. In line with the literature for poro-visco-elastic models in biological applications, we will adopt the following constitutive equations: where u is the displacement of the solid component, p is the Darcy pressure, I is the identity tensor and (w) = (∇w + (∇w) T )/2 is the symmetric part of the gradient of the vector field w. The elastic parameters λ e and µ e and the viscoelastic parameters λ v and µ v are assumed to be given positive constants. We remark that structural viscoelasticity is embodied in the term σ v defined in (2e). If structural viscoelasticity is not present, then σ v = 0 and the set of equations (2a)-(2g) defines a poro-elastic medium. In addition, the permeability K may be a positive constant, as in [6], or may depend nonlinearly on the solution, for example on the structural dilation, namely K = K(∇ · u), as in [39,7]. Finally, notice that (2f) corresponds to (1b) consistently with the assumption of incompressibility of mixture constituents.
Most of the theoretical studies focused on the poro-elastic case without accounting for structural viscoelasticity. In the case of constant permeability, the coupling between the elastic and fluid subproblems is linear and the wellposedness and regularity of solutions have been studied by several authors [40,26,35]. In the case of non-constant permeability, the coupling between the two subproblems becomes nonlinear and only few theoretical results have been obtained. In [36], Showalter utilized monotone operator theory techniques in order to provide well-posedness of solutions in the case where the permeability is a nonlinear function of pressure. To the best of our knowledge, Cao et al in [39] were the first to consider the permeability as a nonlinear function of dilation and provide existence of weak solutions for this nonlinear poro-elastic case. However, the analysis in [39] is performed upon assuming homogeneous boundary conditions for both pressure and elastic displacement, which is often not the case from the viewpoint of applications.
Our recent paper in [7] extends the works mentioned above by considering poro-elastic and poro-visco-elastic models with dilation-dependent permeability, non-zero volumetric sources of mass and momentum and non-homogeneous, mixed Dirichlet-Neumann boundary conditions. More precisely, in [7] we assumed that the boundary of Ω can be decomposed as ∂Ω = Γ N ∪ Γ D , with Γ D = Γ D,p ∪ Γ D,v , where the following boundary conditions are imposed: Our analysis showed that the data time regularity requirements and the smoothness of solutions significantly differ depending on whether the model is poroelastic or poro-visco-elastic. In particular, in the visco-elastic case, if the source of linear momentum F is L 2 in time and space, namely F ∈ L 2 (0, T ; (L 2 (Ω)) d ), and the source of boundary traction g is L 2 in time and H 1/2 on the boundary, namely g ∈ L 2 (0, T ; (H 1/2 (Γ N )) d ), then there exists a weak solution to the system, with elastic displacement u ∈ H 1 in time and space, namely u ∈ H 1 (0, T ; (H 1 (Ω)) d ), and pressure p ∈ L 2 (0, T ; H 1 (Ω)). In comparison, in the poro-elastic case, one needs both F and g to be H 1 in time in order to obtain a solution where both u and p are L 2 in time and H 1 in space. This proves that the system dynamics fundamentally changes as visco-elastic effects fade away. Moreover, the numerical simulations performed in [7] hinted at blowup of fluid energy functional in the poro-elastic case with non-smooth sources F and g. In the following sections, we are going to further elaborate these concepts by means of a 1D poro-visco-elastic model for which we can obtain an explicit solution and assess the effect of structural viscoelasticity on the response of the deformable porous medium to sudden changes in mechanical load. The 1D model is described next.

1D poro-visco-elastic model: description and formulation
Let the space domain Ω be given by the interval (0, L) and the volumetric sources of linear momentum and fluid mass be zero. In this case, the balance of linear momentum for the mixture and the balance of mass for the fluid component simplify to The associated constitutive equations are given by where we have set µ = λ e + 2µ e , η = λ v + 2µ v and we have introduced the symbol σ 0 to denote the contribution of the solid component to the total stress tensor σ. We emphasize that the permeability may be a function of space and dilation, namely We complete the system with the following boundary and initial conditions: where g(t) = −P (t) is the boundary traction. For the ease of reference, the boundary conditions are schematized in Fig. 2. If η = 0 (i.e., no viscoelasticity is present) and under the assumption that the permeability K and the boundary forcing term P are given positive constants, system (6) degenerates into the poro-elastic model studied experimentally and analytically in [37]. Specifically, the boundary and initial conditions (8a)-(8c) correspond to the creep experimental test performed in confined compression, where the boundary condition (8a) represents an impermeable interface at the bottom of the specimen whereas the boundary condition (8b) represents a free-draining permeable piston at the top of the specimen.
System (6a)-(8c) can be rewritten solely in terms of the displacement. Indeed, integration of (6b) with respect to x yields where A (t) is arbitrary. Thus, taking (8a) into account, we obtain that A (t) = 0. From (7c), we then have that On the other hand, using (7b) and (6a), we derive that Therefore the system (6a)-(8c) reduces to the following initial boundary value problem in terms of the sole elastic displacement u: Subsequently, we can recover the solid part of the stress tensor, discharge velocity, pressure and total stress tensor on (0, L) × (0, T ), respectively, as Remark 1. For later reference, we write below the explicit form of the purely elastic problem which corresponds to setting η = 0 in (9): Remark 2. An important quantity associated with the fluid-solid mixture is the fluid power density P (t), defined as From its definition, it follows that P (t) depends on both discharge velocity and permeability.
3.1. Formal derivation of the solution. Let us further assume that the permeability is constant, i.e.
In this case, system (9) is a linear initial boundary value problem, whose solution can be obtained by Fourier series expansion as described below. Case 1. (η > 0). First of all, we homogenize the boundary condition by setting having introduced the auxiliary function where the star symbol denotes convolution. Thus, w (x, t) satisfies the problem where the prime symbol denotes differentiation for functions of a single variable. The associated eigenvalue problem is The eigenvalues and the corresponding eigenfunctions are given by (14) λ n = (2n + 1) 2 π 2 4L 2 and y n (x) = sin (2n + 1)πx 2L for n = 0, 1, . . .
We seek a solution of the form where the coefficients c n (t) of the above series are determined by the differential equation and the initial condition in (13). Using the Fourier Sine series the uniqueness of the Fourier expansion leads to the family of ordinary differential equations whose solution is given by Therefore, we get In conclusion, after performing integration by parts in the convolution term and using the identity ηU (t) + µU (t) = µP (t), the formal solution to problem (9) is given by Case 2. (η = 0). A direct substitution of expression (15) (with η = 0) into problem (11) shows that (16) u is the formal solution of the purely elastic problem. In particular, in the case P (t) = constant, expression (16) coincides with Eq. (8) of [37].

3.2.
Weak solutions and well-posedness. In this section we prove that the formal solutions (15) and (16) indeed solve the visco-elastic problem (9) and the purely elastic problem (11), respectively, in well-defined functional spaces. Let us begin by introducing the functional framework. We consider the real Hilbert space H = L 2 (0, L) equipped with the scalar product and endowed with the induced norm The orthonormal sequence of eigenfunctions If this is the case, the series expansion of v converges in the mean. Next, let be the real Hilbert space equipped with the scalar product and endowed with the induced norm (due to Poincaré's inequality) More generally, the eigenfunction expansion (17) enables us to define a oneparameter family V s (s ∈ R) of spaces: the distribution (17) belongs to V s if and only if its coefficients c n satisfy ∞ n=0 λ s n |c n | 2 < ∞ .
In particular, we have that V −1 ≡ V (the dual of V ). Case 1. (η > 0). Now we introduce the definition of weak solution for the poro-visco-elastic case.
Remark 3. The initial condition (9d) is satisfied by u (t) ∈ V in the pointwise sense of condition (iii), since it is well known that condition (i) implies that u ∈ C 0 ([0, T ] ; V ). The Dirichlet boundary condition (9c) at x = 0 is included in the requirement that u (t) ∈ V , whereas the boundary condition (9b) at x = L is taken into account in the linear form on the right hand side of (18), where v (L) is the pointwise value of the (continuous) test function v ∈ V .
Theorem 3.2. Suppose P (t) ∈ L 2 (0, T ). Then there is a unique weak solution u to problem (9), in the sense of Definition 3.1. Moreover, u is represented by the series expansion (15).
Proof. For sake of exposition, we rewrite u(x, t), given by (15), as [Regularity]. Firstly, we show that u(x, t) has the desired regularity, namely u ∈ H 1 (0, T ; V ) and u (0) = 0. With this goal in mind, we study the convergence of the series expansions in (19) and (21). For any a > 0 Schwarz inequality implies that Upon applying this estimate to (20) and (22) we see that there are positive constants, generically denoted by C, that are independent on n and t and for which it holds that Thus, u, u ∈ L 2 (0, T ; V s ) for every s < 3/2, thereby implying that conditions (i) and (iii) of Def. 3.1 are satisfied.
[Existence]. In order to show that u satisfies condition (ii), it suffices to consider v = y n as test function, since 2 Lλn y n (x) are a basis of V . Each term on the left hand side of (18) can be computed as Adding the above three terms, we obtain that Thus, from (22) and the fact that y n (L) = (−1) n , it follows that the right hand side is equal to −K 0 P (t) y n (L), and this completes the proof of existence.
[Uniqueness] From the linearity of the equation, it suffices to show that P = 0 implies u = 0. Let P = 0 and let u ∈ H 1 (0, T ; V ) be a solution of If we choose v = u (t) as a test function, then we obtain that Integrating with respect to time and using the initial condition u (0) = 0, we get (23), the fact that |y n (x)| ≤ 1, and the standard Weierstrass Test, it follows that the series (19) converges absolutely and uniformly in the space-time rectangle Q = [0, L]×[0, T ] and its sum is continuous, i.e. u (x, t) ∈ C 0 (Q). In a similar way, under the stronger assumption P (t) ∈ L ∞ (0, T ), estimate (24) implies |u n (t)| ≤ Cλ −1 n P L ∞ (0,T ) so that ∂u ∂t (x, t) ∈ C 0 (Q) also holds. Thus, in the physically realistic case of bounded boundary traction P (t), discharge velocity and fluid power density are continuous and bounded both in space and time, and this is true even though the traction has jump discontinuities (e.g., it is a step pulse).
Hence, in this case, u ∈ L 2 (0, T ; V s ) for every s < 3/2 and u ∈ L 2 (0, T ; V s ) for every s < −1/2. In particular, when the boundary traction is bounded, we have u (x, t) ∈ C 0 (Q) whereas the discharge velocity is a distribution. In conclusion, we state the following definition and theorem (without proof): where the brackets ·, · denote pairing between V and V ; (iii): u (0) = 0.

Dimensionless Problem
The goal of this section is to rewrite problem (9) in dimensionless form so that we can identify combinations of geometrical and physical parameters that most influence the solution properties. Dimensional analysis relies on the choice of a set of characteristic values that can be used to scale all the problem variables. Let us use the hat symbol to indicate dimensionless (or scaled) variables and the square brackets to indicate the characteristic value of that quantity. Then, for the problem at hand we would write: It is important to emphasize that there is no trivial choice for the characteristic values and, in general, this choice is not unique. In this particular case, though, we will leverage our knowledge of the forcing terms and the explicit formulas we obtained for the solution to guide us in the choice of some of these values. Since the problem is driven by the boundary condition on the traction with the given function P (t), we set whereas the expressions for v(x, t) in (10b) and P(t) in (12) suggest to choose (30) [ respectively. Using the above scalings, we obtain the following dimensionless problem: Remark 5. As already mentioned above, the choice for the characteristic values is not unique. In this regard, it is worth noticing that our choice for [v], namely K 0 P ref /L, differs from that utilized in [37], which was K 0 µ/L instead. This difference actually follows from a different choice for the pressure scaling, namely P ref in our case, as opposed to µ in [37]. Indeed, the two scalings are equivalent if P ref and µ are of the same order of magnitude. In general, though, the scaling K 0 µ/L represents the characteristic velocity of wave propagation within the medium, whereas the scaling K 0 P ref /L represents the characteristic velocity induced by an external load of magnitude P ref . Since our analysis aims at assessing the influence of the external load on the biomechanical response of the medium, we chose [P ] = P ref . However, the analysis could be easily adapted to other choices, should the scope of the investigation differ from that of this paper.

The caseP t = step pulse
Let the boundary tractionP t be the dimensionless unit step starting at t = 0 defined as In this case, the dimensionless solid displacement (33), discharge velocity (34) 0t P (t ) 1 Figure 3. Schematic representation of the dimensionless step pulseP (t ) defined in (36). Here, the signal is discontinuous at the switch on time.
and power density (35) (hereon denoted with the subscriptη for later reference) become, respectively, The solution in the purely elastic case can be obtained by settingη = 0 in (37).

Remark 6.
If the unit step is shifted att = α > 0, i.e.P t = H t − α , then the solution (33) iŝ The space-time behavior ofûη,vη andPη is reported in Figures 4, 5 and 6, respectively, in the case ofη = 0 andη = 1. We remark thatvη andPη are reported in logarithmic scale to highlight the presence of a blow-up atx = 1 andt = 0 in the purely elastic caseη = 0. Indeed, by settingη =t = 0 in (37c), we can verify that the power density satisfieŝ Proceeding analogously in the case of the discharge velocity, the Fourier expansion in the purely elastic case is given bŷ  From the physical viewpoint, this means that, at the switch on time of the driving term, here set att = 0, the discharge velocity contains high-frequency components whose superposition exhibits the following behaviors depending on x: (i) it relaxes towards zero asx → 0; (ii) it tends to infinity asx → 1, thereby exhibiting a blow-up; (iii) it looks like an amorphous blob for 0 <x < 1.
On the other hand, if viscoelastic effects are present, then the dimensionless discharge velocity att = 0 is given bŷ Here, the n−th coefficient is decreasing as 1/n 2 and the discharge velocity is continuous att = 0. These concepts are illustrated in Fig. 7, which shows the spatial distribution ofvη att = 0 for several values ofη in the interval [0, 1]. The log 10 plot ofvη is reported in the figure to highlight the blow-up of the discharge velocity in the purely elastic case atx = 1.
In order to further investigate this blow-up and its dependence on the structural viscoelasticity, we observe that the maximum value ofvη is attained at x = 1 andt = 0 and can be written as a function ofη as follows: v max (η) = max . Dimensionless discharge velocityvη(x, 0) as a function ofx. In order to highlight the velocity blow-up atx = 1, t = 0, in the caseη = 0, the log 10 plot of |v η | is plotted.
The above series may be summed ( [19], formula no. 1.4212) and the final result is Similarly, the dimensionless power density (37c) is decreasing in time and its maximum is attained att = 0: The behaviors ofv max (η) andP max (η) with respect toη are visualized in Fig. 8.
Clearly, the dimensionless parameterη is the sole determinant of the blowup inv max (η) andP max (η), with smaller values ofη leading to higher values ofv max (η) andP max (η). However, it is important to recall thatη does not depend solely on the viscoelasticity coefficient. The definitionη = ηK 0 /L 2 implies that also small values of the permeability K 0 and/or large values of the domain dimension L would reduce the value ofη. Precisely, going back to dimensional units, we can use (30) and (31) to obtain the following expressions: We can see that the magnitude P ref of the driving term appears as a multiplying constant in the expressions above, meaning that larger values of P ref

The caseP (t) = trapezoidal pulse
Let us now consider the case of a driving term given by a trapezoidal pulse, where the signal switch on and switch off are characterized by linear ramps. Thus, letP t be the dimensionless trapezoidal pulse of unit amplitude and letε,τ > 0. Here we consider the following expression forP t : whereε is the pulse rise/fall time andτ is the pulse duration, as depicted in Fig. 10. Thus, the trapezoidal pulse defined in (41) reduces to a rectangular step pulse asε → 0.
0t P (t) 1ετ +ετ + 2ε Figure 10. Schematic representation of the dimensionless trapezoidal pulseP (t) defined in (41). Here, the signal switch on and switch off are characterized by linear ramps.
Let us now compute the dimensionless discharge velocity resulting from the application of the trapezoidal pulse at the boundary, henceforth denoted bŷ Vη x,t to distinguish it from the discharge velocityvη x,t obtained for the step pulse. To this end, we observe that the trapezoidal pulse is actually the linear combination of four ramp pulses of unit slope starting at timest = 0,ε, ε +τ , 2ε +τ : where the function H(t) was defined in (36). In addition, we can see that u x,t = t 0ûη (x, s) ds whereûη x,t is given by (37a), solves problem (32) withP t =tH t (the ramp pulse of unit slope starting at timet = 0). Thus, using (10b), we have that the discharge velocity is given by −ûη x,t , so that, applying the linear superposition principle,Vη x,t is the linear combination of the velocities corresponding to the various components of the pulse, namely: An illustration of the typical form ofVη x,t forη ≥ 0 andε > 0 is reported in Fig. 11. The maximum possible discharge velocity occurs atx = 1,t =ε and can be written asV The behavior ofV max with respect toη andε is reported in Fig. 12. Interestingly, for allη > 0 andε ≥ 0 we havê V max (η,ε) ≤V max (η, 0) andV max (η, 0) =v max (η) since the trapezoidal pulse reduces to a rectangular pulse asε → 0. Thus, no blow-up takes place in the viscoelastic case when the pulse is rectangular. Similarly, for allη ≥ 0 andε > 0 it iŝ hence no blow-up takes place even in the purely elastic case when the pulse is trapezoidal.

Application in biomechanics: role of structural viscoelasticity in confined compression tests for biological tissues
In this section, we utilize the mathematical analysis developed above to study some interesting features of confined compression tests, which are often utilized in biomechanics to characterize the properties of biological tissues. A schematic of the confined compression experimental setting is depicted in Fig. 13, where a compressive load is applied at the chamber top surface while the bottom surface is maintained fixed. Due to confinement, deformation occurs only in the x direction. The analogy with the 1D setting depicted in Fig. 2 is straightforward, with a note of caution that 1D solutions are most representative of the biomechanical status of the material along vertical lines towards the center of the chamber.
The 1D model described in Section 3 allows us to generalize the mathematical analysis carried out in [37] by quantifying the effect of structural viscoelasticity on the tissue response to sudden changes in external pressure during confined compression experiments. In this perspective, L represents the thickness of the undeformed tissue sample and P ref represents the magnitude of the total compressive stress applied at x = L. It is important to notice that the confined compression creep experiment described in [37] is characterized by the same initial and boundary conditions as those given in (8). However, the mathematical description adopted in [37] (based on the theory developed in [23]), does not account for structural viscoelasticity, which amounts to setting η = 0 in (7a). Table 1 summarizes the parameter values corresponding to the experiments on articular cartilage reported in [37].  Table 1. Numerical values of model parameters in the confined compression experiment for articular cartilage reported in [37].
Characteristic velocities for the confined compression experiments reported in [37] are of the order of 0.35 µm/s. However, our analysis showed that velocities in the sample may attain much higher values, particularly near the permeable piston, should viscoelasticity be too low and/or the time profile of the load protocol be too sharp. In order to make these concepts more specific, let us assume that: (i) the load protocol for the confined compression experiment corresponds to the trapezoidal pulse described in Section 6; and that (ii) we should ensure that discharge velocities always remain below the threshold value V th = 0.35 µm/s in order to preserve the sample from microstructural damage. This means thatη andε should be chosen as to guarantee that whereV max (η,ε) can be computed using the expression (44).  Figure 14. Comparison between the dimensionless maximum velocityV max (η,ε) obtained in the case of trapezoidal pulse using expression (44) and the threshold velocityV th = 16.3 typical of confined compression experiments [37].
To better visualize the various regions of interest in the (η,ˆ )−plane, we reported the colormap of the differenceV th −V max as a function ofη andε in Fig. 15, where we clearly marked the curve in the parameter space for whicĥ V th =V max . These results show that, for example, choosingε ≥ 5 · 10 −3 would lead to fluid velocities below the recommended threshold regardless of the presence of structural viscoelasticity. On the other hand, should the particular experimental conditions constrain us to enforce a rectangular pulse corresponding toε = 0, we can still maintain the maximum velocity below the threshold by modifying the structural properties of the sample to have, for example,η ≥ 4 · 10 −3 .
Interestingly, expression (29) definesη = (K 0 /L 2 )η, meaning thatη can be increased by increasing the viscoelastic parameter η, increasing the permeability parameter K 0 , or decreasing the sample thickness L. In the same spirit as [8], a typical value of the viscoelastic parameter η for biological tissues can be estimated by setting Colormap ofV th −V max (η,ε) Figure 15. Colormap of the differenceV th −V max as a function ofη andε. As in Fig. 14,V max (η,ε) is obtained using expression (44) andV th = 16.3 is the typical velocity of confined compression experiments [37]. The curve in the parameter space for whichV th =V max is reported in a thick black mark.
let us set ρ = 1000 Kg m −3 . Substituting (46a) and (46b) in (29) we obtain Replacing the parameter values of Table 1 into (46c) we obtainη = 1.15 · 10 −8 , which is 5 orders of magnitude smaller than what is needed to attain fluid velocities below the recommended threshold in the case of a rectangular pulse in the pressure load. Thus, in order to control the maximum fluid velocity within the tissue sample, one may decide to: (i) vary the properties of the biological sample by manipulating ρ, µ, η, K 0 or L; and/or (ii) vary the time profile of the load experimental protocol. The choice will depend on the particular conditions and constraints of the experiment at hand.

Conclusions and future perspectives
The exact solutions obtained for the 1D models considered in this article allowed us to clearly identify a blow-up in the solution of certain poroelastic problems. The analysis allowed us to identify the main factors that give rise to the blow-up, namely the absence of structural viscoelasticity and the time-discontinuity of the boundary source of traction. It is very important to emphasize that even a small viscoelastic contributions, namelyη 1, will prevent the blow-up. Similarly, a continuous time profile of the applied boundary load will prevent blow-up even if the structure is purely elastic.
These findings actually provide an evidence of a blow-up that was hypothesized by the theoretical work presented in [7]. Interestingly, the a priori estimates derived in [7] were not sufficient to bound the power density if the data were not regular enough. The 1D examples considered in this article show that indeed it is not possible to bound the power density if the structure is purely elastic and the boundary forcing term is discontinuous in time. Moreover, we have shown that this blow-up occurs even in the simple case of constant permeability.
In real situations, we will never see the fluid velocity spiking to infinity as predicted by the mathematical blow-up, since something will break first! For example, if the poroelastic model represents a biological tissue perfused by blood flowing in capillaries, as the maximum velocity becomes too high, capillaries will break letting blood out. Thus, from a practical perspective, it is crucial to identify parameters that can control the maximum value of the fluid velocity within a deformable porous medium in order to avoid microstructure damage. Our analysis showed that the maximum fluid velocity can be limited by: (i) decreasing P ref , thereby reducing the load on the medium; (ii) increasinĝ ε, thereby smoothing the load action on the medium; (iii) increasing η, thereby changing the structural properties of the solid component of the medium; (iv) decreasing K 0 , thereby changing the pore size, geometry and/or fluid viscosity; (v) decreasing L, thereby changing the domain geometry. We have explored these concepts more concretely by applying them to the confined compression tests for biological tissues.
It is worth noting that the extreme sensitivity of biological tissues to the active role of viscoelasticity predicted by our theoretical analysis seems to agree with the conclusions of [33] where experimental data based on magnetic resonance elastography show that viscoelasticity of the brain is a result of structural alteration occurring in the course of physiological aging, this suggesting that cerebral viscoelasticity may provide a sensitive marker for a variety of neurological diseases such as normal pressure hydrocephalus, Alzheimer's disease, or Multiple Sclerosis.
Further extensions of the current work may include: (1) the study of porous deformable media with compressible components. This case mathematically corresponds to considering (1a) instead of (1b), and is of interest in the wider context of applications of the poro-visco-elastic model to problems in geomechanics (cf. the original contribution by Biot in [6] and the more recent works [28,29,30] and [5,4], devoted to computational analysis and theoretical investigations, respectively); (2) the investigation of time discontinuities in the volumetric sources of linear momentum, which were identified by the analysis in [7] as additional possible causes of blow-up. This case is actually extremely relevant in situations where the gravitational acceleration varies abruptly, such as during space flight take off and landing.
Finally, it is worth emphasizing that the 1D analysis presented in this article has made available a series of testable conditions leading to blow-up, and consequent microstructural damage, that could indeed be verified in laboratory experiments. We see the design and implementation of such experiments as the most interesting development of the present work.