# Contribution to the modeling and the mechanical characterization of the subsoil in the LSBB environment

Hamed Bouare, Arnaud Mesgouez, Gaëlle Lefeuve-Mesgouez

► To cite this version:

Hamed Bouare, Arnaud Mesgouez, Gaëlle Lefeuve-Mesgouez. Contribution to the modeling and the mechanical characterization of the subsoil in the LSBB environment. E3S Web of Conferences, 2019, 88, pp.06001. 10.1051/e3sconf/20198806001. hal-02617618

HAL Id: hal-02617618

https://hal.inrae.fr/hal-02617618

Submitted on 25 May 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 Attribution 4.0 International License
# Contribution to the modeling and the mechanical characterization of the subsoil in the LSBB environment

Hamed Bouare¹,*, Arnaud Mesgouez²,**, and Gaëlle Lefeuve-Mesgouez³,***

¹²³ UMR – EMMAH, University d'Avignon et des Pays de Vaucluse – INRA
Campus Jean Henri Fabre, Agroparc, 301 rue Baruch de S pinoza, BP 21239, 84916 Avignon cedex 9–
France.

**Abstract.** The present research work aims at better characterizing the specific underground environment of the LSBB (Low Noise Inter-Disciplinary Underground Science and Technology, Rustrel, France) using mechanical wave propagation information. The LSBB experimental environment is characterized by a system of cylindrical galleries, some of them presenting a concrete layer. In the global project, three steps are considered : firstly the construction of an efficient forward mechanical wave propagation model to calculate the displacement vector and stress tensor components; secondly a sensitivity analysis to extract the pertinent parameters in the configurations and models under study (viscoelastic or poroviscoelastic media with potential anisotropy); and lastly an inversion strategy to recover some of the pertinent parameters. In this proposal, the first step, under progress, is described. The work carried out is in the continuity of the work presented by Yi et al. (2016) [1] who studied the harmonic response of a cylindrical elastic tunnel, impacted by a plane compressional wave, embedded in an infinite elastic ground. The interface between the rock mass and the linen is an imperfect contact modeled with two spring parameters, Achenbach and Zhu (1989) [2]. We choose a semi-analytical approach to calculate the two-dimensional displacement and stress fields in order to get a fast tool, from the numerical point of view. The main steps of the theoretical approach are : use of the Helmholtz decomposition, solving the wave equation based on the separation method and the expansion in Bessel function series in the harmonic domain. The harmonic results are validated by comparison with Yi et al. (2016) [1] and new ones are presented. Moreover, the transient regime case obtained with the use of a Fourier transform on the time variable, is under progress.

**Keywords:** Dynamic response ; wave propagation ; imperfect interface ; circular tunnel ; semi-analytical approach

## Introduction

In geophysical approaches, it is usual to characterize the mechanical behavior of the ground by taking into account only a part of the signal (for instance the arrival times of the waves propagating in the medium). In this research project, the aim is to consider the full waveform

<footnote>*e-mail: hamed.bouare@univ-avignon.fr</footnote>

<footnote>**e-mail: arnaud.mesgouez@univ-avignon.fr</footnote>

<footnote>***e-mail: gaëlle.mesgouez@univ-avignon.fr</footnote>
information to get a better characterization of the medium. Three steps are planed: the construction of a forward model; a sensitivity analysis to identify the parameters of interest in the configuration under study; an inversion strategy to recover some mechanical parameters. The soil under study is located at the LSBB (Low Noise Inter-Disciplinary Underground Science and Technology) in Rustrel, Vaucluse, France, more specifically the area between the anti-blast gallery and the main gallery. In this proposal, the first step is presented: a forward model, using a semi-analytical approach based on the work of Yi et al. 2016 [1], is constructed. In the full waveform method we consider both the body and the surface waves. Even though the area under study is a karstic massif which is anisotropic and heterogenous, in the first approach, both the rock and the tunnel are considered to be isotropic and elastic, but with different mechanical properties. The semi-analytical method has been chosen because of its efficiency in computational time, feature that will be then needed for the sensitivity analysis and inversion steps. A bibliographical survey has been carried out with a special attention to extract articles focusing on configurations close to the LSBB-kind geometry (cylindrical galleries with a concrete layer). Fang and Jin (2016) [3] study the dynamic response of a non-circular tunnel subjected to P and SV waves. They use a visco-elastic interface to model the contact between the rock mass and the tunnel. They conclude that the viscosity has a small influence on the dynamic stress in high frequency areas compared to low frequency areas. Fang et al. (2018) [4] analyze the dynamic stress distribution around a circular tunnel subjected to SH waves. Using an elastic slip interface to represent the contact between the alluvial valley and the lined tunnel, they find that the dynamic stress increases when the elastic modulus of the tunnel lining increases, and that the interface effect decreases when larger dimension tunnels are employed. Zhao et al. (2018) [5] investigate the interaction between strengthening and isolation layers for tunnels in rock subjected to SH waves. Their results show that the strengthening layer can stabilize the rigid displacement and further reduce the stress in the lining. In this proposal, we have developed a semi-analytical forward model to calculate the dimensionless hoop stress (DSCF) but also the radial and orthoradial displacement fields both in the rock environment and in the lining, taking into account potentially an imperfect boundary condition between the rock and the tunnel. For this purpose, the paper of Yi et al. (2016) [1] is a good starting point. They study the effect of imperfect interface on the dynamic response of a circular lined tunnel impacted by plane P-waves. An imperfect interface, taking into account a two-linear-spring model, is used to model the contact between the rock mass and the tunnel.

# 1 Configuration under study and theoretical approach

An infinitely long cylindrical tunnel embedded in an infinite rock mass is considered, as illustrated in figure 1. The outer radius of the tunnel is b and its inner radius is a. A 2D approach is modeled firstly. We consider a harmonically time-varying incident plane compressional wave that propagates parallel to the x-axis and perpendicular to the tunnel cross-section (Figure 1).

We consider the rock mass and the tunnel as homogeneous, isotropic and elastic media with different mechanical properties: Lamé constants $\lambda$ and $\mu$, and density $\rho$. Taking into account these hypotheses, the displacement equations of motion in each medium can be written as follows

$$ \mu_k \nabla^2 \mathbf{u}_k + (\lambda_k + \mu_k) \nabla \nabla \cdot \mathbf{u}_k = \rho_k \ddot{\mathbf{u}}_k, \quad k = 1, 2 \quad (1) $$

where the displacement vector $\mathbf{u}$ with subscript k represents motion in the rock mass ($k = 1$) and in the lined tunnel ($k = 2$).
<figure><img src="image_2.png" /><figcaption>Figure 1. Problem geometry</figcaption></figure>

ü is the second derivative with respect to time.

The displacement field may be separated into scalar $\phi$ and vector $\Psi$ potentials with the usual Helmholtz decomposition, equation (2). It does not depend on the z-coordinate due to the two dimensional configuration. Consequently potential vector $\Psi$ has only one single component $\psi_z$ and for the sake of clarity it will be written as $\psi$. It comes

$$ \mathbf{u}_k = \nabla\phi_k + \nabla \times \Psi_k \quad (2) $$

Substituting equation (2) into (1), two wave equations are obtained

$$ \nabla^2 \phi_k - \frac{1}{c_{pk}^2} \ddot{\phi}_k = 0 \quad (3) $$

$$ \nabla^2 \Psi_k - \frac{1}{c_{sk}^2} \ddot{\Psi}_k = 0 \quad (4) $$

Equations (3) and (4) are relative to P and SV waves, respectively.

$c_{pk} = \sqrt{\frac{\lambda_k + 2\mu_k}{\rho_k}}$ and $c_{sk} = \sqrt{\frac{\mu_k}{\rho_k}}$ are the P and SV wave velocities in the elastic medium under consideration.

Following the usual methodology proposed by Pao and Mow (1973) [6], Wirgin (2015) [7], Yi et al. (2016) [1] among others, in the context of a harmonic regime, the incident scalar potential can be written in polar coordinates as follows

$$ \phi_i = \varphi_0 e^{i(\alpha_1 r \cos \theta - \omega t)} \quad (5) $$

where $\varphi_0$ is the amplitude of scalar potential of incident wave, $i^2 = -1$ is the imaginary complex number, $\omega = 2\pi f$ is the circular frequency of the incident wave, $f$ is the frequency of the incident wave and $t$ is the time variable.

$\alpha_1$ is the incident P wavenumber.

$(r, \theta)$ are the parameters in polar coordinates.

Equation (5) can be expressed in terms of Fourier Bessel series as

$$ \phi_i = \varphi_0 \sum_{n=0}^{\infty} \epsilon_n i^n J_n(\alpha_1 r) \cos(n\theta) e^{-i\omega t} \quad (6) $$
$J_n$ is the Bessel function of the first kind of an integer n, $\epsilon_n = 1$ for n = 0 and $\epsilon_n = 2$ for n > 0. The presence of the interface between the rock mass and the lining results in two reflected waves (P and SV) in the rock mass when the incident P wave meets the interface. These are expressed with Hankel functions of the first kind, which represent outgoing cylindrical waves. The total displacement potentials in the rock mass are then given by

$$ \phi_{1t} = \sum_{n=0}^{\infty} [\varphi_0 \epsilon_n i^n J_n(\alpha_1 r) + E_n H_n^{(1)}(\alpha_1 r)] \cos(n\theta) e^{-i\omega t} \quad (7) $$

$$ \psi_{1t} = \sum_{n=0}^{\infty} F_n H_n^{(1)}(\beta_1 r) \sin(n\theta) e^{-i\omega t} \quad (8) $$

where $E_n$ and $F_n$ are respectively the expansion coefficients of the two reflected waves (P and SV).

$H_n^{(1)}$ is the Hankel function of the first kind of an integer order n and $\beta_1$ is the SV wavenumber in the rock mass.

In the tunnel medium, in contact with the rock mass at the outer circular surface ($r = b$), two incoming cylindrical waves (P and SV) and two outgoing cylindrical waves propagate. The two incoming waves can be modeled by Hankel functions of the second kind ($H_n^{(2)}$). The two outgoing waves can be written in terms of Hankel functions of the first kind. The total displacement potentials in the tunnel are

$$ \phi_{2t} = \sum_{n=0}^{\infty} [L_n H_n^{(2)}(\alpha_2 r) + P_n H_n^{(1)}(\alpha_2 r)] \cos(n\theta) e^{-i\omega t} \quad (9) $$

$$ \psi_{2t} = \sum_{n=0}^{\infty} [Q_n H_n^{(2)}(\beta_2 r) + Z_n H_n^{(1)}(\beta_2 r)] \sin(n\theta) e^{-i\omega t} \quad (10) $$

where $L_n$, $P_n$, $Q_n$, $Z_n$ are respectively the expansion coefficients of the P and SV waves in the tunnel lined.

$\alpha_2$, $\beta_2$ are the P and SV wavenumbers in the lining.

The six unknown expansion coefficients per mode ($E_n$, $F_n$, $L_n$, $P_n$, $Q_n$, $Z_n$) are determined by the boundary and interface conditions.

### Stresses and displacements in the rock mass and lining

Based on the relations between Helmholtz potentials, stresses and displacements, stress and displacement components in polar coordinates can be expressed as

$$ u_{r_k} = \frac{\partial \phi_k}{\partial r} + \frac{1}{r} \frac{\partial \psi_k}{\partial \theta} \quad (11) $$

$$ u_{\theta_k} = \frac{1}{r} \frac{\partial \phi_k}{\partial \theta} - \frac{\partial \psi_k}{\partial r} \quad (12) $$

$$ \varepsilon_{rr_k} = \frac{\partial u_{r_k}}{\partial r} \quad (13) $$

$$ \varepsilon_{\theta\theta_k} = \frac{1}{r} \frac{\partial u_{\theta_k}}{\partial \theta} + \frac{u_{r_k}}{r} \quad (14) $$

$$ \varepsilon_{r\theta_k} = \frac{1}{2} \left( \frac{\partial u_{\theta_k}}{\partial r} - \frac{u_{\theta_k}}{r} + \frac{1}{r} \frac{\partial u_{r_k}}{\partial \theta} \right) \quad (15) $$
$$ \sigma_{rr_k} = \lambda_k(\varepsilon_{rr_k} + \varepsilon_{\theta\theta_k}) + 2\mu_k\varepsilon_{rr_k} \quad (16) $$

$$ \sigma_{\theta\theta_k} = \lambda_k(\varepsilon_{rr_k} + \varepsilon_{\theta\theta_k}) + 2\mu_k\varepsilon_{\theta\theta_k} \quad (17) $$

$$ \sigma_{r\theta_k} = 2\mu_k\varepsilon_{r\theta_k} \quad (18) $$

Replacing equations (7)-(10) into equations (11)-(18), one can get the displacement and/or stress field components both in the rock mass and in the tunnel.

## 2 Boundary and interface conditions

The area under study is a karstic massif with secondary faults, the contact between the rock mass and the tunnel is not perfect. Three models have been found in the bibliographical survey to model the imperfect contact. The first one is a viscoelastic interface in which a linear spring and a linear dashpot are connected in parallel. The second one is an elastic slip interface. They are proposed by Fang and Jin (2016) [3] and Fang et al. (2018) [4]. In this proposal, we use another model developed by Achenbach and Zhu (1989) [2]. It assumes that tractions are continuous but displacements may be discontinuous across the interface. More precisely, jumps in the displacement components are assumed to be proportional to their respective interface traction components : the proportionality coefficient yields an interface parameter characterizing the kind of contact under consideration. Since we consider different proportionality values depending on the radial or hoop displacement components, two parameters are introduced, equation (19). These interface parameters are assumed to be uniform along the circumference of the interface. The interface model is then said to represent a homogeneously imperfect interface. Consequently, at the outer boundary of the lining ($r = b$), the interface conditions are written as

$$ \left\{ \begin{array}{l} \sigma_{rr_1} = \sigma_{rr_2} \\ \sigma_{r\theta_1} = \sigma_{r\theta_2} \\ u_{r_1} - \frac{1}{K_r}\sigma_{rr_1} = u_{r_2} \\ u_{\theta_1} - \frac{1}{K_\theta}\sigma_{r\theta_1} = u_{\theta_2} \end{array} \right. \quad (19) $$

where $K_r$ and $K_\theta$ are the interface normal and transverse spring constants, respectively. Physically when $K_r = K_\theta = 0$, the rock mass and lining are detached and when $K_r = K_\theta = \infty$ there is no displacement jump between the rock mass and the lining that is to say the contact is perfect.

At the inner boundary of the lining ($r = a$) a free stress surface is considered and the relative equations are

$$ \left\{ \begin{array}{l} \sigma_{rr_2} = 0 \\ \sigma_{r\theta_2} = 0 \end{array} \right. \quad (20) $$

Substituting equations [(11), (12), (16), (18)] into equations (19)-(20), and applying the orthogonality conditions, $n$ sets of linear algebra equations can be obtained as

$$ [X]_{n,\{6,6\}}\{C\}_{n,6} = \{Y\}_{\{n,6\}} \quad (21) $$

where $C$ is a six dimensional column vector containing the unknown expansion coefficients per mode, $X$ is a (6x6) square matrix and $Y$ a six dimensional column vector. They are
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Dimensionless Stress Concentration Factor (DSCF) in the lining at $r = a$, at high frequency ($\alpha_1 a = 2$), comparison with results of Yi et al. [1]</figcaption></figure>

composed by the terms resulting from the boundary and interface conditions (see Yi et al. (2016)). Solving the linear systems allows to get the six unknown coefficients per mode and thus the displacement and stress fields.

## 3 Numerical results and analysis

The boundary conditions at the interface ($r = b$) ensure continuity of stresses (radial and tangential) and displacements (radial and hoop). From the solution of system (21) we can calculate both the stress and displacement components. As an example we propose here to study the hoop stress, displacements (radial and hoop).

Parameters provided by Yi et al. (2016) [1] have been used for the validation of the model. The dynamic hoop stress, radial and hoop displacements are respectively non dimensionalized.

We consider the nondimensionalization process as in chapter (III) of Pao and Mow (1971)[6], the dynamic stress concentration factor (DSCF) is introduced as

$$ \sigma_{\theta\theta}^* = \left| \frac{\sigma_{\theta\theta}}{\sigma_0} \right| \qquad (22) $$

where $\sigma_0 = \mu_1 \beta_1^2 \varphi_0$ denotes the stress intensity of the incident wave in the direction of propagation.

And the dimensionless radial and hoop displacement can be written as follows

$$ u_r^* = \frac{u_r a}{\varphi_0} \quad \text{and} \quad u_\theta^* = \frac{u_\theta a}{\varphi_0} \qquad (23) $$

To compute the displacement and stress fields, the infinite summations of Bessel-Fourier series have to be truncated. Because of the oscillating behaviour, this is a complex task and particular attention has to be paid to the numerical procedure. The divergence of numerical
solutions depends on the incident frequency. To get a good convergence of the numerical solutions we chose to truncate the computations when the summation on the $n^{th}$ terms of the series remains stable, i.e when the contribution in modes of the values of stresses and displacements become negligible (typically we tested until 10 modes at low frequency and 15 at high frequency).

<figure><img src="image_4.png" /><figcaption>**Figure 3.** Dimensionless radial displacement in the lining at $r = a$, left subfigure at low frequency ($\alpha_1 a = 0.1$) and right subfigure at high frequency ($\alpha_1 a = 2$)</figcaption></figure>

The results presented are obtained from the first four modes at low frequency, and from the first thirteen modes at high frequency. We used a thousand points of discretization for the plotting of the respective polar diagrams. For the validation of this first approach, we compare our results with those obtained by Yi et al. (2016) [1]. Figure 2 shows the DSCF in the lining $r = a$ at high frequency $\alpha_1 a = 2$ and a good agreement is achieved.

Nevertheless, some discrepancies appear that may be explained by: the number of modes may be different; the software used for the polar representation may be also different; curves may have been smoothed due to the software.

Figures 3-5 present the DSCF, the radial and hoop displacements in the lining $r = a$ for different frequencies.

Some remarks can be underlined:

* Figure 4 -left: at low frequency the maximum stress is located at $\theta = \frac{\pi}{2}$ and the stress decreases when moving from the perfect contact to the imperfect contact.

* Figures 3 and 5 -left: at low frequency the maxima of radial and hoop displacements are respectively obtained at $\theta = 0$ and $\theta = \frac{\pi}{2}$, and their values slightly increase from the perfect contact to imperfect contact.

* Figures 3-5 -right: at high frequency, wide changes occur in the DSCF, radial and hoop displacements depending on values of $K_r$ and $K_{\theta}$. The dynamic response is different for
each case from the perfect to the imperfect contact. This phenomenon is related to the jump
at the interface $r = b$.

<figure><img src="image_5.png" /><figcaption>**Figure 4.** Dimensionless Stress Concentration Factor (DSCF) in the lining at $r = a$, left subfigure at low frequency ($\alpha_1 a = 0.1$) and right subfigure at high frequency ($\alpha_1 a = 2$)</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Figure 5.** Dimensionless hoop displacement in the lining at $r = a$, left subfigure at low frequency ($\alpha_1 a = 0.1$) and right subfigure at high frequency ($\alpha_1 a = 2$)</figcaption></figure>

# 4 Conclusion and perspectives

The dynamic response of a circular tunnel with imperfect surface interaction has been pre-
sented under the hypothesis of two homogeneous, elastic and isotropic materials for both the
rock mass and the lining. A two linear spring approach has been used to model the imperfect
contact between the rock and the lining. The parameters considered and estimated in this study are the Dimensionless (hoop) Stress Concentration Factor (DSCF) and the dimensionless radial and hoop displacements.

The results have indicated that the imperfect contact has a significant influence on the dynamic response of tunnel. At lower frequency, the stress is more impacted than the displacement field, but at higher frequency the displacement strongly depends on the kind of contact. Further works will also focus on the damaged area around the tunnel which is an important part to account for since the excavation of the tunnel weakened a specific zone. To tackle this, a multilayered cylindrical geometry will be developed. In order to get closer to the real area under study, another point of improvement of the model is to take into account the anisotropy of the soil and/or the tunnel, and the porous aspect of the materials.

## Acknowledgement

This work has been carried out under a PhD grant from Conseil Général Provence-Alpes-Côte d'Azur, France. This support is gratefully acknowledged.

## References

[1] C.P. Yi, W.b. Lu, P. Zhang, D. Johansson, U. Nyberg. Effect of imperfect interface on the dynamic response of a circular lined tunnel impacted by plane P-waves. **Tunnelling and Underground Space Technology** 21 68-74 (2016).

[2] J.D. Achenbach, H. Zhu. Effect of interfacial zone on mechanical behavior and failure of fiber reinforced composites. **Journal of the Mechanics and Physics of Solids** 37 381-393 (1989).

[3] X.Q. Fang, H.X. Jin. Visco-elastic imperfect bonding effect on dynamic response of a non circular lined tunnel subjected to P and SV waves. **Soil Dynamics and Earthquake Engineering** 88 1-7 (2016).

[4] X.Q. Fang, T.F. Zhang, H.Y. Li. Elastic-slip interface effect on dynamic response of a lined tunnel in a semi infinite alluvial valley under SH waves. **Tunnelling and Underground Space Technology** 74 96-106 (2018).

[5] W. Zhao, W. Chen, D. Yang. Interaction between strengthening and isolation layers for tunnels in rock subjected to SH waves. **Tunnelling and Underground Space Technology** 79 121-133 (2018).

[6] YH. Pao, C.C Mow. Diffraction of elastic waves and dynamic stress concentrations. Rand, Crane Russak, New York, 1971, 681 pages.

[7] A. Wirgin. Retrieval of the equivalent acoustic constitutive parameters of an inhomogeneous fluid-like object by nonlinear full waveform inversion. **Ultrasonics** 65 353-369 (2016).