# Ondes dans les milieux poroélastiques - Analyse du modèle de Biot

Abdelaâziz Ezziani

► To cite this version:

Abdelaâziz Ezziani. Ondes dans les milieux poroélastiques - Analyse du modèle de Biot. Revue Africaine de Recherche en Informatique et Mathématiques Appliquées, 2006, Volume 5, Special Issue TAM TAM'05, november 2006, pp.95-109. 10.46298/arima.1857. hal-01263455

HAL Id: hal-01263455

https://inria.hal.science/hal-01263455v1

Submitted on 27 Jan 2016

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_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Ondes dans les milieux poroélastiques
## Analyse du modèle de Biot

Abdelaâziz Ezziani

Projet POems
INRIA Rocquencourt
B.P. 105, 78153 Le Chesnay Cedex
FRANCE
Abdelaaziz.Ezziani@inria.fr

**RÉSUMÉ.** Nous nous intéressons à la modélisation de la propagation d'ondes dans les milieux poroélastiques. Nous considérons le modèle bi-phasique de Biot. Ce papier est consacré à l'analyse mathématique de ce modèle : résultats d'existence et d'unicité, décroissance de l'énergie et le calcul d'une solution analytique.

**ABSTRACT.** We are interested in the modeling of wave propagation in poroelastic media. We consider the biphasic Biot's model. This paper is devoted to the mathematical analysis of such model : existence and uniqueness result, energy decay result and the calculation of an analytical solution.

**MOTS-CLÉS** : Ondes poroélastiques, modèle de Biot, existence et unicité, dissipation d'énergie, solution analytique.

**KEYWORDS** : Poroelastic waves, Biot's model, existence and uniqueness, energy dissipation, analytical solution.
# 1. Introduction

De nombreux sous-sols ne peuvent être considérés comme des matériaux exclusive-
ment solides. Ce sont souvent des milieux poreux c'est-à-dire constitués de solides per-
forés par une multitude de petits trous, appelés pores, occupés par un fluide. C'est no-
tamment souvent le cas des réservoirs pétroliers. Il est clair que l'analyse de résultats par
méthodes sismiques de l'exploration de tels milieux doit tenir compte du fait qu'une onde
se propageant dans un tel milieu rencontre une succession de phases solide et fluide :
on parle de milieu poro-élastique. Lorsque la distance moyenne entre les pores est petite
devant la longueur d'onde, plutôt que de considérer un tel milieu comme un milieu com-
plètement hétérogène, il est légitime de faire appel, au moins localement, à la théorie de
l'homogénéisation [10, 17], on aboutit alors au modèle de Biot [5, 6, 7] qui fait intervenir
comme inconnues non seulement le champ de déplacement dans le solide mais aussi la
pression dans le fluide. La principale caractéristique de ce modèle est qu'aux ondes P et
S dans un solide se rajoute une onde "lente", qu'on pourrait aussi qualifier de "fluide".

Ce papier est consacré à l'analyse mathématique du modèle de Biot. Les méthodes nu-
mériques pour l'approximation de ce modèle seront décrites dans un papier en prépara-
tion [2]. Ici, on s'intéresse essentiellement à la question de l'existence et l'unicité de la
solution forte dans le cas 3D et anisotrope et le calcul d'une solution de référence.

Dans la première section, nous montrerons un résultat d'existence et d'unicité de la so-
lution forte obtenue par le biais de la théorie des semi-groupes, ce résultat semble être
original et différent de celui de J. E. Santos [19] qui est basé sur un autre approche : dans
ce papier l'auteur montre l'unicité de la solution faible pour le problème isotrope à l'aide
de la méthode de Galerkin. Dans ce cadre nous citons aussi les travaux de H. Barucq et
al. [1] pour une classe 1D du problème de Biot non linéaire. Dans une deuxième courte
section, afin de montrer que le modèle de Biot conduit à un mécanisme dissipatif, nous
établirons un résultat de la décroissance d'énergie, voir théorème 4.1. Enfin, dans la der-
nière section nous déterminerons une solution analytique du problème modèle dans un
milieu infini en dimension 2. On se place dans le cas d'un milieu homogène isotrope sans
dissipation avec une source de pression ponctuelle en espace et en temps. Nous utiliserons
la méthode de Cagniard-de Hoop [11, 13] pour déterminer la fonction de Green associée
au problème modèle. Le principe de la méthode consiste à appliquer aux inconnues la
transformation de Laplace en temps et celle de Fourier suivant l'une des variables d'es-
pace (x ou y). Une fois le problème est résolu dans l'espace de Laplace-Fourier, nous
revenons en espace physique temps-espace en manipulant les contours complexes. Cette
 étude préliminaire nous sera utile surtout pour l'étude et la validation du schéma numé-
 rique développé dans [2].
## 2. Problème modèle

On considère le système de Biot du second ordre [5, 6, 7] dans $\mathbb{R}^d \times [0, T] (d = 2, 3)$ :

$$ \rho \ddot{u}_s + \rho_f \ddot{w} - \nabla \cdot \sigma = f_u \quad \text{dans } \mathbb{R}^d \times [0, T], \quad (1a) $$

$$ \rho_f \ddot{u}_s + \rho_w \ddot{w} + \mathcal{K} \dot{w} + \nabla p = f_w \quad \text{dans } \mathbb{R}^d \times [0, T], \quad (1b) $$

$$ \sigma = C\varepsilon(\mathbf{u}_s) - \beta p \mathbf{I}_d \quad \text{dans } \mathbb{R}^d \times [0, T], \quad (1c) $$

$$ \frac{1}{m} p + \beta \nabla \cdot \mathbf{u}_s + \nabla \cdot \mathbf{w} = 0 \quad \text{dans } \mathbb{R}^d \times [0, T], \quad (1d) $$

$$ \mathbf{u}_s(x, 0) = \mathbf{u}_0(x), \quad \partial_t \mathbf{u}_s(x, 0) = \mathbf{u}_1(x) \quad \text{dans } \mathbb{R}^d, \quad (1e) $$

$$ \mathbf{w}(x, 0) = \mathbf{w}_0(x), \quad \partial_t \mathbf{w}(x, 0) = \mathbf{w}_1(x) \quad \text{dans } \mathbb{R}^d. \quad (1f) $$

où $\nabla \cdot \sigma$ est l'opérateur divergence défini par : $(\nabla \cdot \sigma)_i = \sum_{j=1}^d \frac{\partial \sigma_{ij}}{\partial x_j} \quad \forall i = 1, d$.

Les inconnues du problèmes sont :

- $\mathbf{u}_s$ le champ de déplacement du solide,

- $\mathbf{w} = \phi(\mathbf{u}_f - \mathbf{u}_s)$ le déplacement du fluide par rapport au solide, $\mathbf{u}_f$ est le champ de déplacement du fluide et $\phi$ la porosité du milieu,

- $p$ la pression dans le fluide,

- $\sigma$ le tenseur des contraintes.

Les paramètres décrivant les propriétés physiques du milieu sont comme suit :

- $\rho = \rho_f\phi + (1 - \phi)\rho_s$ la densité du milieu saturé,

- $\rho_s$ la densité du solide,

- $\rho_f$ la densité du fluide,

- $\rho_w = a\rho_f/\phi$, $a$ est la tortuosité,

- $C$ un tenseur $4 \times 4$, $\mathbf{I}_d$ la matrice unité de $M_d(\mathbb{R})$,

- $\mathcal{K} = \eta/\kappa$, $\kappa$ est la perméabilité absolue et $\eta$ est la viscosité du fluide.

- $m$ et $\beta$ sont des coefficients positifs :

$$ \beta = 1 - K_0/K_s, \quad m = [\phi/K_f + (\beta - \phi)/K_s]^{-1}, $$

- $K_s$ est le module d'incompressibilité du solide,

- $K_f$ est le module d'incompressibilité du fluide,

- $K_0$ est le module d'incompressibilité sec,

et

- $\varepsilon_{ij}(\mathbf{u}) = (\partial_{x_j} u_i + \partial_{x_i} u_j)/2$ est le tenseur des déformations,
- $f_u$ et $f_w$ sont des densités de source,

- $u_0, w_0, u_1$ et $w_1$ sont des données initiales.

Nous ferons en outre les hypothèses suivantes :

1) $\rho, \rho_f, \rho_w, \mathcal{K}, m, \beta$ et $C$ sont mesurables,

2) il existe $c_-$ and $c_+$, deux constantes positives, telles que :

$$0 < c_{-} \le \rho, \rho_{f}, \rho_{w}, \mathcal{K}, m, \beta \le c_{+} < +\infty \quad \text{p.p. } x \in \mathbb{R}^{d}, \quad (2)$$

3) le tenseur $C$ satisfait :

$$C_{ijkl} = C_{jikl} = C_{klij} \quad \forall i, j, k, l = 1, d, \\ \left\{ \begin{array}{l} \exists M_{-}, M_{+} > 0 \text{ tels que :} \\ \forall \sigma \in \mathcal{L}^{\text{sym}}(\mathbb{R}^{d}), 0 < M_{-}|\sigma|^{2} \le C\sigma : \sigma \le M_{+}|\sigma|^{2} \quad \text{p.p. } x \in \mathbb{R}^{d}, \end{array} \right. \quad (3)$$

avec

$$\mathcal{L}^{\text{sym}}(\mathbb{R}^d) = \{ \sigma \in \mathcal{L}(\mathbb{R}^d) / \sigma_{ij} = \sigma_{ji} \quad \forall i, j = 1, d \},$$

$$\mathcal{L}(\mathbb{R}^d) = \{ \sigma : \mathbb{R}^d \to \mathbb{R}^d / \sigma \text{ est linéaire} \},$$

$$\sigma : \tilde{\sigma} = \sum_{ij} \sigma_{ij} \tilde{\sigma}_{ij} \text{ et } |\sigma| = (\sigma : \sigma)^{\frac{1}{2}} \quad \forall (\sigma, \tilde{\sigma}) \in [\mathcal{L}(\mathbb{R}^d)]^2,$$

$$(C \sigma)_{ij} = \sum_{kl} C_{ijkl} \sigma_{kl}, \quad \forall (C, \sigma) \in \mathcal{L}(\mathcal{L}(\mathbb{R}^d), \mathcal{L}(\mathbb{R}^d)) \times \mathcal{L}(\mathbb{R}^d).$$

**Remarque 2.1** Ici, nous avons écrit la loi (1c) dans un cadre assez général, dans le cas isotrope, le tenseur $C$ vérifie : $(C \sigma)_{ij} = \lambda_0 \sigma_{kk} \delta_{ij} + 2\mu \sigma_{ij}$, où $\mu$ est le module de cisaillement, $\lambda_0 = \lambda_f - \beta^2 m$ est le coefficient de Lamé et $\lambda_f$ est le coefficient de Lamé dans le milieu saturé.

### 3. Résultat d'existence et d'unicité

On note $M$ la matrice définie par :

$$M = \begin{bmatrix} \rho & \rho_f \\ \rho_f & \rho_w \end{bmatrix}.$$

**Remarque 3.1** La matrice $M$ est symétrique définie positive grâce aux conditions physiques [8] : $a > 1$ et $\phi < 1$.

On introduit les variables $\tilde{u}_s = \partial_t u_s$ (la vitesse dans le solide) et $\tilde{w} = \partial_t w$ (la vitesse de filtration). En posant $U = (u_s, w, \tilde{u}_s, \tilde{w}, p)^t$, le problème (1) se récrit alors sous la forme d'un système d'évolution du premier ordre :

$$\left\{ \begin{array}{l} \displaystyle \frac{dU}{dt} + \Lambda U = F, \\ U(0) = U_0, \end{array} \right. \quad (4)$$
avec

$$ \mathbf{\Lambda}U = \begin{pmatrix} -\tilde{u}_s \\ -\tilde{w} \\ -M^{-1} \left[ \begin{array}{l} \nabla \cdot (C\varepsilon(\mathbf{u}_s)) - \nabla(\beta p) \\ -K\tilde{w} - \nabla p \end{array} \right] \\ m\beta \nabla \cdot \tilde{u}_s + m \nabla \cdot \tilde{w} \end{pmatrix}, \quad (5) $$

et

$$ U_0 = (\mathbf{u}_0, \mathbf{w}_0, \mathbf{u}_1, \mathbf{w}_1, p_0)^t, \quad F = (0, 0, (\mathbf{f}_u, \mathbf{f}_w)^t M^{-1}, 0)^t, \quad p_0 = -m[\beta \nabla \cdot \mathbf{u}_0 + \nabla \cdot \mathbf{w}_0]. $$

On introduit l'espace de Hilbert $H = [H^1(\mathbb{R}^d)]^d \times [L^2(\mathbb{R}^d)]^d \times ([L^2(\mathbb{R}^d)]^d)^2 \times L^2(\mathbb{R}^d)$, et on considère l'application bilinéaire symétrique $(\cdot, \cdot)_H : H \times H \mapsto \mathbb{R}$ :

$$ \begin{aligned} (U_1, U_2)_H = & \int_{\mathbb{R}^d} \mathbf{u}_1 \cdot \mathbf{u}_2 \, dx + \int_{\mathbb{R}^d} C\varepsilon(\mathbf{u}_1) : \varepsilon(\mathbf{u}_2) \, dx + \int_{\mathbb{R}^d} \mathbf{w}_1 \cdot \mathbf{w}_2 \, dx \\ & + \int_{\mathbb{R}^d} M \begin{bmatrix} \tilde{u}_1 \\ \tilde{w}_1 \end{bmatrix} \cdot \begin{bmatrix} \tilde{u}_2 \\ \tilde{w}_2 \end{bmatrix} \, dx + \int_{\mathbb{R}^d} \frac{1}{m} p_1 p_2 \, dx, \end{aligned} \quad (6) $$

avec $U_1 = (\mathbf{u}_1, \mathbf{w}_1, \tilde{u}_1, \tilde{w}_1, p_1)^t$, $U_2 = (\mathbf{u}_2, \mathbf{w}_2, \tilde{u}_2, \tilde{w}_2, p_2)^t$. D'après les hypothèses (2)-(3), la remarque 2.1 et le lemme de Korn [15, 18] dans $[H^1(\mathbb{R}^d)]^d$, l'application (6) définie bien un produit scalaire sur $H$.

On considère l'opérateur $\mathbf{\Lambda} : D(\mathbf{\Lambda}) \subset H \mapsto H$ définie par (5), où $D(\mathbf{\Lambda})$ est donné par :

$$ D(\mathbf{\Lambda}) = \left\{ \begin{array}{ll} U = (\mathbf{u}, \mathbf{w}, \tilde{\mathbf{u}}, \tilde{\mathbf{w}}, p) \in H / & \tilde{\mathbf{u}} \in [H^1(\mathbb{R}^d)]^d, \nabla \cdot (C\varepsilon(\mathbf{u})) - \nabla(\beta p) \in [L^2(\mathbb{R}^d)]^d, \\ p \in H^1(\mathbb{R}^d), & \tilde{\mathbf{w}} \in H(\mathrm{div}, \mathbb{R}^d) \end{array} \right\}. $$

La preuve de l'existence et l'unicité de la solution forte du problème (1) se fonde sur l'utilisation de la théorie de Hille-Yosida, ce qui nécessite le lemme suivant :

**Lemme 3.1** L'opérateur $\mathbf{\Lambda}+\lambda\mathbf{I}_d$ est maximal monotone pour tout $\lambda > \max(\rho^{-1/2}, \rho_w^{-1/2})/\sqrt{2}$.

Démonstration. - **Monotonie.** Soit $U = (\mathbf{u}, \mathbf{w}, \tilde{\mathbf{u}}, \tilde{\mathbf{w}}, p)^t \in D(\mathbf{\Lambda})$, on a

$$ \begin{aligned} (\mathbf{\Lambda}U, U)_H = & -\int_{\mathbb{R}^d} \mathbf{u} \cdot \tilde{\mathbf{u}} \, dx - \int_{\mathbb{R}^d} C\varepsilon(\mathbf{u}) : \varepsilon(\tilde{\mathbf{u}}) \, dx - \int_{\mathbb{R}^d} \mathbf{w} \cdot \tilde{\mathbf{w}} \, dx \\ & - \int_{\mathbb{R}^d} \nabla \cdot (C\varepsilon(\mathbf{u})) \cdot \tilde{\mathbf{u}} \, dx - \int_{\mathbb{R}^d} \nabla(\beta p) \cdot \tilde{\mathbf{u}} \, dx + \int_{\mathbb{R}^d} K \tilde{\mathbf{w}} \cdot \tilde{\mathbf{w}} \, dx \\ & + \int_{\mathbb{R}^d} \nabla p \cdot \tilde{\mathbf{w}} \, dx + \int_{\mathbb{R}^d} \beta \nabla \cdot \tilde{\mathbf{u}} p \, dx + \int_{\mathbb{R}^d} \nabla \cdot \tilde{\mathbf{w}} p \, dx. \end{aligned} $$
Utilisant la formule de Green pour les trois intégrales :

$$
\begin{gather*}
\int_{\mathbb{R}^d} \nabla \cdot (C\varepsilon(\mathbf{u})) \cdot \tilde{\mathbf{u}} \, dx = - \int_{\mathbb{R}^d} C\varepsilon(\mathbf{u}) : \varepsilon(\tilde{\mathbf{u}}) \, dx, \\
\int_{\mathbb{R}^d} \nabla(\beta p) \cdot \tilde{\mathbf{u}} \, dx = - \int_{\mathbb{R}^d} \beta p \nabla \cdot \tilde{\mathbf{u}} \, dx, \quad \int_{\mathbb{R}^d} \nabla \cdot \tilde{\mathbf{w}} \, p \, dx = - \int_{\mathbb{R}^d} \tilde{\mathbf{w}} \cdot \nabla p \, dx,
\end{gather*}
$$

on obtient

$$
(\Lambda U, U)_H = - \int_{\mathbb{R}^d} \mathbf{u} \cdot \tilde{\mathbf{u}} \, dx - \int_{\mathbb{R}^d} \mathbf{w} \cdot \tilde{\mathbf{w}} \, dx + \int_{\mathbb{R}^d} \mathcal{K} \tilde{\mathbf{w}} \cdot \tilde{\mathbf{w}} \, dx.
$$

D'autre part :

$$
\|U\|_H^2 = \int_{\mathbb{R}^d} |\mathbf{u}|^2 dx + \int_{\mathbb{R}^d} \varepsilon(\mathbf{u}) : \varepsilon(\mathbf{u}) dx + \int_{\mathbb{R}^d} |\mathbf{w}|^2 dx + \int_{\mathbb{R}^d} \rho |\tilde{\mathbf{u}}|^2 dx \\
+ \int_d \rho_w |\tilde{\mathbf{w}}|^2 dx + 2 \int_{\mathbb{R}^d} \rho_f \tilde{\mathbf{u}} \cdot \tilde{\mathbf{w}} dx.
$$

On en déduit que :

$$
(\Lambda U, U)_H + \lambda \|U\|_H^2 \geq \int_d [\lambda |\mathbf{u}|^2 + \lambda |\mathbf{w}|^2 + \lambda \rho |\tilde{\mathbf{u}}|^2 + \lambda \rho_w |\tilde{\mathbf{w}}|^2 + 2\lambda \rho_f \tilde{\mathbf{u}} \cdot \tilde{\mathbf{w}} - \mathbf{u} \cdot \tilde{\mathbf{u}} - \mathbf{w} \cdot \tilde{\mathbf{w}}] \, dx.
$$

Ce qui montre que $\mathbf{\Lambda} + \lambda \mathbf{I}$ est monotone dès que $\lambda > \max(\rho^{-1/2}, \rho_w^{-1/2}) / \sqrt{2}$.

- Surjectivité. montrons que Λ + νΙ est surjective pour tout ν > 0. Ceci équivaut à
montrer que pour tout F = (f_u, f_w, f_û, f_ŵ, f_p) t ∈ H, il existe U = (u, w, ũ, w̃, p) t ∈
D(Λ) solution du système :

$$
\nu \mathbf{u} - \tilde{\mathbf{u}} = \mathbf{f}_u, \tag{7a}
$$

$$
\nu \mathbf{w} - \tilde{\mathbf{w}} = \mathbf{f}_w, \tag{7b}
$$

$$
\nu \begin{bmatrix} \tilde{\mathbf{u}} \\ \tilde{\mathbf{w}} \end{bmatrix} - M^{-1} \begin{bmatrix} \nabla \cdot (C\varepsilon(\mathbf{u})) - \nabla(\beta p) \\ -\mathcal{K} \tilde{\mathbf{w}} - \nabla p \end{bmatrix} = \begin{bmatrix} \mathbf{f}_{\tilde{\mathbf{u}}} \\ \mathbf{f}_{\tilde{\mathbf{w}}} \end{bmatrix}, \quad (7c)
$$

$$
m\beta \nabla \cdot \tilde{\mathbf{u}} + m \nabla \cdot \tilde{\mathbf{w}} + \nu p = f_p. \quad (7d)
$$

Si le système (7) a une solution, il est facile d'éliminer *u*, *w* et *p* et de voir que *û* et *ŵ* doivent vérifier l'équation :

$$
\begin{equation}
\begin{split}
& - \left[ \nabla \cdot (C\varepsilon(\tilde{\mathbf{u}})) + \nabla(m\beta^2 \nabla \cdot \tilde{\mathbf{u}}) + \nabla(m\beta \nabla \cdot \tilde{\mathbf{w}}) \right] + \nu^2 M \begin{bmatrix} \tilde{\mathbf{u}} \\ \tilde{\mathbf{w}} \end{bmatrix} \\
& \quad \nabla[m\beta \nabla \cdot \tilde{\mathbf{u}} + m \nabla \cdot \tilde{\mathbf{w}}] \\
& \quad + \frac{\nu}{\mathcal{K}} \begin{bmatrix} 0 \\ \tilde{\mathbf{w}} \end{bmatrix} = \nu M \begin{bmatrix} \tilde{\mathbf{f}}_{\tilde{\mathbf{u}}} \\ \tilde{\mathbf{f}}_{\tilde{\mathbf{w}}} \end{bmatrix} + \begin{bmatrix} \nabla \cdot (C\varepsilon(\mathbf{f}_u)) - \nabla(\beta f_p) \\ -\nabla f_p \end{bmatrix}
\end{split}
\tag{8}
\end{equation}
$$
La formulation variationnelle associée à (8) s'écrit :

$$
\left\{
\begin{array}{l}
\text{Trouver } \tilde{U} = (\tilde{u}, \tilde{w}) \in \tilde{H} \text{ tel que :} \\
\quad a(\tilde{U}, V) = l(V), \quad \forall V = (v, s) \in \tilde{H},
\end{array}
\right.
\qquad (9)
$$

où $\tilde{H}$ est l'espace de Hilbert : $\tilde{H} = [H^1(\mathbb{R}^d)]^d \times H(div, \mathbb{R}^d)$, muni du produit scalaire :

$$
(\tilde{U}, V)_{\tilde{H}} = \int_{\mathbb{R}^d} \tilde{u} \cdot v \, dx + \int_{\mathbb{R}^d} \nabla \tilde{u} \cdot \nabla v \, dx + \int_{\mathbb{R}^d} \tilde{w} \cdot s \, dx + \int_{\mathbb{R}^d} \nabla \cdot \tilde{w} \nabla \cdot s \, dx,
$$

$a : \tilde{H} \times \tilde{H} \mapsto \mathbb{R}$ est la forme bilinéaire :

$$
\begin{align*}
a(\tilde{U}, V) &= \int_{\mathbb{R}^d} [C\varepsilon(\tilde{u}) : \varepsilon(v) + m\beta^2 \nabla \cdot \tilde{u} \nabla \cdot v + m\beta \nabla \cdot \tilde{w} \nabla \cdot v + m\beta \nabla \cdot \tilde{u} \nabla \cdot s] \, dx \\
&\quad + \int_{\mathbb{R}^d} m \nabla \cdot \tilde{w} \nabla \cdot s \, dx + \nu^2 \int_{\mathbb{R}^d} [\rho \tilde{u} \cdot v + \rho_f \tilde{w} \cdot v \rho_f \tilde{u} \cdot s \rho_w \tilde{w} \cdot s] \, dx + \nu \int_{\mathbb{R}^d} K \tilde{w} \cdot s \, dx,
\end{align*}
$$

et $l : \tilde{H} \to \mathbb{R}$ est la forme linéaire :

$$
\begin{align*}
l(V) &= \nu \int_{\mathbb{R}^d} [\rho f_{\tilde{u}} \cdot v + \rho_f f_{\tilde{w}} \cdot v + \rho_f f_{\tilde{u}} \cdot s + \rho_w f_{\tilde{w}} \cdot s] \, dx \\
&\quad - \int_{\mathbb{R}^d} [C\varepsilon(f_u) : \varepsilon(v) - \beta f_p \nabla \cdot v - f_p \nabla \cdot s] \, dx.
\end{align*}
$$

D'après (2), (3) et l'inégalité de Korn dans $[H^1(\mathbb{R}^d)]^d$, la forme linéaire $l(\cdot)$ est continue dans $\tilde{H}$ et la forme bilinéaire $a(\cdot, \cdot)$ est continue coercive dans $\tilde{H} \times \tilde{H}$ pour tout $\nu > 0$. Le théorème de Lax-Milgram permet alors d'affirmer que le problème (9) admet une solution unique $\tilde{U} = (\tilde{u}, \tilde{w})^t$ dans l'espace $[H^1(\mathbb{R}^d)]^d \times H(div, \mathbb{R}^d)$.

Les données $(f_u, f_w, f_p)$ étant supposées appartenir à $[H^1(\mathbb{R}^d)]^d \times [L^2(\mathbb{R}^d)]^d \times [L^2(\mathbb{R}^d)]^d$, on obtient alors l'existence de $u \in [H^1(\mathbb{R}^d)]^d$ à partir de l'équation (7a), $w \in [L^2(\mathbb{R}^d)]^d$ à partir de (7b) et $p \in L^2(\mathbb{R}^d)$ à partir de (7d). Enfin en utilisant l'équation (7c), on voit facilement que $u$ et $p$ vérifient :

$$
\nabla \cdot (C\varepsilon(\mathbf{u}) - \beta p \mathbf{I}_d) \in [L^2(\mathbb{R}^d)]^d \text{ et } p \in H^1(\mathbb{R}^d).
$$

Nous avons donc montré que l'opérateur $\mathbf{\Lambda} + \nu\mathbf{I}$ était surjectif, pour tout $\nu > 0$. Pour finir la démonstration du lemme, il suffit de raisonner avec $\nu = \lambda + 1$. $\blacksquare$

Maintenant, on peut énoncer le théorème d'existence et d'unicité de la solution forte :
**Théorème 3.1** Pour toutes conditions initiales $(\mathbf{u}_0, \mathbf{w}_0, \mathbf{u}_1, \mathbf{w}_1, p_0) \in D(\mathbf{\Lambda})$ et tout $(\mathbf{f}_u, \mathbf{f}_w) \in C^1(0, T; [L^2(\mathbb{R}^d)]^d) \times C^1(0, T; [L^2(\mathbb{R}^d)]^d)$, le problème modèle (1) admet une unique solution $(\mathbf{u}_s, \mathbf{w}, p)$ qui vérifie :

$$ egin{cases} 
\mathbf{u}_s \in C^1(0, T; [H^1(\mathbb{R}^d)]^d) \cap C^2(0, T; [L^2(\mathbb{R}^d)]^d), \\
\mathbf{w} \in C^1(0, T; H(div, \mathbb{R}^d)) \cap C^2(0, T; [L^2(\mathbb{R}^d)]^d), \\
p \in C^0(0, T; H^1(\mathbb{R}^d)) \cap C^1(0, T; L^2(\mathbb{R}^d)).
\end{cases} 
\quad (10) $$

Démonstration. Sous l'hypothèse $U_0 \in D(\mathbf{\Lambda})$ et grâce au théorème de Hille-Yosida [9] et le lemme 3.1, on déduit que le problème (4) admet une unique solution $U \in C^0(0, T; D(\mathbf{\Lambda})) \cap C^1(0, T; H)$, ce qui entraîne que la solution $(\mathbf{u}_s, \mathbf{w}, p)$ du problème (1) satisfait (10). ■

# 4. Dissipation de l'énergie

Afin de prouver que le modèle décrit un milieu dissipatif, nous déterminerons une quantité décroissante qui s'appellera l'énergie du modèle :

**Définition 4.1** Soit $(\mathbf{u}_s, \mathbf{w}, p)$ la solution forte du problème (1), On appelle énergie de $(\mathbf{u}_s, \mathbf{w}, p)$ à l'instant t la quantité :

$$ E(t) = \frac{1}{2} \int_{\mathbb{R}^d} \left[ \rho |\dot{\mathbf{u}}_s|^2 + C \varepsilon(\mathbf{u}_s) : \varepsilon(\mathbf{u}_s) + \rho_w |\dot{\mathbf{w}}|^2 + \frac{1}{m} |p|^2 + 2\rho_f \dot{\mathbf{u}}_s \cdot \dot{\mathbf{w}} \right] dx. \quad (11) $$

La quantité E définit bien une énergie positive grâce à la remarque 3.1. On a le résultat de la décroissance d'énergie :

**Théorème 4.1** L'énergie $E(t)$ vérifie l'identité :

$$ \frac{dE}{dt} = - \int_{\mathbb{R}^d} K |\dot{\mathbf{w}}|^2 dx + \int_{\mathbb{R}^d} \mathbf{f}_u \cdot \dot{\mathbf{u}}_s dx + \int_{\mathbb{R}^d} \mathbf{f}_w \cdot \dot{\mathbf{w}} dx. \quad (12) $$

De plus, elle décroît en l'absence des termes sources.

Démonstration. • On applique le produit scalaire dans $\mathbb{R}^d$ à (1a) par $\dot{\mathbf{u}}_s$ et à (1b) par $\dot{\mathbf{w}}$, après une intégration sur $\mathbb{R}^d$, on obtient alors :

$$ \frac{1}{2} \frac{d}{dt} \int_{\mathbb{R}^d} \rho |\dot{\mathbf{u}}_s|^2 dx + \int_{\mathbb{R}^d} \rho_f \ddot{\mathbf{w}} \cdot \dot{\mathbf{u}}_s dx - \int_{\mathbb{R}^d} \nabla \cdot \sigma \cdot \dot{\mathbf{u}}_s dx = \int_{\mathbb{R}^d} \mathbf{f}_u \cdot \dot{\mathbf{u}}_s dx, \quad (13a) $$

$$ \int_{\mathbb{R}^d} \rho_f \ddot{\mathbf{u}}_s \cdot \dot{\mathbf{w}} dx + \frac{1}{2} \frac{d}{dt} \int_{\mathbb{R}^d} \rho |\dot{\mathbf{w}}|^2 dx + \int_{\mathbb{R}^d} K |\dot{\mathbf{w}}|^2 dx + \int_{\mathbb{R}^d} \nabla p \cdot \dot{\mathbf{w}} dx = \quad (13b) $$

$$ \int_{\mathbb{R}^d} \mathbf{f}_w \cdot \dot{\mathbf{w}} dx. $$
En intégrant par parties les derniers termes du premier membre des deux équations, on aura :

$$ \frac{1}{2} \frac{d}{dt} \int_{\mathbb{R}^d} \rho |\dot{\mathbf{u}}_s|^2 dx + \int_{\mathbb{R}^d} \rho_f \ddot{\mathbf{w}} \cdot \dot{\mathbf{u}}_s dx + \int_{\mathbb{R}^d} \sigma : \varepsilon(\dot{\mathbf{u}}_s) dx = \int_{\mathbb{R}^d} \mathbf{f}_u \cdot \dot{\mathbf{u}}_s dx, \quad (14a) $$

$$ \int_{\mathbb{R}^d} \rho_f \ddot{\mathbf{u}}_s \cdot \dot{\mathbf{w}} \, dx + \frac{1}{2} \frac{d}{dt} \int_{\mathbb{R}^d} \rho |\dot{\mathbf{w}}|^2 \, dx + \int_{\mathbb{R}^d} K |\dot{\mathbf{w}}|^2 \, dx - \int_{\mathbb{R}^d} p \nabla \cdot \dot{\mathbf{w}} \, dx = \quad (14b) $$

$$ \int_{\mathbb{R}^d} \mathbf{f}_w \cdot \dot{\mathbf{w}} \, dx. $$

En utilisant l'équation (1c), on peut exprimer $\sigma$ dans (14a) en fonction de $\mathbf{u}_s$ et $p$, ce qui nous donne :

$$ \frac{1}{2} \frac{d}{dt} \int_{\mathbb{R}^d} \rho |\dot{\mathbf{u}}_s|^2 dx + \int_{\mathbb{R}^d} \rho_f \ddot{\mathbf{w}} \cdot \dot{\mathbf{u}}_s dx + \frac{1}{2} \frac{d}{dt} \int_{\mathbb{R}^d} C \varepsilon(\mathbf{u}_s) : \varepsilon(\mathbf{u}_s) dx \\ - \int_{\mathbb{R}^d} \beta p \nabla \cdot \mathbf{u}_s dx = \int_{\mathbb{R}^d} \mathbf{f}_u \cdot \dot{\mathbf{u}}_s dx. \quad (15) $$

En dérivant l'équation (1d) par rapport au temps, on obtient : $\nabla \cdot \dot{\mathbf{w}} = -\frac{1}{m} \dot{p} - \beta \nabla \cdot \dot{\mathbf{u}}_s$, remplaçant $\nabla \cdot \dot{\mathbf{w}}$ dans (14b), on a :

$$ \int_{\mathbb{R}^d} \rho_f \ddot{\mathbf{u}}_s \cdot \dot{\mathbf{w}} \, dx + \frac{1}{2} \frac{d}{dt} \int_{\mathbb{R}^d} \rho |\dot{\mathbf{w}}|^2 \, dx + \int_{\mathbb{R}^d} K |\dot{\mathbf{w}}|^2 \, dx + \frac{1}{2} \frac{d}{dt} \int_{\mathbb{R}^d} \frac{1}{m} |p|^2 \, dx \\ + \int_{\mathbb{R}^d} \beta p \nabla \cdot \mathbf{u}_s \, dx = 0. \quad (16) $$

ce qui nous donne (12), en sommant (15) et (16).

• On en déduit la décroissance de l'énergie pour $\mathbf{f}_u = \mathbf{f}_w = 0$. ■

## 5. Solution analytique dans un milieu homogène infini

L'objectif de ce paragraphe est de calculer la solution analytique du problème modèle dans un milieu infini en dimension 2. On se place dans le cas d'un milieu homogène isotrope sans dissipation, de viscosité $\eta = 0$, avec une source de pression $f_p(x, y, t) = \delta(x)\delta(y)\delta(t)$ ponctuelle en espace et en temps. Nous utilisons la méthode de Cagniard-de Hoop [11, 13] pour déterminer la solution exacte du problème modèle. Nous renvoyons le lecteur à la thèse de J. Diaz [14] pour plus de détails sur la méthode et son application
aux ondes acoustiques et élastiques. En remplaçant $w$ par sa valeur en fonction de $u_f$ et $u_s$ et en éliminant $\sigma$, le système (1) se réécrit en fonction de $u_s$, $u_f$ et $p$ sous la forme :

$$ egin{cases} 
\rho_{11}\ddot{u}_s + \rho_{12}\ddot{u}_f - (\lambda_0 + 2\mu)\nabla(\nabla \cdot \mathbf{u}_s) + \mu \nabla \times (\nabla \times \mathbf{u}_s) + (\beta - \phi)\nabla p = 0, \\
\rho_{12}\ddot{u}_s + \rho_{22}\ddot{u}_f + \phi \nabla p = 0, \\
p + m(\beta - \phi)\nabla \cdot \mathbf{u}_s + m\phi \nabla \cdot \mathbf{u}_f = m f_p.
\end{cases} 
\quad (17) $$

où les coefficients $\rho_{11}$, $\rho_{12}$ et $\rho_{22}$ vérifient : $\rho_{11} = \rho + \phi\rho_f(a - 2)$, $\rho_{12} = \phi\rho_f(1 - a)$, $\rho_{22} = a\phi\rho_f$. Nous simplifions le système (17), en décomposant les champs de déplacement $u_s$ et $u_f$ en champs irrotationnels et isovolumiques (ondes **P** et **S**). On pose $u_s = \nabla\Phi_s + \nabla \times \Psi_s$, $u_f = \nabla\Phi_f + \nabla \times \Psi_f$. Le système (17) devient alors :

$$ A \ddot{\Phi} - B\Delta\Phi = F, \quad (18a) $$

$$ V_s^{-2} \ddot{\Psi}_s - \Delta \Psi_s = 0, \quad (18b) $$

$$ \ddot{\Psi}_f = -\rho_{22}^{-1} \rho_{12} \ddot{\Psi}_s. \quad (18c) $$

où $\Phi$ et $F$ sont deux vecteurs : $(\Phi_s, \Phi_f)^t$, $F = -m\delta(x)\delta(y)\delta(t)(\beta - \phi, \phi)^t$, $A$ et $B$

sont deux matrices symétriques définies positives : $A = \begin{pmatrix} \rho_{11} & \rho_{12} \\ \rho_{12} & \rho_{22} \end{pmatrix}$, $B = \begin{pmatrix} S & R \\ R & T \end{pmatrix}$, $V_s = (\mu\rho_{22}/(\rho_{11}\rho_{22} - \rho_{12}^2))^{1/2}$ est la vitesse des ondes **S** et

$$ S = \lambda + 2\mu, \quad \lambda = \lambda_0 + m(\beta - \phi)^2, \quad R = m\phi(\beta - \phi), \quad T = m\phi^2. $$

La matrice **B** est inversible, nous multiplions l'équation (18a) par $B^{-1}$, le système (18) se réécrit alors sous la forme :

$$ \begin{cases} B^{-1} A \ddot{\Phi} - \Delta \Phi = B^{-1} F, \\ V_s^{-2} \ddot{\Psi}_s - \Delta \Psi_s = 0, \quad \ddot{\Psi}_f = -\rho_{22}^{-1} \rho_{12} \ddot{\Psi}_s. \end{cases} \quad (19) $$

La matrice $B^{-1}A$ est diagonalisable $B^{-1}A = \mathcal{P}D\mathcal{P}^{-1}$, où $\mathcal{P}$ est la matrice de passage, $D = \text{diag}(V_{pf}^{-2}, V_{ps}^{-2})$ est la matrice diagonale semblable à $B^{-1}A$ et $V_{pf}$ et $V_{ps}$ sont respectivement la vitesse de l'onde rapide (fast wave) et de l'onde lente (slow wave) associées aux ondes **P**. En faisant le changement des variables $\Phi^* = \mathcal{P}^{-1}\Phi$ et $F^* = (BP)^{-1}F = (f_1^*, f_2^*)^t$, le système devient :

$$ \begin{cases} D \ddot{\Phi}^* - \Delta \Phi^* = F^*, \\ V_s^{-2} \ddot{\Psi}_s - \Delta \Psi_s = 0, \quad \ddot{\Psi}_f = -\rho_{22}^{-1} \rho_{12} \ddot{\Psi}_s. \end{cases} \quad (20) $$
Si on décompose les opérateurs différentiels en espace, le système (18) se réécrit alors
sous la forme équivalente :

$$
\left\{
\begin{array}{l}
V_{pf}^{-2} \ddot{\Phi}_{p1}^* - (\partial_{xx}^2 \Phi_{p1}^* + \partial_{yy}^2 \Phi_{p1}^*) = f_1^* \delta(x) \delta(y) \delta(t), \\
V_{ps}^{-2} \ddot{\Phi}_{p2}^* - (\partial_{xx}^2 \Phi_{p2}^* + \partial_{yy}^2 \Phi_{p2}^*) = f_2^* \delta(x) \delta(y) \delta(t), \\
V_s^{-2} \ddot{\Psi}_s - (\partial_{xx}^2 \Psi_s + \partial_{yy}^2 \Psi_s) = 0, \\
\ddot{\Psi}_f = -\rho_{22}^{-1} \rho_{12} \ddot{\Psi}_s, \quad (\Phi_s, \Phi_f)^t = \mathcal{P} \Phi^*, \\
\mathbf{u}_s = \nabla \Phi_s + \nabla \times \Psi_s, \quad \mathbf{u}_f = \nabla \Phi_f + \nabla \times \Psi_f.
\end{array}
\right.
\qquad (21)
$$

Appliquons la transformée de Laplace en temps et la transformée de Fourier suivant $x$ à ce système, nous obtenons :

$$
\left\{
\begin{array}{l}
(s^2/V_{pf}^2 + k^2)\hat{\Phi}_{p1}^* - \partial_{yy}^2 \hat{\Phi}_{p1}^* = f_1^* \delta(y), \\
(s^2/V_{ps}^2 + k^2)\hat{\Phi}_{p2}^* - \partial_{yy}^2 \hat{\Phi}_{p2}^* = f_2^* \delta(y), \\
(s^2/V_s^2 + k^2)\hat{\Psi}_s - \partial_{yy}^2 \hat{\Psi}_s = 0, \\
\hat{\Psi}_f = -\rho_{22}^{-1} \rho_{12} \hat{\Psi}_s, \quad (\hat{\Phi}_s, \hat{\Phi}_f)^t = \mathcal{P} \hat{\Phi}^*, \\
\hat{u}_s^x = ik\hat{\Phi}_s + \partial_y \hat{\Psi}_s, \quad \hat{u}_s^y = \partial_y \hat{\Phi}_s - ik\hat{\Psi}_s, \\
\hat{u}_f^x = ik\hat{\Phi}_f + \partial_y \hat{\Psi}_f, \quad \hat{u}_f^y = \partial_y \hat{\Phi}_f - ik\hat{\Psi}_f.
\end{array}
\right.
\qquad (22)
$$

La solution de ce système dans un milieu infini est donnée par :

$$
\left\{
\begin{array}{l}
\hat{\Phi}_{p1}^*(k, y, s) = \frac{f_1^*}{2z_{p1}} e^{-|y|z_{p1}}, \quad \hat{\Phi}_{p2}^*(k, y, s) = \frac{f_1^*}{2z_{p2}} e^{-|y|z_{p2}}, \\
\hat{\Psi}_s = 0, \quad \hat{\Psi}_f = 0, \quad (\hat{\Phi}_s, \hat{\Phi}_f)^t = \mathcal{P} \hat{\Phi}^*, \\
\hat{u}_s^x = ik\hat{\Phi}_s, \quad \hat{u}_s^y = \partial_y \hat{\Phi}_s, \quad \hat{u}_f^x = ik\hat{\Phi}_f, \quad \hat{u}_f^y = \partial_y \hat{\Phi}_f,
\end{array}
\right.
\qquad (23)
$$

avec $z_{p1}^2 = s^2/V_{pf}^2 + k^2$, $z_{p2}^2 = s^2/V_{ps}^2 + k^2$, $\mathbf{u}_s = (\mathbf{u}_s^x, \mathbf{u}_s^y)^t$ et $\mathbf{u}_f = (\mathbf{u}_f^x, \mathbf{u}_f^y)^t$.

A partir du dernier système (23), on en déduit la forme des solutions ($\operatorname{sign}(x) = x/|x|$) :

$$
\left\{
\begin{aligned}
\hat{u}_s^x &= ik \left[ \frac{P_{11} f_1^*}{2z_{p1}} e^{-|y|z_{p1}} + \frac{P_{12} f_2^*}{2z_{p2}} e^{-|y|z_{p2}} \right], \\
\hat{u}_s^y &= -\operatorname{sign}(y)/2 \left[ P_{11} f_1^* e^{-|y|z_{p1}} + P_{12} f_2^* e^{-|y|z_{p2}} \right], \\
\hat{u}_f^x &= ik \left[ \frac{P_{21} f_1^*}{2z_{p1}} e^{-|y|z_{p1}} + \frac{P_{22} f_2^*}{2z_{p2}} e^{-|y|z_{p2}} \right], \\
\hat{u}_f^y &= -\operatorname{sign}(y)/2 \left[ P_{21} f_1^* e^{-|y|z_{p1}} + P_{22} f_2^* e^{-|y|z_{p2}} \right],
\end{aligned}
\right.
\qquad (24)
$$
Pour calculer les solutions $u_s$ et $u_f$, il suffit de déterminer la transformée inverse de Laplace en $s$ et la transformée inverse de Fourier en $k$. On pose alors $\hat{u}(k, y, s) = ik e^{-|y|z}/(2z)$, $\hat{v}(k, y, s) = \text{sign}(y) e^{-|y|z}/2$, avec $z = (k^2+s^2/c^2)^{1/2}$ et $c$ une constante positive. En appliquant la transformée de Fourier inverse en $k$ à $\hat{u}$, on obtient :

$$ \hat{u}(x, y, s) = \frac{1}{4\pi} \int_{-\infty}^{+\infty} \frac{ik}{z} e^{-(|y|z+ikx)} dk. $$

En faisant le changement de variable $k = ps/c$, on aura :

$$ \hat{u}(x, y, s) = \frac{1}{4c\pi} \int_{-\infty}^{+\infty} ips(1+p^2)^{-1/2} e^{-\frac{a}{c}(|y|(1+p^2)^{1/2}+ipx)} dp = \int_{-\infty}^{+\infty} f(x, y, p, s) dp. $$

Pour retourner au domaine temporel, on cherche un chemin $\Gamma$ dans le plan complexe, vérifiant $|y|(1 + p^2)^{1/2} + ipx = ct \quad \forall p \in \Gamma$ avec $t \in \mathbb{R}^+$ et $\int_{-\infty}^{+\infty} f(x, y, p, s) dp = \int_{\Gamma} f(x, y, p, s) dp$. Ceci nous permet par unicité de la transformée de Laplace et en utilisant le changement de variable $p = \gamma(t)$ de calculer $u(x, y, t) = h(x, y, t)$ qui vérifie $\hat{u}(x, y, s) = \int_{0}^{+\infty} h(x, y, t) e^{-st} dt$.

On pose $|y|(1 + p^2)^{1/2} + ipx = ct$, $\Gamma = \{p \in \mathbb{C} / |y|(1 + p^2)^{1/2} + ipx \in \mathbb{R}^+\}$. En posant $x = r \cos \theta$ et $y = r \sin \theta$, on trouve $\Gamma = \Gamma^+ \cup \Gamma^- \cup \gamma^+ \cup \gamma^-$, avec

$$ \Gamma^{\pm} = \left\{ -i \frac{ct}{r} \cos \theta \pm |\sin \theta| (c^2 t^2 / r^2 - 1)^{1/2}, r/c < t \right\}, $$

$$ \gamma^{\pm} = \left\{ -i \left( \frac{ct}{r} \cos \theta \pm |\sin \theta| (1 - c^2 t^2 / r^2)^{1/2} \right), 0 \le t < r/c \right\}. $$

**Remarque 5.1** Puisque on s'intéresse à la propagation d'ondes dans un milieu homogène et infini, il n'aura pas des ondes de tête et par conséquent on n'utilisera pas les contours $\gamma$.

En utilisant cette méthode, on calcule :

$$ \begin{align*} \hat{u}(x, y, s) &= \int_{+\infty}^{-\infty} f(x, y, p, s) dp = \int_{\Gamma^{+} \cup \Gamma^{-}} f(x, y, p, s) dp \\ &= -\frac{1}{4\pi c} \int_{\Gamma^{+} \cup \Gamma^{-}} ips(1+p^2)^{-1/2} e^{-\frac{a}{c}(|y|(1+p^2)^{1/2}+ipx)} dp \\ &= \frac{\cos \theta}{2\pi r} \int_{r/c}^{+\infty} t(t^2 - r^2/c^2)^{-1/2} s e^{-st} dt = \int_{r/c}^{+\infty} \partial_t h_c(x, y, t) e^{-st} dt, \end{align*} $$

avec $h_c(x, y, t) = \frac{xt}{2\pi r^2 \sqrt{t^2 - r^2/c^2}}$. En utilisant l'unicité de la transformée de Laplace, on a alors $u(x, y, t) = H(t - r/c) \partial_t h_c(x, y, t)$, où H est la fonction de Heaviside.
On utilise les mêmes démarches pour $\hat{v}(k, y, s) = sign(y) e^{-|y|/2}$, on obtient dans [16], $v(x, y, t) = H(t - r/c) \partial_t g_c(x, y, t)$, $g_c(x, y, t) = h_c(y, x, t)$.

Finalement, nous obtenons :

$$
\left\{
\begin{aligned}
\boldsymbol{u}_s^x(x, y, t) &= \mathcal{P}_{11}^* H(t - r/c_1) \partial_t h_{c_1}(x, y, t) + \mathcal{P}_{12}^* H(t - r/c_2) \partial_t h_{c_2}(x, y, t), \\
\boldsymbol{u}_s^y(x, y, t) &= \mathcal{P}_{11}^* H(t - r/c_1) \partial_t g_{c_1}(x, y, t) + \mathcal{P}_{12}^* H(t - r/c_2) \partial_t g_{c_2}(x, y, t), \\
\boldsymbol{u}_f^x(x, y, t) &= \mathcal{P}_{21}^* H(t - r/c_1) \partial_t h_{c_1}(x, y, t) + \mathcal{P}_{22}^* H(t - r/c_2) \partial_t h_{c_2}(x, y, t), \\
\boldsymbol{u}_f^y(x, y, t) &= \mathcal{P}_{21}^* H(t - r/c_1) \partial_t g_{c_1}(x, y, t) + \mathcal{P}_{22}^* H(t - r/c_2) \partial_t g_{c_2}(x, y, t), \\
c_1 &= V_{pf}, \quad c_2 = V_{ps}, \quad \mathcal{P}_{ij}^* = \mathcal{P}_{ij} f_j^* \quad \forall i, j = 1, 2.
\end{aligned}
\right.
$$

Dans le cas d'une source quelconque en temps et ponctuelle en espace $\delta(x)\delta(y) f(t)$, la solution de ce problème est obtenue à l'aide de la convolution de la fonction source avec la fonction de Green :

$$ \boldsymbol{u}_s(x, y, t) = \int_0^t \boldsymbol{u}_s^\delta(x, y, \tau) f(t - \tau) \, d\tau, \quad (25) $$

où $\boldsymbol{u}_s^\delta$ est le champ de déplacement dans le solide associé à une source ponctuelle en temps. Pour ne pas alourdir et simplifier cette étude, on n'a traité que le cas d'une source de pression. On peut généraliser cette démarche à tout type de source (voir les exemples qui vont suivre).

On termine cette section par des exemples de propagation avec différents type de source. On se place dans un milieu infini occupé par un matériau poroélastique homogène, iso-trope, non dissipatif ($K=0$) et caractérisé par les données physiques : $\rho = 1.8 \text{ kg/m}^3$, $\rho_f = 1 \text{ kg/m}^3$, $\rho_w = 7.5 \text{ kg/m}^3$, $\mu = 4 \text{ Pa}$, $\lambda_0 = 5.93 \text{ Pa}$, $m = 10 \text{ Pa}$, $\beta = 0.295$, ce qui correspond aux vitesses $V_{pf} = 2.93 \text{ m/s}$, $V_{ps} = 1.19 \text{ m/s}$, $V_s = 1.95 \text{ m/s}$. Dans un premier temps, on considère une source de pression $f_p(x, y, t) = \delta(x)\delta(y) f(t)$ ponctuelle en espace, localisée au point $(0, 0)$, le signal temporel est une gaussienne en temps : $f(t) = \exp(-\pi^2 f_0^2 (t - t_0)^2)$, $t_0 = 1/f_0$, et de fréquence $f_0 = 5.95 \text{ Hz}$. La solution analytique est calculée à l'aide d'un code MATLAB en utilisant la convolution en temps de la source et la fonction de Green à l'aide de l'équation (25). Nous présentons sur la figure 1 la restriction de la première composante de la vitesse du solide sur le carré $[-3, 3] \times [-3, 3]$ dans des différents instants, $t = 0.6s, 0.8s, 1s$ et $1.2s$. Nous observons bien que la source de pression a généré deux types d'ondes $\mathbf{P}$ : une onde lente de vitesse $V_{ps}$ et une onde rapide de vitesse $V_{pf}$. Dans un deuxième temps, nous considérons le cas d'une source de cisaillement $f_u(x, y, t) = f(t) \nabla \times (\delta(x)\delta(y) (1, 1)^t)$, nous suivons les mêmes démarches que dans le cas d'une source de pression et nous montrons que la solution est une onde de cisaillement de vitesse $V_s$, voir la figure 2. Enfin, dans le dernier test, nous considérons une source de compression $f_u(x, y, t) = f(t) \nabla (\delta(x)\delta(y))$, nous présentons sur la figure 3 les instantanés de la norme de la vitesse dans le solide et sur la figure 4 la coupe de la norme de la vitesse sur la ligne y=0.
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Instantanés de la première composante de la vitesse (source de pression)</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>**Figure 2.** Instantanés de la première composante de la vitesse (source de cisaillement)</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>**Figure 3.** Instantanés de la norme de la vitesse (source de compression)</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Figure 4.** Coupe de la norme de la vitesse sur la ligne $y = 0$</figcaption></figure>

A partir de ces tests, nous remarquons que chaque type de source, de pression ou de cisaillement, ne génère que des ondes de même type, nous en déduisons alors que les ondes de compression et de cisaillement dans les milieux poroélastiques homogènes sont découpées.
## 6. Bibliographie

[1] H. BARUCQ AND M. MADAUNE-TORT AND P. SAINT-MACARY, Some existence-uniqueness results for a class of one-dimensional nonlinear Biot models. *Nonlinear Analysis*, 61 :591-612, 2005.

[2] E. BÉCACHE AND A. EZZIANI AND P. JOLY, A higher order mixed finite element method for wave propagation in poroelastic time-domain. *en préparation*.

[3] J. G. BERRYMAN AND S. R. PRIDE, Linear dynamics of double-porosity dual-permeability materials I. Governing equations and acoustic attenuation. *Physical Review E* 68, 2003.

[4] J. G. BERRYMAN AND H. F. WANG, Elastic wave propagation and attenuation in a double-porosity dual-permeability medium. *Int. J. Rock Mech. Min. Sci.*, 37 :63-78, 2000.

[5] M. A. BIOT, Theory of propagation of elastic waves in a fluid-saturated porous solid I. Low-frequency range. *J. Acoust. Soc. Am.*, 28 :168-178, 1956.

[6] M. A. BIOT, Theory of propagation of elastic waves in a fluid-saturated porous solid II. Higher frequency range. *J. Acoust. Soc. Am.*, 28 :179-191, 1956.

[7] M. A. BIOT, Mechanics of deformation and acoustic propagation in porous media. *J. Appl. Phys.*, 33 :1482-1498, 1962.

[8] T. BOURBIÉ AND O. COUSSY AND B. ZINSZNER, *Acoustique des milieux poreux*. Publications de l'Institut Français du Pétrole, 1986.

[9] H. BREZIS, *Analyse fonctionnelle-Théorie et applications*. Masson, 1983.

[10] R. BURRIDGE AND J. B. KELLER, Biot's Poroelasticity Equations by Homogenization. *In Macroscopic Properties of Disordered Media*, pages 51-57, Springer, 1982.

[11] L. CAGNIARD, *Réflexion et Réfraction des Ondes Sismiques Progressives*. Gauthier-Villard, 1939.

[12] J. M. CARCIONE, *Wave Fields in Real Media: Wave propagation in Anisotropic, Anelastic and Porous Media A*. Pergamon, 2001.

[13] A. T. DE HOOP, The surface line source problem. *Appl. Sci. Res. B.*, 8 :349-356, 1959.

[14] J. DIAZ, Approches analytiques et numériques de problèmes de transmission en propagation d'ondes en régime transitoire. Application au couplage fluide-structure et aux méthodes de couches parfaitement adaptées. PhD thesis, Univresité Paris 6, 2005.

[15] G. DUVAUT AND J. L. LIONS, *Inéquations en mécanique et en physique*. Dunod, 1972.

[16] A. EZZIANI, Modélisation mathématique et numérique de la propagation d'ondes dans les milieux viscoélastiques et poroélastiques. PhD thesis, Univresité Paris 9, 2005.

[17] U. HORNUNG, Homogenization and porous media, *volume 6 of Interdisciplinary Applied Mathematics*, Springer, 1997.

[18] J. A. NITSCHÉ, On Korn's second inequality. *RAIRO Numerical Analysis*, 15 :237-248, 1981.

[19] J. E. SANTOS, Elastic wave propagation in fluid-saturated porous medium. I. The existence and uniqueness theorems. *M2AN*, 20 :113-128, 1986.