symbol | value | units |
L | 0.81⋅10−3 | m |
T | 105,104 | s |
θ | 0.97⋅106 | Nm−2 |
K0 | 2.9⋅10−16 | m4N−1s−1 |
Pref | 6⋅104 | Nm−2 |
c0 | 1.67⋅10−5 | m2N−1 |
η | 4.85⋅109,4.85⋅108 | Nsm−2 |
Citation: Lorena Bociu, Giovanna Guidoboni, Riccardo Sacco, Maurizio Verri. On the role of compressibility in poroviscoelastic models[J]. Mathematical Biosciences and Engineering, 2019, 16(5): 6167-6208. doi: 10.3934/mbe.2019308
[1] | Maurizio Verri, Giovanna Guidoboni, Lorena Bociu, Riccardo Sacco . The role of structural viscoelasticity in deformable porous media with incompressible constituents: Applications in biomechanics. Mathematical Biosciences and Engineering, 2018, 15(4): 933-959. doi: 10.3934/mbe.2018042 |
[2] | Xiawei Yang, Wenya Li, Yaxin Xu, Xiurong Dong, Kaiwei Hu, Yangfan Zou . Performance of two different constitutive models and microstructural evolution of GH4169 superalloy. Mathematical Biosciences and Engineering, 2019, 16(2): 1034-1055. doi: 10.3934/mbe.2019049 |
[3] | Matthias Ngwa, Ephraim Agyingi . A mathematical model of the compression of a spinal disc. Mathematical Biosciences and Engineering, 2011, 8(4): 1061-1083. doi: 10.3934/mbe.2011.8.1061 |
[4] | Pengfei Liu, Yantao Luo, Zhidong Teng . Role of media coverage in a SVEIR-I epidemic model with nonlinear incidence and spatial heterogeneous environment. Mathematical Biosciences and Engineering, 2023, 20(9): 15641-15671. doi: 10.3934/mbe.2023698 |
[5] | Jinyi Tai, Chang Liu, Xing Wu, Jianwei Yang . Bearing fault diagnosis based on wavelet sparse convolutional network and acoustic emission compression signals. Mathematical Biosciences and Engineering, 2022, 19(8): 8057-8080. doi: 10.3934/mbe.2022377 |
[6] | Martina Bukač, Sunčica Čanić . Longitudinal displacement in viscoelastic arteries:A novel fluid-structure interaction computational model, and experimental validation. Mathematical Biosciences and Engineering, 2013, 10(2): 295-318. doi: 10.3934/mbe.2013.10.295 |
[7] | Bo Wang, Yabin Li, Jianxiang Zhao, Xue Sui, Xiangwei Kong . JPEG compression history detection based on detail deviation. Mathematical Biosciences and Engineering, 2019, 16(5): 5584-5594. doi: 10.3934/mbe.2019277 |
[8] | Cong Wang, Chang Liu, Mengliang Liao, Qi Yang . An enhanced diagnosis method for weak fault features of bearing acoustic emission signal based on compressed sensing. Mathematical Biosciences and Engineering, 2021, 18(2): 1670-1688. doi: 10.3934/mbe.2021086 |
[9] | Raimund Bürger, Paola Goatin, Daniel Inzunza, Luis Miguel Villada . A non-local pedestrian flow model accounting for anisotropic interactions and domain boundaries. Mathematical Biosciences and Engineering, 2020, 17(5): 5883-5906. doi: 10.3934/mbe.2020314 |
[10] | T. Gayathri Devi, A. Srinivasan, S. Sudha, D. Narasimhan . Web enabled paddy disease detection using Compressed Sensing. Mathematical Biosciences and Engineering, 2019, 16(6): 7719-7733. doi: 10.3934/mbe.2019387 |
Poroviscoelastic models of multi-component mixtures are often utilized in biological applications to describe the flow of fluids within the pores of a deformable solid skeleton, see e.g., [1,2,3,4]. Skeleton viscoelasticity is often due to the complex structures including extracellular matrix, collagen and elastin that are present in biological tissues. Specifically, our work focuses on the bio-fluid-mechanical response of poroviscoelastic media to non-smooth data, since this aspect is crucial in understanding the mechanisms leading to tissue damage in the optic nerve head, and consequent vision loss, associated with glaucoma [5,6,7,8].
In the absence of viscoelasticity, we have recently shown that time irregularities in the volumetric and/or boundary sources of linear momentum lead to a blow-up in the solution of poroelastic models [1,4]. Interestingly, the blow-up can be prevented by including structural viscoelasticity [4]. From the application viewpoint, examples of time-irregularities in the data are discontinuities in intraocular pressure, which acts as a boundary traction for the optic nerve head tissue, or discontinuities in the gravitational force, which acts as a volumetric source of linear momentum. The intraocular pressure exhibits rapid changes every time we rub our eyes or we change posture [9], whereas rapid changes in the gravitational acceleration are experienced by jet pilots and astronauts during flights [10,11]. Since tissue viscoelasticity has been shown to decrease with age and/or disease conditions, the solution blow-up identified by our theoretical work led to hypothesize that rapid changes in intraocular pressure and gravitational acceleration, even if within physiological ranges, could damage the optic nerve head tissue if its viscoelasticity was pathologically reduced.
It is important to emphasize that our previous work was built on the assumption that the poroviscoelastic medium under consideration was made of incompressible components. The incompressibility assumption is quite common in biological applications, since tissues are mostly made of water. However, compressibility is always present in real tissues, and this leads to wonder whether and to what extent compressibility of the mixture components would affect the tissue response to non-smooth data. The present work aims at addressing this question.
Proceeding as in [4], we assume: (i) a one-dimensional (1D) geometry; (ii) negligible inertial terms in the linear momentum balance equation; (iii) a Kelvin-Voigt model for the effective stress tensor and (iv) a constant hydraulic permeability of the porous medium. Then, we express the fluid pressure and the solid stress as functions of the sole solid phase displacement, and we obtain an initial-boundary value problem (IBVP) of pseudo-parabolic type for the solid phase displacement. For this equivalent formulation we are able to construct an analytical solution and prove the well-posedness of weak solutions. Moreover, we recover analytical formulas for fluid pressure and discharge velocity, and discuss their regularity in terms of the regularity of the data. Finally, we analyze the behavior of all the solutions for various continuous and discontinuous boundary loads, which are of particular interest in understanding how changes in intraocular pressure would impact the bio-fluid-mechanics of ocular tissues.
The main conclusion of our analysis is that the compressibility of the components of a poroviscoelastic mixture does not give rise to unbounded responses to non-smooth traction data. Interestingly, compressibility provides an additional mechanism that limits the maximum of the discharge velocity when the imposed boundary traction is irregular in time. This mechanism originates from the capability of the system to store potential energy as its components are elastically compressed, thereby delaying the transmission of irregularities in the linear momentum from the solid to the fluid. As a result, the fluid has the time to accommodate for sudden changes, resulting in bounded velocities.
Our results fit well with other poroviscoelastic studies motivated by geomechanical applications, where viscoelasticity was found to play a crucial role in the response of the medium to impulsive loads. Specifically, the studies focused on evaluating the consequences that different choices for the viscoelastic models would have on the medium response to external loads. For example, in [12], Schanz and Cheng considered a Kelvin-Voigt viscoelastic model and investigated the consequences of adopting it for the bulk compression modulus, the shear modulus and the compression modulus of the solid material. In [2,3], Huang et al. utilized the quasi-linear viscoelastic theory to study the response of articular cartilage under compression and tension experiments. However, these studies did not consider how the mechanical responses of the poroviscoelastic mixture would change depending on the level of compressibility of the mixture components, which is the main focus of the present work.
The outline of this article is as follows. In Section 2 we describe the compressible poro-visco-elastic model under consideration and discuss the finite compressibility of the components of the mixture. Section 3 introduces the energy identity associated with the system. In Section 4 we present an equivalent form of the fluid-solid mixture system, written solely in terms of the solid phase displacement. To simplify the theoretical analysis, in Section 5 we reduce the poroviscoelastic model into dimensionless form. Section 6 studies the well-posedness and regularity of solution for the IBVP introduced in Section 4, and provides analytical formulas for the elastic displacement, fluid pressure and discharge velocity. In Section 7 we compute and analyze the behavior of the solid displacements, the fluid pressures and the discharge velocities associated with different continuous and discontinuous boundary sources. We also display the energies, as well as the dissipation and source terms associated with the various cases. We conclude the article with Section 9, where we discuss our results and draw our final conclusions.
In this paper, we focus on a viscoelastic, compressible Biot model in one spatial dimension. Let x and t denote the spatial and temporal coordinates, respectively. In the case where inertial terms are negligible, displacements are infinitesimal and external sources of linear momentum and mass are absent, the balance equations to be solved in the spatial interval (0,L) and in the temporal observational interval (0,T] can be written as:
∂σ∂x=0, | (2.1a) |
∂ζ∂t+∂v∂x=0, | (2.1b) |
where σ is the total stress, ζ is the fluid content and v is the discharge velocity. Equations (2.1) are complemented with the following constitutive laws:
σ=−αp+σ0, | (2.2a) |
σ0=θ∂u∂x+η∂∂t(∂u∂x), | (2.2b) |
ζ=c0p+α∂u∂x, | (2.2c) |
v=−K∂p∂x. | (2.2d) |
Equation (2.2a) expresses the fact that the total stress is the sum of a contribution due to the interstitial fluid pressure (or pore pressure) p and one due to the effective stress σ0, which is assumed to be characterized by viscoelastic behavior of Kelvin-Voigt type (see Eq (2.2b) where u denotes the solid phase displacement). Equation (2.2c) expresses the fact that the fluid content is altered by changes in fluid pressure and solid deformation. Equation (2.2d) is Darcy's law relating the discharge velocity v with the pressure gradient by means of the hydraulic permeability K. The Biot coefficient α, the storage coefficient c0, the material parameters θ and η and the permeability K are assumed to be given positive constants.
In the classical Biot's theory, the material parameter θ can be expressed as θ=K−(4/3)G, where K is the drained bulk compression modulus and G is the shear modulus. In the following, we will simply refer to θ and η as the elastic and viscoelastic parameters, respectively [12]. In addition, we will use the notation K=K0 to emphasize the fact that the permeability is assumed to be a given constant.
The problem must be equipped with appropriate boundary conditions. In the following, we will assume that the boundary located at x=0 is fixed and impermeable, namely:
u(0,t)=0, | (2.3a) |
v(0,t)=0, | (2.3b) |
and that the boundary located at x=L experiences an external stress P (a compressive stress if P is positive, a tensile stress if P is negative) that is supported entirely by the solid component of the mixture (condition of exposed pores [13]), namely:
p(L,t)=0, | (2.3c) |
σ(L,t)=−P(t). | (2.3d) |
Finally, we complete the formulation of the problem by prescribing the following initial conditions:
u(x,0)=0, | (2.4a) |
p(x,0)=0. | (2.4b) |
We notice that (2.4a) and (2.4b) also imply the following initial conditions for the dilation of the solid material and the discharge velocity, respectively:
∂u∂x(x,0)=0, | (2.5a) |
v(x,0)=0. | (2.5b) |
Remark 1. The case of incompressible mixture components can be obtained by setting c0=0 and α=1 in the model described above, as detailed in [14,15]. The study of analytical solutions for the incompressible model was addressed in [1,4].
The mathematical system described in Section 2 satisfies an energy identity that helps provide a physical interpretation of the solutions, assuming they exist. To this end, let us multiply (2.1a) by ∂u/∂t∈L2(0,L); integrating over (0,L) and using the constitutive laws (2.2a) and (2.2b) and the boundary conditions (2.3a), (2.3c) and (2.3d), we obtain
−θ2ddt‖∂u∂x‖2L2(0,L)−η‖∂2u∂t∂x‖2L2(0,L)+α∫L0p∂2u∂t∂xdx=P(t)∂u∂t(L,t). | (3.6) |
Multiplying (2.1b) by p∈L2(0,L), integrating over (0,L) and using the constitutive laws (2.2c) and (2.2d) and the boundary conditions (2.3b) and (2.3c), we obtain
c02ddt‖p‖2L2(0,L)+K0‖∂p∂x‖2L2(0,L)+α∫L0p∂2u∂t∂xdx=0. | (3.7) |
Subtracting (3.6) from (3.7) we get the following evolution equation for the total energy stored in the poroviscoelastic system
ddtEtot(t)+D(t)=F(t),t∈(0,T], | (3.8) |
where the energy functional Etot=Etot(t), the dissipation functional D=D(t) and the force term F=F(t) are defined as:
Etot(t):=c02‖p(⋅,t)‖2L2(0,L)+θ2‖∂u∂x(⋅,t)‖2L2(0,L), | (3.9a) |
D(t):=K0‖∂p∂x(⋅,t)‖2L2(0,L)+η‖∂2u∂t∂x(⋅,t)‖2L2(0,L), | (3.9b) |
F(t):=−P(t)∂u∂t(L,t). | (3.9c) |
Remark 2. The physical units of Etot are Joules per unit area, namely J m−2. This is due to the fact that we are considering a one-dimensional problem in space and, as a consequence, all the problem variables are assumed to be constant on every plane perpendicular to the chosen direction x. Mathematically, Etot is obtained by integrating in x between 0 and L the energy density εtot defined as
εtot(x,t):=c02p(x,t)2+θ2|∂u∂x(x,t)|2. |
The units of εtot are J m−3. Analogously, the units of D and F are J m−2 s−1.
Since Etot≥0 and D≥0, in absence of forcing terms (i.e., P=0) the energy decreases in time. It is important to emphasize that the terms proportional to the storage coefficient c0 and the elastic parameter θ contribute to the total energy in the form of potential energy, so that we can write
Ec0(t):=c02‖p(⋅,t)‖2L2(0,L),Eθ(t):=θ2‖∂u∂x(⋅,t)‖2L2(0,L)withEtot=Ec0+Eθ. |
Conversely, the terms proportional to K0 and η contribute to dissipate energy via viscous effects within the fluid and solid components. We notice that F does not have a definite sign since it depends on the boundary terms. The energy identity (3.8) will be very useful in interpreting the results presented in Section 7.
Now we express problem (2.1)–(2.4) solely in terms of the solid displacement. Combining (2.1a) and (2.3d), we obtain that the total stress is given by
σ(x,t)=−P(t), for all x∈(0,L] and for all t∈(0,T], |
and, moreover, the fluid pressure p and the discharge velocity v can be written in terms of the solid displacement u as:
p(x,t)=1αP(t)+1ασ0(u(x,t)), | (4.10a) |
v(x,t)=−K0α∂∂xσ0(u(x,t)), | (4.10b) |
where σ0=σ0(u(x,t)) is given by (2.2b). Let us now derive the problem satisfied by u. Integrating Eq (2.1b), where ζ is given by (2.2c), over the space interval [0,x] and taking the boundary conditions (2.3) into account, yields
c0∫x0∂p∂tdy+α∂u∂t+v=0 |
Now we substitute p and v by expressions (4.10) and note that
∫x0σ0dy=θu+η∂u∂t |
by virtue of (2.3a). Then we obtain
c0xP′(t)+c0(θ∂u∂t+η∂2u∂t2)+α2∂u∂t−K0(θ∂2u∂x2+η∂3u∂t∂x2)=0. |
The boundary conditions are (2.3a) and σ0|x=L=−P(t) (coming from (2.3c) and (4.10a)). The initial conditions are (2.4a) and, on using (2.4b), (4.10a), (2.2b) and (2.5a),
η∂2u∂t∂x(x,0)=−P(0) |
or equivalently
∂u∂t(x,0)=−1ηxP(0)+u1 |
u1 being an arbitrary constant. Note that u1=0 if and only if the compatibility condition with (2.3a)
∂u∂t(0,0)=0 | (4.11) |
is satisfied. Summing up:
Find u=u(x,t) such that:
c0η∂2u∂t2+(α2+c0θ)∂u∂t−K0θ∂2u∂x2−K0η∂3u∂t∂x2=−c0xP′(t)in (0,L)×(0,T], | (4.12a) |
u(0,t)=0in (0,T], | (4.12b) |
θ∂u∂x(L,t)+η∂2u∂t∂x(L,t)=−P(t)in (0,T], | (4.12c) |
u(x,0)=0in (0,L), | (4.12d) |
∂u∂t(x,0)=−1ηxP(0)+u1in (0,L). | (4.12e) |
Remark 3. The solution u(x,t) of the problem (4.12) is the sum u1(x,t)+u0(x,t) of the solution u1(x,t) of (4.12) where P=0, and of the solution u0(x,t) of (4.12) where u1=0. In the first case, a perturbation is generated at time zero, with a certain initial velocity of propagation (u1≠0): then u1(x,t) measures how the solid displacement changes through the medium when no disturbance is created at the boundary (P=0). On the other hand, if there is no ''initial impulse" (i.e., u1=0) and a disturbance is generated at the boundary (P≠0), then the corresponding change of the solid displacement is represented by u0(x,t). Therefore, the compatibility condition (4.11) simply means to consider the response of the system to the sole external stress P.
Remark 4. We assume throughout the article that c0>0 and η>0, i.e. the system is characterized by finite compressibility and structural viscoelasticity.
Remark 5. The IBVP (4.12) has a markedly different character compared to the linear poroviscoelastic system studied in [4] because of the second-order time derivative on the left-hand side of (4.12a).
Remark 6. If we define the following quantities:
˜ρ:=1K0c0ηα2+c0θ,˜σ:=1α2+c0θσ0−1K0∫x0∂u∂t(y,t)dy,˜f:=−1K0c0α2+c0θxP′(t), |
then the partial differential equation (4.12a) that describes the dynamics of the solid phase displacement can be written as
˜ρ∂2u∂t2=∂˜σ∂x+˜f. |
Such an equation can be formally interpreted as a linear momentum balance equation for a single phase solid material whose dynamics is equivalent to that of the fluid-solid mixture under the assumptions of Section 2. In particular, we see that the finite compressibility of the mixture components, corresponding to c0>0, provides:
(i) an inertia-like term for the equivalent solid, even though inertial terms were neglected for the original solid phase within the mixture;
(ii) an additional term in the stress tensor introducing nonlocal effects in space;
(iii) a volumetric forcing term that results from the load applied as a boundary condition.
In order to simplify the theoretical analysis, in this section we reduce the 1D model (4.12) into dimensionless form. With this aim, for any variable Y, we define the corresponding non-dimensional variable by
ˆY:=Y[Y], |
where [Y] is a suitably chosen scaling factor that has the same units as Y. The selection of the scaling factor is not unique and, in general, not trivial. In this article we generalize the choice made in [4], by introducing the following scaling factors:
[x]=L, | (5.1a) |
[t]=L2(α2+c0θ)θK0, | (5.1b) |
[η]=θ[t], | (5.1c) |
[σ]=[σ0]=[P]=Pref, | (5.1d) |
[u]=LPrefθ, | (5.1e) |
[u1]=[u][t], | (5.1f) |
[v]=K0PrefαL, | (5.1g) |
[p]=Prefα, | (5.1h) |
[D]=[F]=K0P2refL(α2+c0θ), | (5.1i) |
[Etot]=[Ec0]=[Eθ]=[D][t]=LP2refθ. | (5.1j) |
We also define the non-dimensional quantity
γ:=c0θα2+c0θ. | (5.2) |
Note that 0<γ<1. We recover the same definitions of the scaling factors as in [4] by setting α=1 and c0=0 in (5.1).
Remark 7. For notational simplicity we will drop the 'ˆ⋅' notation for the dimensionless variables and instead use the same symbols adopted for the dimensional variables.
The linear 1D model for a poroviscoelastic mixture in dimensionless form reads:
γη∂2u∂t2+∂u∂t−∂2u∂x2−η∂3u∂t∂x2=−γxP′(t)in (0,1)×(0,T], | (5.3a) |
u(0,t)=0in (0,T], | (5.3b) |
∂u∂x(1,t)+η∂2u∂t∂x(1,t)=−P(t)in (0,T], | (5.3c) |
u(x,0)=0in (0,1), | (5.3d) |
∂u∂t(x,0)=−1ηxP(0)+u1in (0,1). | (5.3e) |
Once u(x,t) is known, the functions σ, p and v can be computed as follows:
σ(x,t)=−P(t), | (5.4a) |
p(x,t)=P(t)+∂u∂x+η∂2u∂t∂x, | (5.4b) |
v(x,t)=−(∂2u∂x2+η∂3u∂t∂x2)=−(γη∂2u∂t2+∂u∂t+γxP′(t)). | (5.4c) |
Lastly, the dimensionless energy equation is written again in the form (3.8) where
Etot(t)=Ec0(t)+Eθ(t), | (5.5a) |
Ec0(t)=12γ1−γ‖p(⋅,t)‖2L2(0,1),Eθ(t)=12‖∂u∂x(⋅,t)‖2L2(0,1), | (5.5b) |
D(t)=11−γ‖∂p∂x(⋅,t)‖2L2(0,1)+η‖∂2u∂t∂x(⋅,t)‖2L2(0,1), | (5.5c) |
F(t)=−P(t)∂u∂t(1,t). | (5.5d) |
We make the following assumption on the boundary traction in (5.3).
Assumption 8. P(t) is a piecewise smooth function on [0,T].
We recall that a function P(t) is piecewise smooth on [0,T] if both P and its derivative P′ are continuous on [0,T], except possibly at a finite number of points in (0,T), where they have simple jump discontinuities.
Define
U(t)=1ηe−t/η∗P(t)=1η∫t0exp(−t−sη)P(s)ds. | (6.6) |
Using Assumption 8 we have that U(t) is absolutely continuous on [0,T] and
ηU′(t)=P(t)−U(t),U(0)=0,U′(0)=P(0)η. | (6.7) |
Now we introduce the change of variable
w(x,t)=u(x,t)+xU(t), | (6.8) |
and note that w solves the following auxiliary IBVP with homogeneous boundary data:
γη∂2w∂t2+∂w∂t−∂2w∂x2−η∂3w∂t∂x2=f(x,t)in (0,1)×(0,T], | (6.9a) |
w(0,t)=0in (0,T], | (6.9b) |
∂w∂x(1,t)+η∂2w∂t∂x(1,t)=0in (0,T], | (6.9c) |
w(x,0)=0in (0,1), | (6.9d) |
∂w∂t(x,0)=φ(x)in (0,1). | (6.9e) |
where the volumetric source and initial datum are given by:
f(x,t)=(1−γ)xU′(t), | (6.10a) |
φ(x)=u1. | (6.10b) |
We shall prove the existence and uniqueness of the solution w(x,t) for a very general class of data f(x,t) and φ(x). For sake of exposition we write H=L2(0,1), and define the real Hilbert space
V={v∈W1,2(0,1):v(0)=0} | (6.11) |
endowed with the equivalent norm ‖v‖V=‖∂v/∂x‖L2(0,1), due to Poincaré's inequality.
Remark 9. Sobolev's Embedding Theorem gives W1,2(0,1)⊂C0[0,1] so that v(0)=0 holds in a strong sense for every v∈V.
Now we make the following assumptions on the functions f(x,t) and φ(x):
f∈L2(0,T;H),φ∈H | (6.12) |
and write w(t)=w(⋅,t), w′(t)=∂w(⋅,t)/∂t, etc. We then define weak solutions of (6.9) as follows.
Definition 10. [Weak solution of problem (6.9)] A function w:[0,T]→V is a weak solution of the auxiliary problem (6.9) if:
D1 w∈W1,2(0,T;V) and w′′∈L2(0,T;V′);
D2 for every v∈V and for t pointwise a.e. in (0,T)
γη⟨w′′(t),v⟩V′×V+(w′(t),v)H+(ηw′(t)+w(t),v)V=(f(t),v)H | (6.13) |
or, equivalently,
ddt(γηw′(t)+w(t),v)H+(ηw′(t)+w(t),v)V=(f(t),v)H∀v∈V; | (6.14) |
D3 w(0)=0 and w′(0)=φ.
Remark 11. Condition [D1] implies that w∈C0([0,T];V) and w′∈C0([0,T];H), and thus condition [D3] is well defined. The Dirichlet boundary condition in (6.9) at x=0 is included in the regularity requirement that w(t)∈V, whereas the boundary condition at x=1 is satisfied in a weak sense.
Lemma 12. Let w be a weak solution of (6.9).Then there are constants Ci's, depending only on γ, η and T, such that the following estimates hold for t pointwise in [0,T]:
‖γηw′(t)+w(t)‖2H≤C1(‖φ‖2H+‖f‖2L2(0,T;H)), | (6.15a) |
‖w(t)‖2V≤C2(‖φ‖2H+‖f‖2L2(0,T;H)), | (6.15b) |
‖w′(t)‖2H≤C3(‖φ‖2H+‖f‖2L2(0,T;H)). | (6.15c) |
Proof. Using (γηw′+w)∈L2(0,T;V) as multiplier in (6.14), we get
ddt‖γηw′+w‖2H+(ηw′+w,γηw′+w)V=(f,γηw′+w)H. |
Integrating in time over (0,t) and using the given initial conditions, we obtain
‖γηw′+w‖2H+η(γ+1)2‖w‖2V+∫t0(‖w‖2V+γη2‖w′‖2V) ds=γ2η2‖φ‖2H+∫t0(f,γηw′+w)H ds. |
Using Cauchy-Schwarz and Young's Inequalities for the last term on the right-hand side, we obtain the following estimate:
‖γηw′+w‖2H+η(γ+1)2‖w‖2V+∫t0(‖w‖2V+γη2‖w′‖2V) ds≤γ2η2‖φ‖2H+12‖f‖2L2(0,T;H)+12∫t0‖γηw′+w‖2H ds. | (6.16) |
Now, dropping the second and third terms on the left-hand side of (6.16) and using Gronwall's Inequality, estimate (6.15a) follows. Similarly, by dropping the first and third terms on the left-hand side and using (6.15a), we get (6.15b).
Lastly, we use triangle inequality and the embedding V↪H to write
γη‖w′‖H≤‖γηw′+w‖H+‖w‖H≤‖γηw′(t)+w(t)‖H+1√2‖w(t)‖V. |
Then, estimate (6.15c) follows from (6.15a) and (6.15b).
Lemma 13. Let w be a weak solution of (6.9).Then there is a constant C, depending only on γ, η and T, such that
‖w‖L∞(0,T;V)+‖w′‖L2(0,T;V)+‖w′′‖L2(0,T;V′)≤C(‖φ‖H+‖f‖L2(0,T;H)). | (6.17) |
Proof. From (6.15b) it immediately follows
‖w‖L∞(0,T;V)≤√C2(‖φ‖H+‖f‖L2(0,T;H)). |
In addition, Eqs (6.16) and (6.15a) give
γη2∫T0‖w′‖2V ds≤γ2η2‖φ‖2H+12‖f‖2L2(0,T;H)+12∫T0‖γηw′+w‖2H ds≤(γ2η2+T2C1)‖φ‖2H+12(1+TC1)‖f‖2L2(0,T;H) |
so that
‖w′‖L2(0,T;V)≤C(‖φ‖H+‖f‖L2(0,T;H)) | (6.18) |
for a suitable C=C(γ,η,T).
Lastly, using the definition of a weak solution, we have that for v∈V
γη|⟨w′′,v⟩V′×V|≤‖w′‖H‖v‖H+‖ηw′+w‖V‖v‖V+‖f‖H‖v‖H≤1√2(‖w′‖V+‖ηw′+w‖V+‖f‖H)‖v‖V. |
This implies that
‖w′′(t)‖V′≤η+1√2(‖w′‖V+‖w‖V+‖f‖H) |
from which, using estimates (6.15b) and (6.15a), we obtain
‖w′′‖L2(0,T;V′)≤C(‖φ‖H+‖f‖L2(0,T;H)) |
for a suitable C=C(γ,η,T).
The following corollary is an immediate consequence of Lemma 13.
Corollary 14. (Uniqueness and continuous dependence on data) The weak solution to problem (6.9) is unique and depends continuously on the data.
The IBVP (6.9) can be solved formally using separation of variables. If we look for solutions of the form w(x,t)=T(t)X(x), then the associated regular Sturm-Liouville Problem is
{X′′+λX=0,0<x<1X(0)=0, X′(1)=0, | (6.19) |
with eigenvalues and eigenfunctions given by:
λn=(nπ+π2)2, n≥0, | (6.20) |
Xn(x)=sin((nπ+π2)x), n≥0. | (6.21) |
Remark 15. The sequence of functions {√2Xn(x)} forms a Hilbert space basis for H, whereas the sequence of functions {√2λnXn(x)} forms a Hilbert space basis for V (see [17]).
The solution w of (6.9) has the expansion
w(x,t)=∞∑n=0Tn(t)Xn(x), | (6.22) |
where Tn(t) can be recovered using the data. Similarly, we use the basis {Xn(x)} to represent φ and f as follows:
φ(x)=∞∑n=0φnXn(x), | (6.23a) |
f(x,t)=∞∑n=0fn(t)Xn(x), | (6.23b) |
where φn and fn(t) are the Fourier coefficients of φ and f(⋅,t) with respect to Xn(x), respectively. Parseval's identity (consequence of Remark 15) provides the following relations:
‖φ‖2H=12∞∑n=0φ2n, | (6.24) |
‖f(⋅,t)‖2H=12∞∑n=0|fn(t)|2, | (6.25) |
‖f‖2L2(0,T;H)=12∫T0∞∑n=0|fn(t)|2dt. | (6.26) |
Note that Tn(t) satisfies the following ordinary differential equation (ODE) for all n≥0:
γηT′′n(t)+(1+ηλn)T′n(t)+λnT(t)=fn(t), | (6.27) |
and initial conditions
Tn(0)=0andT′n(0)=φn. | (6.28) |
The characteristic equation for the homogeneous counterpart of (6.27) is given by
γηΛ2+(1+ηλn)Λ+λn=0. |
Since the discriminant of the equation is
(1+ηλn)2−4γηλn=(1−ηλn)2+4(1−γ)ηλn>0, |
then for each n≥0 the characteristic equation has two real, negative, distinct roots r1n=−Λ1n and r2n=−Λ2n, where:
Λ1n=12γη(1+ηλn−√(1+ηλn)2−4γηλn), and | (6.29a) |
Λ2n=12γη(1+ηλn+√(1+ηλn)2−4γηλn). | (6.29b) |
Remark 16. Λ1n and Λ2n satisfy the following relations:
0<Λ1n<Λ2nΛ1n+Λ2n=1+ηλnγη,Λ1nΛ2n=λnγη,1−γ=−γ(ηΛ1n−1)(ηΛ2n−1)Λ1n=1η−(1−γη2)1λn+O(1λ2n),Λ2n=λnγ+1−γγη+O(1λn)(as n→∞)Λ1n=λn1+ηλn+O(γ),Λ2n=1+ηλnγη−λn1+ηλn+O(γ)(as γ→0) |
Therefore the solution for the homogeneous ODE (6.27) is given by T0n(t)=ane−Λ1nt+bne−Λ2nt. The particular solution is obtained from variation of parameters formula. The Wronskian is given by W(t)=(Λ2n−Λ1n)e−Λ1nte−Λ2nt, so that the particular solution has the following form
Tpn(t)=1γη∫t0e−Λ1n(t−s)−e−Λ2n(t−s)Λ2n−Λ1n f(s) ds. | (6.30) |
We introduce the following notation
Gn(t):=exp(−Λ1nt)−exp(−Λ2nt)Λ2n−Λ1n | (6.31) |
and then write the solution to (6.27) as
Tn(t)=ane−Λ1nt+bne−Λ2nt+1γη(Gn∗fn)(t) | (6.32) |
where we used the following formula for convolution
(Gn∗fn)(t)=∫t0Gn(t−s)fn(s)ds. |
Now we use (6.32) back into (6.22) and impose the initial conditions. We have:
w(x,0)=an+bn=0, andwt(x,0)=−Λ1nan−Λ2nbn=φn |
which yields an=−bn=φn/(Λ2n−Λ1n). In conclusion, we obtain:
Tn(t)=Gn(t)φn+1γη(Gn∗fn)(t), and | (6.33a) |
w(x,t)=∞∑n=0[Gn(t)φn+1γη(Gn∗fn)(t)]Xn(x). | (6.33b) |
We note that the terms Gn(t) given in (6.31) satisfy the following estimates.
Lemma 17. There exists a constant C>0 such that for all n≥0 and for all t∈[0,T] we have:
0≤Gn(t)≤Cλn, | (6.34a) |
|G′n(t)|≤C. | (6.34b) |
Proof. Since 0<Λ1n<Λ2n, then 0<exp(−Λ2nt)≤exp(−Λ1nt)≤1 and thus for all t∈[0,T] we have
0≤Gn(t)<1Λ2n−Λ1n |
and we get (6.34a) since the sequence λnΛ2n−Λ1n=γηλn√(ηλn+1)2−4γηλn→γ as n→∞. Moreover, we have that for all t∈[0,T]
|G′n(t)|=Λ2nexp(−Λ2nt)−Λ1nexp(−Λ1nt)Λ2n−Λ1n≤2Λ2nΛ2n−Λ1n |
and (6.34b) follows since Λ2nΛ2n−Λ1n→1 as n→∞.
Now we can state and prove our well-posedness result.
Theorem 18. (Well-posedness of problem (6.9)) For every f∈L2(0,T;H) and φ∈H there is a unique weak solution of (6.9), in the sense of Definition 10. Moreover, the solution depends continuously on the data.
Proof. Uniqueness and continuous dependence of weak solution have already been proved, see Corollary 14, so that it remains to prove existence, i.e. that w(x,t) given in (6.33b) satisfies conditions (D1)–(D3) of Definition 10.
(D1) (A) First we show that w∈L2(0,T;V). For all t∈[0,T], we have
‖w(t)‖2V=12∞∑n=0λnT2n(t)≤∞∑n=0λn|Gn(t)φn|2+∞∑n=0λn|1γη∫t0Gn(t−s)fn(s)ds|2. |
Using Lemma 17, we obtain*
*In what follows, for notational convenience, C will denote possibly different constants depending only on γ,η,T.
‖w(t)‖2V≤C∞∑n=0|φn|2λn+C∞∑n=01λn(∫t0|fn(s)|ds)2≤C∞∑n=0|φn|2+C∞∑n=0∫T0|fn(s)|2ds≤C(‖φ‖2H+‖f‖2L2(0,T;H)). | (6.35) |
Therefore w∈L∞(0,T;V)⊂L2(0,T;V).
(B) The weak derivative of w with respect to time is given by
w′(t)=∞∑n=0T′n(t)Xn(x)=∞∑n=0[G′n(t)φn+1γη(G′n∗fn)(t)]Xn(x), |
and we obtain the following estimate
‖w′(t)‖2H=12∞∑n=0(T′n(t))2≤∞∑n=0|G′n(t)φn|2+∞∑n=0|∫t0G′n(t−s)fn(s)ds|2. |
Again by Lemma 17, we obtain
‖w′(t)‖2H≤C∞∑n=0|φn|2+C∞∑n=0(∫t0|fn(s)|ds)2≤C∞∑n=0|φn|2+C∞∑n=0∫T0|fn(s)|2ds≤C(‖φ‖2H+‖f‖2L2(0,T;H)). | (6.36) |
Thus w′∈L∞(0,T;H)⊂L2(0,T;H). In order to show that w′∈L2(0,T;V), we first use the Monotone Convergence Theorem to write
‖w′‖2L2(0,T;V)=12∫T0∞∑n=0λn|T′n(s)|2ds=12∞∑n=0∫T0λn|T′n(s)|2ds. | (6.37) |
We use multiplier T′n in (6.27) and the initial conditions (6.28) and obtain the following estimate:
γηT′nT′′n+(1+ηλn)T′2n+λnTnT′n=fnT′n ⇒ |
γη2ddt(T′n)2+(T′n)2+ηλn(T′n)2+λn2ddt(Tn)2=fnT′n ⇒ |
γη2(T′n(t))2−γη2φ2n+∫t0(T′n)2 ds+∫t0ηλn(T′n)2 ds+λn2(Tn)2=∫t0fnT′n ds ⇒ |
∫t0ηλn(T′n)2 ds≤γη2φ2n++12∫T0|fn(s)|2 ds+12∫T0|T′n(s)|2 ds. | (6.38) |
Now we use (6.38) back into (6.37), and take advantage of estimate (6.36) to obtain
‖w′‖2L2(0,T;V)≤C∞∑n=0φ2n+C∞∑n=0∫T0|fn(s)|2ds+C∞∑n=0∫T0|T′n(s)|2ds=C{‖φ‖2H+‖f‖2L2(0,T;H)+‖w′‖2L2(0,T;H)}≤C{‖φ‖2H+‖f‖2L2(0,T;H)}. | (6.39) |
and this gives the desired conclusion that w′∈L2(0,T;V).
(C) Left to show that w′′∈L2(0,T;V′). For any v(x)=∑∞n=0cnXn(x)∈V and for t∈[0,T], we use (6.27) to define the linear functional w′′(t) as follows
⟨w′′(t),v⟩V′×V=12∞∑n=0cnT′′n(t)=−12γη∞∑n=0T′n(t)cn−12γη∞∑n=0λn(ηT′n(t)+Tn(t))cn+12γη∞∑n=0fn(t)cn=−1γη(w′(t),v)H−1γη(ηw′(t)+w(t),v)V+1γη(f(t),v)H. |
Hence, by virtue of the embedding V↪H, we obtain that
|⟨γηw′′(t),v⟩V′×V|≤‖w′(t)‖H‖v‖H+‖ηw′(t)+w(t)‖V‖v‖V+‖f(t)‖H‖v‖H≤C{‖w′(t)‖V+‖w(t)‖V+‖f(t)‖H}‖v‖V. |
This shows that w′′(t)∈V′ and
‖w′′(t)‖V′≤C(‖w′(t)‖V+‖w(t)‖V+‖f(t)‖H) |
or, equivalently,
‖w′′(t)‖2V′≤C(‖w′(t)‖2V+‖w(t)‖2V+‖f(t)‖2H). |
In conclusion, using the estimates in part (A) and part (B), we obtain that w′′∈L2(0,T;V′).
(D2) Now we show that w satisfies condition (D2). Since the {√2λnXn} is a basis in V, it suffices to consider the test function v=Xn. For t pointwise almost everywhere in (0,T), we use (6.27) and obtain
γη⟨w′′(t),v⟩V′×V+(w′(t),v)H+(ηw′(t)+w(t),v)V=γη2T′′n(t)+12T′n(t)+ηλn2T′n(t)+λn2Tn(t)=12fn(t)=(f(t),v)H. |
(D3) As stated in Remark 11, condition (D1) implies that w∈C0([0,T];V) and w′∈C0([0,T];H). Moreover, solution w(x,t) satisfies the initial conditions (D3) of Definition 10 by virtue of (6.28). This concludes the proof of our well-posedness theorem.
We have established the existence and uniqueness of the weak solution w to the auxiliary problem (6.9). Now we examine its regularity.
Proposition 19. (Regularity of w and wxx) Under the hypotheses in Theorem 18, the following is true:
w(x,t)∈C0(Q), where Q=[0,1]×[0,T], | (6.40a) |
wxx∈L∞(0,T;H). | (6.40b) |
Proof. (A) The result follows from Lemma 17, the fact that |Xn(x)|≤1 and Weierstrass Test. Indeed, for all x∈[0,1] and t∈[0,T], and every n≥0, we have
|Tn(t)Xn|≤|Gn(t)φnXn(x)|+1γη|(Gn∗fn)(t)Xn(x)|. |
Now, due to Lemma 17, we have
∞∑n=0|Gn(t)φnXn(x)|≤∞∑n=0|Gn(t)φn|≤C∞∑n=01λn|φn|≤C∞∑n=0(1λ2n+φ2n)=C(∞∑n=01λ2n+‖φ‖2H)<∞, |
and
∞∑n=0|(Gn∗fn)(t)Xn(x)|≤∞∑n=0|(Gn∗fn)(t)|≤∞∑n=0∫t0|Gn(t−s)fn(s)|ds≤C∞∑n=0∫T01λn|fn(s)|ds≤C∞∑n=0∫T0(1λ2n+|fn(s)|2)ds≤C(∞∑n=01λ2n+‖f‖2L2(0,T;H))<∞. |
Applying Weierstrass Test, we obtain that the series ∑∞n=0 Tn(t)Xn(x) converges absolutely and uniformly in Q=[0,1]×[0,T]. Moreover, we note that Gn(t)φnXn(x)+1γη(Gn∗fn)(t)Xn(x)∈C0(Q), for every n≥0, and thus w(x,t)∈C0(Q).
(B) The second order weak derivative in space of w is given by
∂2w∂x2=−∞∑n=0Gn(t)φnλnXn(x)−∞∑n=0(Gn∗fn)(t)λnXn(x). |
Then we use estimate (6.34a) in Lemma 17 and obtain
‖∂2w∂x2(⋅,t)‖2H≤∞∑n=0|Gn(t)φnλn|2+∞∑n=0|(Gn∗fn)(t)λn|2≤C∞∑n=0φ2n+C∞∑n=0|∫t0|fn(s)|ds|2≤C∞∑n=0φ2n+C∞∑n=0∫t0|fn(s)|2ds=C(‖φ‖2H+‖f‖2L2(0,T;H)) |
so that wxx∈L∞(0,T;H).
If the data are more regular with respect to x the weak solution w(x,t) enjoys stronger regularity properties.
Proposition 20. (Regularity of wt, wxx and wtxx) In addition to the hypotheses in Theorem 18, we make the following ones:
f∈L2(0,T;V)andφ∈V. |
Then the following is true:
wt(x,t)∈C0(Q), where Q=[0,1]×[0,T], | (6.41a) |
wxx∈L∞(0,T;V), | (6.41b) |
wtxx∈L2(0,T;H). | (6.41c) |
Proof. (A) For all x∈[0,1] and t∈[0,T], and every n≥0, we have
|T′n(t)Xn|≤|G′n(t)φnXn(x)|+1γη|(G′n∗fn)(t)Xn(x)|. |
Now, due to Lemma 17, we have
∞∑n=0|G′n(t)φnXn(x)|≤∞∑n=0|G′n(t)φn|≤C∞∑n=0|φn|≤C∞∑n=0(1λn+λnφ2n)=C(∞∑n=01λn+‖φ‖2V)<∞, |
and
∞∑n=0|(G′n∗fn)(t)Xn(x)|≤∞∑n=0|(G′n∗fn)(t)|≤∞∑n=0∫t0|G′n(t−s)fn(s)|ds≤C∞∑n=0∫T0|fn(s)|ds≤C∞∑n=0∫T0(1λn+λn|fn(s)|2)ds≤C(∞∑n=01λn+‖f‖2L2(0,T;V))<∞. |
Applying Weierstrass Test, we obtain that the series ∑∞n=0 T′n(t)Xn(x) converges absolutely and uniformly in Q=[0,1]×[0,T] hence wt(x,t)∈C0(Q).
(B) Like in the proof of (6.40a), we have
‖∂2w∂x2(⋅,t)‖2V≤∞∑n=0λn|Gn(t)φnλn|2+∞∑n=0λn|(Gn∗fn)(t)λn|2≤C∞∑n=0λnφ2n+C∞∑n=0λn|∫t0|fn(s)|ds|2≤C∞∑n=0λnφ2n+C∞∑n=0∫t0λn|fn(s)|2ds=C(‖φ‖2V+‖f‖2L2(0,T;V)) |
so that wxx∈L∞(0,T;V).
(C) Since the second order weak derivative in space of w′ is given by
∂3w∂t∂x2(x,t)=−∞∑n=0T′n(t)λnXn(x) |
then
‖∂3w∂t∂x2‖2L2(0,T;H)=12∫T0∞∑n=0λ2n|T′n(s)|2ds | (6.42) |
In order to estimate the right hand side, we use multiplier λnT′n on the ODE (6.27) that Tn solves to obtain
ηλ2nT′2n=−λnT′2n−12(γηλnT′2n+λ2nT2n)′+λnfn(t)T′n≤−12(γηλnT′2n+λ2nT2n)′+λnfn(t)T′n. |
Now we integrate from 0 to t and use the initial conditions (6.28), to get
η∫t0λ2n|T′n(s)|2ds≤−12(γηλn|T′n(t)|2+λ2n|Tn(t)|2)+12γηλnφ2n+∫t0λnfn(s)T′n(s)ds≤12(γηλnφ2n+∫t0λn|fn(s)|2ds+∫t0λn|T′n(s)|2ds). | (6.43) |
Thus, using (6.43) back into (6.42), we obtain
‖∂3w∂t∂x2‖2L2(0,T;H)≤C(∞∑n=0λnφ2n+∫T0∞∑n=0λn|fn(t)|2dt+∫T0∞∑n=0λn|T′n(t)|2dt)=C(‖φ‖2V+‖f‖2L2(0,T;V)+‖w′‖2L2(0,T;V)) |
and the assertion follows from estimate (6.17).
Now we are in a position to return to our original linear 1D problem (5.3). The solid displacement solution u(x,t) of (5.3) is the sum
u(x,t)=−xU(t)+w(x,t) | (6.44) |
where we recall that
U(t)=1ηe−t/η∗P(t) |
and w(x,t) is the unique weak solution of the auxiliary problem (6.9) with special data (6.10). From the Fourier expansions
∞∑n=01√λnXn(x)=12,∞∑n=0(−1)nλnXn(x)=x2,0≤x≤1, |
as well as the Fourier series (6.23) of φ and f with respect to the basis {Xn}, their Fourier coefficients are given by:
φn=2u1√λn, | (6.45a) |
fn(t)=2(−1)nλn(1−γ)U′(t). | (6.45b) |
To summarize, we can write the solution u as a sum
u(x,t)=u1(x,t)+u0(x,t), | (6.46) |
where:
u1(x,t)=u1∞∑n=02√λnGn(t)Xn(x), | (6.47a) |
u0(x,t)=−xU(t)+1−γγη∞∑n=02(−1)nλnGn(t)∗U′(t)Xn(x). | (6.47b) |
Remark 21. The solid displacement solution u(x,t)∈C0(Q) since the physical data (6.10) satisfy f∈ L2(0,T;V) and φ∈H. Notice, however, that our special φ belongs to V if and only if u1=0, i.e. φ=0 and u1(x,t)≡0.
By virtue of (6.46), the discharge velocity (5.4c) is given by
v(x,t)=v1(x,t)+v0(x,t) | (6.48) |
where:
v1(x,t)=u1∞∑n=02√λn(Gn(t)+ηG′n(t))Xn(x), | (6.49a) |
v0(x,t)=1−γγη∞∑n=02(−1)n[Gn(t)∗U′(t)+η(Gn(t)∗U′(t))′]Xn(x). | (6.49b) |
Remark 22. It should be stressed that the Fourier coefficients of v1(x,t) are O(n), hence the component of the discharge velocity due to the presence of an initial impulse u1 always exibits a blow-up whatever boundary traction P(t) is considered.
In this section, we analyze the behavior of the solutions obtained for model (5.3) in the case where the compatibility conditions are satisfied (i.e., u1=0), and the boundary traction P(t) is characterized by continuous or discontinuous waveforms. To graphically represent the model solutions we proceed as follows: (1) we define the numerical values of model parameters (in dimensional form) consistently with the experimental data illustrated in [16]; (2) we perform the non-dimensionalization of the model equations according to the procedure described in Section 5; (3) we compute the model solutions in non-dimensional form; (4) we multiply the non-dimensional input data and model solutions by the scaling factors introduced in (5.1) and plot the obtained results for subsequent analysis. For each considered waveform of P(t), we provide the non-dimensional expressions of the solutions u, p and v. These latter expressions allow us to compute the non-dimensional energies Eθ, Ec0 and Etot, the non-dimensional dissipation functional D and the non-dimensional force term F using the expressions (5.5). The resulting non-dimensional energies are then multiplied by the scaling factor [Etot] introduced in (5.1j) and the same is done for the non-dimensional dissipation and force terms that are multiplied by the scaling factor [D] introduced in (5.1i).
Table 1 reports the numerical values of model parameters (in dimensional form) that are used in the next sections. The values for L, T, θ, K0 and Pref have been chosen consistently with the experimental data illustrated in [16]. To observe the asymptotic behavior of the modeled system, T has been taken equal to 105s in the example illustrated in Section 7.1. In this example, the corresponding value of η is 4.85⋅109Nsm−2. The value of the Biot coefficient α has been set equal to 0.9, whereas in [16] it was equal to 1 since both mixture components were assumed to be incompressible. The values of the viscoelastic parameter η and of the compressibility parameter c0 (which were not present in the model studied in [16]) have been set equal to θτe and 1/Pref, respectively, where τe is an elastic time constant that has been set equal to T/20.
symbol | value | units |
L | 0.81⋅10−3 | m |
T | 105,104 | s |
θ | 0.97⋅106 | Nm−2 |
K0 | 2.9⋅10−16 | m4N−1s−1 |
Pref | 6⋅104 | Nm−2 |
c0 | 1.67⋅10−5 | m2N−1 |
η | 4.85⋅109,4.85⋅108 | Nsm−2 |
The numerical values of the parameters (in dimensionless form) that are used in the computations discussed in the next sections are reported in Table 2. These values have been obtained by applying the scaling procedure described in Section 5 to the values in Table 1. With a slight abuse of notation, the symbols used to denote the dimensionless parameters are the same ones that we used for the corresponding parameters expressed with their physical units.
symbol | value |
L | 1 |
T | 2.5, 0.25 |
η | 1.26⋅10−1, 1.26⋅10−2 |
γ | 0.95 |
For convenience, we recall here the formulas that we use to compute the solid displacement, the fluid pressure and the discharge velocity in dimensionless form:
u0(x,t)=−xU(t)+1−γγη∞∑n=02(−1)nλnGn(t)∗U′(t)Xn(x), | (7.1a) |
p0(x,t)=∞∑n=02(−1)nλnBn(t)X′n(x), | (7.1b) |
v0(x,t)=∞∑n=02(−1)nBn(t)Xn(x), | (7.1c) |
where
Bn(t)=1−γγη[(Gn∗U′)(t)+η((Gn∗U′)(t))′]. | (7.1d) |
Recalling formula (6.31) and using the fact that (Gn∗U′)′=G′n∗U′, we obtain the following simplification for the coefficients Bn(t)
Bn(t)=1−γγη(1−ηΛ1nΛ2n−Λ1ne−Λ1nt∗U′−1−ηΛ2nΛ2n−Λ1ne−Λ2nt∗U′). | (7.1e) |
Let t∗∈[0,T). We consider the following boundary source
P(t)=H(t−t∗)={0ift<t∗1ift≥t∗. | (7.2) |
Figure 1 gives a graphical representation of P(t) for t∗∈(0,T).
Replacing (7.2) into (6.6) we obtain
U(t)={0,t<t∗1−e−t−t∗η,t≥t∗ and U′(t)={0,t<t∗1ηe−t−t∗η,t≥t∗. | (7.3) |
We can now compute the convolution needed in (7.1a)
Gn(t)∗U′(t)={0,t<t∗,1η(Λ2n−Λ1n)[e−(t−t∗)/η−e−Λ1n(t−t∗)Λ1n−1/η−e−(t−t∗)/η−e−Λ2n(t−t∗)Λ2n−1/η],t≥t∗. | (7.4) |
We observe that for t≥t∗,
(Gn(t)∗U′(t))′=1η(Gn(t−t∗)−(Gn∗U′)(t)), |
and therefore the coefficients (7.1d) present in the expansions of pressure and discharge velocity become
Bn(t)=1−γγηGn(t−t∗), for any t≥t∗. |
Then displacement, pressure and discharge velocity are all zero for t<t∗, and have the following representations for t≥t∗:
u(x,t)=−x−∞∑n=02(−1)nλne−Λ2n(t−t∗)(ηΛ1n−1)−e−Λ1n(t−t∗)(ηΛ2n−1)η(Λ2n−Λ1n)Xn(x), | (7.5a) |
p(x,t)=1−γγη∞∑n=02(−1)nλnGn(t−t∗)X′n(x), | (7.5b) |
v(x,t)=1−γγη∞∑n=02(−1)nGn(t−t∗)Xn(x). | (7.5c) |
Note that all three series in the expressions (7.5) converge absolutely and uniformly on [0,1]×[t∗,T], and therefore p,v∈C0([0,1]×[t∗,T]). Moreover, we observe that p(x,t)→0 and v(x,t)→0, as t→t∗.
Remark 23. Note that if t∗=0, then this is the case of a continuous constant boundary source P(t)=1, for all t≥0. Then the solid displacement, the fluid pressure and the discharge velocity have the following representations for t≥0:
u(x,t)=−x−∞∑n=02(−1)nλne−Λ2nt(ηΛ1n−1)−e−Λ1nt(ηΛ2n−1)η(Λ2n−Λ1n)Xn(x), | (7.6a) |
p(x,t)=1−γγη∞∑n=02(−1)nλnGn(t)X′n(x), | (7.6b) |
v(x,t)=1−γγη∞∑n=02(−1)nGn(t)Xn(x). | (7.6c) |
All three series in (7.6) converge absolutely and uniformly on [0,1]×[0,T], and therefore p,v∈C0([0,1]×[0,T]).
Figure 2 illustrates displacement, fluid pressure and discharge velocity (left, middle and right panel, respectively) as a function of t. Displacement and velocity are evaluated at x=1 (right boundary) whereas the pressure is evaluated at x=0 (left boundary). Within the observational time interval, the displacement decreases monotonically till it reaches an asymptotic value of approximately 50μm. Conversely, pressure and discharge velocity exhibit a nonmonotonic behaviour, characterized by a rapid increase, attaining a maximum value of approximately 2.3kPa and 0.002μms−1, respectively, followed by a monotonic decrease that approaches zero asymptotically The behavior of the simulated solutions are consistent with those reported in [12,16].
Figure 3 illustrates displacement, fluid pressure and discharge velocity (left, middle and right panel, respectively) as a function of x and t. We see that, after an initial transient due to structural viscoelasticity, the solid displacement tends to a linear behavior along the domain length, being maximum (in absolute value) at x=L. Mathematically, this behavior is due to the fact that, for long times, the viscoelastic terms in the series in (7.6a) become negligible with respect to the linear elastic component. Both pressure and discharge velocity distributions decrease as time →+∞, in agreement with the fact that definition (6.31) implies that limt→+∞Gn(t)=0 for all n≥0.
Figure 4 illustrates the energies per unit area Eθ, Ec0 and Etot (left, middle and right panels, respectively) as a function of t. In agreement with the previous analysis for the displacement and pressure profiles, we see that Eθ tends to a constant value as time increases because the deformation of the mixture tends to become uniform in all the domain. At the same time, Ec0 decreases in time following the decrease of the pressure, being significant only during the initial transient. As a result, the total potential energy Etot of the mixture tends to coincide with Eθ, as demonstrated by the right panel of Figure 4.
Figure 5 illustrates the dissipation D and force term F (left and right panel, respectively) as a function of t. The temporal profiles of both functions rapidly decay to zero as time increases. This is consistent with the fact that the domain of the mixture tends to deform in a uniform manner as time gets large so that its time variation becomes rapidly negligible.
We assume that the external traction applied at the right boundary x=L has a jump discontinuity at t∗=0.25T.
Figure 6 illustrates the computed displacement, fluid pressure and discharge velocity (left, middle and right panel, respectively) as a function of t. Displacement and velocity are evaluated at x=L (right boundary) whereas the pressure is evaluated at x=0 (left boundary). We notice that the three graphs are the translation of the corresponding graphs in Figure 2. In particular, we see that u, p and v are continuous at t=t∗, where their value is equal to zero.
Figure 7 illustrates the displacement, fluid pressure and discharge velocity (left, middle and right panel, respectively) as a function of x and t.
Figure 8 illustrates the energies per unit area Eθ, Ec0 and Etot (left, middle and right panels, respectively) as a function of t.
Figure 9 illustrates the dissipation D and force term F (left and right panel, respectively) as a function of t. We notice that both D and F are discontinuous at t=t∗ where they experience a finite jump.
Let P(t) be the dimensionless unbounded ramp pulse of unit-slope starting at t=0, represented in Figure 10
P(t)=tH(t)=t, t≥0. | (7.7) |
In this case, we have
U(t)=t−η(1−e−tη) and U′(t)=1−e−tη,t≥0, | (7.8) |
so that:
Gn(t)∗1=1(Λ2n−Λ1n)[1−e−Λ1ntΛ1n−1−e−Λ2ntΛ2n], t≥0, | (7.9a) |
Gn(t)∗e−tη=1(Λ2n−Λ1n)[e−t/η−e−Λ1ntΛ1n−1/η−e−t/η−e−Λ2ntΛ2n−1/η], t≥0. | (7.9b) |
Summing the two above expressions we obtain
Gn(t)∗U′(t)=1(Λ2n−Λ1n)[1−e−Λ1ntΛ1n−1−e−Λ2ntΛ2n−e−t/η−e−Λ1ntΛ1n−1/η+e−t/η−e−Λ2ntΛ2n−1/η], t≥0. | (7.10) |
Note that the above expression can be rewritten as
Gn(t)∗U′(t)=1Λ1nΛ2n+1(Λ2n−Λ1n)(e−Λ1ntΛ1n(ηΛ1n−1)−e−Λ2ntΛ2n(ηΛ2n−1))+η2γ1−γe−tη. |
The time derivative of expression (7.10) is given by: for t≥0,
(Gn(t)∗U′(t))′=1(Λ2n−Λ1n)(e−Λ1nt−e−Λ2nt+(1/η)e−t/η−Λ1ne−Λ1ntΛ1n−1/η−(1/η)e−t/η−Λ2ne−Λ2ntΛ2n−1/η)=1(Λ2n−Λ1n)(−e−Λ1ntηΛ1n−1+e−Λ2ntηΛ2n−1)−ηγ1−γe−tη. | (7.11) |
Note that
η(Gn(t)∗U′(t))′=Gn(t)∗exp(−t/η)=Gn(t)∗1−Gn(t)∗U′(t), |
and the coefficients Bn(t) become
Bn(t)=1−γγηGn(t)∗1=1−γγη∫t0Gn(s)ds. | (7.12) |
Then the displacement, the pressure and the discharge velocity have the following representations for t≥0:
u0(x,t)=−x(t−η)+1−γγη∞∑n=02(−1)nλn[1Λ1nΛ2n+1Λ2n−Λ1n(e−Λ1ntΛ1n(ηΛ1n−1)−e−Λ2ntΛ2n(ηΛ2n−1))]Xn(x), | (7.13a) |
p0(x,t)=1−γγη∞∑n=02(−1)nλn1(Λ2n−Λ1n)[1−e−Λ1ntΛ1n−1−e−Λ2ntΛ2n]X′n(x), | (7.13b) |
v0(x,t)=1−γγη∞∑n=02(−1)n1(Λ2n−Λ1n)[1−e−Λ1ntΛ1n−1−e−Λ2ntΛ2n]Xn(x). | (7.13c) |
Figure 11 illustrates the displacement, fluid pressure and discharge velocity (left, middle and right panel, respectively) as a function of t. Displacement and velocity are evaluated at x=L (right boundary) whereas the pressure is evaluated at x=0 (left boundary). Unlike the case presented in Section 7.1, here the external pressure load continues to increase linearly with time, thereby inducing a continuous increase in displacement, pressure and velocity as time goes by. Over the observational time interval, the magnitude of the external load is smaller than what considered in Section 7.1 and this leads to a smaller (absolute) value of the maximum displacement, which is about 10μm in this case.
Figure 12 illustrates the displacement, fluid pressure and discharge velocity (left, middle and right panel, respectively) as a function of x and t. The displacement profile exhibits an approximately bilinear variation with respect to temporal and spatial coordinates, whereas pressure and velocity display a nonlinear behavior in the space-time domain. All dependent variables tend to increase in magnitude as a function of time at any spatial position of the mixture, the discharge velocity being closer to reach a stationary condition than the pressure.
Figure 13 illustrates the energies per unit area Eθ, Ec0 and Etot (left, middle and right panels, respectively) as a function of t. Interestingly, Eθ exhibits a nonlinear increase with respect to time in accordance with the fact that deformation is not constant in space. Ec0 follows a similar pattern because of the nonlinear trend of the pressure, but has a much smaller value, so that the Etot almost coincides with Eθ for all times.
Figure 14 illustrates the dissipation D and forcing term F (left and right panel, respectively) as a function of t. Dissipation increase with time is characterized by two markedly different slopes. In a first time interval, approximately equal to T/5, dissipation increases rapidly and is mainly determined by the fluid component of the mixture. In the remaining part of the observational time interval, dissipation increases less rapidly and is mainly determined by the structural viscoelasticity of the mixture. The forcing term increases linearly with time in accordance with the trend of the temporal variation of solid displacement at the right boundary of the domain, as shown in Figure 12.
Let P(ε)(t) be the dimensionless ramp pulse of unit amplitude and finite rise time (=ε>0) starting at t=0
P(ε)(t)=1ε[tH(t)−(t−ε)H(t−ε)]={0ift≤0t/εif0<t≤ε1ift>ε, | (7.14) |
represented graphically in Figure 15.
Using the linear superposition principle, the solid displacement, fluid pressure, and discharge velocity are given by:
uε0(x,t)={1εu0(x,t)if0≤t≤ε1ε[u0(x,t)−u0(x,t−ε)]ift>ε | (7.15a) |
pε0(x,t)={1εp0(x,t)if0≤t≤ε1ε[p0(x,t)−p0(x,t−ε)]ift>ε. | (7.15b) |
vε0(x,t)={1εv0(x,t)if0≤t≤ε1ε[v0(x,t)−v0(x,t−ε)]ift>ε. | (7.15c) |
where the expressions of u0(x,t), p0(x,t) and v0(x,t) are given by (7.13).
Notice that pε0(x,t),vε0(x,t)∈C0(Q) since p0(x,t),v0(x,t)∈C0(Q) and p0(x,0)=v0(x,0)=0.
Figure 16 illustrates the displacement, fluid pressure and discharge velocity (left, middle and right panel, respectively) as a function of t. Displacement and velocity are evaluated at x=L (right boundary) whereas the pressure is evaluated at x=0 (left boundary). We see that the displacement increases in magnitude almost linearly during the increase in time of the externally applied pressure. Then, it rapidly tends to stationary conditions once the external applied pressure becomes constant. The asymptotic value of the displacement is the same as in the case of the step pulse illustrated in Section 7.1. The profile of fluid pressure increases in time and the externally applied pressure increases; once the solid deformation attains stationary conditions, the fluid pressure decreases. A similar trend is shown by the discharge velocity. The maximum value of pressure is the same as in the case of a step pulse external pressure whereas the maximum value of the discharge velocity is slightly smaller and coincides with the value of the velocity at t=T/2 that is obtained in the case of an unbounded ramp external pressure (cf. Figure 11, right panel).
Figure 17 illustrates the displacement, fluid pressure and discharge velocity (left, middle and right panel, respectively) as a function of x and t. The spatial variation of the displacement becomes linear after the external pressure ceases to increase, whereas, in the same time interval, both fluid pressure and discharge velocity exhibit a spatial decrease.
Figure 18 illustrates the energies per unit area Eθ, Ec0 and Etot (left, middle and right panels, respectively) as a function of t. Both Eθ and Ec0 increase in time until the external pressure increases. Then, Eθ becomes constant since deformation is constant in time, whereas Ec0 decreases because of the decrease of fluid pressure. As in all previously considered examples, the contribution of Ec0 is much smaller than that of Eθ so that the Etot almost coincides with Eθ.
Figure 19 illustrates the dissipation D and forcing term F (left and right panel, respectively) as a function of t. Both terms exhibit an increase with respect to time during the increase of the externally applied pressure. Then, they both experience a sudden decay once the externally applied pressure becomes constant. Dissipation tends to an asymptotic value that is much smaller than the value attained at t=T/2 whereas the force term tends to zero since the boundary displacement is almost constant in time after t=T/2.
In this section we investigate the dependence of the solution of model (4.12) on the compressibility parameter c0. In terms of the dimensionless equation system (5.3), this amounts to analyzing the solutions as a function of the quantity γ defined in (5.2). We denote by c0,ref=1.67⋅10−5m2N−1 the reference value of the compressibility parameter. This value has been used in all the computations illustrated in Section 7. Then, we let c0 assume the following values:
c0=[0.001,0.01,0.1,1,10100,1000]c0,ref |
and we compute the solid displacement and discharge velocity at x=L, the fluid pressure at x=0 and the energies Eθ, Ec0 and Etot, as functions of time in the interval [0,T], T=10000s being the width of the observational window considered in Section 7. All the other model parameters have been set equal to the values adopted in Section 7. The applied boundary traction is a step pulse of amplitude Pref at t∗=0.
Figure 20 shows the profiles of solid displacement (left panel), fluid pressure (middle panel) and discharge velocity (right panel) as functions of time and of the compressibility parameter c0. We notice that compressibility has a significant impact on the quantitative values attained by the solution variables. In particular we see that increasing c0 gives rise to:
● an increase of the magnitude of the solid displacement;
● a decrease of pressure and velocity;
● a right shift of the peak of pressure and velocity.
These behaviors are indicative of the fact that the compressibility of the mixture components allows the body to deform more under the same pressure load, thereby reducing the internal level of fluid pressure and limiting the impact on the fluid velocity. We also notice that decreasing c0 has a much more significant impact than increasing c0 on solution range variation. From the theoretical viewpoint, this is due to the fact that γ→1− for large values of c0 and, as a consequence, the ratio (1−γ)/(γη), that characterizes the mathematical form of the solution expressions (7.5), tends to 0. From the physical viewpoint, decreasing the value of c0 implies that the body cannot significantly deform upon applying a pressure load, thereby inducing higher levels of fluid pressure whose gradients lead to larger fluid velocities.
Figure 21 illustrates the discharge velocity computed by the model studied in this work and the discharge velocity computed by the model studied in [4]. The principal difference between the two models is that no compressibility was included in [4]. We see that:
● the velocity predicted by the model of [4] is a (very sharp) upper bound for all the velocities predicted by the model studied in the present work, when c0→0+;
● for increasing values of c0, the velocities predicted by the present model differ substantially from the upper bound velocity yielded by the model in [4]. In particular, we see that increasing c0 gives rise to a decrease and to a right shift of the peak in the velocity profile.
Figure 22 shows the profiles of Eθ, Ec0 and Etot as a function of t and of c0. We see that increasing c0 has the effect of increasing substantially and monotonically the magnitude of Eθ because the deformation profile gets larger as the components become more compressible, in accordance with Figure 20, left panel. On the contrary, Ec0 exhibits a nonmonotonic dependence on c0. Specifically, for c0<c0,ref we see that increasing c0 gives rise to an increase of Ec0 with a right shift of its peak; then, for c0≥c0,ref, we see that increasing c0 leads to a significant monotonic decrease of Ec0. For the theoretical viewpoint, this is due to the fact that the ratio γ/(1−γ)→+∞ as c0→+∞ and the pressure profile tends to decrease in accordance with Figure 20, middle panel. The physical meaning of these results can be appreciated by observing that the predicted total energy is mainly determined by Eθ for large values of c0, since larger deformations can occur for more compressible components, whereas the contribution to the total energy given by Ec0 becomes more significant for smaller values of c0, since larger pressures develop for less compressible components.
The analysis presented in this work shows that the solutions of a poroviscoelastic model with compressible components remain bounded even in the case when the imposed boundary traction is irregular in time. In particular, given a certain functional form for the boundary traction and a certain level of structural viscoelasticity, the discharge velocity attains a maximum value that is lower when compressibility is higher. By investigating the dynamics of the energy functionals characterizing the system, we showed that this limiting effect is due to the capability of the system to store potential energy as its components are elastically compressed, thereby delaying the transmission of irregularities in the linear momentum from the solid to the fluid. As a result, the fluid has the time to accommodate for sudden changes, resulting in bounded velocities. This mechanism is very different from that provided by structural viscoelasticity, whose limiting effect on the discharge velocity is due to increased viscous dissipation, as shown in [4].
Ultimately, this work elucidates the specific role that compressibility plays in the control of fluid flow through complex deformable porous structure, which finds numerous applications in science and engineering. The work presented here offers many future directions of research. For example, it would be very interesting to investigate how the findings concerning the role of compressibility would translate to a more realistic three-dimensional setting. Furthermore, different viscoelastic models could be considered and the specific roles of fluid and solid viscosities could be investigated and compared.
Dr. Bociu has been partially supported by NSF CAREER DMS-1555062. Dr. Guidoboni has been partially supported by the award NSF DMS-1853222. Dr. Sacco has been partially supported by Micron Semiconductor Italia S.r.l., SOW nr. 4505462139.
The authors declare there is no conflict of interest.
[1] | L. Bociu, G. Guidoboni, R. Sacco, et al., Analysis of nonlinear poro-elastic and poro-viscoelastic models. Arch. Rational Mech. Anal., 222 (2016), 1445–1519. |
[2] | C.-Y. Huang, V. C. Mow and G. A. Ateshian, The role of flow-independent viscoelasticity in the biphasic tensile and compressive responses of articular cartilage. J. Biomech. Eng., 123 (2001), 410–417. |
[3] | C.-Y. Huang, M. A. Soltz, M. Kopacz, et al., Experimental verification of the roles of intrinsic matrix viscoelasticity and tension-compression nonlinearity in the biphasic response of cartilage. J. Biomech. Eng., 125 (2003), 84–93. |
[4] | M. Verri, G. Guidoboni, L. Bociu, et al., The role of structural viscoelasticity in deformable porous media with incompressible constituents: Applications in biomechanics. Math. Biosci. Eng., 15 (2018), 933. |
[5] | D. Prada, A. Harris, G. Guidoboni, et al., Autoregulation and neurovascular coupling in the optic nerve head. Surv. Ophthalmol., 61 (2016), 164–186. |
[6] | P. Causin, G. Guidoboni, A. Harris, et al., A poroelastic model for the perfusion of the lamina cribrosa in the optic nerve head. Math. Biosci., 257 (2014), 33–41. |
[7] | J. C. Gross, A. Harris, B. A. Siesky, et al., Mathematical modeling for novel treatment approaches to open-angle glaucoma. Expert Rev. Ophthalmol., 12 (2017), 443–455. |
[8] | A. Harris, G. Guidoboni, J. C. Arciero, et al., Ocular hemodynamics and glaucoma: the role of mathematical modeling. Eur. J. Ophthalmol., 23 (2013), 139–146. |
[9] | J. H. Kim and J. Caprioli,. Intraocular pressure fluctuation: Is it important? J. Ophthalmic Vis. Res., 13 (2018), 170. |
[10] | G. Guidoboni, F. Salerni, A. Harris, et al., Ocular and cerebral hemo-fluid dynamics in microgravity: a mathematical model. Invest. Ophth. Vis. Sci., 58 (2017), 3036–3036. |
[11] | A. G. Lee, T. H. Mader, C. R. Gibson, et al., Space flight-associated neuro-ocular syndrome (sans). Eye, 32 (2018), 1164–1167. |
[12] | M. Schanz and A.-D. Cheng. Dynamic analysis of a one-dimensional poroviscoelastic column. J. Appl. Mech., 68 (2001), 192–198. |
[13] | R. E. Showalter. Diffusion in poro-elastic media. J. Math. Anal. Appl., 251 (2000), 310–340. |
[14] | M. Biot. General theory of three-dimensional consolidation. J. Appl. Phys., 12 (1941), 155–164. |
[15] | E. Detournay and A.-D. Cheng. Fundamentals of poroelasticity. In J. A. Hudson, editor, Comprehensive Rock Eng., 2 (1993), 113–171. Pergamon. |
[16] | M. A. Soltz and G. A. Ateshian. Experimental verification and theoretical prediction of cartilage interstitial fluid pressurization at an impermeable contact interface in confined compression. J. Biomech., 31 (1998), 927–934. |
[17] | F. Treves. Basic Linear Partial Differential Equations. Academic Press, 1975. |
1. | Giovanna Guidoboni, Riccardo Sacco, Marcela Szopos, Lorenzo Sala, Alice Chandra Verticchio Vercellin, Brent Siesky, Alon Harris, Neurodegenerative Disorders of the Eye and of the Brain: A Perspective on Their Fluid-Dynamical Connections and the Potential of Mechanism-Driven Modeling, 2020, 14, 1662-453X, 10.3389/fnins.2020.566428 | |
2. | Alon Harris, Giovanna Guidoboni, Brent Siesky, Sunu Mathew, Alice C. Verticchio Vercellin, Lucas Rowe, Julia Arciero, Ocular blood flow as a clinical observation: Value, limitations and data analysis, 2020, 78, 13509462, 100841, 10.1016/j.preteyeres.2020.100841 | |
3. | Maia M. Svanadze, Potential Method in the Coupled Theory of Viscoelasticity of Porous Materials, 2021, 144, 0374-3535, 119, 10.1007/s10659-021-09830-y | |
4. | Lorena Bociu, Sarah Strikwerda, Optimal control in poroelasticity, 2022, 101, 0003-6811, 1774, 10.1080/00036811.2021.2008372 | |
5. | Lorena Bociu, Giovanna Guidoboni, Riccardo Sacco, Daniele Prada, Numerical simulation and analysis of multiscale interface coupling between a poroelastic medium and a lumped hydraulic circuit: Comparison between functional iteration and operator splitting methods, 2022, 466, 00219991, 111379, 10.1016/j.jcp.2022.111379 | |
6. | Lorena Bociu, Sunčica Canic, Boris Muha, Justin T. Webster, Multilayered Poroelasticity Interacting with Stokes Flow, 2021, 53, 0036-1410, 6243, 10.1137/20M1382520 | |
7. | Lorena Bociu, Sarah Strikwerda, 2022, Chapter 5, 978-3-031-04495-3, 103, 10.1007/978-3-031-04496-0_5 | |
8. | Lorena Bociu, Boris Muha, Justin T. Webster, Weak solutions in nonlinear poroelasticity with incompressible constituents, 2022, 67, 14681218, 103563, 10.1016/j.nonrwa.2022.103563 | |
9. | Lorena Bociu, Justin T. Webster, Nonlinear quasi-static poroelasticity, 2021, 296, 00220396, 242, 10.1016/j.jde.2021.05.060 | |
10. | Xiong Tang, Pei Zheng, Liangwei Zhong, Keming Zhang, A Poroviscoelastic Model at Finite Strains for a Viscous Compressible Solid Skeleton, 2024, 16, 1758-8251, 10.1142/S1758825124500583 | |
11. | Lorena Bociu, Boris Muha, Justin T. Webster, Mathematical effects of linear visco-elasticity in quasi-static Biot models, 2023, 527, 0022247X, 127462, 10.1016/j.jmaa.2023.127462 | |
12. | Lorenzo Sala, Christophe Prud'homme, Giovanna Guidoboni, Marcela Szopos, Alon Harris, The ocular mathematical virtual simulator: A validated multiscale model for hemodynamics and biomechanics in the human eye, 2024, 40, 2040-7939, 10.1002/cnm.3791 | |
13. | Maia M. Svanadze, On the coupled linear theory of thermoviscoelasticity of porous materials, 2024, 235, 0001-5970, 4253, 10.1007/s00707-024-03937-8 | |
14. | Lorena Bociu, Matthew Broussard, Giovanna Guidoboni, Sarah Strikwerda, Analysis of a Multiscale Interface Problem Based on the Coupling of Partial and Ordinary Differential Equations to Model Tissue Perfusion, 2025, 23, 1540-3459, 1, 10.1137/23M1628073 | |
15. | L. Bociu, M. Broussard, G. Guidoboni, D. Prada, S. Strikwerda, Comparing Interface Conditions for a 3D–0D Multiscale Interface Coupling With Applications in Tissue Perfusion, 2025, 41, 2040-7939, 10.1002/cnm.70017 |
symbol | value | units |
L | 0.81⋅10−3 | m |
T | 105,104 | s |
θ | 0.97⋅106 | Nm−2 |
K0 | 2.9⋅10−16 | m4N−1s−1 |
Pref | 6⋅104 | Nm−2 |
c0 | 1.67⋅10−5 | m2N−1 |
η | 4.85⋅109,4.85⋅108 | Nsm−2 |
symbol | value |
L | 1 |
T | 2.5, 0.25 |
η | 1.26⋅10−1, 1.26⋅10−2 |
γ | 0.95 |