# Hybrid modal reduction for poroelastic materials

Cédric Batifol, Mohamed Ichchou, Marie-Annick Galland

► To cite this version:

Cédric Batifol, Mohamed Ichchou, Marie-Annick Galland. Hybrid modal reduction for poroelastic materials. Comptes Rendus. Mécanique, 2008, 336, pp.757-765. 10.1016/j.crme.2008.09.005. hal-00442506

HAL Id: hal-00442506

https://hal.science/hal-00442506v1

Submitted on 23 Oct 2020

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
Hybrid modal reduction for poroelastic materials

Cédric Batifol<sup>a,b</sup>, Mohamed N. Ichchou<sup>a,*</sup>, Marie-Annick Galland<sup>b</sup>

<sup>a</sup> Laboratoire de tribologie et de dynamique des systèmes, CNRS-UMR 5513, École centrale de Lyon, 36, avenue Guy-de-Collongue, 69134 Ecully cedex, France

<sup>b</sup> Laboratoire de mécanique des fluides et d'acoustique, CNRS-UMR 5509, École centrale de Lyon, 36, avenue Guy-de-Collongue, 69134 Ecully cedex, France

**Abstract**

A modal-like projection method for poroelastic materials is proposed and implemented for finite element calculations. Non-physical Dirichlet conditions are imposed at the junction interface, involving constrained fluid displacements and free solid displacements. The $(u^s, U^f)$ formulation is used. The resulting frequency-dependent eigenproblem is solved without simplification using the non-linear Arnoldi algorithm. The projection subspace is spanned by calculated dynamic modes and fluid static boundary functions. A convergence study is performed and results are compared to classical Craig and Bampton and MacNeal approaches. The hybrid basis proves to be efficient.

**Résumé**

**Réduction modale hybride appliquée aux matériaux poroélastiques.** Une méthode de réduction modale pour les matériaux poroélastiques est proposée. Cette procédure de réduction est implémentée lors de calculs élémentiels finis. Des conditions de Dirichlet non-physiques sont appliquées à l'interface. La phase fluide est ainsi encastrée, contrairement à la phase solide qui est libre. La formulation $(u^s, U^f)$ est utilisée. Le problème spectral, fréquentiellement dépendant, est résolu sans approximation par l'emploi de l'algorithme d'Arnoldi non-linéaire. Le sous-espace de projection est généré par les modes dynamiques calculés et le relèvement statique fluide. Une étude de convergence est menée, et les résultats sont comparés aux approches classiques de type Craig et Bampton et MacNeal. La base hybride apparaît efficace.

**Keywords:** Poroelastic materials; Dynamic substructuring; Modal synthesis

**Mots-clés:** Matériaux poroélastiques ; Sous-structuration dynamique ; Synthèse modale

<footnote>* Corresponding author.</footnote>

<footnote>E-mail addresses: cedric.batifol@ec-lyon.fr (C. Batifol), mohamed.ichchou@ec-lyon.fr (M.N. Ichchou), marie-annick.galland@ec-lyon.fr (M.-A. Galland).</footnote>
## Version française abrégée

La modélisation éléments finis (EF) des milieux poroélastiques reste une tâche coûteuse en temps de calcul. En effet, le modèle biphasique de Biot nécessite 6 degrés de liberté (dds) par noeud dans sa version ($\mathbf{u}^s, \mathbf{U}^f$) [1], ou 4 dds en utilisant la formulation ($\mathbf{u}^s, p$) [2]. De plus, les paramètres mécaniques nécessaires à ce modèle sont à valeurs complexes et dépendent de la fréquence. Le modèle EF de la sous-structure est réduit à l'aide d'une projection de type modale. Pour le calcul des modes dynamiques, des conditions de Dirichlet non-physiques sont appliquées à l'interface de jonction. La phase fluide est bloquée, au contraire de la phase solide qui est libre. Cet artefact permet de rapprocher les comportements résonants des phases fluide et solide. Ce raisonnement est déjà appliqué dans les problèmes classiques d'interaction fluide/structure, le fluide (généralement de l'air) étant traditionnellement moins rigide que le solide [5]. Le problème spectral est résolu à l'aide de l'algorithme itératif d'Arnoldi non-linéaire [8]. L'algorithme est initialisé par les valeurs et vecteurs propres correspondant à la limite basse fréquence du problème spectral. Au final, les modes obtenus sont réorthogonalisés puis normalisés afin de former une matrice unitaire. Cette base de projection est complétée par le relèvement statique dû aux dds fluides présents à l'interface de jonction. Le passage de la base réduite à la base physique est exprimée par l'équation (25). La méthode est validée à l'aide de deux cas test, l'un monodimensionnel et l'autre bidimensionnel. Dans le premier cas, les réponses harmoniques dûes à un déplacement et une pression imposés sont calculées. Les propriétés de convergence sont comparées aux résultats obtenus à l'aide des méthodes classiques de Craig et Bampton (C-B) [10] et de MacNeal (M-N) [11]. La pertinence de la méthode est démontrée. L'erreur relative, entre résultat direct et résultat obtenu après réduction, est calculée. Le calcul de 11 modes suffit à obtenir une erreur de l'ordre de 1%. Enfin, l'évolution de l'erreur, en fonction du retrait d'un des 11 modes de la base, est calculée. En utilisant la base hybride, les erreurs des phases fluide et solide présentent la même évolution, contrairement aux cas où les bases (C-B) et (M-N) sont utilisées. Le second cas test correspond à la réponse d'un échantillon 2D soumis à un déplacement imposé. L'étude est menée jusqu'à une fréquence de 2300 Hertz. Le comportement de la base hybride est intermédiaire à ceux des bases (C-B) et (M-N). L'erreur résiduelle reste inférieure à 2 pourcents pour une fréquence d'étude égale à la moitié de la plus grande fréquence propre calculée.

## 1. Introduction

Finite element modelling of poroelastic materials is still a complicated and time-consuming task. Indeed, the classical biphasic Biot model involves 6 degrees of freedom (dofs) for the ($\mathbf{u}^s, \mathbf{U}^f$) formulation [1] or 4 dofs for the ($\mathbf{u}^s, p$) formulation [2]. Moreover, materials' parameters are complex-valued and frequency-dependent. This prevents from using such formulations in optimisation procedures. Reduced models are needed, and some have already been proposed. Recently, Jaouen [3] developed a plate-like model and Dazel [4] extended the complex modes theory to the poroelastic case. Approximations still remain in these models and the accuracy for a wide range of materials and configurations is not obvious.

In this Note, a hybrid clamped-free (C-F) component mode synthesis (CMS) method is developed. Coupled dynamic modes are calculated for simultaneously fixed fluid dofs and free solid dofs at the junction boundary. Using such boundary conditions allows to shift closer both fluid and solid phases resonances. This assumption is widely used in more classical fluid/structure interaction problems, where the fluid (mainly air) is softer than the solid part [5]. In this contribution, eigenmodes are calculated, using the ($\mathbf{u}^s, \mathbf{U}^f$) formulation, solving the frequency-dependent eigenproblem without simplifying it. Then, the projection subspace is spanned by a dynamic basis and fluid static boundary functions. Firstly, the complete derivation of the projection basis is detailed. Then, the last part of this paper presents a convergence study where results are compared to classical Craig and Bampton (C-B) and MacNeal (M-N) approaches, showing the effectiveness of the method.

## 2. The Clamped-Free poroelastic eigenproblem

The configuration studied is depicted on Fig. 1. A poroelastic domain $\Omega$ is subjected to harmonic forcing on its Neumann boundary $\partial\Omega_N$. Fixing conditions involve linear constraints on its Dirichlet boundary $\partial\Omega_D$. Finally, this sample is coupled to other layers on its junction boundary $\partial\Omega_J$. The whole boundary $\partial\Omega$ admits the direct summation of Eq. (1).

$$ \partial\Omega = \partial\Omega_D \oplus \partial\Omega_N \oplus \partial\Omega_J \quad (1) $$

The associated eigenproblem, using the solid displacement $\mathbf{u}^s$ and the fluid displacement $\mathbf{U}^f$ as kinematic unknowns, with hybrid (C-F) conditions on $\partial\Omega_J$, states as below (Eqs. (2) and (3)), with boundary equations (4)–(6).
<figure><img src="image_2.png" /><figcaption>Fig. 1. Poroelastic domain $\Omega$ and its boundary $\partial\Omega$.</figcaption></figure>

<caption>Fig. 1. Domaine poroélastique $\Omega$ et sa frontière $\partial\Omega$.</caption>

Find $\omega$ and $(\mathbf{u}^s, \mathbf{U}^f) \neq (\mathbf{0}, \mathbf{0})$ solution of

$$ \nabla \cdot \mathbf{T}^s = -\omega^2(\tilde{\rho}^{11}\mathbf{u}^s + \tilde{\rho}^{12}\mathbf{U}^f) \qquad (2) $$

$$ \nabla \cdot \mathbf{T}^f = -\omega^2(\tilde{\rho}^{22}\mathbf{U}^f + \tilde{\rho}^{12}\mathbf{u}^s) \qquad (3) $$

$$ \mathbf{u}^s = \mathbf{u}_{\text{fix}} \quad \text{and} \quad \mathbf{U}^f \cdot \mathbf{n} = \mathbf{u}_{\text{fix}} \cdot \mathbf{n} \quad \text{on } \partial\Omega_D \qquad (4) $$

$$ \mathbf{T}^s \cdot \mathbf{n} = \mathbf{0} \quad \text{and} \quad \mathbf{T}^f \cdot \mathbf{n} = \mathbf{0} \quad \text{on } \partial\Omega_N \qquad (5) $$

$$ \mathbf{T}^s \cdot \mathbf{n} = \mathbf{0} \quad \text{and} \quad \mathbf{U}^f = \mathbf{0} \quad \text{on } \partial\Omega_J \qquad (6) $$

where $\mathbf{n}$, $\mathbf{T}^s$, $\mathbf{T}^f$ and $\mathbb{E}^s$ stands for the unit normal vector, the solid stress tensor, the fluid stress tensor and the solid Green tensor, respectively. $\mathbf{u}_{\text{fix}}$ stands for the displacement vector imposed by the fixing conditions. The stress–strain relations are described in Eqs. (7) and (8) as follows.

$$ \mathbf{T}^s = \tilde{\lambda}^s \nabla \cdot \mathbf{u}^s \mathbf{I} + 2\mu_s \mathbb{E}^s + \tilde{\lambda}^{sf} \nabla \cdot \mathbf{U}^f \mathbf{I} \qquad (7) $$

$$ \mathbf{T}^f = \tilde{\lambda}^f \nabla \cdot \mathbf{U}^f \mathbf{I} + \tilde{\lambda}^{sf} \nabla \cdot \mathbf{u}^s \mathbf{I} \qquad (8) $$

where $\mathbf{I}$ is the identity matrix. In above equations, $\tilde{\lambda}^s$ ans $\mu^s$ stand for the two Lamé-like coefficients of the solid phase. $\tilde{\lambda}^f$ and $\tilde{\lambda}^{sf}$ are the fluid phase and coupling Lamé-like coefficients. Finally, $\tilde{\rho}^{11}$, $\tilde{\rho}^{22}$ and $\tilde{\rho}^{12}$ are the solid, fluid and coupling relative densities. All these coefficients, except $\mu^s$, are complex-valued and frequency-dependent parameters. For the sake of clarity, the derivation of all these parameters is not recalled. The reader is referred to [6] for additional informations.

Let $\mathbf{C}_{\mathbf{u}^s}(\Omega)$ and $\mathbf{C}_{\mathbf{U}^f}(\Omega)$ be the admissible spaces for the solid and fluid displacements, respectively. Multiplying Eqs. (2) and (3) by admissible displacement variations $\delta\mathbf{u}^s$ and $\delta\mathbf{U}^f$, applying the Green formula, the stress–strain relations, and considering boundary conditions lead to the variational formulation of the eigenproblem.

Find $(\mathbf{u}^s, \mathbf{U}^f, \omega) \in \mathbf{C}_{\mathbf{u}^s}(\Omega) \times \mathbf{C}_{\mathbf{U}^f}(\Omega) \times \mathbb{C}$ satisfying Eqs. (9) and (10) $\forall(\delta\mathbf{u}^s, \delta\mathbf{U}^f) \in \mathbf{C}_{\mathbf{u}^s}(\Omega) \times \mathbf{C}_{\mathbf{U}^f}(\Omega)$.

$$ -\left(\tilde{\lambda}^s \int_{\Omega} \nabla \cdot \mathbf{u}^s \nabla \cdot \delta\mathbf{u}^s \, d\Omega + 2\mu_s \int_{\Omega} \mathbb{E}^s : \delta\mathbb{E}^s \, d\Omega + \tilde{\lambda}^{sf} \int_{\Omega} \nabla \cdot \mathbf{U}^f \nabla \cdot \delta\mathbf{u}^s \, d\Omega\right) + \omega^2 \tilde{\rho}^{11} \int_{\Omega} \mathbf{u}^s \cdot \delta\mathbf{u}^s \, d\Omega \\ + \omega^2 \tilde{\rho}^{12} \int_{\Omega} \mathbf{U}^f \cdot \delta\mathbf{u}^s \, d\Omega = 0 \qquad (9) $$

$$ -\left(\tilde{\lambda}^f \int_{\Omega} \nabla \cdot \mathbf{U}^f \nabla \cdot \delta\mathbf{U}^f \, d\Omega + \tilde{\lambda}^{sf} \int_{\Omega} \nabla \cdot \mathbf{u}^s \nabla \cdot \delta\mathbf{U}^f \, d\Omega\right) + \omega^2 \tilde{\rho}^{22} \int_{\Omega} \mathbf{U}^f \cdot \delta\mathbf{U}^f \, d\Omega \\ + \omega^2 \tilde{\rho}^{12} \int_{\Omega} \mathbf{u}^s \cdot \delta\mathbf{U}^f \, d\Omega = 0 \qquad (10) $$
These formulations are discretized using the finite element method. The numerical eigenproblem is then written as in Eq. (11).

$$ \begin{bmatrix} \tilde{\lambda}^s K_{\text{int},1}^s + 2\mu_s K_{\text{int},2}^s - \omega^2 \tilde{\rho}^{11} M_{\text{int},1}^s & \tilde{\lambda}^{sf} K_{\text{int},c}^s - \omega^2 \tilde{\rho}^{12} M_{\text{int},c}^s \\ \tilde{\lambda}^{sf} K_{\text{int},c}^f - \omega^2 \tilde{\rho}^{12} M_{\text{int},c}^f & \tilde{\lambda}^f K_{\text{int},1}^f - \omega^2 \tilde{\rho}^{22} M_{\text{int},1}^f \end{bmatrix} \begin{bmatrix} \mathbf{u}^s \\ \mathbf{U}^f \end{bmatrix} = \begin{bmatrix} \mathbf{0} \\ \mathbf{0} \end{bmatrix} \quad (11) $$

where $\mathbf{u}^s$ and $\mathbf{U}^f$ are vectors of solid and fluid dofs. The above dimensionless matrices are discretized forms of domain integrals. As an example, the meaning of $K_{\text{int},1}^s$ is detailed in Eq. (12).

$$ \int_{\Omega} \nabla \cdot \mathbf{u}^s \nabla \cdot \delta \mathbf{u}^s \, d\Omega \Rightarrow \delta \mathbf{u}^s T K_{\text{int},1}^s \mathbf{u}^s \quad (12) $$

where $T$ is the transpose operator.

Finally, linear constraints are applied thanks to Lagrange multipliers. An incompatibles meshes linking procedure [7] can be applied.

### 3. The non-linear Arnoldi eigensolver

The non-linear Arnoldi eigensolver is an iterative solver for sparse problems. The version used for solving the poroelastic eigenproblem is adapted from the one originally proposed by Voss [8]. It relies on three main steps. First, approximated eigenvectors and eigenvalues are calculated. It provides accurate first vectors forming a projection basis. Then the non-linear projected problem is solved using a dense solver. Finally, the error between calculated and exact eigenvector is estimated. If convergency is not reached, a new relevant vector is added to the search space.

#### 3.1. Approximated problem

For the poroelastic case, the approximated eigenpairs are calculated by solving the frequency-independent eigenproblem (13), with the constant approximated poroelastic parameters expressed in Eqs. (14), (15).

$$ \begin{bmatrix} \lambda^{s0} K_{\text{int},1}^s + 2\mu_s K_{\text{int},2}^s - \omega^2 \rho^{11} M_{\text{int},1}^s & \lambda^{sf0} K_{\text{int},c}^s - \omega^2 \rho^{12} M_{\text{int},c}^s \\ \lambda^{sf0} K_{\text{int},c}^f - \omega^2 \rho^{12} M_{\text{int},c}^f & \lambda^f0 K_{\text{int},1}^f - \omega^2 \rho^{22} M_{\text{int},1}^f \end{bmatrix} \begin{bmatrix} \mathbf{u}^s \\ \mathbf{U}^f \end{bmatrix} = \begin{bmatrix} \mathbf{0} \\ \mathbf{0} \end{bmatrix} \quad (13) $$

$$ \rho^{12} = -\phi\rho_0(\alpha_{\infty} - 1), \quad \rho^{11} = (1 - \phi)\rho_s - \rho^{12}, \quad \rho^{22} = \phi\rho_0 - \rho^{12} \quad (14) $$

$$ \lambda_0^f = \phi P_0, \quad \lambda_0^{sf} = (1 - \phi) P_0, \quad \lambda_0^s = \frac{(1 - \phi)^2}{\phi} P_0 + 2\mu_s \frac{\nu}{1 - 2\nu} \quad (15) $$

A few characteristic coefficients have to be introduced. $\rho_0$, $\rho_s$, $\phi$, $\alpha_\infty$, $P_0$, $\nu$ are the fluid density, the density of the material of the skeleton, the porosity, the tortuosity, the static pressure, and the Poisson ratio, respectively. Problem (13) involves symmetric matrices and is solved using standard sparse eigensolvers. The calculated eigenvalues and eigenvectors are denoted $\tau_j$ and $\mathbf{B}_j$.

#### 3.2. Resolution of the projected problem

The eigenpair $(\omega_j, \mathbf{V}_j)$ now needs to be searched for. Denoting $\mathcal{M}(\omega)$ the matrix involved in problem (11), a first order development is written as in Eq. (16).

$$ \mathcal{M}(\omega_j) \approx \mathcal{M}(\tau_j) - \theta \mathcal{M}'(\tau_j) \quad (16) $$

The method of successive linear problems introduced by Ruhe [9] allows to search for an eigenpair near a pole $\xi_j$. Obviously, the relevance of such a search depends on the chosen pole. For the poroelastic problem, one may choose $\xi_j = 0.9\tau_j$. If the pole is too close to the exact eigenvalue, matrices could be ill-conditioned. This solver being not directly applicable to large sparse matrices, a projection of problem (11) is to be solved. A projection basis $\mathbf{B}$ is spanned by appropriate vectors. The first of them is the approximated eigenvector $\mathbf{B}_j$. The resulting linear eigenproblem is presented in Eq. (17). The first approximation of the seek eigenvalue is denoted $r_0$ and is equal to $\tau_j$. At the lth-step, the retained value of $\theta$ is such that $\|\xi_j - (r_l - \theta)\|$ is minimal. The (l + 1)th eigenvalue approximation, $r_{l+1}$, is then equal to $r_l - \theta$. Finally, the algorithm stops when the convergence criterion (18) is reached.
$$ \mathbb{B}^H \mathbf{M}(r_l) \mathbb{B} \mathbf{Y} = \theta \mathbb{B}^H \mathbf{M}'(r_l) \mathbb{B} \mathbf{Y}, \quad \text{with } \mathbb{B}^H \mathbb{B} = \mathbf{I} \qquad (17) $$

$$ \frac{\|\xi_j - (r_{l+1} - \theta)\|}{\|\xi_j - (r_l - \theta)\|} - 1 < \epsilon \qquad (18) $$

$\mathbb{B}^H$ is the adjoint operator. $\epsilon$ is a small scalar parameter.

### 3.3. Error estimate and basis updating

The preceding step gives the estimated eigenpair $(\hat{\omega}_j, \hat{\mathbf{U}}_j)$, with $\hat{\mathbf{U}}_j = \mathbb{B}\mathbf{Y}$. The error between estimated and exact eigenpair is expressed via the calculation of the residue $\mathbf{R}$ expressed in Eq. (19).

$$ \mathbf{R} = \mathbf{M}(\hat{\omega}_j)\hat{\mathbf{U}}_j \qquad (19) $$

Convergence is reached if the norm of $\mathbf{R}$ is sufficiently smaller than the norm of $\hat{\mathbf{U}}_j$. If this is not the case, the basis $\mathbb{B}$ has to be expanded with the new vector $\mathbf{b}_n$ satisfying conditions (20).

$$ \mathbf{M}(\xi_j)\mathbf{b}_n = \mathbf{R}, \quad \mathbb{B}^H\mathbf{b}_n = 0, \quad \|\mathbf{b}_n\| = 1 \qquad (20) $$

## 4. Construction of the hybrid (C-F) reduced poroelastic model

We consider now the response of a poroelastic sample subjected to non-zero prescribed loads and/or displacements. The admissible search space for such forced problem is denoted $C_{\mathbf{u}^s}^{\text{for}} \times C_{\mathbf{U}^f}^{\text{for}}$. The hybrid (C-F) reduced poroelastic model relies on the direct decomposition (21).

$$ \left\{ \begin{matrix} C_{\mathbf{u}^s}^{\text{for}} \\ \times \\ C_{\mathbf{U}^f}^{\text{for}} \end{matrix} \right\} = \left\{ \begin{matrix} 0 \\ \times \\ C_{\mathbf{U}^f}^{\text{stat}} \end{matrix} \right\} \oplus C_{\mathbf{u}^s, \mathbf{U}^f}^{\text{dyn}} \qquad (21) $$

where $C_{\mathbf{U}^f}^{\text{stat}}$ is the space of fluid static boundary functions, and $C_{\mathbf{u}^s, \mathbf{U}^f}^{\text{dyn}}$ is the dynamic (C-F) space.

For the discrete case, this dynamic (C-F) space is spanned by a dynamic basis $\mathbb{P}^d$ calculated thanks to dynamic modes. Indeed, the dynamic eigenvectors, calculated using the non-linear Arnoldi solver, do not form a unitary matrix. Due to the frequency-dependence of the eigenproblem, eigenvectors have to be orthogonalised using the Gramm-Schmidt process. Calculated eigenvectors $\hat{\mathbf{U}}_j$ are first classified in ascending order in regards to eigenvalues. At iterate $j$, the vector $\Phi_j$ is then calculated thanks to the procedure described in Eq. (22). The dynamic basis $\mathbb{P}^d$, formed by vectors $\Phi_j$, is then a modal-like basis.

$$ \mathbf{D}_j = \hat{\mathbf{U}}_j - \mathbb{P}_{j-1}^d \mathbb{P}_{j-1}^{d\mathbf{H}} \hat{\mathbf{U}}_j, \quad \Phi_j = \frac{\mathbf{D}_j}{\sqrt{\mathbf{D}_j^H \mathbf{D}_j}} \qquad (22) $$

Next, the fluid static functions matrix, $\mathbb{P}^{\text{stat}}$, is formed by static fluid functions calculated for each junction dofs. As an example, the $l$th static fluid function is obtained by solving the discretized linear problem expressed by Eqs. (23) and (24).

$$ \begin{bmatrix} \lambda_0^s \mathbb{K}_{\text{int},1}^s + 2\mu_s \mathbb{K}_{\text{int},2}^s & \lambda_0^{sf} \mathbb{K}_{\text{int},c}^s & \lambda_0^{sf} \mathbb{K}_{\text{int},c}^f \\ \lambda_0^{sf} \mathbb{K}_{\text{int},c}^f & \lambda_0^f \mathbb{K}_{\text{int},1}^s & \lambda_0^f \mathbb{K}_{\text{int},1}^f \end{bmatrix} \begin{bmatrix} \mathbf{u}^s \\ \mathbf{U}^f \end{bmatrix} = \begin{bmatrix} \mathbf{0} \\ \mathbf{0} \end{bmatrix} \qquad (23) $$

$$ \mathbf{U}^f \cdot \mathbf{n} = 0 \quad \text{on } \partial\Omega_J \setminus l, \quad \mathbf{U}_l^f \cdot \mathbf{n} = 1 \qquad (24) $$

where $\partial\Omega_J \setminus l$ is the junction boundary without the $l$th fluid junction dof $\mathbf{U}_l^f$. Remaining Dirichlet conditions expressed in Eq. (4) have also to be fulfilled. Finally, the reduced model is written as in Eqs. (25).

$$ \begin{bmatrix} \mathbf{u}^s \\ \mathbf{U}^f \\ \mathbf{U}^f \end{bmatrix}_{|\partial\Omega_J} \approx \mathbb{M}_{C-F} \mathbf{Z} = \begin{bmatrix} \mathbf{0} \\ \mathbb{I}_{|\partial\Omega_J} \\ \mathbb{P}_{f|\partial\Omega_J}^{\text{stat}} \end{bmatrix}_{|\partial\Omega_J} \begin{bmatrix} \mathbf{U}^f \\ \mathbf{Z}_\alpha \end{bmatrix}_{|\partial\Omega_J} \qquad (25) $$

where
$Z_α$ is the vector of generalised dofs,

$P_{f|\Omega\setminus\partial\Omega_j}^{\text{stat}}$ is the fluid static functions matrix restricted to fluid dofs belonging to $\Omega \setminus \partial\Omega_j$,

$P_{f|\Omega\setminus\partial\Omega_j}^d$ is the dynamic (C-F) matrix restricted to fluid dofs belonging to $\Omega \setminus \partial\Omega_j$,

$P_s^d$ is the dynamic (C-F) matrix restricted to solid dofs.

The relevance of the projected model depends obviously on the number of dynamic modes retained.

## 5. Convergence study for a 1D sample

The test case involves a one-dimensional felt sample, 12 centimetres long. This poroelastic material is clamped at one side and is subjected to a harmonic loading or an imposed displacement on its other side. Reference results are obtained using a standard finite element formulation involving 962 dofs, and quadratic Lagrange shape functions. The hybrid (C-F) projection basis is compared to classical Craig and Bampton (C-B) [10] and MacNeal (M-N) [11] approaches. Errors on calculated fluid and solid displacements are estimated using formula (26).

$$ \text{err} = \frac{|\hat{\mathbf{u}} - \mathbf{u}|}{|\hat{\mathbf{u}}|} \qquad (26) $$

where $\hat{\mathbf{u}}$ and $\mathbf{u}$ are reference and projected solid or fluid displacements, respectively.

Results presented are calculated for a 2000 Hertz forcing frequency. Eleven dynamic modes are calculated. The projection basis is then spanned by calculated modes, their complex conjugates and relevant static functions. The number of generalised coordinates is equal to 22, 23 and 24 for (M-N), (C-F) and (C-B) projection basis, respectively. Remaining relative errors are listed on Table 1 for both phases and excitation types. The (C-F) basis proves to be the most efficient. Indeed, residual errors are small for both fluid and solid phases, which is not the case for the (C-B) and (M-N) basis.

On Fig. 2 are plotted error for all three basis, when one of the eleventh dynamic modes is not retained. Modes around the fifth one are shown to be the most relevant. Indeed, the fifth eigenfrequency is near 2000 Hertz for all

<caption>Table 1
Error estimate at 2000 Hertz for an imposed displacement (a) and an imposed pressure (b).</caption>

<table><caption>Tableau 1
Eritreur résiduelle à 2000 Hertz pour un déplacement (a) et une pression (b) imposée.</caption><thead><tr><th>u<sub>imp</sub></th><th>C-F</th><th>C-B</th><th>M-N</th><th>p<sub>imp</sub></th><th>C-F</th><th>C-B</th><th>M-N</th></tr></thead><tbody><tr><td>Solid phase err (%)</td><td>0.9</td><td>4</td><td>35</td><td>Solid phase err (%)</td><td>0.8</td><td>3.4</td><td>1.7</td></tr><tr><td>Fluid phase err (%)</td><td>1.7</td><td>0.4</td><td>23</td><td>Fluid phase err (%)</td><td>0.1</td><td>2.3</td><td>8.2</td></tr></tbody></table>

<figure><img src="image_3.png" /><figcaption>Fig. 2. Error evolution versus non-retained modes at frequency 2000 Hertz for imposed displacement (a) and pressure (b). ○ and ● stand for (C-F) solid and fluid errors. ▷ and ▶ stand for (C-B) solid and fluid errors. □ and ■ stand for (M-N) solid and fluid errors.</figcaption></figure>

<caption>Fig. 2. Évolution de l'erreur à 2000 Hertz en fonction des modes non-retenus. (a) : déplacement imposé. (b) pression imposée. ○ et ● : erreurs solide et fluide avec la base (C-F). ▷ et ▶ : erreurs solide et fluide avec la base (C-B). □ et ■ : erreurs solide et fluide avec la base (M-N).</caption>
basis. For both (C-B) and (M-N) basis, solid and fluid phases have opposite behaviours around this mode. For the (C-B) case, the fluid error continue to diminish whereas the solid one increases. The free interface basis exhibits clearly two different peaks for the two different phases. The fluid displacement converges faster than the solid one due to the use of similar interface boundary conditions for both phases. In contrast, (C-F) errors behave the same way for both phases. Each calculated eigenmode is then relevant for both phases displacements reconstruction. This is not the case for the (C-B) and (M-N) basis.

## 6. Convergence study for a 2D sample

The test case involves a 2D sample, 10 centimetres wide and 1 centimetre deep. The material is subjected to a prescribed displacement on its junction boundary $\partial\Omega_J$. In the cases of the (C-F) and (M-N) basis, the constraint is corrected by using residual attachment modes. Finally, remaining boundaries are subjected to Dirichlet conditions as written in Eq. (27). A sketch of the studied configuration can be found on Fig. 3.

$$ \mathbf{u}^s \cdot \mathbf{n} = \mathbf{U}^f \mathbf{n} = 0 \quad \text{on } \partial\Omega_{D1}, \quad \mathbf{u}^s = \mathbf{0} \quad \text{on } \partial\Omega_{D2}, \quad \mathbf{U}^f \cdot \mathbf{n} = 0 \quad \text{on } \partial\Omega_{D2} \qquad (27) $$

The forcing frequency varies from 100 to 2300 Hertz. For each frequency, the relative error between direct and reduced results is calculated following Eq. (26). Without reduction, 3600 dofs are needed to accurately model the sample behaviour. Sizes of the (C-B), (C-F) and (M-N) basis are 303 (30 eigenmodes), 131 (40 eigenmodes) and 36 (36 eigenmodes), respectively. For all three basis, only eigenpairs with eigenfrequency less or equal to 2300 are retained. In Fig. 4 are plotted errors with respect to the frequency for both fluid and solid phases. For the (M-N) basis, both fluid and solid errors are always greater than 2 percents. More eigenmodes should have been calculated. In the cases of

<figure><img src="image_4.png" /><figcaption>Fig. 3. Sketch of the 2D test case. The shape of the imposed normal displacement is plotted on the left.</figcaption></figure>

<caption>Fig. 3. Configuration étudiée pour le cas test 2D. L'allure du déplacement imposé est visible à gauche.</caption>
<figure><img src="image_5.png" /><figcaption>Fig. 4. X-displacement relative errors. ○ and ●: Solid and fluid errors, basis (C-F). ▷ and ▶: Solid and fluid errors, basis (C-B). ☐ and ■: Solid and fluid errors, basis (M-N).</figcaption></figure>

<caption>Fig. 4. ○ et Erreurs relatives pour le déplacement selon x. ● : Erreurs solide et fluide, base (C-F). ▷ et ▶ : Erreurs solide et fluide, base (C-B). ☐ et ■ : Erreurs solide et fluide, base (M-N).</caption>

(C-B) and (C-F) basis, both fluid and solid errors are less than 2 percents until the frequency 1000 Hertz is reached. One has to note that first eigenvalues for the (C-B) and (C-F) basis are equal to 1300 and 850 Hertz, respectively. The great behaviour of the (C-B) basis is then mainly due to the use of both static and fluid functions. Another remark has to be done on the slopes of fluid and solid errors. In the case of the (C-F) basis, both errors behave the same way, whereas for the (C-B) basis, the solid error increases while the fluid error stay almost constant. Calculated eigenmodes, using the (C-B) boundary conditions, are very efficient only for the fluid phase response. Finally, the (C-F) basis leads to reasonable results, for both phases, until a frequency equal to half the highest calculated eigenfrequency, with almost three times less dofs than in a Craig and Bampton approach.

## 7. Conclusion

A reduction procedure for poroelastic finite element models has been proposed. The projection subspace is spanned by hybrid clamped-free dynamic modes and static fluid functions. The poroelastic eigenproblem is solved without simplification using an iterative process. The fluid phase is fixed to the junction boundary whereas solid dofs are not constrained. The whole process of building such a basis has been detailed. The convergence properties of the (C-F) basis have been studied on both 1D and 2D test cases. For the 1D test case, the reconstructed response is better converged, using the clamped-free basis, than with Craig and Bampton or MacNeal approaches. Each calculated hybrid mode is relevant to describe both fluid and solid dynamic behaviours. This is not the case for classical substructures' basis. In the 2D test case, the (C-F) basis clearly appears as a good compromise between accuracy and reduction ratio. Both phases errors are less than 2 percents until half the highest eigenfrequency, with almost three times less dofs than in the (C-B) basis. Next steps consist in studying convergence properties for a 3D sample, formulating coupling with others substructures and deriving a simple criterion for choosing modes to be retained.

## References

[1] M.A. Biot, Theory of propagation of elastic waves in a fluid-saturated porous solid, J. Acoust. Soc. Am. 28 (1956) 168-178.

[2] N. Atalla, R. Panneton, P. Debergue, A mixed displacement-pressure formulation for poroelastic materials, J. Acoust. Soc. Am. 104 (1998) 1444-1452.

[3] L. Jaouen, B. Brouard, N. Atalla, C. Langlois, A simplified numerical model for a plate backed by a thin foam layer in the low frequency range, J. Sound Vib. 280 (2005) 681-698.

[4] O. Dazel, F. Sgard, C.-H. Lamarque, N. Atalla, An extension of complex modes for the resolution of finite-element poroelastic problems, J. Sound Vib. 253 (2002) 421-445.
[5] R. Ohayon, Reduced symmetric models for modal analysis of internal structural-acoustic and hydroelastic-shloshing problems, Comput. Meth. Appl. Mech. Engrg. 190 (2001) 3009-3019.

[6] J.F. Allard, Propagation of Sound in Porous Media. Modelling Sound Absorbing Materials, Elsevier Science, 1993.

[7] R. Ohayon, R. Sampaio, C. Soize, Dynamic substructuring of damped structures using singular value decomposition, J. Appl. Mech. 64 (2) (1997) 292-298.

[8] H. Voss, An Arnoldi method for non-linear eigenvalue problems, Tech. report, Dept math., TU Hamburd-Harburg, 2003.

[9] A. Ruhe, Rational Krylov, a practical algorithm for large sparse nonsymmetric matrix pencils, SIAM J. Sci. Comput. 19 (1998) 1535-1551.

[10] R.R. Craig, M.C.C. Bampton, Coupling of substructures for dynamic analyses, AIAA J. 6 (7) (1968) 1313-1319.

[11] R.H. MacNeal, A hybrid method of component mode synthesis, J. Computers Structures 1 (4) (1971) 581-601.