1. 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 [9, 11] and fluid flow inside cartilages, bones and engineered tissue scaffolds [10, 13, 21, 31, 37]. 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 [22, 23, 26, 30]. 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 [5] 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 [8, 12, 14, 15, 16, 25, 27, 28, 29, 34, 35, 36, 37, 40], 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 [6]. 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 [2] 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 [24], osteoporosis [1], renal disease [20] and glaucoma [17].

Interestingly, the study in [6] 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 [39, 23, 33]. 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 $\zeta$ and volumetric dilation $\nabla \cdot \mathbf{u}$, $\mathbf{u}$ denoting the solid displacement

$\begin{eqnarray}\zeta = c_0 p + \alpha \nabla \cdot \mathbf{u}, \end{eqnarray}$

(1a)

where $c_0$ is the storage coefficient and $\alpha$ is the Biot coefficient (see [16]). In the case of incompressible solid and fluid components, we have $c_0=0$ and $\alpha=1$ (see [16]) so that (1a) becomes

$\begin{eqnarray}\zeta = \nabla \cdot \mathbf{u}. \end{eqnarray}$

(1b)

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 $\nabla \cdot \mathbf{u}$ is positive for extension while $\zeta$ is positive for a "gain" of fluid by the porous solid.

Building upon our theoretical and numerical results presented in [6], 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 [6, 8] 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 [6], let $\Omega \subset \mathbb 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

$\begin{align} \nabla \cdot \boldsymbol{\sigma} + \mathbf{F} = \mathbf 0 && \mbox{in}\; \Omega\times (0,T)\end{align}$

(2a)

$\begin{align} \dfrac{\partial \zeta}{\partial t} + \nabla \cdot \mathbf v = S && \mbox{in}\; \Omega\times (0,T) \end{align}$

(2b)

respectively, where $\boldsymbol{\sigma}$ is the total stress tensor, $\zeta$ is the fluid content, $\mathbf v$ is the discharge velocity and $\mathbf{F}$ and $S$ are volumetric sources of linear momentum and fluid mass, respectively. The system must be completed by constitutive equations for $\boldsymbol{\sigma}$, $\zeta$ and $\mathbf v$. In line with the literature for poro-visco-elastic models in biological applications, we will adopt the following constitutive equations:

$\begin{align} \boldsymbol{\sigma} &= \boldsymbol{\sigma}_{e} + \boldsymbol{\sigma}_{v} - p \mathbf{I}\end{align}$

(2c)

$\begin{align} \boldsymbol{\sigma}_{e} &= \lambda_e (\nabla \cdot \mathbf u )~ \mathbf I + 2\mu_e \boldsymbol{\epsilon}(\mathbf u)\end{align}$

(2d)

$\begin{align} \boldsymbol{\sigma}_{v} &= \lambda_v \left(\nabla \cdot \dfrac{\partial \mathbf u}{\partial t} \right)~ \mathbf I + 2\mu_v \boldsymbol{\epsilon}\left(\dfrac{\partial \mathbf u}{\partial t} \right)\end{align}$

(2e)

$\begin{align} \zeta & = \nabla \cdot \mathbf u \end{align}$

(2f)

$\begin{align} \mathbf v &= -K \nabla p \end{align}$

(2g)

where $\mathbf u$ is the displacement of the solid component, $p$ is the Darcy pressure, $\mathbf{I}$ is the identity tensor and $\boldsymbol \epsilon (\mathbf w)=(\nabla \mathbf w + (\nabla \mathbf w)^T)/2$ is the symmetric part of the gradient of the vector field $\mathbf w$. The elastic parameters $\lambda_e$ and $\mu_e$ and the viscoelastic parameters $\lambda_v$ and $\mu_v$ are assumed to be given positive constants. We remark that structural viscoelasticity is embodied in the term $\boldsymbol{\sigma}_v$ defined in (2e). If structural viscoelasticity is not present, then $\boldsymbol{\sigma}_v = \mathbf 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 [5], or may depend nonlinearly on the solution, for example on the structural dilation, namely $K=K(\nabla \cdot \mathbf u)$, as in [6, 8]. 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 well-posedness and regularity of solutions have been studied by several authors [25, 35, 40]. 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 [8] 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 [8] 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 [6] 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 [6] we assumed that the boundary of $\Omega$ can be decomposed as $\partial \Omega = \Gamma_N \cup \Gamma_D$, with $\Gamma_D = \Gamma_{D,p}\cup \Gamma_{D,v}$, where the following boundary conditions are imposed:

$\begin{align} \boldsymbol\sigma \cdot \mathbf n = \mathbf g, \quad \mathbf v \cdot \mathbf n = 0 && \mbox{on}\; \Gamma_N\times (0,T) \end{align}$

(3)

$\begin{align} \mathbf u = \mathbf 0, \quad p = 0 && \mbox{on}\; \Gamma_{D_p}\times (0,T) \end{align}$

(4)

$\begin{align} \mathbf u = \mathbf 0, \quad \mathbf v \cdot \mathbf n = \psi && \mbox{on}\; \Gamma_{D_v}\times (0,T).\end{align}$

(5)

Our analysis showed that the data time regularity requirements and the smoothness of solutions significantly differ depending on whether the model is poro-elastic or poro-visco-elastic. In particular, in the visco-elastic case, if the source of linear momentum $\mathbf F$ is $L^2$ in time and space, namely $\mathbf F\in L^2(0,T;(L^2(\Omega))^d)$, and the source of boundary traction $\mathbf g$ is $L^2$ in time and $H^{1/2}$ on the boundary, namely $\mathbf g \in L^2(0,T;(H^{1/2}(\Gamma_N))^d)$, then there exists a weak solution to the system, with elastic displacement $\mathbf u \in H^1$ in time and space, namely $\mathbf u \in H^1(0,T;(H^1(\Omega))^d)$, and pressure $p \in L^2(0,T;H^1(\Omega))$. In comparison, in the poro-elastic case, one needs both $\mathbf F$ and $\mathbf g$ to be $H^1$ in time in order to obtain a solution where both $\mathbf 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 [6] hinted at blow-up of fluid energy functional in the poro-elastic case with non-smooth sources $\mathbf F$ and $\mathbf 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.

---

3. 1D poro-visco-elastic model: Description and formulation

Let the space domain $\Omega$ 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

$\begin{align} \dfrac{\partial \sigma }{\partial x}=0& & & \mbox{in}\;(0,L)\times (0,T) \end{align}$

(6a)

$\begin{align} \dfrac{\partial ^{2}u}{\partial t\partial x}+\dfrac{\partial v}{\partial x} =0& & & \mbox{in}\;(0,L)\times (0,T). \end{align}$

(6b)

The associated constitutive equations are given by

$\begin{align} \sigma _{0}=\mu \dfrac{\partial u}{\partial x}+\eta \dfrac{\partial ^{2}u}{ \partial t\partial x}& & & \mbox{in}\;(0,L)\times (0,T) \end{align}$

(7a)

$\begin{align} \sigma =\sigma _{0}-p& & & \mbox{in}\;(0,L)\times (0,T) \end{align}$

(7b)

$\begin{align} v=-K\dfrac{\partial p}{\partial x}& & & \mbox{in}\;(0,L)\times (0,T) \end{align}$

(7c)

where we have set $\mu =\lambda _{e}+2\mu _{e}$, $\eta =\lambda _{v}+2\mu _{v}$ and we have introduced the symbol $\sigma_0$ to denote the contribution of the solid component to the total stress tensor $\sigma$. We emphasize that the permeability may be a function of space and dilation, namely

$\begin{eqnarray*} K=K\left(x,\frac{\partial u}{\partial x}\right). \end{eqnarray*}$

We complete the system with the following boundary and initial conditions:

$\begin{align} u\left( 0,t\right) =v\left( 0,t\right) =0& & & \mbox{for}\;0<t<T \end{align}$

(8a)

$\begin{align} p\left( L,t\right) =0\ ;\ \sigma \left( L,t\right) =-P\left( t\right) & & & \mbox{for}\;0<t<T \end{align}$

(8b)

$\begin{align} u\left( x,0\right) =0& & & \mbox{for}\;0<x<L \end{align}$

(8c)

where $g(t) = -P(t)$ is the boundary traction. For the ease of reference, the boundary conditions are schematized in Fig. 2.

If $\eta =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

$\begin{eqnarray*} \dfrac{\partial u}{\partial t}(x,t)+v(x,t)=A(t) \end{eqnarray*}$

where $A\left( t\right)$ is arbitrary. Thus, taking (8a) into account, we obtain that $A\left( t\right) =0$. From (7c), we then have that

$\begin{eqnarray*} \dfrac{\partial u}{\partial t}-K\dfrac{\partial p}{\partial x}=0\quad \text{ in}\ (0,L)\times (0,T)~. \end{eqnarray*}$

On the other hand, using (7b) and (6a), we derive that

$\begin{eqnarray*} \dfrac{\partial p}{\partial x}=\dfrac{\partial \sigma _{0}}{\partial x}~. \end{eqnarray*}$

Therefore the system (6a)-(8c) reduces to the following initial boundary value problem in terms of the sole elastic displacement $u$:

$\begin{align} \dfrac{\partial u}{\partial t}-K\mu \dfrac{\partial ^{2}u}{\partial x^{2}} -K\eta \dfrac{\partial ^{3}u}{\partial t\partial x^{2}}=0 & & & \mbox{in}\ (0,L)\times (0,T) \end{align}$

(9a)

$\begin{align} \mu \dfrac{\partial u}{\partial x}\left( L,t\right) +\eta \dfrac{\partial ^{2}u}{\partial t\partial x}\left( L,t\right) =-P\left( t\right) & & & \mbox{for}\;0<t<T \end{align}$

(9b)

$\begin{align} u\left( 0,t\right) =0 & & & \mbox{for}\;0<t<T \end{align}$

(9c)

$\begin{align} u\left( x,0\right) =0 & & & \mbox{for}\;0<x<L ~.\end{align}$

(9d)

Subsequently, we can recover the solid part of the stress tensor, discharge velocity, pressure and total stress tensor on $(0,L)\times (0,T)$, respectively, as

$\begin{eqnarray} \sigma _{0} &=&\mu \dfrac{\partial u}{\partial x}+\eta \dfrac{\partial ^{2}u }{\partial t\partial x} \end{eqnarray}$

(10a)

$\begin{eqnarray} v &=&-\dfrac{\partial u}{\partial t}=-K\dfrac{\partial \sigma _{0}}{\partial x} \end{eqnarray}$

(10b)

$\begin{eqnarray} p &=&\sigma _{0}+P\left( t\right) \end{eqnarray}$

(10c)

$\begin{eqnarray} \sigma &=&-P\left( t\right) ~. \end{eqnarray}$

(10d)

Remark 1.  For later reference, we write below the explicit form of the purely elastic problem which corresponds to setting $\eta =0$ in (9):

$\begin{align} \dfrac{\partial u}{\partial t}-K\mu \dfrac{\partial ^{2}u}{\partial x^{2}} =0 & & & \mbox{in}\ (0,L)\times (0,T) \end{align}$

(11a)

$\begin{align} \mu \dfrac{\partial u}{\partial x}\left( L,t\right) =-P\left( t\right) & & & \mbox{for}\;0<t<T \end{align}$

(11b)

$\begin{align} u\left( 0,t\right) =0 & & & \mbox{for}\;0<t<T \end{align}$

(11c)

$\begin{align} u\left( x,0\right) =0 & & & \mbox{for}\;0<x<L~. \end{align}$

(11d)

Remark 2.  An important quantity associated with the fluid-solid mixture is the fluid power density $\mathcal{P}\left( t\right)$, defined as

$\begin{eqnarray} \mathcal{P}\left( t\right) =\int_{0}^{L}\frac{1}{K}\left\vert v\left( x,t\right) \right\vert ^{2}dx. \end{eqnarray}$

(12)

From its definition, it follows that $\mathcal{P}\left( t\right)$ 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.

$\begin{eqnarray*} K=K\left(x,\frac{\partial u}{\partial x}\right)\equiv K_{0}=\mathrm{constant} >0. \end{eqnarray*}$

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. ($\eta >0$). First of all, we homogenize the boundary condition by setting

$\begin{eqnarray*} w\left( x,t\right) =u\left( x,t\right) +\frac{U\left( t\right) }{\mu }x \end{eqnarray*}$

having introduced the auxiliary function

$\begin{eqnarray*} U\left( t\right) =\frac{\mu }{\eta }\int_{0}^{t}\exp \left( -\frac{\mu }{ \eta }\left( t-s\right) \right) P\left( s\right) ds=\frac{\mu }{\eta }\exp \left( -\frac{\mu }{\eta }t\right) \ast P\left( t\right) \end{eqnarray*}$

where the star symbol denotes convolution. Thus, $w\left( x,t\right)$ satisfies the problem

$\begin{eqnarray} \begin{array}{rll} \dfrac{\partial w}{\partial t}-K_{0}\mu \dfrac{\partial ^{2}w}{\partial x^{2} }-K_{0}\eta \dfrac{\partial ^{3}w}{\partial t\partial x^{2}}=\dfrac{x}{\mu } U^{\prime }\left( t\right) & \qquad & \mbox{in}\ (0,L)\times (0,T) \\ \mu \dfrac{\partial w}{\partial x}\left( L,t\right) +\eta \dfrac{\partial ^{2}w}{\partial t\partial x}\left( L,t\right) =0 & \qquad & \mbox{for} \;0<t<T \\ w\left( 0,t\right) =0 & \qquad & \mbox{for}\;0<t<T \\ w\left( x,0\right) =0 & \qquad & \mbox{for}\;0<x<L \end{array} \end{eqnarray}$

(13)

where the prime symbol denotes differentiation for functions of a single variable. The associated eigenvalue problem is

$\begin{eqnarray*} Find~y=y\left( x\right) \textit{,}~ 0<x<~L \textit{,}~such~that \end{eqnarray*}$

$\begin{eqnarray*} y^{\prime \prime }+\lambda y=0,\qquad y\left( 0\right) =y^{\prime }\left( L\right) =0~. \end{eqnarray*}$

The eigenvalues and the corresponding eigenfunctions are given by

$\begin{eqnarray} \lambda _{n}=\frac{(2n+1)^{2}\pi ^{2}}{4L^{2}}\quad \mbox{and}\quad y_{n}(x)=\sin \frac{(2n+1)\pi x}{2L} \quad \mbox{for}\quad n=0,1,\dots \end{eqnarray}$

(14)

We seek a solution of the form

$\begin{eqnarray*} w\left( x,t\right) =\sum_{n=0}^{\infty }c_{n}\left( t\right) y_{n}\left( x\right) \end{eqnarray*}$

where the coefficients $c_{n}\left( t\right)$ of the above series are determined by the differential equation and the initial condition in (13 ). Using the Fourier Sine series expansion

$\begin{eqnarray*} x=\frac{2}{L}\sum_{n=0}^{\infty }\frac{\left( -1\right) ^{n}}{\lambda _{n}} y_{n}\left( x\right) ,\qquad 0\leq x\leq L \end{eqnarray*}$

the uniqueness of the Fourier expansion leads to the family of ordinary differential equations

$\begin{eqnarray*} \left( 1+K_{0}\eta \lambda _{n}\right) c_{n}^{\prime }\left( t\right) +K_{0}\mu \lambda _{n}c_{n}\left( t\right) =\frac{2\left( -1\right) ^{n}}{ \mu L\lambda _{n}}U^{\prime }\left( t\right) ,\qquad c_{n}\left( 0\right) =0 \end{eqnarray*}$

whose solution is given by

$\begin{eqnarray*} c_{n}\left( t\right) =\frac{2\left( -1\right) ^{n}}{\mu L\lambda _{n}\left( 1+K_{0}\eta \lambda _{n}\right) }\exp \left( -\frac{K_{0}\mu \lambda _{n}}{ 1+K_{0}\eta \lambda _{n}}t\right) \ast U^{\prime }\left( t\right)~. \end{eqnarray*}$

Therefore, we get

$\begin{eqnarray*} u\left( x,t\right) =-\frac{U\left( t\right) }{\mu }x+\frac{2}{\mu L} \sum_{n=0}^{\infty }\frac{\left( -1\right) ^{n}y_{n}\left( x\right) }{ \lambda _{n}\left( 1+K_{0}\eta \lambda _{n}\right) }\exp \left( -\frac{ K_{0}\mu \lambda _{n}}{1+K_{0}\eta \lambda _{n}}t\right) \ast U^{\prime }\left( t\right) . \end{eqnarray*}$

In conclusion, after performing integration by parts in the convolution term and using the identity $\eta U^{\prime }\left( t\right) +\mu U\left( t\right) =\mu P\left( t\right)$, the formal solution to problem (9) is given by

$\begin{eqnarray} u\left( x,t\right) =-\frac{2K_{0}}{L}\sum_{n=0}^{\infty }\frac{\left( -1\right) ^{n}y_{n}\left( x\right) }{1+K_{0}\eta \lambda _{n}}\exp \left( - \frac{K_{0}\mu \lambda _{n}}{1+K_{0}\eta \lambda _{n}}t\right) \ast P\left( t\right) ~. \end{eqnarray}$

(15)

Case 2. ($\eta =0$). A direct substitution of expression (15) (with $\eta=0$) into problem (11) shows that

$\begin{eqnarray} u\left( x,t\right) =-\frac{2K_{0}}{L}\sum_{n=0}^{\infty } \left(-1\right) ^{n}y_{n}\left( x\right) \exp \left(-K_{0}\mu \lambda _{n} t\right) \ast P\left( t\right) \end{eqnarray}$

(16)

is the formal solution of the purely elastic problem. In particular, in the case $P(t)=\mathrm{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}\left( 0,L\right)$ equipped with the scalar product

$\begin{eqnarray*} (f,g)_{H}=\int_{0}^{L}f(x)g(x)~dx\qquad \forall f,g\in H, \end{eqnarray*}$

and endowed with the induced norm

$\begin{eqnarray*} \left\Vert f\right\Vert _{H}=\sqrt{(f,f)_{H}}. \end{eqnarray*}$

The orthonormal sequence of eigenfunctions $\left\{ \sqrt{\frac{2}{L}} y_{n}\left( x\right) \right\}$ forms a Hilbert space basis for $H$. A distribution $v$ in $\left( 0,L\right)$ belongs to $H$ if and only if it has a series expansion (converging in $\mathcal{D}^{\prime }\left( 0,L\right)$)

$\begin{eqnarray} v\left( x\right) =\sum_{n=0}^{\infty }c_{n}y_{n}\left( x\right), \end{eqnarray}$

(17)

with coefficients $c_{n}$ satisfying

$\begin{eqnarray*} \sum_{n=0}^{\infty }\left\vert c_{n}\right\vert ^{2}<\infty. \end{eqnarray*}$

If this is the case, the series expansion of $v$ converges in the mean. Next, let

$\begin{eqnarray*} V=\left\{ v\in H:v^{\prime }\in H,v\left( 0\right) =0\right\} \end{eqnarray*}$

be the real Hilbert space equipped with the scalar product

$\begin{eqnarray*} \left( v,w\right) _{V}=\left( v^{\prime },w^{\prime }\right) _{H}\qquad \forall v,w\in V, \end{eqnarray*}$

and endowed with the induced norm (due to Poincaré's inequality)

$\begin{eqnarray*} \left\Vert v\right\Vert _{V}=\left\Vert v^{\prime }\right\Vert _{H}. \end{eqnarray*}$

Sobolev's Embedding Theorem ensures that $V\subset C^{0}\left( \left[ 0,L\right] \right)$, so that the boundary value $v\left( 0\right) =0$ for $v\in V$ is assumed in a strong sense. The sequence of functions $\left\{ \sqrt{\frac{2}{L\lambda _{n}}}y_{n}\left( x\right) \right\}$ forms a Hilbert space basis for $V$. A function $v\in H$ belongs to $V$ if and only if the Fourier coefficients of the series expansion (17) satisfy

$\begin{eqnarray*} \sum_{n=0}^{\infty }\lambda _{n}\left\vert c_{n}\right\vert ^{2}<\infty. \end{eqnarray*}$

More generally, the eigenfunction expansion (17) enables us to define a one-parameter family $V^{s}$ ($s\in \mathbb{R}$) of spaces: the distribution (17) belongs to $V^{s}$ if and only if its coefficients $c_{n}$ satisfy

$\begin{eqnarray*} \sum_{n=0}^{\infty }\lambda _{n}^{s}\left\vert c_{n}\right\vert ^{2} <\infty ~. \end{eqnarray*}$

In particular, we have that $V^{-1}\equiv V^{\prime }$ (the dual of $V$).

Case 1. ($\eta >0$). Now we introduce the definition of weak solution for the poro-visco-elastic case.

Definition 3.1.  A function $u:\left[ 0,T\right] \rightarrow V$ is a weak solution to problem (9) if:

(ⅰ) $u\in H^{1}\left( 0,T;V\right)$, implying that $u,u^{\prime }\in L^{2}\left( 0,T;V\right)$;

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

$\begin{eqnarray} \left( u^{\prime }\left( t\right) ,v\right) _{H}+K_{0}\left( \mu u\left( t\right) +\eta u^{\prime }\left( t\right) ,v\right) _{V}=-K_{0}P\left( t\right) v\left( L\right) ; \end{eqnarray}$

(18)

(ⅲ) $u\left( 0\right) =0,$

where we used the notation $u\left( t\right) =u\left( \cdot ,t\right)$ and $u^{\prime }\left( t\right) = \dfrac{\partial u}{\partial t}\left( \cdot ,t\right)$.

Remark 3.  The initial condition (9d) is satisfied by $u\left( t\right) \in V$ in the pointwise sense of condition (ⅲ), since it is well known that condition (ⅰ) implies that $u\in C^{0}\left( \left[ 0,T\right] ;V\right)$. The Dirichlet boundary condition (9c) at $x=0$ is included in the requirement that $u\left( t\right) \in 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\left( L\right)$ is the pointwise value of the (continuous) test function $v\in V$.

Theorem 3.2.  Suppose $P(t)\in 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

$\begin{eqnarray} u\left( x,t\right) =\sum_{n=0}^{\infty }u_{n}\left( t\right) y_{n}\left( x\right) \end{eqnarray}$

(19)

where

$\begin{eqnarray} u_{n}\left( t\right) =-\frac{2K_{0}}{L}\frac{\left( -1\right) ^{n}}{ 1+K_{0}\eta \lambda _{n}}\exp \left( -\frac{K_{0}\mu \lambda _{n}}{ 1+K_{0}\eta \lambda _{n}}t\right) \ast P\left( t\right) \end{eqnarray}$

(20)

so that

$\begin{eqnarray} \frac{\partial u}{\partial t}\left( x,t\right) =\sum_{n=0}^{\infty }u_{n}^{\prime }\left( t\right) y_{n}\left( x\right) \end{eqnarray}$

(21)

where

$\begin{eqnarray} u_{n}^{\prime }\left( t\right) =-\frac{2K_{0}}{L}\frac{\left( -1\right) ^{n} }{1+K_{0}\eta \lambda _{n}}P\left( t\right) -\frac{K_{0}\mu \lambda _{n}}{ 1+K_{0}\eta \lambda _{n}}u_{n}\left( t\right) . \end{eqnarray}$

(22)

[Regularity]. Firstly, we show that $u(x,t)$ has the desired regularity, namely $u\in H^{1}\left( 0,T;V\right)$ and $u\left( 0\right) =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

$\begin{eqnarray*} \left\vert e^{-at}\ast P\left( t\right) \right\vert \leq \frac{1}{\sqrt{2a}} \left\Vert P\right\Vert _{L^{2}\left( 0,T\right) }. \end{eqnarray*}$

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

$\begin{eqnarray} \left\vert u_{n}\left( t\right) \right\vert \leq \frac{C}{\lambda _{n}} \left\Vert P\right\Vert _{L^{2}\left( 0,T\right) } \end{eqnarray}$

(23)

and

$\begin{eqnarray} \left\vert u_{n}^{\prime }\left( t\right) \right\vert \leq \frac{C}{\lambda _{n}}\left( \left\vert P\left( t\right) \right\vert +\left\Vert P\right\Vert _{L^{2}\left( 0,T\right) }\right) .\end{eqnarray}$

(24)

Thus, $u,u^{\prime }\in L^{2}\left( 0,T;V^{s}\right)$ for every $s<3/2$, thereby implying that conditions (ⅰ) and (ⅲ) of Def. 3.1 are satisfied.

[Existence]. In order to show that $u$ satisfies condition (ⅱ), it suffices to consider $v=y_{n}$ as test function, since $\left\{ \sqrt{\frac{2 }{L\lambda _{n}}}y_{n}\left( x\right) \right\}$ are a basis of $V$. Each term on the left hand side of (18) can be computed as

$\begin{eqnarray*} \left( u^{\prime }\left( t\right) ,y_{n}\right) _{H}=\frac{L}{2} u_{n}^{\prime }\left( t\right) \end{eqnarray*}$

$\begin{eqnarray*} \left( u\left( t\right) ,y_{n}\right) _{V}=\left( \dfrac{\partial u}{ \partial x}\left( \cdot ,t\right) ,\dfrac{dy_{n}}{dx}\left( \cdot \right) \right) _{H}=\frac{L}{2}\lambda _{n}u_{n}\left( t\right) \end{eqnarray*}$

$\begin{eqnarray*} \left( u^{\prime }\left( t\right) ,y_{n}\right) _{V}=\left( \dfrac{\partial ^{2}u}{\partial t\partial x}\left( \cdot ,t\right) ,\dfrac{dy_{n}}{dx}\left( \cdot \right) \right) _{H}=\frac{L}{2}\lambda _{n}u_{n}^{\prime }\left( t\right) ~. \end{eqnarray*}$

Adding the above three terms, we obtain that

$\begin{eqnarray*} \left( u^{\prime }\left( t\right) ,y_{n}\right) _{H}+K_{0}\left( \mu u\left( t\right) +\eta u^{\prime }\left( t\right) ,y_{n}\right) _{V}=\frac{L}{2} \left[ \left( 1+K_{0}\eta \lambda _{n}\right) u_{n}^{\prime }\left( t\right) +K_{0}\mu \lambda _{n}u_{n}\left( t\right) \right]~. \end{eqnarray*}$

Thus, from (22) and the fact that $y_{n}\left( L\right) =\left( -1\right) ^{n}$, it follows that the right hand side is equal to $-K_{0}P\left( t\right) y_{n}\left( L\right)$, 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\in H^{1}\left( 0,T;V\right)$ be a solution of

$\begin{eqnarray*} \left( u^{\prime }\left( t\right) ,v\right) _{H}+K_{0}\left( \mu u\left( t\right) +\eta u^{\prime }\left( t\right) ,v\right) _{V}=0,\qquad v\in V~. \end{eqnarray*}$

If we choose $v=u\left( t\right)$ as a test function, then we obtain that

$\begin{eqnarray*} \left( u^{\prime }\left( t\right) ,u\left( t\right) \right) _{H}+K_{0}\mu \left\Vert u\left( t\right) \right\Vert _{V}^{2}+K_{0}\eta \left( u^{\prime }\left( t\right) ,u\left( t\right) \right) _{V}=0 \end{eqnarray*}$

implying that

$\begin{eqnarray*} \frac{1}{2}\frac{d}{dt}\left( \left\Vert u\left( t\right) \right\Vert _{H}^{2}+K_{0}\eta \left\Vert u\left( t\right) \right\Vert _{V}^{2}\right) =-K_{0}\mu \left\Vert u\left( t\right) \right\Vert _{V}^{2}\leq 0~. \end{eqnarray*}$

Integrating with respect to time and using the initial condition $u\left( 0\right) =0$, we get

$\begin{eqnarray*} \left\Vert u\left( t\right) \right\Vert _{H}^{2}+K_{0}\eta \left\Vert u\left( t\right) \right\Vert _{V}^{2}\leq 0 \end{eqnarray*}$

Hence $u\left( t\right) =0$ for every $t$ in $\left[ 0,T\right]$.

Remark 4.  From (23), the fact that $\left\vert y_{n}\left( x\right) \right\vert \leq 1$, and the standard Weierstrass Test, it follows that the series (19) converges absolutely and uniformly in the space-time rectangle $Q=\left[ 0,L\right] \times \left[ 0,T\right]$ and its sum is continuous, i.e. $u\left( x,t\right) \in C^{0}\left( Q\right)$. In a similar way, under the stronger assumption $P\left( t\right) \in L^{\infty }\left( 0,T\right)$, estimate (24) implies $\left\vert u_{n}^{\prime }\left( t\right) \right\vert \leq C\lambda _{n}^{-1}\left\Vert P\right\Vert _{L^{\infty }\left( 0,T\right) }$ so that $\frac{ \partial u}{\partial t}\left( x,t\right) \in C^{0}\left( Q\right)$ also holds. Thus, in the physically realistic case of bounded boundary traction $P\left( t\right)$, 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).

Case 2. ($\eta =0$). In the purely elastic problem, the formal solution (16) is again written as given in (19) where now

$\begin{eqnarray*} u_{n}\left( t\right) =-\frac{2K_{0}}{L}\left( -1\right) ^{n}\exp \left( -K_{0}\mu \lambda _{n}t\right) \ast P\left( t\right). \end{eqnarray*}$

Then, estimates (23) and (24) become, respectively,

$\begin{eqnarray*} \left\vert u_{n}\left( t\right) \right\vert \leq \frac{C}{\sqrt{\lambda _{n}} }\left\Vert P\right\Vert _{L^{2}\left( 0,T\right) } \end{eqnarray*}$

and

$\begin{eqnarray*} \left\vert u_{n}^{\prime }\left( t\right) \right\vert \leq C\sqrt{\lambda _{n}}\left( \left\vert P\left( t\right) \right\vert +\left\Vert P\right\Vert _{L^{2}\left( 0,T\right) }\right). \end{eqnarray*}$

As a consequence, it can only be asserted that $u\in L^{2}\left( 0,T;V^{s}\right)$ for every $s<1/2$ and $u^{\prime }\in L^{2}\left( 0,T;V^{s}\right)$ for every $s<-3/2$. On the other hand, if one assumes $P\left( t\right) \in L^{\infty }\left( 0,T\right)$, we find that

$\begin{eqnarray*} \left\vert u_{n}\left( t\right) \right\vert \leq \frac{C}{\lambda _{n}} \left\Vert P\right\Vert _{L^{\infty }\left( 0,T\right) } \end{eqnarray*}$

and

$\begin{eqnarray*} \left\vert u_{n}^{\prime }\left( t\right) \right\vert \leq C\left\Vert P\right\Vert _{L^{\infty }\left( 0,T\right) }. \end{eqnarray*}$

Hence, in this case, $u\in L^{2}\left( 0,T;V^{s}\right)$ for every $s<3/2$ and $u^{\prime }\in L^{2}\left( 0,T;V^{s}\right)$ for every $s<-1/2$. In particular, when the boundary traction is bounded, we have $u\left( x,t\right) \in C^{0}\left( Q\right)$ whereas the discharge velocity is a distribution. In conclusion, we state the following definition and theorem (without proof):

Definition 3.3.  A function $u:\left[ 0,T\right] \rightarrow V$ is a weak solution to problem (11) if:

(ⅰ) $u\in L^{2}\left( 0,T;V\right)$, $u^{\prime }\in L^{2}\left( 0,T;V^{\prime }\right)$;

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

$\begin{eqnarray} \left\langle u^{\prime }\left( t\right) ,v\right\rangle +K_{0}\mu \left( u\left( t\right) ,v\right) _{V}=-K_{0}P\left( t\right) v\left( L\right) \end{eqnarray}$

(25)

where the brackets $\left\langle \cdot ,\cdot \right\rangle$ denote pairing between $V^{\prime }$ and $V$;

(ⅲ) $u\left( 0\right) =0.$

Theorem 3.4.  Suppose $P(t)\in L^{\infty }(0,T)$. Then there exists a unique weak solution $u$ to problem (11), in the sense of Definition 3.3. Moreover, $u$ is represented by the series expansion (16).

---

4. 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:

$\begin{eqnarray} \begin{array}{c} \displaystyle\hat{x} =\frac{x}{[x]},\qquad \hat{t}=\frac{t}{[t]},\qquad \hat{\eta}=\frac{\eta}{[\eta]}, \qquad \hat{\lambda} _{n} =\frac{\lambda _{n}}{[\lambda _{n}]}, \\[.1in] \displaystyle\hat{P} = \frac{P}{[P]}, \qquad \hat{u} =\frac{u}{[u]},\qquad \hat{v}=\frac{v}{[v]}, \qquad \hat{\mathcal{P}} =\frac{\mathcal P}{[\mathcal{P}]} ~. \end{array} \end{eqnarray}$

(26)

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

$\begin{eqnarray} [P]=P_{\text{ref}} \end{eqnarray}$

(27)

where $P_{\text{ref}}$ is a reference value, for example the mean value, of the given function $P(t)$. The expression for $\lambda_n$ in (14) suggests to choose

$\begin{eqnarray} [\lambda_n]=\frac{1}{L^2} \quad \mbox{and} \quad [x]=L \end{eqnarray}$

(28)

and, consequently, the expression for $u(x,t)$ in (15) suggests to choose

$\begin{eqnarray}[\eta] = \frac{1}{K_0[\lambda_n]} = \frac{L^2}{K_0}, \quad [t] = \frac{1}{K_0\mu[\lambda_n]} = \frac{L^2}{K_0\mu}, \quad [u] = \frac{K_0}{L}[P][t] = \frac{P_{\text{ref}} L}{\mu} ~, \end{eqnarray}$

(29)

whereas the expressions for $v(x,t)$ in (10b) and $\mathcal P(t)$ in (12) suggest to choose

$\begin{eqnarray}[v] = \frac{[u]}{[t]} = \frac{K_0}{L}P_{\text{ref}} \end{eqnarray}$

(30)

and

$\begin{eqnarray}[\mathcal P] = \frac{L}{K_0}[v]^2 = \frac{K_0}{L}P_{\text{ref}}^2 \end{eqnarray}$

(31)

respectively. Using the above scalings, we obtain the following dimensionless problem:

$\begin{align} \dfrac{\partial \hat{u}}{\partial \hat{t}}- \dfrac{\partial ^{2} \hat{u}}{\partial \hat{x} ^{2}} -\hat{\eta} \dfrac{\partial ^{3}\hat{u}}{\partial \hat{t}\partial \hat{x} ^{2}}=0 & & & \mbox{in}\ (0,1)\times (0,\hat{T}) \end{align}$

(32a)

$\begin{align} \dfrac{\partial \hat{u}}{\partial \hat{x}}\left( 1,\hat{t}\right) +\hat{\eta} \dfrac{\partial ^{2}\hat{u}}{\partial \hat{t}\partial \hat{x}}\left( 1, \hat{t}\right) =-\hat{P}\left( \hat{t}\right) & & & \mbox{for}\;0<\hat{t}<\hat{T} \end{align}$

(32b)

$\begin{align} \hat{u}\left( 0,\hat{t}\right) =0 & & & \mbox{for}\;0<\hat{t}<\hat{T} \end{align}$

(32c)

$\begin{align} \hat{u}\left( \hat{x},0\right) =0 & & & \mbox{for}\;0<\hat{x}<1 \end{align}$

(32d)

where $\hat{T} = T/[t]$ and $\hat{\eta} \geq 0$. Recalling again the expressions (15), (10b) and (12), we can now write solid displacement, discharge velocity and power density in dimensionless form as

$\begin{eqnarray} \hat{u}\left( \hat{x},\hat{t}\right) =-2\sum_{n=0}^{\infty }\frac{\left( -1\right) ^{n}y_{n}\left( \hat{x}\right) }{1+\hat{\eta} \hat{\lambda} _{n}}\exp \left( - \frac{ \hat{\lambda} _{n}}{1+\hat{\eta} \hat{\lambda} _{n}}\hat{t}\right) \ast \hat{P}\left( \hat{t}\right) \end{eqnarray}$

(33)

$\begin{eqnarray} \hat{v}\left( \hat{x},\hat{t}~\right) =2\sum_{n=0}^{\infty }\frac{\left( -1\right) ^{n}y_{n}\left( \hat{x}\right) }{1+\hat{\eta} \hat{\lambda} _{n}} \left\{ \hat{P}\left( \hat{t}~\right) -\frac{ \hat{\lambda} _{n}}{1+\hat{\eta} \hat{\lambda} _{n}}\exp \left( - \frac{ \hat{\lambda} _{n}}{1+\hat{\eta} \hat{\lambda} _{n}}\hat{t}\right) \ast \hat{P}\left( \hat{t}\right)\right\} \end{eqnarray}$

(34)

$\begin{eqnarray}\hat{\mathcal{P}}\left( \hat{t}\right) =2\sum_{n=0}^{\infty }\frac{1}{\left( 1+\hat{\eta} \hat{\lambda} _{n}\right) ^{2}}\left\{ \hat{P}\left( \hat{t}~\right) -\frac{\hat{\lambda} _{n} }{1+\hat{\eta} \hat{\lambda} _{n}}\exp \left( - \frac{ \hat{\lambda} _{n}}{1+\hat{\eta} \hat{\lambda} _{n}}\hat{t}\right) \ast \hat{P}\left( \hat{t}\right) \right\} ^{2}~. \end{eqnarray}$

(35)

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_{\text{ref}}/L$, differs from that utilized in [37], which was $K_0 \mu/L$ instead. This difference actually follows from a different choice for the pressure scaling, namely $P_{\text{ref}}$ in our case, as opposed to $\mu$ in [37]. Indeed, the two scalings are equivalent if $P_{\text{ref}}$ and $\mu$ are of the same order of magnitude. In general, though, the scaling $K_0 \mu/L$ represents the characteristic velocity of wave propagation within the medium, whereas the scaling $K_0 P_{\text{ref}}/L$ represents the characteristic velocity induced by an external load of magnitude $P_{\text{ref}}$. Since our analysis aims at assessing the influence of the external load on the biomechanical response of the medium, we chose $[P]=P_{\text{ref}}$. However, the analysis could be easily adapted to other choices, should the scope of the investigation differ from that of this paper.

---

5. The case $\hat{P}\left( \hat{t}~\right) =\mathrm{step}$ $\mathrm{pulse}$

Let the boundary traction $\hat{P}\left( \hat{t}~\right)$ be the dimensionless unit step starting at $\hat{t}=0$ defined as

$\begin{eqnarray}\hat{P}\left( \hat{t}~\right) =H\left( \hat{t}~\right) =\left\{ \begin{array}{lll} 0 & \text{if} & \hat{t}<0 \\ 1 & \text{if} & \hat{t}\geq 0 \end{array} \right. ~. \end{eqnarray}$

(36)

In this case, the dimensionless solid displacement (33), discharge velocity (34) and power density (35) (hereon denoted with the subscript $\hat{\eta}$ for later reference) become, respectively,

$\begin{align} \hat{u}_{\hat{\eta}}\left( \hat{x},\hat{t}~\right) & =-2\sum_{n=0}^{\infty }\frac{\left( -1\right) ^{n}}{\hat{\lambda}_{n}}\left\{ 1-\exp \left( -\frac{\hat{\lambda}_{n}\hat{t}}{ 1+\hat{\eta}\hat{\lambda}_{n}}\right) \right\} y_{n}\left( \hat{x}\right) & & \nonumber\\ & =-\hat{x}+2\sum_{n=0}^{\infty }\frac{\left( -1\right) ^{n}}{\hat{\lambda}_{n}}\exp \left( -\frac{\hat{\lambda}_{n}\hat{t}}{1+\hat{\eta}\hat{\lambda}_{n}}\right) y_{n}\left( \hat{x} \right) & & \end{align}$

(37a)

$\begin{align} \hat{v}_{\hat{\eta}}\left( \hat{x},\hat{t}\right) & =2\sum_{n=0}^{\infty }\frac{\left( -1\right) ^{n}}{1+\hat{\eta}\hat{\lambda}_{n}}\exp \left( -\frac{\hat{\lambda}_{n}\hat{t}}{1+ \hat{\eta}\hat{\lambda}_{n}}\right) y_{n}\left( \hat{x}\right) & & \end{align}$

(37b)

$\begin{align} \hat{\mathcal{P}}_{\hat{\eta}}\left( \hat{t}\right) & =2\sum_{n=0}^{\infty }\frac{1 }{\left( 1+\hat{\eta}\hat{\lambda}_{n}\right) ^{2}}\exp \left( -\frac{2\hat{\lambda}_{n} \hat{t}}{1+\hat{\eta}\hat{\lambda}_{n}}\right) ~. & & \end{align}$

(37c)

The solution in the purely elastic case can be obtained by setting $\hat{\eta} =0$ in (37).

Remark 6.  If the unit step is shifted at $\hat{t}=\alpha >0$, i.e. $\hat{P}\left( \hat{ t}\right) =H\left( \hat{t}-\alpha \right)$, then the solution (33) is

$\begin{eqnarray*} \hat{u}\left( \hat{x},\hat{t}\right) =\hat{u}_{\hat{\eta}}\left( \hat{x}, \hat{t}-\alpha \right) H\left( \hat{t}-\alpha \right) =\left\{ \begin{array}{lll} 0 & \text{if} & 0\leq \hat{t}<\alpha \\ \hat{u}_{\hat{\eta}}\left( \hat{x},\hat{t}-\alpha \right) & \text{if} & \hat{t}\geq \alpha. \end{array} \right. \end{eqnarray*}$

The space-time behavior of $\hat{u}_{\hat{\eta}}$, $\hat{v}_{\hat{\eta}}$ and $\hat{\mathcal P}_{\hat{\eta}}$ is reported in Figures 4, 5 and 6, respectively, in the case of $\hat{\eta}=0$ and $\hat{\eta}=1$. We remark that $\hat{v}_{\hat{\eta}}$ and $\hat{\mathcal P}_{\hat{\eta}}$ are reported in logarithmic scale to highlight the presence of a blow-up at $\hat{x}=1$ and $\hat{t}=0$ in the purely elastic case $\hat{\eta}=0$. Indeed, by setting $\hat{\eta}=\hat{t}=0$ in (37c), we can verify that the power density satisfies

$\begin{eqnarray*} \hat{\mathcal{P}}_{0}\left( 0 \right) =2\sum_{n=0}^{\infty } 1 = +\infty. \end{eqnarray*}$

Proceeding analogously in the case of the discharge velocity, the Fourier expansion in the purely elastic case is given by

$\begin{eqnarray*} \hat{v}_{0}\left( \hat{x},0\right) =2\sum_{n=0}^{\infty }\left( -1\right) ^{n}y_{n}\left( \hat{x}\right) \end{eqnarray*}$

which clearly lacks pointwise convergence for any $\hat{x}\neq 0$.

From the physical viewpoint, this means that, at the switch on time of the driving term, here set at $\hat{t}=0$, the discharge velocity contains high-frequency components whose superposition exhibits the following behaviors depending on $\hat{x}$: (ⅰ) it relaxes towards zero as $\hat{x}\rightarrow 0$; (ⅱ) it tends to infinity as $\hat{x} \rightarrow 1$, thereby exhibiting a blow-up; (ⅲ) it looks like an amorphous blob for $0<\hat{x}<1$. On the other hand, if viscoelastic effects are present, then the dimensionless discharge velocity at $\hat{t}=0$ is given by

$\begin{eqnarray*} \hat{v}_{\hat{\eta}}\left( \hat{x},0\right) =2\sum_{n=0}^{\infty }\frac{ \left( -1\right) ^{n}}{1+\hat{\eta}\hat{\lambda}_{n}}y_{n}\left( \hat{x} \right)~. \end{eqnarray*}$

Here, the $n-$th coefficient is decreasing as $1/n^{2}$ and the discharge velocity is continuous at $\hat{t}=0$. These concepts are illustrated in Fig. 7, which shows the spatial distribution of $\hat{v}_{\hat{\eta}}$ at $\hat{t}=0$ for several values of $\hat{\eta}$ in the interval $[0,1]$. The $\log_{10}$ plot of $\hat{v}_{\hat{\eta}}$ is reported in the figure to highlight the blow-up of the discharge velocity in the purely elastic case at $\hat{x}=1$.

In order to further investigate this blow-up and its dependence on the structural viscoelasticity, we observe that the maximum value of $\hat{v}_{\hat{\eta}}$ is attained at $\hat{x}=1$ and $\hat{t}=0$ and can be written as a function of $\hat{\eta}$ as follows:

$\begin{eqnarray*} \hat{v}_{\max }\left( \hat{\eta}\right) =\max_{\substack{ 0\leq \hat{x}\leq 1 \nonumber\nonumber\nonumber\nonumber\nonumber\nonumber\\ \hat{t} \geq 0}}\left\vert \hat{v}_{\hat{\eta} }\left( \hat{x},\hat{t}\right) \right\vert =\hat{v}_{\hat{\eta} }\left( 1,0\right) =2\sum_{n=0}^{\infty }\frac{1}{1+\hat{\eta} \hat{\lambda} _{n}}~. \end{eqnarray*}$

The above series may be summed ([19], formula no. 1.4212) and the final result is

$\begin{eqnarray} \hat{v}_{\max }\left( \hat{\eta}\right) =\frac{1}{\sqrt{\hat{\eta}}}\tanh \left( \frac{1}{ \sqrt{\hat{\eta}}}\right) ~. \end{eqnarray}$

(38)

Similarly, the dimensionless power density (37c) is decreasing in time and its maximum is attained at $\hat{t}=0$:

$\begin{eqnarray*} \hat{\mathcal{P}}_{\max }\left( \hat{\eta}\right) &=&\max_{\hat{t}\geq 0} \hat{\mathcal{P}}_{\hat{\eta}}\left( \hat{t}\right) =\hat{\mathcal{P}}_{\hat{ \eta}}\left( 0\right) =2\sum_{n=0}^{\infty }\frac{1}{\left( 1+\hat{\eta}\hat{ \lambda}_{n}\right) ^{2}} \\ &=&\frac{1}{2\hat{\eta}}\left( \tanh ^{2}\frac{1}{\sqrt{\hat{\eta}}}+\sqrt{ \hat{\eta}}\tanh \frac{1}{\sqrt{\hat{\eta}}}-1\right)~. \end{eqnarray*}$

The behaviors of $\hat{v}_{\max }\left( \hat{\eta}\right)$ and $\hat{\mathcal{P}}_{\max }\left( \hat{\eta}\right)$ with respect to $\hat{\eta}$ are visualized in Fig. 8.

Clearly, the dimensionless parameter $\hat{\eta}$ is the sole determinant of the blow-up in $\hat{v}_{\max }\left( \hat{\eta}\right)$ and $\hat{\mathcal{P}}_{\max }\left( \hat{\eta}\right)$, with smaller values of $\hat{\eta}$ leading to higher values of $\hat{v}_{\max }\left( \hat{\eta}\right)$ and $\hat{\mathcal{P}}_{\max }\left( \hat{\eta}\right)$. However, it is important to recall that $\hat{\eta}$ does not depend solely on the viscoelasticity coefficient. The definition $\hat{\eta}=\eta 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 $\hat{\eta}$. Precisely, going back to dimensional units, we can use (30) and (31) to obtain the following expressions:

$\begin{align} & v_{\max}\left( \hat{\eta} \right) = \frac{P_{\text{ref}} K_0}{L} \hat{v}_{\max }\left( \hat{\eta}\right) & \end{align}$

(39)

$\begin{align} & \mathcal{P}_{\max} \left( \hat{\eta} \right) = \frac{P_{\text{ref}}^2 K_0}{L} \hat{\mathcal{P}}_{\max} \left( \hat{\eta} \right). & \end{align}$

(40)

We can see that the magnitude $P_{\text{ref}}$ of the driving term appears as a multiplying constant in the expressions above, meaning that larger values of $P_{\text{ref}}$ will lead to larger values of $v_{\max}$ and $\mathcal{P}_{\max }$. However, the factor that controls the blow-up still remains the dimensionless parameter $\hat{\eta} = \eta K_0/L^2$. These concepts are illustrated in Fig. 9, where we set $K_0/L = 1 ~\mathrm{m^2 s ~Kg^{-1}}$ and we consider increasing values of $P_{\text{ref}}$ in the range $[10^{-3}, 10^3] ~ \mathrm{N m^{-2}}$.

---

6. The case $\hat{P}(\hat{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, let $\hat{P}\left( \hat{t}\right)$ be the dimensionless trapezoidal pulse of unit amplitude and let $\hat{\varepsilon}$, $\hat{\tau}>0$. Here we consider the following expression for $\hat{P}\left( \hat{t}\right)$:

$\begin{eqnarray}\hat{P}\left( \hat{t}\right) =\left\{ \begin{array}{lll} 0 & \text{if} & \hat{t}<0 \\ \dfrac{\hat{t}}{\hat{\varepsilon}} & \text{if} & 0\leq \hat{t}<\hat{ \varepsilon} \\ 1 & \text{if} & \hat{\varepsilon}\leq \hat{t}<\hat{\varepsilon}+\hat{\tau} \\ \dfrac{\hat{\tau}-\hat{t}}{\hat{\varepsilon}}+2 & \text{if} & \hat{ \varepsilon}+\hat{\tau}\leq \hat{t}<2\hat{\varepsilon}+\hat{\tau} \\ 0 & \text{if} & \hat{t}\geq 2\hat{\varepsilon}+\hat{\tau} \end{array} \right. \end{eqnarray}$

(41)

where $\hat{\varepsilon}$ is the pulse rise/fall time and $\hat{\tau}$ is the pulse duration, as depicted in Fig. 10. Thus, the trapezoidal pulse defined in (41) reduces to a rectangular step pulse as $\hat{\varepsilon}\rightarrow 0$.

Let us now compute the dimensionless discharge velocity resulting from the application of the trapezoidal pulse at the boundary, henceforth denoted by $\hat{V}_{\hat{\eta}}\left( \hat{x},\hat{t} \right)$ to distinguish it from the discharge velocity $\hat{v}_{\hat{\eta}}\left( \hat{x},\hat{t} \right)$ 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 times $\hat{t}=0$, $\hat{\varepsilon}$ , $\hat{\varepsilon}+\hat{\tau}$,$\ 2\hat{\varepsilon}+\hat{\tau}$:

$\begin{eqnarray} \hat{P}\left( \hat{t}\right) &=&\frac{1}{\hat{\varepsilon}}\left\{ \hat{t} H\left( \hat{t}\right) -\left( \hat{t}-\hat{\varepsilon}\right) H\left( \hat{ t}-\hat{\varepsilon}\right) \right. \end{eqnarray}$

(42)

$\begin{eqnarray} &&\left. -\left( \hat{t}-\hat{\tau}-\hat{\varepsilon}\right) H\left( \hat{t}- \hat{\tau}-\hat{\varepsilon}\right) +\left( \hat{t}-\hat{\tau}-2\hat{ \varepsilon}\right) H\left( \hat{t}-\hat{\tau}-2\hat{\varepsilon}\right) \right\} \nonumber \end{eqnarray}$

where the function $H(\hat{t})$ was defined in (36). In addition, we can see that

$\begin{eqnarray*} \hat{u}\left( \hat{x},\hat{t}\right) =\int_{0}^{\hat{t}}\hat{u}_{\hat{\eta} }\left( \hat{x},s\right) ds \end{eqnarray*}$

solves problem (32) with $\hat{P}\left( \hat{t}\right) =\hat{t}H\left( \hat{t}\right)$ (the ramp pulse of unit slope starting at time $\hat{t}=0$), where $\hat{u}_{\hat{\eta}}\left( \hat{x},\hat{t}\right)$ is given by (37a). Thus, using (10b), we have that the discharge velocity is given by $-\hat{u}_{\hat{\eta}}\left( \hat{x},\hat{t}\right)$, so that, applying the linear superposition principle, $\hat{V}_{\hat{\eta}}\left( \hat{x},\hat{t}\right)$ is the linear combination of the velocities corresponding to the various components of the pulse, namely:

$\begin{eqnarray} \hat{V}_{\hat{\eta}}\left( \hat{x},\hat{t}\right) &=&\frac{1}{\hat{ \varepsilon}}\left\{ -\hat{u}_{\hat{\eta}}\left( \hat{x},\hat{t}\right) H\left( \hat{t}\right) +\hat{u}_{\hat{\eta}}\left( \hat{x},\hat{t}-\hat{ \varepsilon}\right) H\left( \hat{t}-\hat{\varepsilon}\right) \right. \nonumber\\ &&+\hat{u}_{\hat{\eta}}\left( \hat{x},\hat{t}-\hat{\tau}-\hat{ \varepsilon}\right) H\left( \hat{t}-\hat{\tau}-\hat{\varepsilon}\right)\\ &&\left. - \hat{u}_{\hat{\eta}}\left( \hat{x},\hat{t}-\hat{\tau}-2\hat{\varepsilon} \right) H\left( \hat{t}-\hat{\tau}-2\hat{\varepsilon}\right) \right\}~. \nonumber \end{eqnarray}$

(43)

An illustration of the typical form of $\hat{V}_{\hat{\eta}}\left( \hat{x},\hat{t}\right)$ for $\hat{\eta} \geq 0$ and $\hat{\varepsilon}>0$ is reported in Fig. 11.

The maximum possible discharge velocity occurs at $\hat{x}=1$, $\hat{t}=\hat{\varepsilon}$ and can be written as

$\begin{align} \hat{V}_{\max }\left( \hat{\eta},\hat{\varepsilon}\right) &=\max_{\substack{ 0\leq \hat{x}\leq 1 \nonumber\\ \hat{t}\geq 0}}\left\vert \hat{V}_{\hat{\eta}}\left( \hat{x},\hat{t}\right) \right\vert =\hat{V}_{\hat{\eta}}\left( 1,\hat{ \varepsilon}\right) \nonumber\\ &=\frac{2}{\hat{\varepsilon}}\sum_{n=0}^{\infty }\frac{1}{ \hat{\lambda}_{n}}\left\{ 1-\exp \left( -\frac{\hat{\lambda}_{n}\hat{ \varepsilon}}{1+\hat{\eta}\hat{\lambda}_{n}}\right) \right\} .\end{align}$

(44)

The behavior of $\hat{V}_{\max }$ with respect to $\hat{\eta}$ and $\hat{\varepsilon}$ is reported in Fig. 12. Interestingly, for all $\hat{\eta}>0$ and $\hat{\varepsilon}\geq 0$ we have

$\begin{eqnarray*} \hat{V}_{\max }\left( \hat{\eta},\hat{\varepsilon}\right) \leq \hat{V}_{\max }\left( \hat{\eta},0\right) \end{eqnarray*}$

and

$\begin{eqnarray*} \hat{V}_{\max }\left( \hat{\eta},0\right) =\hat{v}_{\max }\left( \hat{\eta} \right) \end{eqnarray*}$

since the trapezoidal pulse reduces to a rectangular pulse as $\hat{ \varepsilon}\rightarrow 0$. Thus, no blow-up takes place in the viscoelastic case when the pulse is rectangular. Similarly, for all $\hat{ \eta}\geq 0$ and $\hat{\varepsilon}>0$ we have

$\begin{eqnarray*} \hat{V}_{\max }\left( \hat{\eta},\hat{\varepsilon}\right) \leq \hat{V}_{\max }\left( 0,\hat{\varepsilon}\right) =\frac{2}{\hat{\varepsilon}} \sum_{n=0}^{\infty }\frac{1-\exp \left( -\hat{\lambda}_{n}\hat{\varepsilon} \right) }{\hat{\lambda}_{n}}, \end{eqnarray*}$

hence no blow-up takes place even in the purely elastic case when the pulse is trapezoidal.

---

7. 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_{\text{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 [22]), does not account for structural viscoelasticity, which amounts to setting $\eta=0$ in (7a). Table 1 summarizes the parameter values corresponding to the experiments on articular cartilage reported in [37].

Characteristic velocities for the confined compression experiments reported in [37] are of the order of 0.35 $\mu$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: (ⅰ) the load protocol for the confined compression experiment corresponds to the trapezoidal pulse described in Section 6; and that (ⅱ) we should ensure that discharge velocities always remain below the threshold value $V_{\text{th}}=0.35$ $\mu$m/s in order to preserve the sample from microstructural damage. This means that $\hat{\eta}$ and $\hat{\varepsilon}$ should be chosen as to guarantee that

$\begin{eqnarray} \hat{V}_{\max } (\hat{\eta},\hat{\varepsilon}) < \hat{V}_{\text{th}} = \frac{L}{P_{\text{ref}} K_0}~ V_{\text{th}} = 16.3~, \end{eqnarray}$

(45)

where $\hat{V}_{\max } (\hat{\eta},\hat{\varepsilon})$ can be computed using the expression (44). Fig. 14 compares $\hat{V}_{\text{th}}$ with the values of $\hat{V}_{\max }$ obtained for different values of $\hat{\eta}$ and $\hat{\varepsilon}$. Indeed, we can identify a whole region in the $(\hat{\eta},\hat{\varepsilon}) -$plane for which $\hat{V}_{\max } (\hat{\eta},\hat{\varepsilon}) < \hat{V}_{\text{th}}$.

To better visualize the various regions of interest in the $(\hat{\eta},\hat{\epsilon}) -$plane, we reported the colormap of the difference $\hat{V}_{\text{th}}-\hat{V}_{\max }$ as a function of $\hat{\eta}$ and $\hat{\varepsilon}$ in Fig. 15, where we clearly marked the curve in the parameter space for which $\hat{V}_{\text{th}}=\hat{V}_{\max }$. These results show that, for example, choosing $\hat{\varepsilon} \geq 5 \cdot 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 $\hat{\varepsilon}=0$, we can still maintain the maximum velocity below the threshold by modifying the structural properties of the sample to have, for example, $\hat{\eta}\geq4\cdot 10^{-3}$.

Interestingly, expression (29) defines $\hat{\eta}= (K_0/L^2)\eta$, meaning that $\hat{\eta}$ can be increased by increasing the viscoelastic parameter $\eta$, increasing the permeability parameter $K_0$, or decreasing the sample thickness $L$. In the same spirit as [7], a typical value of the viscoelastic parameter $\eta$ for biological tissues can be estimated by setting

$\begin{align} & \eta= \mu \tau_e & \end{align}$

(46a)

where

$\begin{align} & \tau_e = L \sqrt{\frac{\rho}{\mu}} & \end{align}$

(46b)

is the characteristic elastic time constant of the porous material under compression and $\rho$ is the mass density of the fluid component of the mixture. Since the densities of most of the biological fluids are very similar to that of water, let us set $\rho=1000 \mathrm{~Kg~ m^{-3}}$. Substituting (46a) and (46b) in (29) we obtain

$\begin{align} & \hat{\eta} = \frac{\tau_e}{[t]}= \frac{K_0}{L} \sqrt{\rho \mu}. & \end{align}$

(46c)

Replacing the parameter values of Table 1 into (46c) we obtain $\hat{\eta} = 1.15 \cdot 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: (ⅰ) vary the properties of the biological sample by manipulating $\rho$, $\mu$, $\eta$, $K_0$ or $L$; and/or (ⅱ) vary the time profile of the load experimental protocol. The choice will depend on the particular conditions and constraints of the experiment at hand.

---

8. 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. Compared to the analysis carried out in [6], in which singularities in the solution were demonstrated by numerical evidences when the hypotheses of the existence theorem were not satisfied, the 1D study conducted here allowed us, via analytical solution, to prove that singularities may arise even in the simple case of constant permeability. 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 contribution, namely $\hat{\eta}\ll 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 [6]. Interestingly, the a priori estimates derived in [6] 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. We expect that similar results can be obtained for the 3D fluid-solid mixture described in (2a)-(2g) with boundary conditions (3)-(5). In particular, if one considers the coupling driven by a source of linear momentum and a source of boundary traction that are not time differentiable, then one can only obtain local in time existence of weak solutions for the system (in comparison to [6]), and prove blow-up of solutions in finite time (in the norm of the fluid energy). We plan to provide this analysis in a subsequent paper.

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: (ⅰ) decreasing $P_{\text{ref}}$, thereby reducing the load on the medium; (ⅱ) increasing $\hat{\varepsilon}$, thereby smoothing the load action on the medium; (ⅲ) increasing $\eta$, thereby changing the structural properties of the solid component of the medium; (ⅳ) decreasing $K_0$, thereby changing the pore size, geometry and/or fluid viscosity; (ⅴ) 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 [32] 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 [5] and the more recent works [27, 28, 29] and [3, 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 [6] 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.

---

Acknowledgments

Dr. Bociu has been partially supported by NSF CAREER DMS-1555062. Dr. Guidoboni has been partially supported by the award NSF DMS-1224195, the Chair Gutenberg funds of the Cercle Gutenberg (France) and the LabEx IRMIA (University of Strasbourg, France). Dr. Sacco has been partially supported by Micron Semiconductor Italia S.r.l., SOW nr. 4505462139.

---