# Transient Propagation of Longitudinal and Transverse Waves in Cancellous Bone: Application of Biot Theory and Fractional Calculus

Djihane Benmorsli, Zine El Abidine Fellah, Djema Belgroune, Nicholas Ongwen, Erick Ogam, Claude L. Depollier, Mohamed Fellah

<figure><img src="image_1.png" /></figure>

## ► To cite this version:

Djihane Benmorsli, Zine El Abidine Fellah, Djema Belgroune, Nicholas Ongwen, Erick Ogam, et al.. Transient Propagation of Longitudinal and Transverse Waves in Cancellous Bone: Application of Biot Theory and Fractional Calculus. Symmetry, 2022, 14, pp.1971. 10.3390/sym14101971. hal-03808042

HAL Id: hal-03808042

https://hal.science/hal-03808042

Submitted on 10 Oct 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>
Article

# Transient Propagation of Longitudinal and Transverse Waves in Cancellous Bone: Application of Biot Theory and Fractional Calculus

Djihane Benmorsli<fnref n="1" />, Zine El Abidine Fellah <fnref n="2" />* ID, Djema Belgroune<fnref n="3" />, Nicholas O. Ongwen <fnref n="4" />, Erick Ogam <fnref n="2" /> ID,
Claude Depollier <fnref n="5" /> ID and Mohamed Fellah <fnref n="1" />

<figure><img src="image_3.png" /></figure>

check for
updates

**Citation:** Benmorsli, D; Fellah, Z.E.A.;
Belgroune, D.; Ongwen, N.O.;
Ogam, E.; Depollier, C; Felllah, M.
Transient Propagation of
Longitudinal and Transverse Waves
in Cancellous Bone: Application of
Biot Theory and Fractional Calculus.
*Symmetry* **2022**, *14*, 1971. https://
doi.org/10.3390/sym14101971

Academic Editors: Serkan Araci and
Sergei D. Odintsov

Received: 1 August 2022
Accepted: 14 September 2022
Published: 21 September 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

<figure><img src="image_4.png" /></figure>

**Copyright:** © 2022 by the authors.
Licensee MDPI, Basel, Switzerland.
This article is an open access article
distributed under the terms and
conditions of the Creative Commons
Attribution (CC BY) license (https://
creativecommons.org/licenses/by/
4.0/).

**Abstract:** In this paper, the influence of the transverse wave on sound propagation in a porous medium with a flexible structure is considered. The study is carried out in the time domain using the modified Biot theory obtained by the symmetry of the Lagrangian (invariance by translation and rotation). The viscous exchanges between the fluid and the structure are described by fractional calculus. When a sound pulse arrives at normal incidence on a porous material with a flexible structure, the transverse waves interfere with the longitudinal waves during propagation because of the viscous interactions that appear between the fluid and the structure. By performing a calculation in the Laplace domain, the reflection and transmission operators are derived. Their time domain expressions depend on the Green functions of the longitudinal and transverse waves. In order to study the effects of the transverse wave on the transmitted longitudinal waves, numerical simulations of the transmitted waves in the time domain by varying the characteristic parameters of the medium are realized whether the transverse wave is considered or not.

**Keywords:** porous materials; fractional calculus; Biot theory; transverse wave; transmission

## 1. Introduction

Osteoporosis is defined as a sickness that affects the bones, and to understand it, an overview of the structure of bone tissues is necessary. Bones are in perpetual renewal; this phenomenon is called bone remodeling: constant construction of new tissue and destruction of old tissue. In osteoporosis, biochemical and hormonal changes occur and affect this balance. The bone is, therefore, considerably weakened by reductions in its tissue and mass [1], a fragility that is increased tenfold when the contribution of bone minerals is reduced [2–4]. Therefore, there is a growing need for improved diagnosis and management of osteoporosis. Indeed, the sooner the disease is detected, the more effective the care can be.

As it is known, bone is an inhomogeneous porous material [5–9] and the determination of its characteristics by means of ultrasound techniques can be a game-changer in the diagnosis of osteoporosis since its elastic properties affect acoustic propagation. Within this context, several studies have been conducted on cancellous and trabecular bone [10–13]. There are ultrasonic methods that rely on the speed of sound and the broadband attenuation of a single wave to estimate the different stages of osteoporosis, but the poroelastic nature of the bone generates two types of waves during the propagation of the acoustic wave

<footnote n="1">¹ Laboratory of Theoretical Physics, Faculty of Physics, University of Science and Technology Houari-Boumediene, BP 32 El Alia, Bab Ezzouar 16111, Algeria</footnote>

<footnote n="2">² Aix Marseille Univ, CNRS, Centrale Marseille, LMA UMR 7031, Laboratory of Mechanics and Acoustics, CEDEX 20, 13402 Marseille, France</footnote>

<footnote n="3">³ Materials Physics Laboratory, Faculty of Physics, University of Science and Technology Houari-Boumediene, BP 32 El Alia, Bab Ezzouar 16111, Algeria</footnote>

<footnote n="4">⁴ Department of Physics and Materials Science, Maseno University, Maseno 40105, Kenya</footnote>

<footnote n="5">⁵ Laboratoire d'Acoustique de l'Université du Mans (LAUM), UMR 6613, Institut d'Acoustique, CNRS, Le Mans Université, Avenue O. Messiaen, CEDEX 09, 72085 Le Mans, France</footnote>

<footnote>* Correspondence: fellah@lma.cnrs-mrs.fr</footnote>
in the bone. The most appropriate theory for the study of acoustic wave propagation in porous media such as trabecular bone is the Biot theory [14–27]. This theory [14–27] has been proven to be a very advantageous model for the description of wave propagation in porous media saturated with fluids; it takes into account the movements of the solid and fluid phases, and through the mass, elastic, and viscous couplings, the interaction between these two phases is considered. Its main success lies in the prediction of the existence of three waves in a porous material subjected to external stress: two longitudinal waves and a transverse wave. The first compressional wave has the highest speed and is called the “fast wave”, while the second compressional wave has a lower speed and is called the “slow wave”.

Because of its very general and rather fundamental character, the Biot theory has been used widely in the field of petroleum exploration and geoscientific testing, but many have also applied it to the study of bone structures [14, 15, 28–31]. McKelvie [28,31] qualitatively predicted that, for cancelous bone, attenuation is related to ultrasound frequency. The second compressional wave with a lower velocity was reported in the studies by Lakes, Yoon, and Katz [32] on wet human and bovine cortical bone. The connection between porosity and longitudinal wave velocities and the relevance of the presence of these two waves to bone characterization were demonstrated by Hosokawa and Otani [14] in their studies of cancelous bone. By implementing the concepts of dynamic tortuosity, permeability, and viscous characteristic length, in 1987, Johnson et al. [33] made a major input to the study of porous media.

The modified Biot theory has been used to study the transmission and reflection of acoustic waves for different media in the frequency domain [34–42]. In order to offer an alternative to the frequency methods, Fellah and Depollier [43,44] have proposed a temporal approach [43–52] that allows the phenomena that follow the propagation of acoustic waves in porous media to be seen from a different angle. The viscous fluid/structure interactions are described using fractional calculus, which has been proven to be a well-adapted tool to describe the dispersion effects that develop when a wave passes in a saturated porous medium. Fractional derivatives and integrals [53–65] have long been the focus of research in mathematics, but their use outside mathematics is only recent. Indeed, the fractional derivative was only used towards the end of the 1960s by Caputo [57–59] to describe dissipation phenomena in solids and, later, by Bagley and Torvik [60] to treat viscoelastic materials, but the most significant contribution in the field of acoustics remains the one made by Matignon [61]. The reflection and transmission operators have been established as a result to external solicitation of a saturated porous medium by a viscous fluid in the frequency domain for normal incidence by considering the three Biot waves instead of only the longitudinal waves, by Hodaei, Maghoul, and Popplewell [66].

The purpose of this article is to study, from a new perspective, the behavior of a porous medium subjected to an external load when the three Biot waves (i.e., the transverse wave in addition to the longitudinal waves) are considered. For that, we establish the expressions of the diffusion operators in the time domain while including the effect of the transverse wave; the fractional calculus as well as the Biot model modified by Johnson et al. [33] will be used. Numerical simulations are given for the transmitted signals, with shear being taken into account or not while applying a variation of the characteristic parameters of the porous medium. The effect of transverse waves as well as the sensitivity of the parameters are studied. The advantage of this model is that it presents a temporal modeling of the propagation by taking into account both longitudinal and transverse waves in normal incidence.
## 2. Modified Biot Theory

The displacements of the solid phase **u** and the fluid phase **U** are modeled via the Biot theory, where the equations of motion are given by the following:

$$ 
\begin{aligned}
\rho_{11} \frac{\partial^2 \mathbf{u}}{\partial t^2} + \rho_{12} \frac{\partial^2 \mathbf{U}}{\partial t^2} &= P \nabla (\nabla \cdot \mathbf{u}) + Q \nabla (\nabla \cdot \mathbf{U}) - N \nabla \wedge (\nabla \wedge \mathbf{u}), \\
\rho_{22} \frac{\partial^2 \mathbf{U}}{\partial t^2} + \rho_{12} \frac{\partial^2 \mathbf{u}}{\partial t^2} &= R \nabla (\nabla \cdot \mathbf{U}) + Q \nabla (\nabla \cdot \mathbf{u}).
\end{aligned}
\quad (1) $$

$N$ is the shear modulus of the solid structure. $P$, $Q$, and $R$ are the elastic coefficients of Biot; they are functions of the bulk modulus of the pore fluid $K_f$, the bulk modulus of the elastic solid $K_s$, the bulk modulus of the porous skeletal frame $K_b$ as well as of the porosity $\phi$ of the porous medium and are defined by the relations below:

$$ P = \frac{(1 - \phi) \left(1 - \phi - \frac{K_b}{K_s}\right) K_s + \phi \frac{K_s}{K_f}}{1 - \phi - \frac{K_b}{K_s} + \phi \frac{K_s}{K_f}} + \frac{4}{3} N, \quad (2) $$

$$ Q = \frac{\left(1 - \phi - \frac{K_b}{K_s}\right) \phi K_s}{1 - \phi - \frac{K_b}{K_s} + \phi \frac{K_s}{K_f}}, $$

$$ R = \frac{\phi^2 K_s}{1 - \phi - \frac{K_b}{K_s} + \phi \frac{K_s}{K_f}}, $$

$K_f$, $K_s$, and $K_b$ depend on the classical mechanical parameters, which are Young's moduli, Poisson's coefficients of the solid $E_s$ and $v_s$, and Poisson's coefficients of the skeleton $E_b$ and $v_b$ as follows:

$$ K_s = \frac{E_s}{3(1 - 2v_s)}, \quad K_b = \frac{E_b}{3(1 - 2v_b)}, \quad N = \frac{E_b}{2(1 + v_b)}. \quad (3) $$

The coefficients $\rho_{ij}$ involved in the motion equations are mass coefficients. They have a density dimension and are linked to the skeleton-constituting solid density $\rho_s$, which constitutes the skeleton, and to the fluid density $\rho_f$ by:

$$ \rho_{11} = (1 - \phi)\rho_s - \rho_{12}, \quad (4) $$

$$ \rho_{22} = \phi \rho_f - \rho_{12}, $$

The term $\rho_{12}$ expresses the fluid-structure inertial coupling. It appears as a contribution to the densities of the fluid and the solid because of the interactions between the fluid and the structure.

$$ \rho_{12} = -\phi \rho_f (\alpha_{\infty} - 1), \quad (5) $$

$\alpha_{\infty}$ is the tortuosity of the medium. When an acoustic wave crosses a fluid-saturated porous medium, it is attenuated. This attenuation is mainly due to the viscous exchanges that occur between the structure and the saturating fluid because of relative fluid/solid displacements, leading to a modification of the densities. The tortuosity becomes a function of the frequency and is called dynamic tortuosity $\alpha(\omega)$:

$$ 
\begin{aligned}
\rho_{12}(\omega) &= -\phi \rho_f (\alpha(\omega) - 1), \\
\rho_{11}(\omega) &= (1 - \phi) \rho_s + \phi \rho_f (\alpha(\omega) - 1), \\
\rho_{22}(\omega) &= \phi \rho_f + \phi \rho_f (\alpha_{\infty} - 1).
\end{aligned}
\quad (6) $$

In order to estimate the fluid domain affected by these exchanges, the thickness of the viscous boundary layer $\delta = \left(\frac{2\eta}{\omega\rho_0}\right)^{1/2}$ ($\eta$ fluid viscosity and $\omega$ angular frequency) is
compared to a microscopic dimension characteristic of the medium such as the pore size “a”. In the high-frequency mode, the viscous boundary layer thickness is small compared to the pore size “a” i.e., if ($\frac{\delta}{a} \ll 1$). In this mode, the viscous effects take place only in a small fluid thickness, at the walls where the velocities are affected by the fluid/solid displacements. In the rest of the fluid volume, the fluid behaves like a perfect fluid.

In this case, the dynamic tortuosity takes the following form:

$$ \alpha(\omega) \approx \alpha_{\infty} \left( 1 + \frac{2}{\Lambda} \sqrt{\frac{\eta}{j\omega\rho_f}} \right), \qquad (7) $$

This model was proposed by Johnson et al. [33]. It involves the tortuosity $\alpha_{\infty}$ and the viscous characteristic length $\Lambda$, a geometric parameter indicative of small pores, in the high-frequency limit and represents the thickness where viscous effects are important.

By replacing the new expressions for the mass coefficients in (6) where $\alpha(\omega)$ is given by (7) in the system of Equation (1), the motion equations becomes the following:

$$ \begin{aligned} \rho_{11}(\omega) \frac{\partial^2 \mathbf{u}}{\partial t^2} + \rho_{12}(\omega) \frac{\partial^2 \mathbf{U}(t)}{\partial t^2} &= P \nabla (\nabla \cdot \mathbf{u}) + Q \nabla (\nabla \cdot \mathbf{U}) - N \nabla \wedge (\nabla \wedge \mathbf{u}), \\ \rho_{22}(\omega) \frac{\partial^2 \mathbf{U}}{\partial t^2} + \rho_{12}(\omega) \frac{\partial^2 \mathbf{u}}{\partial t^2} &= R \nabla (\nabla \cdot \mathbf{U}) + Q \nabla (\nabla \cdot \mathbf{u}). \end{aligned} \qquad (8) $$

This system is written neither in the time domain nor in the frequency domain, since the densities depend on the frequency and the derivatives are in the time domain. In the following section, we try to describe the above problem in the time domain.

### 3. Temporal Formulation of the Modified Biot Theory

In order to study the waves transmitted in the time domain through a porous medium, we focus on the temporal formulation of the Biot theory. To achieve this, we introduce the concept of the fractional derivative of order $\nu$, the definition of which [56] is given by the following:

$$ D^{\nu}[x(t)] = \frac{1}{\Gamma(-\nu)} \int_{0}^{t} (t-u)^{-\nu-1}x(u)du, \qquad (9) $$

where $\nu$ is a real number and $\Gamma$ is the Gamma function.

The high-frequency expression of the dynamic tortuosity includes the term $\sqrt{\frac{1}{j\omega}}$, of which the time equivalent is a fractional derivative of order 1/2. Following the definition in (9), the expression of the response factor $\alpha$ takes the following expression in the time domain [43,49]:

$$ \tilde{\alpha}(t) = \alpha_{\infty} \left( \delta(t) + \frac{2}{\Lambda} \sqrt{\frac{\eta}{\pi \rho_f}} t^{-\frac{1}{2}} \right), \qquad (10) $$

where $\delta(t)$ is the Dirac function. With this new expression of the dynamic tortuosity $\tilde{\alpha}(t)$ in the system of Equation (1), the equations of motion in the time domain take the following shape:

$$ \begin{aligned} & \int_{0}^{t} \tilde{\rho}_{11}(t-t') \frac{\partial^2 \mathbf{u}(t')}{\partial t'^2} dt' + \int_{0}^{t} \tilde{\rho}_{12}(t-t') \frac{\partial^2 \mathbf{U}(t')}{\partial t'^2} dt' \\ & = P \nabla (\nabla \cdot \mathbf{u}(t)) + Q \nabla (\nabla \cdot \mathbf{U}(t)) - N \nabla \wedge (\nabla \wedge \mathbf{u}(t)), \end{aligned} \qquad (11) $$

$$ \int_{0}^{t} \tilde{\rho}_{22}(t-t') \frac{\partial^2 \mathbf{U}(t')}{\partial t'^2} dt' + \int_{0}^{t} \tilde{\rho}_{12}(t-t') \frac{\partial^2 \mathbf{u}(t')}{\partial t'^2} dt' = R \nabla (\nabla \cdot \mathbf{U}(t)) + Q \nabla (\nabla \cdot \mathbf{u}(t)), $$

where the temporal operators $\tilde{\rho}_{11}$, $\tilde{\rho}_{12}$, and $\tilde{\rho}_{22}$ are obtained by replacing $\alpha(\omega)$ in Equation (6) by its temporal equivalent $\tilde{\alpha}$ given by Equation (10).
### 3.1. Solution of the Motion Equations

According to Helmholtz's theorem, a solution to the system in (12) has the following form:

$$ \mathbf{u}(x,t) = \nabla\Phi_s(x,t) + \nabla \wedge \mathbf{\Psi}_s(x,t), \quad (12) $$

$$ \mathbf{U}(x,t) = \nabla\Phi_f(x,t) + \nabla \wedge \mathbf{\Psi}_f(x,t), $$

where $\Phi_s$ and $\Phi_f$ are the scalar displacement potentials of the solid $\mathbf{u}$ and fluid $\mathbf{U}$, respectively, and $\mathbf{\Psi}_s$ and $\mathbf{\Psi}_f$ are the vector potentials of $\mathbf{u}$ and $\mathbf{U}$, respectively. The Biot equations of motion are linear, and it is possible to look for the scalar and vector solutions separately.

#### 3.1.1. Longitudinal Waves

In the case of the longitudinal wave, the contribution of scalar displacement potentials are obtained by substituting $\mathbf{u}(x, t) = \nabla\Phi_s(x, t)$ and $\mathbf{U} = \nabla\Phi_f(x, t)$ in the equations of motion, which gives

$$ \begin{pmatrix} \rho_{11} & \rho_{12} \\ \rho_{12} & \rho_{22} \end{pmatrix} \frac{\partial^2}{\partial t^2} \begin{pmatrix} \Phi_s(x,t) \\ \Phi_f(x,t) \end{pmatrix} + a \begin{pmatrix} 1 & -1 \\ -1 & 1 \end{pmatrix} \frac{\partial^3}{\partial t^3} \begin{pmatrix} \Phi_s(x,t) \\ \Phi_f(x,t) \end{pmatrix} = \begin{pmatrix} P & Q \\ Q & R \end{pmatrix} \Delta \begin{pmatrix} \Phi_s(x,t) \\ \Phi_f(x,t) \end{pmatrix}, \quad (13) $$

where $a = 2 \frac{\phi \rho_f \alpha_\infty}{\Lambda} \sqrt[2]{\frac{\eta}{\rho_f}}$ $\Delta$ is the Laplacian. The viscous effects are included in the fractional derivative term $\frac{\partial^3}{\partial t^3}$. To solve this system of equations, we use the Laplace transform. The Biot equations are then written as follows:

$$ \Delta \begin{pmatrix} \tilde{\Phi}_s(x,s) \\ \tilde{\Phi}_f(x,s) \end{pmatrix} = M \begin{pmatrix} \tilde{\Phi}_s(x,s) \\ \tilde{\Phi}_f(x,s) \end{pmatrix}, \quad (14) $$

where $\tilde{\Phi}_s$ and $\tilde{\Phi}_f$ are, respectively, the Laplace transform of $\Phi_s$ and $\Phi_f$. We note that $\tilde{\Phi}_s(x, s)$ is the Laplace transform of $\Phi_s(x, t)$ defined by the following:

$$ \tilde{\Phi}_s(s) = L[\Phi_s(t)] = \int_0^{+\infty} \exp(-st) \Phi_s(t) dt. \quad (15) $$

The matrix $M$ is given by the following:

$$ M = \begin{pmatrix} (R' \rho_{11} - Q' \rho_{12})s^2 + a s^{\frac{3}{2}}(R' + Q') & (R' \rho_{12} - Q' \rho_{22})s^2 - a s^{\frac{3}{2}}(R' + Q') \\ (P' \rho_{12} - Q' \rho_{11})s^2 - a s^{\frac{3}{2}}(P' + Q') & (P' \rho_{22} - Q' \rho_{12})s^2 + a s^{\frac{3}{2}}(P' + Q') \end{pmatrix}, \quad (16) $$

with

$$ \begin{pmatrix} R' & -Q' \\ -Q' & P' \end{pmatrix} = \begin{pmatrix} P & Q \\ Q & R \end{pmatrix}^{-1} = \begin{pmatrix} \frac{R}{PR-Q^2} & -\frac{Q}{PR-Q^2} \\ -\frac{Q}{PR-Q^2} & \frac{P}{PR-Q^2} \end{pmatrix}, $$

$$ R' = \frac{R}{PR - Q^2}, \quad Q' = \frac{Q}{PR - Q^2}, \quad P' = \frac{P}{PR - Q^2}. $$

The system is solved by finding the eigenvalues of the matrix $M$:

$$ \Delta \begin{pmatrix} \Phi_1(x,s) \\ \Phi_2(x,s) \end{pmatrix} = \begin{pmatrix} \lambda_1(s) & 0 \\ 0 & \lambda_2(s) \end{pmatrix} \begin{pmatrix} \Phi_1(x,s) \\ \Phi_2(x,s) \end{pmatrix}. \quad (17) $$
This corresponds to the following equations:

$$ \frac{\partial^2 \Phi_1(x, s)}{\partial x^2} = \lambda_1(s) \Phi_1(x, s), \qquad (18) $$

$$ \frac{\partial^2 \Phi_2(x, s)}{\partial x^2} = \lambda_2(s) \Phi_2(x, s), $$

The solutions to the above are in the following form:

$$ \Phi_1(x, s) = \Phi_{11}(s)e^{-x\sqrt{\lambda_1(s)}} + \Phi_{12}(s)e^{x\sqrt{\lambda_1(s)}}, \qquad (19) $$

$$ \Phi_2(x, s) = \Phi_{21}(s)e^{-x\sqrt{\lambda_2(s)}} + \Phi_{22}(s)e^{x\sqrt{\lambda_2(s)}}, \qquad (20) $$

The displacements of the solid and fluid are the linear combinations of the displacements due to the fast wave (first-type wave) characterized by the scalar potential $\Phi_1(x, s)$ and the slow wave (second-type wave) characterized by the scalar potential $\Phi_2(x, s)$.

We associate the eigenvalues $\lambda_1(s)$ and $\lambda_2(s)$ with eigenvectors $(1, V_1(s))$ and $(1, V_2(s))$, which link the fast and slow waves $\Phi_1(x, s)$ and $\Phi_2(x, s)$ to the solid and fluid waves $\tilde{\Phi}_s(x, s)$ and $\tilde{\Phi}_f(x, s)$.

$$ \begin{pmatrix} \tilde{\Phi}_s(x, s) \\ \tilde{\Phi}_f(x, s) \end{pmatrix} = \begin{pmatrix} 1 & 1 \\ V_1(s) & V_2(s) \end{pmatrix} \begin{pmatrix} \Phi_1(x, s) \\ \Phi_2(x, s) \end{pmatrix}. \qquad (21) $$

The expressions for the eigenvalues $\lambda_1(s)$ and $\lambda_2(s)$ and the components $V_1(s)$ and $V_2(s)$ of the eigenvectors are given in Appendix A.

### 3.1.2. Transverse Wave

For the wave, the displacements of the solid and fluid phases can be written as follows:

$$ \mathbf{u}(t) = \nabla \wedge \mathbf{\Psi}_s(t), \qquad (22) $$

$$ \mathbf{U}(t) = \nabla \wedge \mathbf{\Psi}_f(t). $$

The equations of motion in (12) are then

$$ \tilde{\rho}_{11}(t) \frac{\partial^2 \mathbf{\Psi}_s(x, t)}{\partial t^2} + \tilde{\rho}_{12}(t) \frac{\partial^2 \mathbf{\Psi}_f(x, t)}{\partial t^2} = -N \nabla \wedge (\nabla \wedge \mathbf{\Psi}_s(t)), \qquad (23) $$

$$ \tilde{\rho}_{22}(t) \frac{\partial^2 \mathbf{\Psi}_f(x, t)}{\partial t^2} + \tilde{\rho}_{12}(t) \frac{\partial^2 \mathbf{\Psi}_s(x, t)}{\partial t^2} = 0. \qquad (24) $$

A relation of proportionality between the solid and fluid displacement potentials is obtained from Equation (24) by

$$ \mathbf{\Psi}_f(x, t) = V_3(t) \mathbf{\Psi}_s(x, t), \qquad (25) $$

where

$$ V_3(t) = - \frac{\widetilde{\rho_{12}(t)}}{\widetilde{\rho_{22}(t)}}. \qquad (26) $$

This last relation allows us to write the system of Equation (12) in the following form:

$$ \Delta \mathbf{\Psi}_s(x, t) - \frac{1}{N} \left( \widetilde{\rho_{11}(t)} - \frac{\widetilde{\rho_{12}^2(t)}}{\widetilde{\rho_{22}(t)}} \right) \frac{\partial^2 \mathbf{\Psi}_s(x, t)}{\partial t^2} = 0. \qquad (27) $$
Using the Laplace transform, where $\tilde{\Psi}_s(x,s)$ is the Laplace transform of $\Psi_s(x,t)$ as defined previously in (15) gives the following:

$$ \Delta\Psi_s(x,s) - \frac{1}{N} \left( \widetilde{\rho_{11}(s)} - \frac{\widetilde{\rho_{12}^2(s)}}{\widetilde{\rho_{22}(s)}} \right) s^2 \Psi_s(x,s) = 0. \quad (28) $$

Posing

$$ \begin{aligned} \lambda_3(s) &= \frac{1}{N} \left( \widetilde{\rho_{11}(s)} - \frac{\widetilde{\rho_{12}^2(s)}}{\widetilde{\rho_{22}(s)}} \right) s^2, \\ &= C_3 s^2 + D_3 s^{\frac{3}{2}} + G_3 s, \end{aligned} $$

where

$$ C_3 = \frac{1}{N} \left( \rho_{11} - \frac{\rho_{12}^2}{\rho_{22}} \right), D_3 = \frac{A}{N} \left( 1 + \frac{\rho_{12}^2}{\rho_{22}^2} + 2 \frac{\rho_{12}}{\rho_{22}} \right), G_3 = \frac{A}{N \rho_{22}} \left( 1 + \frac{\rho_{12}^2}{\rho_{22}^2} + 2 \frac{\rho_{12}}{\rho_{22}} \right). $$

The solution in the Laplace domain takes the following form

$$ \tilde{\Psi}_s(x,s) = \tilde{\Psi}_1(s)e^{-x\sqrt{\lambda_3(s)}} + \tilde{\Psi}_2(s)e^{x\sqrt{\lambda_3(s)}}. \quad (29) $$

The proportionality relation in (25) in the Laplace domain becomes

$$ \Psi_f(x,s) = V_3(s) \Psi_s(s,t), \quad (30) $$

with

$$ V_3(s) = - \frac{\widetilde{\rho_{12}(s)}}{\widetilde{\rho_{22}(s)}} = A_3 + \frac{B_3}{\sqrt{s}}, \quad (31) $$

and

$$ A_3 = - \frac{\rho_{12}}{\rho_{22}}, \quad B_3 = A \left( \frac{\rho_{12}}{\rho_{22}^2} + \frac{1}{\rho_{22}} \right). $$

The following section focuses on the determination of the diffusion operators in the time domain; the expressions of the longitudinal and transverse waves in (19), (20), and (29); the expressions of the pressure and velocity fields in the different media; as well as the boundary conditions at the interfaces.

## 4. The Scattering Operators

In a two-dimensional problem (X-Z), a poroelastic material of thickness L, which is supposed to be homogeneous and isotropic, is saturated by a viscous fluid and found to be disturbed by a high-frequency sound wave that arrives in normal incidence on the left boundary of the material ($x = 0$); the geometry of the problem is represented in the Figure 1.
<figure><img src="image_5.png" /><figcaption>Figure 1. Geometry of the problem.</figcaption></figure>

A part of the incident wave is reflected in the fluid and the other one is transmitted through the porous medium in the shape of the three waves of Biot: the two compressional waves and a shear wave. By pointing out the fact that the porous medium occupies the region (0 ≤ x ≤ L), we can see that the expression of the pressure field in the region (I), where x < 0, is the sum of the incident p^i and the reflected p^r fields.

$$p_1(x,t) = p^i \left( t - \frac{x}{c_0} \right) + p^r \left( t + \frac{x}{c_0} \right), \quad x < 0,$$

where $p_1(x,t)$ is the total pressure field at the left of the medium. In the limit $x > L$, the pressure field results from the transmission of the incident wave through the porous material and is given by the following:

$$p_3(x,t) = p^t \left( t - \frac{x-L}{c_0} \right), \quad x > L,$$

where $p^t$ is the pressure of the transmitted wave.

The expressions of the reflected and transmitted pressure fields can be established by integrating the product of the scattering operators ($R$ and $T$) by the incident wave, which gives the following:

$$p^r(x,t) = \int_0^t \tilde{R}(\tau)p^i\left(t - \tau + \frac{x}{c_0}\right)d\tau \quad (32)$$

$$p^t(x,t) = \int_0^t \tilde{T}(\tau)p^i\left(t - \tau - \frac{L}{c} - \frac{(x-L)}{c_0}\right)d\tau. \quad (33)$$

The fact that the lower limit of the integrals of these expressions is zero indicates that it is exactly at $t = 0$ that the incident wavefront hits the material for the first time.

Accordingly, $p_1(x,t)$ and $p_3(x,t)$ are expressed by the following equations:

$$p_1(x,t) = \int_0^t \left[ \delta\left(\tau - \frac{x}{c_0}\right) + \int_0^\tau \tilde{R}(\tau')\delta\left(\tau - \tau' + \frac{x}{c_0}\right)d\tau' \right] p^i(t - \tau)d\tau$$

$$p_3(x,t) = \int_0^t \left[ \left( \int_0^\tau \tilde{T}_1(\tau')p^i(x, \tau - \tau')d\tau' \right) \delta\left( t - \tau \frac{L}{c_1} - \frac{x-L}{c_0} \right) \right] d\tau \\ + \int_0^t \left[ \left( \int_0^\tau \tilde{T}_2(\tau')p^i(x, \tau - \tau')d\tau' \right) \delta\left( t - \tau \frac{L}{c_2} - \frac{x-L}{c_0} \right) \right] d\tau.$$
$T_1(t)$ and $c_1$ are the transmission operator and the propagation speed of the fast wave, respectively. Similarly, $T_2(t)$ and $c_2$ are the transmission operator and the propagation speed of the slow wave, respectively.

To simplify the analysis, as previously achieved, we use the Laplace transform, which is more adequate for our problem. We note that, following the definition in (15), $P(x, s)$ is the Laplace transform of $p(x, t)$. The pressure fields on both sides of the medium in the Laplace domain are given by the following:

$$P_1(x, s) = (e^{-\frac{x}{c_0}s} + R(s)e^{\frac{x}{c_0}s})P(s), \qquad (34)$$

$$P_3(x, s) = T(s)P(s)e^{-\frac{x-L}{c_0}s}, \qquad (35)$$

with

$$T(s) = \left( T_1(s)e^{-\frac{L}{c_1}s} + T_2(s)e^{-\frac{L}{c_2}s} \right).$$

$P_1(x, s)$ and $P_3(x, s)$ are, respectively, the Laplace transform of the field at the left and right of the material; $P(s)$ is the Laplace transform of the incident field; and finally, $R(s)$ and $T(s)$ are the Laplace transforms of the reflection and transmission operators.

The stress-strain equations in the porous medium are given by

$$\begin{aligned} \sigma_{ij}^{s} &= ((P - 2N)\nabla \cdot \mathbf{u} + Q\nabla \cdot \mathbf{U})\delta_{ij} + N(u_{i,j} + u_{j,i}), \\ \sigma_{ij}^{f} &= (R\nabla \cdot \mathbf{U} + Q\nabla \cdot \mathbf{u}) = -p_f \delta_{ij}, \end{aligned} \qquad (36)$$

where $\sigma_{ij}^s$ are the components of the stress tensor that act on the solid (a = s) or fluid (a = f) phase and $\epsilon_{ij} = \frac{1}{2}(u_{i,j} + u_{j,i})$ is the skeletal strain tensor.

Since the pressure and stress fields at the boundaries of the medium are continuous [39], the relations between the pressure field and the stresses in the Laplace domain at the interfaces $x = 0$ and $x = L$ are given by the following:

$$\begin{aligned} \sigma^s(0^+, s) &= -(1-\phi)P_1(0^-, s), \\ \sigma^f(0^+, s) &= -\phi P_1(0^-, s), \\ \tilde{\sigma}^s(0^+, s) &= 0 \\ \sigma^s(L^-, s) &= -(1-\phi)P_3(L^+, s), \\ \sigma^f(L^-, s) &= -\phi P_3(L^+, s), \\ \tilde{\sigma}^s(L^-, s) &= 0 \end{aligned} \qquad (37)$$

$\sigma^s$ and $\sigma^f$ are the normal stresses in the skeleton and fluid, respectively, and $\tilde{\sigma}^s$ is the shear stress in the skeleton; their expressions in the Laplace domain take the following form:

$$\begin{aligned} \sigma_{xx}^{s}(x, s) &= (P - 2N) \frac{\partial^2 \Phi_s(x, s)}{\partial x^2} + Q \frac{\partial^2 \Phi_f(x, s)}{\partial x^2} + 2N \frac{\partial^2 \Phi_s(x, s)}{\partial x^2}, \\ \sigma_{xx}^{f}(x, s) &= Q \frac{\partial^2 \Phi_s(x, s)}{\partial x^2} + R \frac{\partial^2 \Phi_f(x, s)}{\partial x^2}, \\ \tilde{\sigma}^{s}(x, s) &= 2N \left( \frac{\partial^2 \Phi_s(x, s)}{\partial x^2} + \frac{\partial^2 \tilde{\Psi}_s(x, s)}{\partial x^2} \right), \end{aligned}$$

where $\Phi_s$, $\Phi_f$ and $\tilde{\Psi}_s$ are the Laplace transforms of $\phi_s$, $\phi_f$, and $\Psi_s$.

Hence, using the boundary conditions in (38), the expressions of the scalar functions $\Phi_{11}(s)$, $\Phi_{12}(s)$, $\Phi_{21}(s)$, $\Phi_{22}(s)$, $\tilde{\Psi}_1(s)$ and $\tilde{\Psi}_1(s)$ are calculated in Appendix B.
To determine the reflection and transmission coefficients, we consider the conservation conditions of the flow at the limits of the medium ($x = 0$ and $x = L$):

$$V_1(0^+, s) = (1 - \phi)V_s(L^-, s) + \phi V_f(L^-, s), \quad (38)$$

$$V_3(L^-, s) = (1 - \phi)V_s(L^-, s) + \phi V_f(L^-, s).$$

In those expressions, $V_1(x, s)$ and $V_3(x, s)$ are, respectively, the acoustic velocity fields in the regions $x < 0$ and $x > L$, while $V_s(x, s)$ and $V_f(x, s)$ are the velocities of the solid and fluid in the porous medium, respectively.

The acoustic velocity fields are obtained by considering the field pressure of the surrounding fluid in (34) and (35), and the Euler equation:

$$\rho_f s V_i(x, s) = - \frac{\partial P_i(x, s)}{\partial x}, \quad i = 1, 2, 3.$$

The expressions of $V_s$ and $V_f$ are found using the following expressions:

$$V_s(x, s) = s(\nabla\Phi_s + \nabla \wedge \Psi_s), \quad (39)$$

$$V_f(x, s) = s(\nabla\Phi_f + \nabla \wedge \Psi_f).$$

Using the boundary conditions in (38) and the expressions for the fluid and solid velocity fields, $V_s(x, s)$, $V_f(x, s)$, $V_1(x, s)$, and $V_3(x, s)$, we obtain the expressions for the transmission and reflection coefficients in the Laplace domain:

$$R(s) = \frac{s^2(F_5^2(s) - F_4^2(s)) + 1}{(sF_4(s) - 1)^2 - s^2F_5^2(s)}, \quad (40)$$

$$T(s) = \frac{-2sF_5(s)}{(sF_4(s) - 1)^2 - s^2F_5^2(s)}, \quad (41)$$

where

$$F_4(s) = F_1(s) \cosh(l\sqrt{\lambda_1(s)}) + F_2(s) \cosh(l\sqrt{\lambda_2(s)}) + F_3(s) \cosh(l\sqrt{\lambda_3(s)}),$$

$$F_5(s) = F_1(s) + F_2(s) + F_3(s),$$

$$F_i(s) = \rho_f c_0 [1 + \phi(V_i(s) - 1)] \sqrt{\lambda_1(s)} \frac{2\Psi_i''(s)}{\sinh(l\sqrt{\lambda_i(s)})\Psi''(s)},$$

$$V_i(s) = A_i + \frac{B_i}{\sqrt{s}}, \dots, i = 1, 2,$$

$$\lambda_i(s) = C_i s^2 + D_i s^2 + G_i s, \dots, i = 1, 2$$

$$V_3(s) = A_3 + \frac{B_3}{\sqrt{s}},$$

$$\lambda_3(s) = C_3 s^2 + D_3 s^2 + G_3 s.$$

The coefficients $\Psi_1''(s)$, $\Psi_2''(s)$, $\Psi_3''(s)$ and $\Psi''(s)$ are given by

$$\Psi_1''(s) = \phi Z_2(s) - (1 - \phi) Z_4(s),$$

$$\Psi_2''(s) = (1 - \phi) Z_3(s) - \phi Z_1(s),$$

$$\Psi_3''(s) = (1 - \phi) Z_8(s) - \phi Z_9(s),$$

$$\Psi''(s) = 2(Z_1(s)Z_4(s) - Z_2(s)Z_3(s)).$$
The coefficients $Z_i(s)$ are given by

$$
\begin{aligned}
Z_1(s) &= (V_1 Q + P)\lambda_1(s), Z_2(s) = (V_2 Q + P)\lambda_2(s), \\
Z_3(s) &= (V_1 R + Q)\lambda_1(s), Z_4 = (V_2 R + Q)\lambda_2(s), \\
Z_5(s) &= 2N\lambda_1(s), Z_6(s) = 2N\lambda_2(s), Z_7 = 2N\lambda_3(s), \\
Z_8(s) &= \frac{Z_6(s)Z_3(s) - Z_5(s)Z_4(s)}{Z_7(s)}, \\
Z_9(s) &= \frac{Z_6(s)Z_1(s) - Z_5(s)Z_2(s)}{Z_7(s)}.
\end{aligned}
$$

The scattering operators are obtained in the Laplace domain. The purpose of our work is to determine these operators in the time domain, but using the inverse Laplace transform to go directly to the time domain is a very difficult task to perform because of the complexity of these operators. For this reason, modifications and simplifications are applied (Appendix C) to the relations in (40) and (41), which allows us to find the time equivalent of the reflection and transmission operators:

$$R(t) = r(t) + \tilde{R}(t), \qquad (42)$$

where

$$r(t) = \frac{1+U}{1-U}\delta(t) + \frac{2a}{1-U}\left[-\frac{1}{\sqrt{\pi t}} + a \exp(a^2 t) \text{Erfc}(a\sqrt{t})\right], \qquad (43)$$

and

$$\tilde{R}(t) = \frac{4}{(1-U)^2} [x_1 G_1(t, 2L) + x_2 G_2(t, 2L) + x_3 G_3(t, 2L)] \\ + \frac{4}{(1-U)^2} [P_1(t) * G_1(t, 2L) + P_2(t) * G_2(t, 2L) + P_3(t) * G_3(t, 2L)]. \qquad (44)$$

$$T(t) = - \frac{4}{(1-U)^2} [x_1 G_1(t, L) + x_2 G_2(t, L) + x_3 G_3(t, L)] \\ - \frac{4}{(1-U)^2} [P_1(t) * G_1(t, L) + P_2(t) * G_2(t, L) + P_3(t) * G_3(t, L)]. \qquad (45)$$

The expressions of $P_i(t)$ and $G_i(t, L)$ are given by

$$P_i(t) = \frac{y_i - 2ax}{\sqrt{\pi t}} + 2a^2(y_i - ax_i)\sqrt{\frac{t}{\pi}} + [a^2x_i(3 + 2a^2t) - 2ay_i(1 + a^2t)] \\ \times \exp(a^2t)\text{Erfc}(a\sqrt{t}), i = 1, 2, 3 \qquad (46)$$

$$G_i(t, jL) = L^{-1}\left[\exp(-jL\sqrt{\lambda_1(s)})\right], i = 1, 2, 3, j = 1, 2 \qquad (47)$$

$G_i(t), i = 1, 2, 3$, are the Green functions of compressional (first and second species) and shear waves, respectively; they describe the propagation of each of these waves inside the porous medium, and their expressions are given in Appendix D.

$G_i(t), i = 1, 2, 3$, are the Green functions of the longitudinal (fast and slow) and transverse waves, respectively. They describe the propagation of each of these waves inside the porous medium, and their expressions are given in Appendix D.

Knowing the strong attenuation that the acoustic wave undergoes in a porous media saturated by a viscous fluid makes multiple reflections negligible; thus, we take only the first reflections at the interfaces $x = 0$ and $x = L$ in the expressions of the operators in (42) and (45).
The first term of the second member of Equation (43) corresponds to the instantaneous reflection of the wave by the first interface $x = 0$. The reflected wave does not have time to propagate in the medium, which gives it the advantage of not being dispersive but only attenuated by the term $\frac{1+U}{1-U}$.

The second term of $r(t)$ expresses the attenuation and dispersion process since it depends on the fractional derivative operator, which takes into account the viscous fluid-structure interaction effects during the wave propagation. This memory term expresses the phenomenon of relaxation of the wave reflected by the first interface.

The term of the second member of Equation (45) corresponds to the reflection by the second interface $x = L$. This term depends on the Green function of the longitudinal and transverse waves that describe the propagation and dispersion of these waves inside the porous material just as $\tilde{R}(t)$ the temporal expression of the transmission operator in (45) depends on the Green functions of the three Biot waves, highlighting the contribution of each of these waves to the total transmitted wave.

From the temporal expressions of the scattering operators obtained in (42) and (45), the effect of the shear wave plus the effect of the compressional waves (fast and slow) on the porous medium's response to an external load can be examined. It is also possible to study the influence of the characteristic parameters of the porous medium on the amplitudes and speed of these waves.

## 5. Numerical Validation

In order to test the validity of our analytical development and to show the influence that the shear waves can have on the response of bone-type porous materials, we proceed to compare transmitted signals obtained during this study, the shear waves being taken into account, and the transmitted signals, where only the longitudinal waves are observed because the shear waves are very weak in normal incidence; however, whether or not they are taken into account is visible on the longitudinal waveforms, as shown in the following study. The characteristics of the bone samples used in this study were taken from work carried out previously by Fellah [67] and are shown in Table 1.

<table><caption>**Table 1.** Biot's model parameters of cancellous bone sample (S1).</caption><thead><tr><th>Parameters</th><th>L (mm)</th><th>φ</th><th>α<sub>∞</sub></th><th>Λ (μm)</th><th>ρ<sub>s</sub> Kg/m<sup>3</sup></th><th>K<sub>s</sub> (GPa)</th><th>K<sub>b</sub> (GPa)</th><th>N (GPa)</th></tr></thead><tbody><tr><td>Bone (S1)</td><td>0.7</td><td>0.85</td><td>1.3</td><td>8</td><td>1960</td><td>2.45</td><td>1</td><td>6.7</td></tr></tbody></table>

The transmitted signals are computed by performing the convolution product of the transmission operator given by (45) with the incident signal, as in the relation in (33).

### 5.1. Porosity Variation

The two waves observed in the numerical simulations (Figure 2) correspond to the slow and fast longitudinal waves. The shear wave is not noticeable because of its low amplitude at normal incidence; however, its effect is visible on the velocities and amplitudes of the longitudinal waves. This is the subject of the following sensitivity study.
<figure><img src="image_6.png" /><figcaption>**Figure 2.** Comparison between simulated transmitted signals by considering or not considering the transverse effects for two values of the porosity: $\phi = 0.85$ and $\phi = 0.2$.</figcaption></figure>

The transmitted signals are simulated by taking or not taking into account the transverse wave for two distinct porosities: $\phi = 0.2$ and $\phi = 0.85$. Figure 2 depicts the impact of the shear wave on different values of the porosity. When the porosity value is equal to 0.2, the shear effect is not felt in the slow longitudinal wave, and its amplitude and speed do not change, unlike the amplitude and speed of the fast wave. For the value of the porosity of 0.85, we notice that, whether the shear wave is considered is much more important than in the previous case (when the porosity was equal to 0.2) for the amplitudes of the fast and slow longitudinal waves. We can therefore conclude at this stage that the effect of shear is more visible as the value of porosity is high and at the level of the amplitudes of both longitudinal waves. The velocities of the longitudinal waves do not seem to be affected by whether or not shear is taken into account.

## 5.2. Tortuosity Variation

Figure 3 represents a comparison of the simulated transmitted signals when only the effect of longitudinal waves is considered (black dotted line) and when the effect of the transverse waves is combined with that of the longitudinal waves (red line) for two different values of the tortuosity $a_\infty = 1$ and $a_\infty = 1.3$.

<figure><img src="image_7.png" /><figcaption>**Figure 3.** Comparison between simulated transmitted signals by considering or not considering the transverse effects for two values of tortuosity $a_\infty = 1$ and $a_\infty = 1.3$.</figcaption></figure>

Readers can notice that the amplitudes of the two compressional waves (the fast and slow) are different for the two values of tortuosity, even though it is more apparent for $a_\infty = 1.3$ than it is for $a_\infty = 1$. The velocities of both compressional waves remain unchanged, regardless of whether the effect of the transverse waves is included or not.

Tortuosity is an important parameter to describe the inertial effects, which occur between the fluid and solid parts of a porous material because it provides information on the geometry of the pores of the porous medium, which are rarely straight and have
variable sections. Like porosity variations, tortuosity variations cause changes in the fast and slow waves and their velocities.

Comparing the figures for two values of the tortuosity, we notice that an increase in the tortuosity induces an increase in the amplitude of the fast wave and a decrease in the amplitude of the slow-wave. A temporal shift of the two Biot waves occurs (decrease in the velocities), which is more distinct on the slow wave than it is on the fast wave as the tortuosity increases.

### 5.3. Variation of the Viscous Characteristic Length $\Lambda$

The effect of whether transverse waves are considered along with longitudinal waves on the transmitted signals for two different viscous characteristic lengths is illustrated in Figure 4. In the comparison of the simulated transmitted signals, the effect of shear can be observed on the amplitudes of both the fast and slow waves. This effect leads to an increase in the amplitudes of the fast and slow waves regardless of the value of the viscous characteristic length. The speeds of the fast and slow waves remain unchanged in all cases. The comparison between the figures corresponding to two values of the viscous characteristic length provides information on the effect of the viscous characteristic length on the transmitted signals. This parameter was first introduced by Johnson et al. [33] in 1987 and is considered to be a geometrical parameter indicative of small pore sizes; it represents the thickness where viscous effects are important.

<figure><img src="image_8.png" /><figcaption>**Figure 4.** Comparison Please add a space between number and unit. between simulated transmitted signals by considering or not considering the transverse effects for two values of the viscous characteristic length $\Lambda = 5$ µm and $\Lambda = 9$ µm.</figcaption></figure>

Through our comparison, we can see that, by increasing the viscous characteristic length from 5 µm to 9 µm, the amplitude of the slow wave increases while the velocities of the slow and fast waves as well as the amplitude of the fast wave remain practically the same.

The viscous characteristic length plays a less important role than the tortuosity, since the variation in its value barely affects the transmitted signals whether the shear wave is considered or not.

### 5.4. Solid Density $\rho_s$ Variation

The difference between the transmitted signals when shear effects are considered and not considered along with longitudinal waves for two distinct values of the solid density $\rho_s = 1000 \text{ kg/m}^3$ and $\rho_s = 2400 \text{ kg/m}^3$ is given in Figure 5. Unlike for the previous case of the effect of the viscous characteristic length, we notice that the only significant difference is for the amplitudes of the fast waves that increase when the shear waves are considered for both values of the solid density; the slow waves are practically the same for both signals in Figure 5. It is worth mentioning that, for a high value of the density $\rho_s = 2400 \text{ kg/m}^3$, the fast wave with shear taken into account arrives faster than the one with shear not taken
into account, which means that its speed is greater without than with shear taken into account, which is not observed for the slow wave.

<figure><img src="image_9.png" /><figcaption>**Figure 5.** Comparison between simulated transmitted signals by considering or not considering the transverse effects for two values of the solid density $\rho_s = 1000 \text{ kg/m}^3$ and $\rho_s = 2400 \text{ kg/m}^3$.</figcaption></figure>

We can also notice from these figures that the fast waves (amplitude and speed) are sensitive to the variation of the solid density; indeed, a decrease in the amplitudes of the fast waves as well as a decrease in their speed appear when the solid density goes from $1000 (\text{kg/m}^3)$ to $2400 (\text{kg/m}^3)$.

### 5.5. Bulk Modulus of the Elastic Solid $K_s$ Variation

The simulated transmitted signals with and without considering the effect of transverse waves in addition to longitudinal waves are compared in Figure 6 for two values of the bulk modulus of the elastic solid $K_s$. We notice that the effect of the shear wave is visible mainly on the amplitude of the fast wave for the value of $K_s = 2.45 \text{ GPa}$, and the amplitude of the slow wave varies but in a lesser way. For the value of $K_s = 10 \text{ GPa}$, it is the amplitude of the slow wave that is most sensitive to the incorporation of shear in the model, with the change in the amplitude of the fast wave being, however, smaller for this value of $K_s$. Shear affects the amplitudes of the fast and slow waves differently depending on the values of $K_s$.

<figure><img src="image_10.png" /><figcaption>**Figure 6.** Comparison between simulated transmitted signals by considering or not considering the transverse effects (TE) for two values of the bulk modulus of the porous skeletal frame $K_s = 2.45 \text{ GPa}$ and $K_s = 10 \text{ GPa}$.</figcaption></figure>

The speed of the fast wave also seems to be sensitive to the shear; in fact, when decreasing the value of $K_s$, the fast wave arrives delayed when the shear is not taken into account. In other words, taking into account the shear increases the speed of the fast wave, which makes it arrive faster in comparison with the fast wave obtained without shear. This
phenomenon is not observable for the speed of the fast wave when $K_s = 10$ GPa or for the speeds of the slow waves, regardless of the value of $K_s$.

### 5.6. Bulk Modulus of the Porous Skeletal Frame $K_b$ Variation

The sensitivity study of the bulk modulus of the porous skeletal frame $K_b$ taking or not taking into account the shear shows quite interesting results compared with the previous cases (see Figure 7). For example, we notice as an example that, for the high value of $K_b = 3.6$ GPa, the amplitudes of the slow and fast waves are very affected by whether shear is considered. The amplitude of the slow wave in this case can be doubled when the shear is taken into account, while the amplitude of the slow wave decreases sharply. For the value of $K_b = 1.1$ GPa, the variation in the slow wave amplitude is very small while that of the fast wave is much larger. We also notice generally that the velocities of the waves do not practically change with the taking into account of the shear, except for that of the fast wave for $K_b = 1.1$ GPa, where we see that taking shear into account of the shear decreases the velocity of the fast wave, which slows down compared with the case where the shear is taken into account.

<figure><img src="image_11.png" /><figcaption>**Figure 7.** Comparison between simulated transmitted signals by considering or not considering the transverse effects (TE) for two values of the bulk modulus of the porous skeletal frame $K_b = 1.1$ (GPa) and $K_b = 3.6$ (GPa).</figcaption></figure>

### 5.7. The Shear Modulus $N$ Variation

The simulated transmitted signals including and excluding the effect of the shear wave in addition to the longitudinal waves for two values of the shear modulus are shown in Figure 8. The comparison between the two signals for each value of $N$ allows us to say that the inclusion of the shear wave in the model affects the fast wave more than the slow wave and more particularly for large values of $N$ ($N = 8.9$ GPa). However, the effect of shear is more visible on the speed of the fast wave for low values of $N$ ($N = 5.8$ GPa). However, the slow wave velocities remain unchanged whether shear is considered or not.

<figure><img src="image_12.png" /><figcaption>**Figure 8.** Comparison between simulated transmitted signals by considering or not the transverse effects (TE) for two values of the bulk modulus of the porous skeletal frame $N = 5.8$ (GPa) and $N = 8.9$ (GPa).</figcaption></figure>
**6. Discussion and Conclusions**

This study allows us to examine the influence of the shear wave in the modeling of ultrasound propagation in porous media. This effect is studied by comparing the transmitted waveforms by varying each physical parameter of the porous material (porosity, tortuosity, viscous characteristic length, density of the solid, bulk modulus of the elastic solid, bulk modulus of the porous skeletal frame, and shear modulus). The amplitudes and velocities of the transmitted longitudinal waves (slow and fast) are observed and compared in order to see the impact of taking into account the shear wave in the model. In our previous work [19,62], we did not take into account the shear wave at normal incidence, but recent studies [62,67] have shown that the shear wave must be taken into account because of the viscous exchanges between fluid and structure. At normal incidence, the shear wave is not as visible as the fast and slow longitudinal waves because its amplitude is very low; however, its effect is felt in the amplitudes and speeds of longitudinal waves, and this effect varies depending on the values of the physical parameters of the porous medium, as shown in this study.

The originality of our approach is that it makes it possible not only to model in the time domain acoustic propagation in a porous medium using fractional calculus but also to obtain analytical solutions to the response (reflection and transmission operators) of the porous medium to an acoustic excitation (incident wave). This analytical solution in the time domain expressed by the Green functions of longitudinal (fast and slow) and transverse waves using fractional calculation is the strong point (original) of our work. Indeed, the existing studies do not provide analytical solutions but only propose numerical solutions.

The study of the sensitivity of the physical parameters of the longitudinal waves by taking or not taking into account the shear wave in the model allows us to see and quantify the importance of the shear wave in the modeling and its impact on the longitudinal waves. It appears that the variation in the values of some parameters such as the bulk modulus of the porous skeletal frame reflects a significant change in the amplitude of the fast and slow longitudinal waves. The velocity of the fast wave is also modified, contrary to that of the fast wave, which is unchanged, and thus independent of whether shear is taken into account. The variation in other parameters (such as porosity, tortuosity, viscous characteristic length, solid density, and shear modulus) modifies the amplitude of the fast wave more than that of the slow wave. Practically for all the parameters, their variation has no impact on the speed of the slow wave, only on that of the fast wave.

The results of this study clearly show that the inclusion of the shear wave in the model has a considerable impact on the amplitude of the slow and fast waves as well as on the speed of the fast wave. The fractional calculus have proven to be very useful in obtaining the analytical solutions in the time domain. We can then conclude that the shear wave at normal incidence should no longer be ignored as it was the case in the previous studies and that it is necessary to consider the shear effects, as this will have consequences for the acoustic and mechanical characterization of porous materials by solving the inverse problems using experimental data of transient waves.

**Author Contributions:** Writing—original draft preparation, D.B. (Djihane Benmorsli), D.B. (Djema Belgroune), Z.E.A.F., M.F., E.O., N.O.O. and C.D.; analysis, analytical calculations, and processing of results, D.B. (Djihane Benmorsli), D.B. (Djema Belgroune), Z.E.A.F., M.F., E.O. and C.D.; writing—review and editing, D.B. (Djihane Benmorsli), D.B. (Djema Belgroune), Z.E.A.F., M.F., E.O. and C.D.; resources, D.B. (Djihane Benmorsli), D.B. (Djema Belgroune), Z.E.A.F., M.F., E.O. and C.D. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.
**Conflicts of Interest:** The authors declare no conflicts of interest.

## Appendix A. The Eigenvalues and Eigenvectors of the Matrix M

The matrix M is given in the Laplace domain by the following:

$$M = \begin{pmatrix} (R'\rho_{11} - Q'\rho_{12})s^2 + As^3(R' + Q') & (R'\rho_{12} - Q'\rho_{22})s^2 - As^3(R' + Q') \\ (P'\rho_{12} - Q'\rho_{11})s^2 - As^3(P' + Q') & (P'\rho_{22} - Q'\rho_{12})s^2 + As^3(P' + Q') \end{pmatrix},$$

with

$$\begin{pmatrix} R' & -Q' \\ -Q' & P' \end{pmatrix} = \begin{pmatrix} P & Q \\ Q & R \end{pmatrix}^{-1} = \begin{pmatrix} \frac{R}{PR-Q^2} & -\frac{Q}{PR-Q^2} \\ -\frac{Q}{PR-Q^2} & \frac{P}{PR-Q^2} \end{pmatrix},$$

and

$$A = \frac{2\phi\alpha_{\infty}\rho_f}{\Lambda} \left(\frac{\eta}{\rho_f}\right)^{\frac{1}{2}}.$$

The eigenvalues of the matrix M are written as follows:

$$\lambda_{1,2}(s) = \frac{-\mathrm{tr}(M) \pm \sqrt{\mathrm{tr}(M)^2 - 4 \det(M)}}{2},$$

where Tra(M) and Det(M) are, respectively, the trace and determinant of this same matrix:

$$\mathrm{Tra}(M) = (R'\rho_{11} - 2Q'\rho_{12} + P'\rho_{22})s^2 + As^3(R' + 2Q' + P'),$$

$$\mathrm{Det}(M) = (P'R' - Q'^2) \left[ (\rho_{11}\rho_{12} - \rho_{12}^2)s^4 + As^2(\rho_{11} + \rho_{22} - 2\rho_{12}) \right].$$

Posing

$$\tau_1 = R'\rho_{11} - 2Q'\rho_{12} + P'\rho_{22},$$

$$\tau_2 = A(R' + 2Q' + P'),$$

$$\tau_3 = (P'R' - Q'^2)(\rho_{11}\rho_{12} - \rho_{12}^2),$$

$$\tau_4 = (P'R' - Q'^2)(\rho_{11} + \rho_{22} - 2\rho_{12})A,$$

and

$$C_i = \frac{1}{2} \left( \tau_1 + (-1)^i \sqrt{(\tau_1^2 - 4\tau_3)} \right),$$

$$D_i = \frac{1}{2} \left( \tau_2 + (-1)^i \frac{\tau_1 \tau_2 - 2\tau_4}{\sqrt{(\tau_1^2 - 4\tau_3)}} \right),$$

$$G_i = (-1)^i \frac{1}{4} \left( \frac{\tau_2^2}{\sqrt{(\tau_1^2 - 4\tau_3)}} - \frac{(\tau_1 \tau_2 - 2\tau_4)^2}{2(\tau_1^2 - 4\tau_3)^{\frac{3}{2}}} \right),$$

we obtain the expressions of the eigenvalues:

$$\lambda_i(s) = C_i s^2 + D_i s^3 + G_i s, i = 1, 2.$$

The expressions of the components $V_{1,2}$ of the eigenvectors are given by

$$V_{1,2}(s) = \frac{\lambda_{1,2}(s) - \left[ (R'\rho_{11} - Q'\rho_{12})s^2 + As^3(R' + Q') \right]}{(R'\rho_{12} - Q'\rho_{22})s^2 - As^3(R' + Q')}.$$
By posing

$$ \tau_5 = (R'\rho_{11} - Q'\rho_{12}), $$

$$ \tau_6 = A(R' + Q'), $$

$$ \tau_7 = (R'\rho_{12} - Q'\rho_{22}), $$

$$ A_i = \frac{\tau_1 - 2\tau_5 + (-1)^i \sqrt{(\tau_1^2 - 4\tau_3)}}{2\tau_7}, $$

$$ B_i = \frac{1}{4\tau_7^2} \left[ \left( \tau_2 - 2\tau_6 + (-1)^i \frac{(\tau_1\tau_2 - 2\tau_4)}{(\tau_1^2 - 4\tau_3)^{\frac{1}{2}}} \right) 2\tau_7 + \left( \tau_1 - 2\tau_5 - (\tau_1^2 - 4\tau_3)^{\frac{1}{2}} 2\tau_6 \right) \right], i = 1, 2 $$

we obtain

$$ V_i(s) = A_i + \frac{B_i}{\sqrt{s}}, \quad i = 1, 2. \qquad (A1) $$

## Appendix B. Scalar Functions

The components of the stress tensor $\sigma_{ij}^s$ that act on the solid ($a = s$) or fluid ($a = f$) phase are given by

$$ \sigma_{ij}^s = (P - 2N)\vec{\nabla} \cdot \vec{u} + Q\vec{\nabla} \cdot \vec{U} \delta_{ij} + N(u_{i,j} + u_{j,i}), $$

$$ \sigma_{ij}^f = (R\vec{\nabla} \cdot \vec{U} + Q\vec{\nabla} \cdot \vec{u}) = -p_f \delta_{ij}. $$

The stress–strain relationship for the fluid in the Laplace domain is in the following form:

$$ \tilde{\sigma}_{xx}^f(x, s) = (R\vec{\nabla} \cdot \vec{U} + Q\vec{\nabla} \cdot \vec{u}), $$

with

$$ u(x, s) = \nabla \Phi_s(x, s), \qquad (A2) $$

$$ U(x, s) = \nabla \Phi_f(x, s), $$

$$ \sigma_{xx}^f(x, s) = R \frac{\partial^2 \Phi_s(x, s)}{\partial x^2} + Q \frac{\partial^2 \Phi_f(x, s)}{\partial x^2}, $$

where

$$ \Phi_s(x, s) = \Phi_1(x, s) + \Phi_2(x, s), $$

$$ \Phi_f(x, s) = V_1(s)\Phi_1(x, s) + V_2(s)\Phi_2(x, s), $$

and we find

$$ \sigma_{xx}^f(x, s) = (V_1(s)R + Q) \frac{\partial^2 \Phi_1(x, s)}{\partial x^2} + (V_2(s)R + Q) \frac{\partial^2 \Phi_2(x, s)}{\partial x^2}. $$

By replacing $\Phi_1(x, s)$ and $\Phi_2(x, s)$ by their expressions in (19) and (20), we find

$$ \sigma_{xx}^f(x, s) = Z_3(s)(\Phi_{11} + \Phi_{12}) \cosh(x\sqrt{\lambda_1(s)}) - Z_3(\Phi_{11} - \Phi_{12}) \sinh(x\sqrt{\lambda_1(s)}) + Z_4(\Phi_{21} + \Phi_{22}) \cosh(x\sqrt{\lambda_2}) - Z_4(\Phi_{21} - \Phi_{22}) \sinh(x\sqrt{\lambda_2}) $$
with:

$$Z_3(s) = (V_1 R + Q) \lambda_1(s),$$

$$Z_4(s) = (V_2 R + Q) \lambda_2(s),$$

The stress-strain relation for the solid in the Laplace domain is in the following form:

$$\sigma_{xx}^{s}(x,s) = (P - 2N)\vec{\nabla} \cdot \vec{u}(x,s) + Q\vec{\nabla} \cdot \vec{U}(x,s) + 2N \frac{\partial u(x,s)}{\partial x},$$

where the solid displacement vectors are given by

$$u(x,s) = \nabla \Phi_s(x,s), \qquad (A3)$$

$$U(x,s) = \nabla \Phi_f(x,s).$$

The normal component of the stress tensor acting on the solid part becomes the following:

$$\sigma_{xx}^{s}(x,s) = (P - 2N)\vec{\nabla} \cdot \vec{u}(x,s) + Q\vec{\nabla} \cdot \vec{U}(x,s) + 2N \frac{\partial u(x,s)}{\partial x},$$

$$\sigma_{xx}^{s}(x,s) = P \frac{\partial^2 \Phi_s(x,s)}{\partial x^2} + Q \frac{\partial^2 \Phi_f(x,s)}{\partial x^2}.$$

Considering the expressions for $\Phi_1(x,s)$ and $\Phi_2(x,s)$, we obtain the following:

$$\sigma_{xx}^{s}(x,s) = Z_1(s)(\Phi_{11} + \Phi_{12}) \cosh(x\sqrt{\lambda_1(s)}) - Z_1(s)(\Phi_{11} - \Phi_{12}) \sinh(x\sqrt{\lambda_1(s)}) \\ + Z_2(s)(\Phi_{21} + \Phi_{22}) \cosh(x\sqrt{\lambda_2(s)}) - Z_2(s)(\Phi_{21} - \Phi_{22}) \sinh(x\sqrt{\lambda_2(s)}),$$

with

$$Z_1(s) = (V_1(s)Q + P)\lambda_1(s),$$

$$Z_2(s) = (V_2(s)Q + P)\lambda_2(s).$$

For shear wave, the stress-strain relationship is written as follows:

$$\tilde{\sigma}^{s}(x,s) = N(u_{i,j} + u_{j,i})$$

$$\tilde{\sigma}^{s}(x,s) = 2N \frac{\partial u(x,s)}{\partial x},$$

Using the fact that the displacement vector of the solid is

$$u = \nabla \Phi_s + \nabla \wedge \Psi_s,$$

the shear stress becomes

$$\tau_{xx}(x,s) = \left[ Z_5(s)(\Phi_{11} + \Phi_{12}) \cosh(x\sqrt{\lambda_1(s)}) - Z_5(s)(\Phi_{11} - \Phi_{12}) \sinh(x\sqrt{\lambda_1(s)}) \right. \\ \left. + Z_6(s)(\Phi_{21} + \Phi_{22}) \cosh(x\sqrt{\lambda_2(s)}) - Z_6(s)(\Phi_{21} - \Phi_{22}) \sinh(x\sqrt{\lambda_2(s)}) \right. \\ \left. + Z_7(s)(\Psi_1 + \Psi_2) \cosh(x\sqrt{\chi(s)}) - Z_7(s)(\Psi_1 - \Psi_2) \sinh(x\sqrt{\chi(s)}) \right],$$

with

$$Z_5 = 2N\lambda_1, Z_6 = 2N\lambda_2, Z_7 = 2N\chi.$$
The continuity conditions for solid and fluid stresses at $x = 0$ and $x = L$ are as follows:

$$ \sigma_{xx}^{s}(0^{+}, s) = -(1 - \phi)P_{1}(0^{-}, s), \qquad (A4) $$

$$ \sigma_{xx}^{f}(0^{+}, s) = -\phi P_{1}(0^{-}, s), $$

$$ \tau_{xx}(0^{+}, s) = 0. $$

$$ \sigma_{xx}^{s}(L^{-}, s) = -(1 - \phi)P_{3}(L^{+}, s), \qquad (A5) $$

$$ \sigma_{xx}^{f}(L^{-}, s) = -\phi P_{3}(L^{+}, s), $$

$$ \tau_{xx}(L^{-}, s) = 0. $$

By using the expressions of the constraints and those of the pressure fields in (34) and (35), we obtain

$$ \begin{aligned} \Phi_{11}(s) &= P(s) \frac{\Psi_1''}{\Psi''} \left[ (1 + R(s)) \left( 1 + \frac{\cosh(l\sqrt{\lambda_1(s)})}{\sinh(l\sqrt{\lambda_1(s)})} \right) - \frac{T(s)}{\sinh(l\sqrt{\lambda_1(s)})} \right], \\ \Phi_{12}(s) &= P(s) \frac{\Psi_1''}{\Psi''} \left[ (1 + R(s)) \left( 1 - \frac{\cosh(l\sqrt{\lambda_1(s)})}{\sinh(l\sqrt{\lambda_1(s)})} \right) + \frac{T(s)}{\sinh(l\sqrt{\lambda_1(s)})} \right], \\ \Phi_{21}(s) &= P(s) \frac{\Psi_2''}{\Psi''} \left[ (1 + R(s)) \left( 1 + \frac{\cosh(l\sqrt{\lambda_2(s)})}{\sinh(l\sqrt{\lambda_2(s)})} \right) - \frac{T(s)}{\sinh(l\sqrt{\lambda_2(s)})} \right], \end{aligned} \qquad (A6) $$

$$ \begin{aligned} \Phi_{11}(s) &= P(s) \frac{\Psi_2''}{\Psi''} \left[ (1 + R(s)) \left( 1 - \frac{\cosh(l\sqrt{\lambda_2(s)})}{\sinh(l\sqrt{\lambda_2(s)})} \right) + \frac{T(s)}{\sinh(l\sqrt{\lambda_2(s)})} \right], \\ \Psi_1(s) &= P(s) \frac{\Psi_3''}{\Psi''} \left[ (1 + R(s)) \left( 1 + \frac{\cosh(l\sqrt{\chi(s)})}{\sinh(l\sqrt{\chi(s)})} \right) - \frac{T(s)}{\sinh(l\sqrt{\chi(s)})} \right], \end{aligned} $$

$$ \Psi_2(s) = P(s) \frac{\Psi_3''}{\Psi''} \left[ (1 + R(s)) \left( 1 - \frac{\cosh(l\sqrt{\chi(s)})}{\sinh(l\sqrt{\chi(s)})} \right) + \frac{T(s)}{\sinh(l\sqrt{\chi(s)})} \right], $$

with

$$ \Psi_{1}''(s) = \phi Z_{2}(s) - (1 - \phi) Z_{4}(s), \qquad (A7) $$

$$ \Psi_{2}''(s) = (1 - \phi) Z_{3}(s) - \phi Z_{1}(s), $$

$$ \Psi_{3}''(s) = (1 - \phi) Z_{8}(s) - \phi Z_{9}(s), $$

$$ Z_{8}(s) = \frac{Z_{5}(s)Z_{4}(s) - Z_{6}(s)Z_{3}(s)}{Z_{7}(s)}, \quad Z_{8}(s) = \frac{Z_{6}(s)Z_{1}(s) - Z_{5}(s)Z_{2}(s)}{Z_{7}(s)}, $$

$$ \Psi''(s) = 2(Z_{1}(s)Z_{4}(s) - Z_{2}(s)Z_{3}(s)). $$

## Appendix C. The Scattering Operators in the Time Domain

The expressions for the scattering operators in the Laplace domain are given by the following:

$$ R(s) = \frac{s^2(F_5^2(s) - F_4^2(s)) + 1}{(sF_4(s) - 1)^2 - s^2F_5^2(s)}, \qquad (A8) $$

$$ T(s) = \frac{-2sF_5(s)}{(sF_4(s) - 1)^2 - s^2F_5^2(s)}. \qquad (A9) $$
These expressions can be decomposed into simple elements:

$$R(s) = -1 + \frac{1}{(1 - s(F_4(s) - F_5(s)))} + \frac{1}{(1 - s(F_4(s) + F_5(s)))}, \quad (A10)$$

$$T(s) = \frac{1}{(1 - s(F_4(s) - F_5(s)))} - \frac{1}{(1 - s(F_4(s) + F_5(s)))}. \quad (A11)$$

Using the expressions of $F_4(s)$ and $F_5(s)$ given by (42), the calculus of $(F_4(s) + F_5(s))$ as well as that of $(F_4(s) - F_5(s))$ gives the following:

$$F_4(s) + F_5(s) = 2\rho_f c_0 \left[ [1 + \phi(V_1(s) - 1)] \sqrt{\lambda_1(s)} \frac{\Psi_1''(s)}{\Psi''(s)} \coth\left(\frac{l}{2} \sqrt{\lambda_1(s)}\right) \right. \\ \left. + [1 + \phi(V_2(s) - 1)] \sqrt{\lambda_2(s)} \frac{\Psi_2''(s)}{\Psi''(s)} \coth\left(\frac{l}{2} \sqrt{\lambda_2(s)}\right) \right. \\ \left. + [1 + \phi(V_3(s) - 1)] \sqrt{\lambda_3(s)} \frac{\Psi_3''(s)}{\Psi''(s)} \coth\left(\frac{l}{2} \sqrt{\lambda_3(s)}\right) \right],$$

$$F_4(s) - F_5(s) = 2\rho_f c_0 \left[ [1 + \phi(V_1(s) - 1)] \sqrt{\lambda_1(s)} \frac{\Psi_1''(s)}{\Psi''(s)} \tanh\left(\frac{l}{2} \sqrt{\lambda_1(s)}\right) \right. \\ \left. + [1 + \phi(V_2(s) - 1)] \sqrt{\lambda_2(s)} \frac{\Psi_2''(s)}{\Psi''(s)} \tanh\left(\frac{l}{2} \sqrt{\lambda_2(s)}\right) \right. \\ \left. + [1 + \phi(V_3(s) - 1)] \sqrt{\lambda_3(s)} \frac{\Psi_3''(s)}{\Psi''(s)} \tanh\left(\frac{l}{2} \sqrt{\lambda_3(s)}\right) \right].$$

Posing

$$X(s) = 2\rho_f c_0 [1 + \phi(V_1(s) - 1)] \sqrt{\lambda_1(s)} \frac{\Psi_1''(s)}{\Psi''(s)},$$

$$Y(s) = 2\rho_f c_0 [1 + \phi(V_2(s) - 1)] \sqrt{\lambda_2(s)} \frac{\Psi_2''(s)}{\Psi''(s)},$$

$$Z(s) = 2\rho_f c_0 [1 + \phi(V_3(s) - 1)] \sqrt{\lambda_3(s)} \frac{\Psi_3''(s)}{\Psi''(s)},$$

and

$$\tanh\left(\frac{l}{2} \sqrt{\lambda_i(s)}\right) = 1 - \frac{2e^{-l\sqrt{\lambda_i(s)}}}{1 + e^{-l\sqrt{\lambda_i(s)}}}, \quad (A12)$$

$$\coth\left(\frac{l}{2} \sqrt{\lambda_i(s)}\right) = 1 + \frac{2e^{-l\sqrt{\lambda_i(s)}}}{1 - e^{-l\sqrt{\lambda_i(s)}}},$$

we obtain

$$R(s) = -1 + \frac{1}{1 - X(s) - Y(s) - Z(s) + \frac{2X(s)W_1(s)}{1 + W_1(s)} + \frac{2Y(s)W_2(s)}{1 + W_2(s)} + \frac{2Z(s)W_3(s)}{1 + W_3(s)}} \\ + \frac{1}{1 - X(s) - Y(s) - Z(s) - \frac{2X(s)W_1(s)}{1 - W_1(s)} - \frac{2Y(s)W_2(s)}{1 - W_2(s)} - \frac{2Z(s)W_3(s)}{1 - W_3(s)}}.$$
$$T(s) = \frac{1}{1 - X(s) - Y(s) - Z(s) + \frac{2X(s)W_1(s)}{1+W_1(s)} + \frac{2Y(s)W_2(s)}{1+W_2(s)} + \frac{2Z(s)W_3(s)}{1+W_3(s)}} - \frac{1}{1 - X(s) - Y(s) - Z(s) - \frac{2X(s)W_1(s)}{1-W_1(s)} - \frac{2Y(s)W_2(s)}{1-W_2(s)} - \frac{2Z(s)W_3(s)}{1-W_3(s)}}$$

where

$$W_i(s) = e^{-l\sqrt{\lambda_i(s)}} \quad i = 1, 2, 3. \qquad (A13)$$

$e^{-l\sqrt{\lambda_i(s)}}$ are assumed to be very small for large values of $s$ because of the high-frequency condition assumed in our study, which allows us to use the following development:

$$R(s) = -1 + \frac{1}{1 - X(s) - Y(s) - Z(s)} \sum_{n \ge 0} \frac{2^n}{(1 - X(s) - Y(s) - Z(s))^n} \left[ \left( \frac{X(s)W_1(s)}{1 - W_1(s)} + \frac{Y(s)W_2(s)}{1 - W_2(s)} + \frac{Z(s)W_3(s)}{1 - W_3(s)} \right)^n + (-1)^n \left( \frac{X(s)W_1(s)}{1 + W_1(s)} + \frac{Y(s)W_2(s)}{1 + W_2(s)} + \frac{Z(s)W_3(s)}{1 + W_3(s)} \right)^n \right],$$

$$T(s) = \frac{1}{1 - X(s) - Y(s) - Z(s)} \sum_{n \ge 0} \frac{2^n}{(1 - X(s) - Y(s) - Z(s))^n} \left[ - \left( \frac{X(s)W_1(s)}{1 - W_1(s)} + \frac{Y(s)W_2(s)}{1 - W_2(s)} + \frac{Z(s)W_3(s)}{1 - W_3(s)} \right)^n + (-1)^n \left( \frac{X(s)W_1(s)}{1 + W_1(s)} + \frac{Y(s)W_2(s)}{1 + W_2(s)} + \frac{Z(s)W_3(s)}{1 + W_3(s)} \right)^n \right].$$

By considering only the small powers of $W_i(s)$, we obtain the following approximations:

$$R(s) = -1 + \frac{2}{(1 - X(s) - Y(s) - Z(s))} + \frac{4X(s)W_1^2(s)}{(1 - X(s) - Y(s) - Z(s))^2(1 - W_1^2(s))} + \frac{4Y(s)W_2^2(s)}{(1 - X(s) - Y(s) - Z(s))^2(1 - W_2^2(s))} + \frac{4Z(s)W_3^2(s)}{(1 - X(s) - Y(s) - Z(s))^2(1 - W_3^2(s))}.$$

$$T(s) = - \frac{4X(s)W_1(s)}{(1 - X(s) - Y(s) - Z(s))^2(1 - W_1^2(s))} - \frac{4Y(s)W_2(s)}{(1 - X(s) - Y(s) - Z(s))^2(1 - W_2^2(s))} - \frac{4Z(s)W_3(s)}{(1 - X(s) - Y(s) - Z(s))^2(1 - W_3^2(s))}.$$

For large values of $s$, we obtain the expressions of the reflection and transmission operators in the Laplace domain:
$$
\begin{align}
R(s) = & -1 + \frac{2}{(1 - X(s) - Y(s) - Z(s))} \nonumber \\
& + \frac{4X(s)W_1^2(s)}{(1 - X(s) - Y(s) - Z(s))^2} + \frac{4Y(s)W_2^2(s)}{(1 - X(s) - Y(s) - Z(s))^2} \nonumber \\
& + \frac{4Z(s)W_3^2(s)}{(1 - X(s) - Y(s) - Z(s))^2}. \tag{A14} \\
T(s) = & -\frac{4X(s)W_1(s)}{(1 - X(s) - Y(s) - Z(s))^2} - \frac{4Y(s)W_2(s)}{(1 - X(s) - Y(s) - Z(s))^2} \nonumber \\
& - \frac{4Z(s)W_3(s)}{(1 - X(s) - Y(s) - Z(s))^2}. \nonumber
\end{align}
$$

To obtain the Laplace inverse transform for these operators, it would be easier to
explicitly define X(s), Y(s) and Z(s). Knowing that

$$
\begin{align*}
\sqrt{\lambda_1(s)} \frac{\Psi_1''(s)}{\Psi''(s)} &= \frac{\phi(P+Q) - Q + V_2(s)(\phi(P+Q) - R)}{2(Q^2 - PR)(V_1(s) - V_2(s))} \frac{1}{\sqrt{\lambda_1(s)}}, \\
\sqrt{\lambda_2(s)} \frac{\Psi_2''(s)}{\Psi''(s)} &= \frac{\phi(P+Q) - Q + V_1(s)(\phi(P+Q) - R)}{2(Q^2 - PR)(V_2(s) - V_1(s))} \frac{1}{\sqrt{\lambda_2(s)}}, \\
\sqrt{\lambda_3(s)} \frac{\Psi_3''(s)}{\Psi''(s)} &= \frac{\phi(P+Q) - R}{2(Q^2 - PR)} \frac{1}{\sqrt{\lambda_3(s)}},
\end{align*}
$$

we find

$$
\frac{X(s)}{\rho_f c_0} = \frac{[1 + \phi(V_1(s) - 1)][\phi(P + Q) - Q + V_2(s)(\phi(P + Q) - R)]}{(Q^2 - PR)(V_1(s) - V_2(s))} \frac{s}{\sqrt{\lambda_1(s)}}.
$$

By replacing $V_i(s)$ with their expressions,

$$
\frac{X(s)}{\rho_f c_0} = \frac{\left[1 + \phi(A_1 - 1) + \frac{\phi B_1}{\sqrt{s}}\right] \left[\phi(P+Q) - Q + A_2(\phi(P+Q) - R) + \frac{B_1}{\sqrt{s}}(\phi(P+Q) - R)\right]}{(Q^2 - PR)(A_1 - A_2 + \frac{B_1 - B_2}{\sqrt{s}})} \frac{s}{\sqrt{\lambda_1(s)}}. \quad (A15)
$$

Using the high-frequency approximation,

$$
\frac{s}{\sqrt{\lambda_1(s)}} = \frac{1}{\sqrt{C_1}} \left( 1 - \frac{1}{2} \frac{D_1}{C_1} \frac{1}{\sqrt{s}} \right),
$$

we obtain

$$
\frac{X(s)}{\rho_f c_0} = \frac{1}{\sqrt{C_1(Q^2 - PR)(A_1 - A_2)}} \left[ \alpha_1 \gamma_2 + \frac{1}{\sqrt{s}} \left[ \phi \gamma_2 B_1 + \alpha_1 \beta_2 - \alpha_1 \gamma_2 \left( \frac{B_1 - B_2}{A_1 - A_2} + \frac{1}{2} \frac{D_1}{C_1} \right) \right] \right].
$$

The same steps for Y(s) and Z(s) lead us to the following relations,

$$
\begin{align*}
\frac{Y(s)}{\rho_f c_0} &= \frac{1}{\sqrt{C_2(Q^2 - PR)(A_1 - A_2)}} \left[ \alpha_2 \gamma_1 + \frac{1}{\sqrt{s}} \left[ \phi \gamma_2 B_2 + \alpha_2 \beta_1 - \alpha_2 \gamma_1 \left( \frac{B_2 - B_1}{A_2 - A_1} + \frac{1}{2} \frac{D_2}{C_2} \right) \right] \right], \\
\frac{Z(s)}{\rho_f c_0} &= \frac{1}{\sqrt{C_3(Q^2 - PR)}} \left[ \alpha_3 \beta_3 + \frac{1}{\sqrt{s}} \left[ B_3 \beta_3 \phi - \frac{1}{2} \frac{D_3}{C_3} \alpha_3 \beta_3 \right] \right],
\end{align*}
$$
with the following notations:

$$
\begin{align*}
\alpha_1 &= 1 + \phi(A_1 - 1), \alpha_2 = 1 + \phi(A_2 - 1), \alpha_3 = 1 + \phi(A_3 - 1), \\
\beta_1 &= B_1(\phi(P + Q) - R), \beta_2 = B_2(\phi(P + Q) - R), \beta_3 = (\phi(Q + R) - R), \\
\gamma_1 &= \phi(P + Q) - Q + A_1(\phi(P + Q) - R), \gamma_2 = \phi(P + Q) - Q + A_2(\phi(P + Q) - R),
\end{align*}
$$

and by posing

$$
X = x_1 + \frac{y_1}{\sqrt{s}}, Y = x_2 + \frac{y_2}{\sqrt{s}}, Z = x_3 + \frac{y_3}{\sqrt{s}}, \quad (A16)
$$

where

$$
x_1 = \frac{\alpha_1 \gamma_2 \rho_f c_0}{\sqrt{C_1}(Q^2 - PR)(A_1 - A_2)}, x_1 = \frac{\alpha_2 \gamma_1 \rho_f c_0}{\sqrt{C_2}(Q^2 - PR)(A_1 - A_2)}, x_3 = \frac{\alpha_3 \beta_3 \rho_f c_0}{\sqrt{C_3}(Q^2 - PR)},
$$

$$
y_1 = \frac{\rho_f c_0}{\sqrt{C_1}(Q^2 - PR)(A_1 - A_2)} \left[ \phi \gamma_2 B_2 + \alpha_2 \beta_1 - \alpha_2 \gamma_1 \left( \frac{B_2 - B_1}{A_2 - A_1} + \frac{1}{2} \frac{D_2}{C_2} \right) \right],
$$

$$
y_2 = \frac{\rho_f c_0}{\sqrt{C_2}(Q^2 - PR)(A_1 - A_2)} \left[ \phi \gamma_2 B_2 + \alpha_2 \beta_1 - \alpha_2 \gamma_1 \left( \frac{B_2 - B_1}{A_2 - A_1} + \frac{1}{2} \frac{D_2}{C_2} \right) \right],
$$

$$
y_3 = \frac{\rho_f c_0}{\sqrt{C_3}(Q^2 - PR)} \left[ B_3 \beta_3 \phi - \frac{1}{2} \frac{D_3}{C_3} \alpha_3 \beta_3 \right].
$$

Then, by putting

$$
U = x_1 + x_2 + x_3, \quad V = y_1 + y_2 + y_3,
$$

we obtain

$$
\frac{1}{1 - X - Y - Z} = \frac{1}{1 - U - \frac{V}{\sqrt{s}}} = \frac{1}{1 - U} \left( 1 - \frac{a}{\sqrt{s} + a} \right),
$$

with

$$
a = \frac{V}{(U-1)}.
$$

This gives us

$$
\frac{X(s)}{(1 - X(s) - Y(s) - Z(s))^2} = \frac{1}{(1 - U)^2} \left[ x_1 - \frac{2ax_1}{\sqrt{s} + a} + \frac{a^2 x_1}{(\sqrt{s} + a)^2} + \frac{y_1}{\sqrt{s}} - \frac{2ay_1}{\sqrt{s}(\sqrt{s} + a)} + \frac{a^2 y_1}{\sqrt{s}(\sqrt{s} + a)^2} \right]
$$

$$
\frac{Y(s)}{(1 - X(s) - Y(s) - Z(s))^2} = \frac{1}{(1 - U)^2} \left[ x_2 - \frac{2ax_1}{\sqrt{s} + a} + \frac{a^2 x_2}{(\sqrt{s} + a)^2} + \frac{y_2}{\sqrt{s}} - \frac{2ay_2}{\sqrt{s}(\sqrt{s} + a)} + \frac{a^2 y_2}{\sqrt{s}(\sqrt{s} + a)^2} \right]
$$

$$
\frac{Z(s)}{(1 - X(s) - Y(s) - Z(s))^2} = \frac{1}{(1 - U)^2} \left[ x_3 - \frac{2ax_3}{\sqrt{s} + a} + \frac{a^2 x_3}{(\sqrt{s} + a)^2} + \frac{y_3}{\sqrt{s}} - \frac{2ay_3}{\sqrt{s}(\sqrt{s} + a)} + \frac{a^2 y_3}{\sqrt{s}(\sqrt{s} + a)^2} \right].
$$
Using the inverse Laplace transforms,

$$L^{-1}\left[\frac{1}{\sqrt{s} + a}\right] = \frac{1}{\sqrt{\pi t}} - a \exp(a^2 t) \text{Erfc}(a\sqrt{t}),$$

$$L^{-1}\left[\frac{1}{(\sqrt{s} + a)^2}\right] = -2a\sqrt{\frac{t}{\pi}} + (1 + 2a^2t) \exp(a^2t) \text{Erfc}(a\sqrt{t}),$$

$$L^{-1}\left[\frac{1}{\sqrt{s}}\right] = \frac{1}{\sqrt{\pi t}}, L^{-1}\left[\frac{1}{\sqrt{s}(\sqrt{s} + a)}\right] = \exp(a^2t) \text{Erfc}(a\sqrt{t}),$$

$$L^{-1}\left[\frac{1}{\sqrt{s}(\sqrt{s} + a)^2}\right] = 2a\sqrt{\frac{t}{\pi}} - 2at \exp(a^2t) \text{Erfc}(a\sqrt{t}),$$

we find the reflection and transmission operators in the time domain:

$$R(t) = r(t) + \tilde{R}(t),$$

where

$$r(t) = \frac{1+U}{1-U}\delta(t) + \frac{2a}{1-U}\left[-\frac{1}{\sqrt{\pi t}} + a \exp(a^2t) \text{Erfc}(a\sqrt{t})\right],$$

and

$$\tilde{R}(t) = \frac{4}{(1-U)^2} [x_1 G_1(t, 2L) + x_2 G_2(t, 2L) + x_3 G_3(t, 2L)] \\ + \frac{4}{(1-U)^2} [P_1(t) * G_1(t, 2L) + P_2(t) * G_2(t, 2L) + P_3(t) * G_3(t, 2L)]. \quad (A17)$$

$$T(t) = -\frac{4}{(1-U)^2} [x_1 G_1(t, L) + x_2 G_2(t, L) + x_3 G_3(t, L)] \\ -\frac{4}{(1-U)^2} [P_1(t) * G_1(t, L) + P_2(t) * G_2(t, L) + P_3(t) * G_3(t, L)], \quad (A18)$$

with

$$P_i(t) = \frac{y_i - 2ax}{\sqrt{\pi t}} + 2a^2(y_i - ax_i)\sqrt{\frac{t}{\pi}} \\ + [a^2x_i(3 + 2a^2t) - 2ay_i(1 + a^2t)] \exp(a^2t) \text{Erfc}(a\sqrt{t}), i = 1, 2, 3, \quad (A19)$$

and

$$G_i(t, L) = L^{-1}\left[\exp\left(-jL\sqrt{\lambda_1(s)}\right)\right], i = 1, 2, 3, j = 1, 2. \quad (A20)$$

## Appendix D. Green Functions of Longitudinal and Transverse Waves

The following relations are used for the Green functions $G_i(t, K)$ of the porous material for the fast wave ($i = 1$), the slow ($i = 2$), and the rotational ($i = 3$) waves.

$$e^{-K\sqrt{\lambda_1(s)}} = e^{-K\sqrt{C_i}} \sqrt{s^2 + b_i s\sqrt{s} + e_i s}$$

$$b_i = \frac{D_i}{C_i}, e_i = \frac{G_i}{C_i}, \Delta_i = (b_i^2 - 4e_i)$$

$$G_i(t, K) = L^{-1}(e^{-K\sqrt{\lambda_i(z)}}) = 0 \quad \text{if } 0 \le t \le K \le \sqrt{c_i} \\ = \Xi_i(t) + \Delta_i \int_0^{t-k\sqrt{c_i}} h_i(t, \zeta) d\zeta \quad \text{if } t \ge K \ge \sqrt{c_i}, \quad i = 1, 2, 3 \text{ with}$$
$$ \Xi_i(t) = \frac{b_i K \sqrt{C_i}}{4 \sqrt{\pi} (t - K \sqrt{C_i})^{2/3}} \exp \left( - \frac{b_i^2 K^2 C_i}{16(t - K \sqrt{C_i})} \right), $$

where $h_i(\tau, \zeta)$ has the following form,

$$ h_i(\tau, \zeta) = -\frac{1}{4\pi^{3/2}} \frac{1}{\sqrt{(\tau - \zeta)^2 - K^2 c_i}} \frac{1}{\zeta^{3/2}} \int_{-1}^{1} \exp\left(-\frac{X_i(\mu, \tau, \zeta)}{2}\right) \times (X_i(\mu, \tau, \zeta) - 1) \times \frac{\mu d\mu}{\sqrt{1 - \mu^2}} a $$

and where

$$ X_i(\mu, \tau, \zeta) = \left( \Delta_{i\mu} \sqrt{(\tau - \zeta)^2 - K^2 C_i} + b_i (\tau - \zeta)^2 \right)^2 / 8\zeta. $$

## References

1. Osterhoff, G.; Morgan, E.F.; Shefelbine, S.J.; Karim, L.; McNamara, L.M.; Augat, P. Bone mechanical properties and changes with osteoporosis. *Injury* **2016**, *47*, S11-S20. [CrossRef]

2. Cummings, S.R.; Browner, W.; Black, D.M.; Nevitt, M.C.; Browner, W.; Genant, H.K.; Cauley, J.; Ensrud, K.; Scott, J.; Vogt, T.M. Bone density at various sites for prediction of hip fractures. *Lancet* **1993**, *341*, 72-75. [CrossRef]

3. Hui, S.L.; Slemenda, C.W.; Johnston, C. Baseline Measurement of Bone Mass Predicts Fracture in White Women. *Ann. Intern. Med.* **1989**, *111*, 355-361. [CrossRef] [PubMed]

4. Stegman, M.; Recker, R.R.; Davies, K.; Ryan, R.; Heaney, R. Fracture risk as determined by prospective and retrospective study designs. *Osteoporos. Int.* **1992**, *2*, 290-297. [CrossRef]

5. Wang, E.; Carcione, J.M.; Cavallini, F. Generalized Thermo-poroelasticity Equations and Wave Simulation. *Surv. Geophys.* **2021**, *42*, 133-157. [CrossRef]

6. Abbas, I.; Hobiny, A. The thermomechanical response of a poroelastic medium with two thermal relaxation times. *Multidiscip. Model. Mater. Struct.* **2021**, *17*, 493-506. [CrossRef]

7. Saeed, T.; Abbas, I.; Marin, M. A GL Model on Thermo-Elastic Interaction in a Poroelastic Material Using Finite Element Method. *Symmetry* **2020**, *12*, 488. [CrossRef]

8. Saeed, T. A Study on Thermoelastic Interaction in a Poroelastic Medium with and without Energy Dissipation. *Mathematics* **2020**, *8*, 1286. [CrossRef]

9. Alzahrani, F.; Abbas, I.A. Generalized Thermoelastic Interactions in a Poroelastic Material Without Energy Dissipations. *Int. J. Thermophys.* **2020**, *41*, 95. [CrossRef]

10. Yousefian, O.; White, R.D.; Karbalaeisadegh, Y.; Banks, H.T.; Muller, M. The effect of pore size and density on ultrasonic attenuation in porous structures with mono-disperse random pore distribution: A two dimensional in-silico Study. *J. Acoust. Soc. Am.* **2018**, *144*, 709-719. [CrossRef]

11. Fry, F.J.; Barger, J.E. Acoustical properties of the human skull. *J. Acoust. Soc. Am.* **1978**, *63*, 1576-1590. [CrossRef]

12. Ashman, R.B.; Corin, J.D.; Turner, C.H. Elastic properties of cancellous bone: Measurement by an ultrasonic technique. *J. Biomech.* **1987**, *20*, 979-989. [CrossRef]

13. Ashman, R.B.; Rho, J.Y. Elastic modulus of trabecular bone material. *J. Biomech.* **1988**, *21*, 177-181. [CrossRef]

14. Hosokawa, A.; Otani, T. Ultrasonic wave propagation in bovine cancellous bone. *J. Acoust. Soc. Am.* **1997**, *101*, 558-562. [CrossRef]

15. Hosokawa, A.; Otani, T. Acoustic anisotropy in bovine cancellous bone. *J. Acoust. Soc. Am.* **1998**, *103*, 2718-2722. [CrossRef]

16. Haire, T.J.; Langton, C.M. Biot theory: A review of its application on ultrasound propagation through cancellous bone. *Bone* **1999**, *24*, 291-295. [CrossRef]

17. Cardoso, L.; Teboul, F.; Sedel, L.; Meunier, A.; Oddou, C. In vitro acoustic waves propagation in human and bovine cancellous bone. *J. Bone Miner. Res.* **2003**, *18*, 1803-1812. [CrossRef]

18. Cardoso, L.; Cowin, S.C. Fabric dependence of quasi-waves in anisotropic porous media. *J. Acoust. Soc. Am.* **2011**, *129*, 3302-3316. [CrossRef]

19. Fellah, Z.E.A.; Chapelon, J.Y.; Berger, S.; Lauriks, W.; Depollier, C. Ultrasonic wave propagation in human cancellous bone: Application of Biot theory. *J. Acoust. Soc. Am.* **2004**, *116*, 61-73. [CrossRef]

20. Sebaa, N.; Fellah, Z.E.A.; Fellah, M.; Ogam, E.; Wirgin, A.; Mitri, F.G.; Depollier, C.; Lauriks, W. Ultrasonic characterization of human cancellous bone using the Biot theory: Inverse problem. *J. Acoust. Soc. Am.* **2006**, *120*, 1816-1824. [CrossRef]

21. Marutyan, K.R.; Holland, M.R.; Miller, J.G. Anomalous negative dispersion in bone can result from the interference of fast and slow waves. *J. Acoust. Soc. Am.* **2006**, *120*, EL55-EL61. [CrossRef]

22. Hughes, E.R.; Leighton, T.G.; White, P.R.; Petley, G.W. Investigation of an anisotropic tortuosity in a Biot model of ultrasonic propagation in cancellous bone. *J. Acoust. Soc. Am.* **2006**, *121*, 568-574. [CrossRef]

23. Pakula, M.; Padilla, F.; Laugier, P.; Kaczmarek, M. Application of Biot's theory to ultrasonic characterization of human cancellous bones: Determination of structural, material, and mechanical properties. *J. Acoust. Soc. Am.* **2008**, *123*, 2415-2423. [CrossRef]
24. Anderson, C.C.; Marutyan, K.R.; Holland, M.R.; Wear, K.A.; Miller, J.G. Interference between wave modes may contribute to the apparent negative dispersion observed in cancellous bone. *J. Acoust. Soc. Am.* **2008**, *124*, 1781–1789. [CrossRef]

25. Mizuno, K.; Matsukawa, M.; Otani, T.; Laugier, P.; Padilla, F. Propagation of two longitudinal waves in human cancellous bone: An in vitro study. *J. Acoust. Soc. Am.* **2009**, *125*, 3460–3466. [CrossRef]

26. Wear, K.A. Cancellous bone analysis with modified least squares Prony's method and chirp filter: Phantom experiments and simulation. *J. Acoust. Soc. Am.* **2010**, *128*, 2191–2203. [CrossRef]

27. Nelson, A.M.; Hoffman, J.J.; Anderson, C.C.; Holland, M.R.; Nagatani, Y.; Mizuno, K.; Matsukawa, M.; Miller, J.G. Determination attenuation properties of interfering fast and slow ultrasonic waves in cancellous bone. *J. Acoust. Soc. Am.* **2011**, *130*, 2233–2239. [CrossRef]

28. McKelvie, M.L. Ultrasonic Propagation in Cancellous Bone. Ph.D. Thesis, University of Hull, Hull, UK, 1988.

29. Williams, J.L. Ultrasonic wave propagation in cancellous bone and cortical bone: Prediction of some experimental results by Biot's theory. *J. Acoust. Soc. Am.* **1992**, *91*, 1106–1112. [CrossRef]

30. Lauriks, W.; Thoen, J.; Asbroeck, I.V.; Lowet, G.; Vanderperre, G. Propagation of ultrasonic pulses through trabecular bone. *J. Phys. Paris Colloq.* **1994**, *4*, 1255–1258. [CrossRef]

31. Mckelvie, M.L.; Palmer, S.B. The interaction of ultrasound with cancellous bone. *Phys. Med. Biol.* **1991**, *10*, 1331–1340. [CrossRef] [PubMed]

32. Lakes, R.; Yoon, H.S.; Katz, J.L. Slow compressional wave propagation in wet human and bovine cortical bone. *Science* **1983**, *220*, 513–515. [CrossRef] [PubMed]

33. Johnson, D.L.; Koplik, J.; Dashen, R. Theory of dynamic permeability and tortuosity in fluid-saturated porous media. *J. Fluid Mech.* **1987**, *176*, 379–402. [CrossRef]

34. Buchanan, J.L.; Gilbert, R.P.; Ou, M.J. Transfer functions for a onedimensional fluid-poroelastic system subject to an ultrasonic pulse. *Nonlinear Anal. Real World Appl.* **2012**, *13*, 1030–1043. [CrossRef]

35. Buchanan, J.L.; Gilbert, R.P.; Ou, M.Y. Recovery of the parameters of cancellous bone by inversion of effective velocities, and transmission and reflection coefficients. *Inverse Probl.* **2011**, *27*, 125006. [CrossRef]

36. Nguyen, V.H.; Naili, S.; Sansalone, V. A closed-form solution for in vitro transient ultrasonic wave propagation in cancellous bone. *Mech. Res. Commun.* **2010**, *37*, 377–383. [CrossRef]

37. Jocker, J.; Smeulders, D. Ultrasonic measurements on poroelastic slabs: Determination of reflection and transmission coefficients and processing for Biot input parameters. *Ultrasonics* **2009**, *49*, 319–330. [CrossRef]

38. Derible, S. Acoustical measurement of the bulk characteristics of a water saturated porous plate obeying Biot's theory. *Acta Acust. Acust.* **2004**, *90*, 85–90.

39. Wu, K.; Xue, Q.; Adler, L. Reflection and transmission of elastic waves from a fluid-saturated porous solid boundary. *J. Acoust. Soc. Am.* **1990**, *87*, 2349–2358. [CrossRef]

40. Santos, J.E.; Corbero, J.M.; Ravazzoli, C.L.; Hensley, J.L. Reflection and transmission coefficients in fluid-saturated porous medium. *J. Acoust. Soc. Am.* **1992**, *91*, 1911–1923. [CrossRef]

41. Johnson, D.L.; Plona, T.J.; Kojima, H. Probing porous media with first and second sound. II. Acoustic properties of water-saturated porous media. *J. Appl. Phys.* **1994**, *76*, 115–125. [CrossRef]

42. Belhocine, F.; Derible, S.; Franklin, H. Transition term method for the analysis of the reflected and the transmitted acoustic signals from watersaturated porous plates. *J. Acoust. Soc. Am.* **2007**, *122*, 1518–1526. [CrossRef]

43. Fellah, Z.E.A.; Depollier, C. Transient wave propagation in rigid porous media: A time domain approach. *J. Acoust. Soc. Am.* **2000**, *107*, 683–688. [CrossRef]

44. Fellah, Z.E.A.; Wirgin, A.; Fellah, M.; Sebaa, N.; Depollier, C.; Lauriks, W. A time-domain model of transient acoustic wave propagation in double-layered porous media. *J. Acoust. Soc. Am.* **2005**, *118*, 661–670. [CrossRef]

45. Umnova, O.; Turo, D. Time domain formulation of the equivalent fluid model for rigid porous media. *J. Acoust. Soc. Am.* **2009**, *125*, 1860–1863. [CrossRef]

46. Szabo, T.L. Time domain wave equations for lossy media obeying a frequency power law. *J. Acoust. Soc. Am.* **1994**, *96*, 491–500. [CrossRef]

47. Caviglia, G.; Morro, A. A closed-form solution for reflection and transmission of transient waves in multilayers. *J. Acoust. Soc. Am.* **2004**, *116*, 643–654. [CrossRef]

48. Norton, G.V.; Novarini, J.C. Including dispersion and attenuation directly in the time domain for wave propagation in isotropic media. *J. Acoust. Soc. Am.* **2003**, *113*, 3024–3031. [CrossRef]

49. Fellah, Z.E.A.; Fellah, M.; Lauriks, W.; Depollier, C. Direct and inverse scattering of transient acoustic waves by a slab of rigid porous material. *J. Acoust. Soc. Am.* **2003**, *113*, 61–71. [CrossRef]

50. Moura, A. Causal analysis of transient viscoelastic wave propagation. *J. Acoust. Soc. Am.* **2006**, *119*, 751–755. [CrossRef]

51. Wilson, D.K.; Ostanhev, V.E.; Collier, S.L. Time domain equations for sound propagation in rigid-frame porous media. *J. Acoust. Soc. Am.* **2004**, *116*, 1889–1892. [CrossRef]

52. Chen, W.; Holm, S. Modified Szabo's wave equation models for lossy media obeying frequency power law. *J. Acoust. Soc. Am.* **2003**, *114*, 2570–2574. [CrossRef] [PubMed]

53. Miller, K.S.; Ross, B. *An Introduction to the Fractional Calculus and Fractional Differential Equations*; John Willey & Sons, Inc.: New York, NY, USA, 1993.
54. Oldham, K.B.; Spanier, J. *The Fractional Calculus*; Academic Press: New York, NY, USA, 1974.

55. Podlubny, I. *Fractional Differential Equations*; Serie: Mathematics in Science and Engineering no 198; Academic Press: San Diego, CA, USA, 1999.

56. Samko, S.G.; Kilbas, A.A.; Marichev, O.I. *Fractional Integrals and Derivatives Theory and Applications*; Gordon and Breach Publishers: Amsterdam, The Netherlands, 1993.

57. Caputo, M. Linear models of dissipation whose Q is almost frequency independant. *Ann. Geofis.* **1966**, *19*, 383–393.

58. Caputo, M. Linear models of dissipation whose Q is almost frequency independant, Part II. *Geophys. J. R. Astron. Soc.* **1967**, *13*, 529–539. [CrossRef]

59. Caputo, M.; Mainardi, F. Linear models of dissipation in anelastic solids. *Riv. Nuovo Cimento (Ser. II)* **1971**, *1*, 161–198. [CrossRef]

60. Bagley, R.L.; Torvik, P.J. On the fractional calculus model of viscoelastic behavior. *J. Rheol.* **1983**, *30*, 133–155. [CrossRef]

61. Matignon, D. Représentations en Variables d'état de Modèles de Guides d'ondes avec Dérivation Fractionnaire. Ph.D. Thesis, l'Université Paris XI, Orsay, France, 1994.

62. Fellah, M.; Fellah, Z.E.A.; Mitri, F.; Ogam, E.; Depollier, C. Transient ultrasound propagation in porous media using Biot theory and fractional calculus: Application to human cancellous bone. *J. Acoust. Soc. Am.* **2013**, *133*, 1867–1881. [CrossRef]

63. Fellah, M.; Fellah, Z.E.A.; Ongwen, N.; Ogam, E.; Depollier, C. Acoustics of fractal porous material and fractional calculus. *Mathematics* **2021**, *9*, 1774. [CrossRef]

64. Fellah, Z.E.A.; Fellah, M.; Roncen, R.; Ongwen, N.O.; Ogam, E.; Depollier, C. Transient Propagation of Spherical Waves in Porous Material: Application of Fractional Calculus. *Symmetry* **2022**, *14*, 233. [CrossRef]

65. Fellah, Z.E.A.; Depollier, C.; Fellah, M. Direct and inverse scattering problem in porous material having a rigid frame by fractional calculus based method. *J. Sound. Vib.* **2001**, *244*, 3659. [CrossRef]

66. Hodaei, M.; Maghoul, P.; Popplewell, N. An overview of the acoustic studies of bone-like porous materials, and the effect of transverse acoustic waves. *Int. Eng. Sci.* **2020**, *147*, 103189. [CrossRef]

67. Ogam, E.; Fellah, Z.E.A.; Ogam, G.; Ongwen, N.O.; Oduor, A.O. Investigation of long acoustic waveguides for the very low frequency characterization of monolayer and stratified air-saturated poroelastic materials. *Appl. Acoust.* **2021**, *182*, 108200. [CrossRef]