# An efficient model to predict guided wave radiation by finite-sized sources in multilayered anisotropic plates with account of caustics

Mathilde Stévenin, Alain Lhémery, Sébastien Grondel

► To cite this version:

Mathilde Stévenin, Alain Lhémery, Sébastien Grondel. An efficient model to predict guided wave radiation by finite-sized sources in multilayered anisotropic plates with account of caustics. Journal of Physics: Conference Series, 2016, 684, pp.012004. 10.1088/1742-6596/684/1/012004. cea-01811899

HAL Id: cea-01811899

https://cea.hal.science/cea-01811899

Submitted on 13 Jun 2024

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
# An efficient model to predict guided wave radiation by finite-sized sources in multilayered anisotropic plates with account of caustics

Mathilde Stévenin, Alain Lhémery, Sébastien Grondel

► To cite this version:

Mathilde Stévenin, Alain Lhémery, Sébastien Grondel. An efficient model to predict guided wave radiation by finite-sized sources in multilayered anisotropic plates with account of caustics. Journal of Physics: Conference Series, 2016, 684, pp.012004. 10.1088/1742-6596/684/1/012004. cea-01811899

HAL Id: cea-01811899

https://cea.hal.science/cea-01811899v1

Submitted on 13 Jun 2024

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
# An efficient model to predict guided wave radiation by finite-sized sources in multilayered anisotropic plates with account of caustics

M Stévenin<fnref n="1" />,<fnref n="2" />, A Lhémery<fnref n="1" /> and S Grondel<fnref n="2" />

E-mail: alain.lhémery@cea.fr

**Abstract.** Elastic guided waves (GW) are used in various non-destructive testing (NDT) methods to inspect plate-like structures, generated by finite-sized transducers. Thanks to GW long range propagation, using a few transducers at permanent positions can provide a full coverage of the plate. Transducer diffraction effects take place, leading to complex radiated fields. Optimizing transducers positioning makes it necessary to accurately predict the GW field radiated by a transducer. Fraunhofer-like approximations applied to GW in isotropic homogeneous plates lead to fast and accurate field computation but can fail when applied to multi-layered anisotropic composite plates, as shown by some examples given.

Here, a model is proposed for composite plates, based on the computation of the approximate Green's tensor describing modal propagation from a source point, with account of caustics typically seen when strong anisotropy is concerned. Modal solutions are otherwise obtained by the Semi-Analytic Finite Element method. Transducer diffraction effects are accounted for by means of an angular integration over the transducer surface as seen from the calculation point, that is, over energy paths involved, which are mode-dependent. The model is validated by comparing its predictions with those computed by means of a full convolution integration of the Green's tensor with the source over transducer surface. Examples given concern disk and rectangular shaped transducers commonly used in NDT.

## 1. Introduction

In most non-destructive testing (NDT) methods used to inspect plate-like structures by means of elastic guided waves (GW), transducers of finite size are used to generate these waves in the structure. GW propagate at long range, so that inspecting a large structure may be made by using a few permanently attached transducers, avoiding any scan to cover the whole structure. However, transducer diffraction effects take place, leading to complex field radiation characteristics. The optimization of transducer positioning for insuring the expected coverage must be carefully addressed and makes it necessary to accurately predict the guided wave field radiated by a transducer. In the literature, Fraunhofer-like approximations [1] have been developed for elastic guided waves in the case of isotropic and homogeneous plates, leading to fast computation of typical diffraction effects. In the case of multi-layered anisotropic plates such as those used in the aircraft industry, this approach fails at predicting the field radiated; some examples are given herein to illustrate this point. To solve the problem of efficient computing of transducer diffraction effect in GW radiation into multilayered anisotropic plates, a new model

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

Content from this work may be used under the terms of the Creative Commons Attribution 3.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.

<footnote n="1">¹ CEA LIST, point courrier 120, bâtiment 565, 91191 Gif-Sur-Yvette Cedex, France</footnote>

<footnote n="2">² IEMN-DOAE, UVHC, Le Mont Houy, 59313 Valenciennes Cedex 9, France</footnote>
# An efficient model to predict guided wave radiation by finite-sized sources in multilayered anisotropic plates with account of caustics

M Stévenin<fnref n="1" />,<fnref n="2" />, A Lhémery<fnref n="1" /> and S Grondel<fnref n="2" />

E-mail: alain.lhémery@cea.fr

**Abstract.** Elastic guided waves (GW) are used in various non-destructive testing (NDT) methods to inspect plate-like structures, generated by finite-sized transducers. Thanks to GW long range propagation, using a few transducers at permanent positions can provide a full coverage of the plate. Transducer diffraction effects take place, leading to complex radiated fields. Optimizing transducers positioning makes it necessary to accurately predict the GW field radiated by a transducer. Fraunhofer-like approximations applied to GW in isotropic homogeneous plates lead to fast and accurate field computation but can fail when applied to multi-layered anisotropic composite plates, as shown by some examples given.

Here, a model is proposed for composite plates, based on the computation of the approximate Green's tensor describing modal propagation from a source point, with account of caustics typically seen when strong anisotropy is concerned. Modal solutions are otherwise obtained by the Semi-Analytic Finite Element method. Transducer diffraction effects are accounted for by means of an angular integration over the transducer surface as seen from the calculation point, that is, over energy paths involved, which are mode-dependent. The model is validated by comparing its predictions with those computed by means of a full convolution integration of the Green's tensor with the source over transducer surface. Examples given concern disk and rectangular shaped transducers commonly used in NDT.

## 1. Introduction

In most non-destructive testing (NDT) methods used to inspect plate-like structures by means of elastic guided waves (GW), transducers of finite size are used to generate these waves in the structure. GW propagate at long range, so that inspecting a large structure may be made by using a few permanently attached transducers, avoiding any scan to cover the whole structure. However, transducer diffraction effects take place, leading to complex field radiation characteristics. The optimization of transducer positioning for insuring the expected coverage must be carefully addressed and makes it necessary to accurately predict the guided wave field radiated by a transducer. In the literature, Fraunhofer-like approximations [1] have been developed for elastic guided waves in the case of isotropic and homogeneous plates, leading to fast computation of typical diffraction effects. In the case of multi-layered anisotropic plates such as those used in the aircraft industry, this approach fails at predicting the field radiated; some examples are given herein to illustrate this point. To solve the problem of efficient computing of transducer diffraction effect in GW radiation into multilayered anisotropic plates, a new model

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

Content from this work may be used under the terms of the Creative Commons Attribution 3.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.

<footnote n="1">¹ CEA LIST, point courrier 120, bâtiment 565, 91191 Gif-Sur-Yvette Cedex, France</footnote>

<footnote n="2">² IEMN-DOAE, UVHC, Le Mont Houy, 59313 Valenciennes Cedex 9, France</footnote>
is proposed. It is based on the computation of the approximate Green's tensor describing modal propagation from a source point, with account of caustics which arise in anisotropic plates. The modal solution itself is obtained thanks to well-established Semi-Analytic Finite Element method [2]. The overall principle of the method is to proceed to an angular integration over the transducer surface as seen from the calculation point, based upon the energy paths involved, which are mode-dependent. The first part of the paper is dedicated to a description of the proposed model. In the second part, the model is validated thanks to comparisons of its predictions to those computed by means of a full convolution integral of the Green's tensor with the source description. Examples given concern disk and rectangular shaped transducers which are the most commonly used sources in the NDT practice.

## 2. Formulation of the 3D Green's function

### 2.1. Geometry of the system

We consider an infinite plate as shown in Fig. 1, of thickness $2d$, which upper plane is at $z = d$ and lower plane at $z = -d$. The plate is assumed to be made of an arbitrary composite material described as multiple parallel layers of arbitrary anisotropy.

<figure><img src="image_3.png" /><figcaption>Figure 1. System geometry.</figcaption></figure>

### 2.2. Far field expression

The field radiated by a finite-sized source can be expressed as a convolution integral of the Green's function $g^{(3)}(x, y, z)$ with a source term $q^{(3)}(x, y)$:

$$u^{(3)}(x, y, z, q^{(3)}) = \int \int_{S} g^{(3)}(x - x', y - y', z)q^{(3)}(x', y')dx'dy'. \quad (1)$$

This Green's function can be written as a sum over the propagating modes [3]:

$$g^{(3)}(x, y, z) = \sum_{m} \sum_{n(\phi)} g_{m,n}^{(3)}(x, y, z). \quad (2)$$

These modal contributions are calculated thanks to the Semi Analytic Finite Element (SAFE) method [2]. With this time harmonic method, arbitrary sequences of multi-layered composite plates can be accounted for, results being given as propagative and evanescent modes described by their wavenumber and modal amplitude of components of the displacement along the plate
is proposed. It is based on the computation of the approximate Green's tensor describing modal propagation from a source point, with account of caustics which arise in anisotropic plates. The modal solution itself is obtained thanks to well-established Semi-Analytic Finite Element method [2]. The overall principle of the method is to proceed to an angular integration over the transducer surface as seen from the calculation point, based upon the energy paths involved, which are mode-dependent. The first part of the paper is dedicated to a description of the proposed model. In the second part, the model is validated thanks to comparisons of its predictions to those computed by means of a full convolution integral of the Green's tensor with the source description. Examples given concern disk and rectangular shaped transducers which are the most commonly used sources in the NDT practice.

## 2. Formulation of the 3D Green's function

### 2.1. Geometry of the system

We consider an infinite plate as shown in Fig. 1, of thickness $2d$, which upper plane is at $z = d$ and lower plane at $z = -d$. The plate is assumed to be made of an arbitrary composite material described as multiple parallel layers of arbitrary anisotropy.

<figure><img src="image_3.png" /><figcaption>**Figure 1.** System geometry.</figcaption></figure>

### 2.2. Far field expression

The field radiated by a finite-sized source can be expressed as a convolution integral of the Green's function $g^{(3)}(x, y, z)$ with a source term $q^{(3)}(x, y)$:

$$ u^{(3)}(x, y, z, q^{(3)}) = \int \int_{S} g^{(3)}(x - x', y - y', z)q^{(3)}(x', y')dx'dy'. \quad (1) $$

This Green's function can be written as a sum over the propagating modes [3]:

$$ g^{(3)}(x, y, z) = \sum_{m} \sum_{n(\phi)} g_{m,n}^{(3)}(x, y, z). \quad (2) $$

These modal contributions are calculated thanks to the Semi Analytic Finite Element (SAFE) method [2]. With this time harmonic method, arbitrary sequences of multi-layered composite plates can be accounted for, results being given as propagative and evanescent modes described by their wavenumber and modal amplitude of components of the displacement along the plate
thickness, for one given phase direction. The calculation is made once for all possible phase directions and results can be easily reused. The *n* subscript is never necessary when using isotropic plates. But, in some anisotropic cases it appears in several observation directions *φ* that for certain modes there are more than one phase contribution. Caustics can also appear for these modes.

To illustrate this, Fig. 2 presents the group velocity in mm.µs⁻¹ of the SHo mode at a frequency of 300 kHz in a plate whose characteristics are given on the section 3.2. For this mode at this frequency we must deal with caustics. On the left, in the considered direction of observation, there are three phase contributions. On the right, in the considered direction of observation, two phase contributions coalesce, leading to a caustic.

<figure><img src="image_4.png" /><figcaption>**Figure 2.** Three phase contributions far from caustics (left) and near a caustic (right).</figcaption></figure>

To calculate the Green's function, one must distinguish propagation directions close to a caustic from the others. Far from caustic the phase term can be developed at the second order approximation whereas near a caustic, a third order developed must be used. The expressions of these Green's functions can be found in literature [3, 4] and are the following:

$$g_{m,n}^{(3)}(x, y, z) = d(x, y, z)a_{m,n}(x, y, z)p_{m,n}(x, y, z), \quad (3)$$

with,

$$d_f(x, y) = \sqrt{\frac{x^2 + y^2}{2\pi}}^{-1/2}, \quad (4)$$

$$a_{fm,n}(x, y, z) = \operatorname{res}[G]_{k=k_{m,n}(\gamma_{m,n})} \alpha \frac{|k_{m,n}(\gamma_{m,n})|}{2\pi|\frac{\partial^2 \varphi_{m,n}(\gamma_{m,n}, \phi)}{\partial \gamma^2}|}, \quad (5)$$

where *f* index stands for "far from caustics" and with *G* being 2D spatial Fourier transform of *g*,

$$\varphi_{m,n}(\gamma_{m,n}, \phi) = k_{m,n}(\gamma_{m,n})\cos(\gamma_{m,n} - \phi),$$

*k*<sub>*m*,*n*</sub> is the wave number of the mode,

and *γ*<sub>*m*,*n*</sub> is the phase direction for a given direction of observation *φ*, and where *p* denotes the following phase term

$$p_{fm,n}(x, y, z) = e^{i \sqrt{\frac{x^2 + y^2}{\gamma_{m,n}(\gamma_{m,n}, \phi)}}} e^{i \frac{\pi}{4} \operatorname{sgn}(\frac{\partial^2 \varphi_{m,n}(\gamma_{m,n}, \phi)}{\partial \gamma^2})} e^{i \frac{\pi}{2}}, \quad (6)$$
thickness, for one given phase direction. The calculation is made once for all possible phase directions and results can be easily reused. The *n* subscript is never necessary when using isotropic plates. But, in some anisotropic cases it appears in several observation directions *φ* that for certain modes there are more than one phase contribution. Caustics can also appear for these modes.

To illustrate this, Fig. 2 presents the group velocity in mm.µs⁻¹ of the SHo mode at a frequency of 300 kHz in a plate whose characteristics are given on the section 3.2. For this mode at this frequency we must deal with caustics. On the left, in the considered direction of observation, there are three phase contributions. On the right, in the considered direction of observation, two phase contributions coalesce, leading to a caustic.

<figure><img src="image_4.png" /><figcaption>**Figure 2.** Three phase contributions far from caustics (left) and near a caustic (right).</figcaption></figure>

To calculate the Green's function, one must distinguish propagation directions close to a caustic from the others. Far from caustic the phase term can be developed at the second order approximation whereas near a caustic, a third order developed must be used. The expressions of these Green's functions can be found in literature [3, 4] and are the following:

$$g_{m,n}^{(3)}(x, y, z) = d(x, y, z)a_{m,n}(x, y, z)p_{m,n}(x, y, z), \quad (3)$$

with,

$$d_f(x, y) = \sqrt{\frac{x^2 + y^2}{2\pi}}^{-1/2}, \quad (4)$$

$$a_{fm,n}(x, y, z) = \operatorname{res}[G]_{k=k_{m,n}(\gamma_{m,n})} \alpha \frac{|k_{m,n}(\gamma_{m,n})|}{2\pi |\frac{\partial^2 \varphi_{m,n}(\gamma_{m,n}, \varphi)}{\partial \gamma^2}|}, \quad (5)$$

where *f* index stands for "far from caustics" and with *G* being 2D spatial Fourier transform of *g*,

$$\varphi_{m,n}(\gamma_{m,n}, \varphi) = k_{m,n}(\gamma_{m,n})\cos(\gamma_{m,n} - \varphi),$$

*k*<sub>*m*,*n*</sub> is the wave number of the mode,

and *γ*<sub>*m*,*n*</sub> is the phase direction for a given direction of observation *φ*, and where *p* denotes the following phase term

$$p_{fm,n}(x, y, z) = e^{i \sqrt{\frac{x^2 + y^2}{\gamma_{m,n}(\gamma_{m,n}, \phi)}}} e^{i \frac{\pi}{4} \operatorname{sgn}(\frac{\partial^2 \varphi_{m,n}(\gamma_{m,n}, \phi)}{\partial \gamma^2})} e^{i \frac{\pi}{2}}, \quad (6)$$
far from caustics as presented on Fig. 2 and

$$d_c(x, y, z) = \sqrt{x^2 + y^2}^{-1/3} \qquad (7)$$

$$a_{cm,(1)+(2)}(x, y, z) = \begin{aligned} & \underset{k=k_{m,1}(\gamma_{m,1})}{\operatorname{res}[G]} \left|k_{m,1}(\gamma_{m,1})\right| \frac{u - 2 \sqrt{S(\gamma_{m,1}, \phi)}}{\frac{\partial^2 \phi_{m,1}(\gamma_{m,1}, \phi)}{\partial \gamma^2}} \\ & + \underset{k=k_{m,2}(\gamma_{m,2})}{\operatorname{res}[G]} \left|k_{m,2}(\gamma_{m,2})\right| \frac{u - 2 \sqrt{S(\gamma_{m,2}, \phi)}}{\frac{\partial^2 \phi_{m,2}(\gamma_{m,2}, \phi)}{\partial \gamma^2}} \end{aligned} \qquad (8)$$

$$p_{cm,(1)+(2)}(x, y, z) = e^{i \sqrt{x^2+y^2}L(\phi)} A i(\kappa) \frac{i^{\frac{\pi}{2}}}{2} \qquad (9)$$

near caustics when the two phase contributions (1) and (2) cannot be distinguished and treated separately as interfering contributions, as seen on Fig. 2, with

$$L(\phi) = \frac{1}{2} (\phi_{m,1}(\gamma_{m,1}, \phi) + \phi_{m,2}(\gamma_{m,2}, \phi)),$$

$$\kappa = -\left(\sqrt{x^2 + y^2}\right)^{\frac{2}{3}} S(\phi),$$

$$S(\phi) = \frac{3}{4} (\phi_{m,1}(\gamma_{m,1}, \phi) + \phi_{m,2}(\gamma_{m,2}, \phi))^{\frac{2}{3}},$$

Ai the Airy's function and where c index stands for "close to caustics".

Finally, we want to calculate the following integral:

$$u^{(3)}(x, y, z, q^{(3)}) = \sum_{m,n(\phi)} \int_{S} g_{m,n}^{(3)}(x - x', y - y', z) q^{(3)}(x', y') dx' dy'. \qquad (10)$$

In the isotropic case, such an integral can be accurately and efficiently evaluated thanks to the Fraunhofer approximation [1]. In the next section, a Fraunhofer-like approximation is proposed for the anisotropic case.

### 3. Fraunhofer-like approximation

#### 3.1. Theory

To start with this approximation, modes having only one phase contribution in each direction are considered. The source is centred on $(x_c, y_c)$ and radiates with a uniform amplitude, so that for all source points, one has:

$$q^{(3)}(x', y') = q^{(3)}(x_c, y_c). \qquad (11)$$

The Fraunhofer approximation consists in approximating the phase term by its first order Taylor expansion and the amplitude term by its zero order Taylor expansion. By doing so, Eq. (3) can be re-written as:

$$\begin{aligned} d_f(x - x', y - y')a_{fm,n}(x - x', y - y', z)q^{(3)}(x', y') \\ = d_f(x - x_c, y - y_c)a_{fm,n}(x - x_c, y - y_c, z)q^{(3)}(x_c, y_c) \\ = \left(\sqrt{(x-x_c)^2 + (y-y_c)^2}\right)^{-1/2} \times \operatorname{res}[G]_{k=k_{m,n}(\gamma_{m,n})} \frac{q\left|\frac{k_{m,n}(\gamma_{m,n})}{2\pi|\frac{\partial^2\phi_{m,n}(\gamma_{m,n},\phi)}{\partial\gamma^2}|}\right|}{q^{(3)}(x_c, y_c)}, \end{aligned} \qquad (12)$$
far from caustics as presented on Fig. 2 and

$$d_c(x, y, z) = \sqrt{x^2 + y^2} - 1/3 \qquad (7)$$

$$a_{cm,(1)+(2)}(x, y, z) = \begin{aligned} & \underset{k=k_{m,1}(\gamma_{m,1})}{\operatorname{res}[G]} \left|k_{m,1}(\gamma_{m,1})\right| \underbrace{u - 2 \frac{\sqrt{S(\gamma_{m,1}, \phi)}}{\partial^2 \phi_{m,1}(\gamma_{m,1}, \phi)}}_{\Rightarrow} \\ & + \underset{k=k_{m,2}(\gamma_{m,2})}{\operatorname{res}[G]} \left|k_{m,2}(\gamma_{m,2})\right| \underbrace{u - \frac{\sqrt{S(\gamma_{m,2}, \phi)}}{\partial^2 \phi_{m,2}(\gamma_{m,2}, \phi)}}_{\Rightarrow} \end{aligned} \qquad (8)$$

$$p_{cm,(1)+(2)}(x, y, z) = e^{i \sqrt{x^2 + y^2 L(\phi)}} A i(\kappa) \frac{i^{\frac{\pi}{2}}}{2} \qquad (9)$$

near caustics when the two phase contributions (1) and (2) cannot be distinguished and treated separately as interfering contributions, as seen on Fig. 2, with

$$L(\phi) = \frac{1}{2} (\phi_{m,1}(\gamma_{m,1}, \phi) + \phi_{m,2}(\gamma_{m,2}, \phi)),$$

$$\kappa = -\left(\sqrt{x^2 + y^2}\right)^{\frac{2}{3}} S(\phi),$$

$$S(\phi) = \frac{3}{4} (\phi_{m,1}(\gamma_{m,1}, \phi) + \phi_{m,2}(\gamma_{m,2}, \phi))^{\frac{2}{3}},$$

Ai the Airy's function and where c index stands for "close to caustics".

Finally, we want to calculate the following integral:

$$u^{(3)}(x, y, z, q^{(3)}) = \sum_{m,n(\phi)} \int_{S} g_{m,n}^{(3)}(x - x', y - y', z) q^{(3)}(x', y') dx' dy'. \qquad (10)$$

In the isotropic case, such an integral can be accurately and efficiently evaluated thanks to the Fraunhofer approximation [1]. In the next section, a Fraunhofer-like approximation is proposed for the anisotropic case.

### 3. Fraunhofer-like approximation

#### 3.1. Theory

To start with this approximation, modes having only one phase contribution in each direction are considered. The source is centred on $(x_c, y_c)$ and radiates with a uniform amplitude, so that for all source points, one has:

$$q^{(3)}(x', y') = q^{(3)}(x_c, y_c). \qquad (11)$$

The Fraunhofer approximation consists in approximating the phase term by its first order Taylor expansion and the amplitude term by its zero order Taylor expansion. By doing so, Eq. (3) can be re-written as:

$$d_f(x - x', y - y')a_{fm,n}(x - x', y - y', z)q^{(3)}(x', y') = d_f(x - x_c, y - y_c)a_{fm,n}(x - x_c, y - y_c, z)q^{(3)}(x_c, y_c) = \left(\sqrt{(x-x_c)^2 + (y-y_c)^2}\right)^{-1/2} \times \operatorname{res}[G]_{k=k_{m,n}(\gamma_{m,n})} \frac{q\left|\frac{k_{m,n}(\gamma_{m,n})}{2\pi|\frac{\partial^2\phi_{m,n}(\gamma_{m,n},\phi)}{\partial\gamma^2}|}\right|}{q^{(3)}(x_c, y_c)} \qquad (12)$$
$$p_{fm,n}(x-x^r, y-y^r, z) = e^{i\varphi_{m,n}(\gamma_{m,n}, \phi_c)} \frac{\sqrt{(x-x_c)^2+(y-y_c)^2} e^{-\frac{i\pi}{4} \operatorname{sgn}(\frac{\partial^2\varphi_{m,n}(\gamma_{m,n}, \phi_c)}{\partial\nu^2})}}{e^{i\varphi_{m,n}(\gamma_{m,n}, \phi_c)} \frac{\sqrt{(x-x_c)^2+(y-y_c)^2} + \sqrt{(\gamma_{m,n}-\phi_c)(\gamma_{m,n}-\phi_c)}}{(x-x_c)^2+(y-y_c)^2}} e^{i\frac{\pi}{2}} \quad (13)$$

Fields radiated by most commonly used disk or rectangular transducers are now considered.

### Disk-shaped transducers

The transducer is of radius $a$. To ease the calculation, a polar coordinate system is now used, so that the radiation integral is now given by:

$$u^{(3)}(x, y, z, q^{(3)}) = \sum_{m} d_f(x-x_c, y-y_c) a_{fm,n} (x-x_c, y-y_c, z) q^{(3)}(x, y) e^{i\frac{\pi}{2}} \times e^{i\frac{\pi}{4} \operatorname{sgn}(\frac{\partial^2\varphi_{m,n}(\gamma_{m,n}, \phi_c)}{\partial\nu^2})} \int_{-\pi}^{\pi} \int_{0}^{a} e^{i\varphi_{m,n}(\gamma_{m,n}, \phi_c)R(1-\cos(\vartheta-\phi_c)\frac{r'}{R})} r dr d\vartheta, \quad (14)$$

where $R = \sqrt{(x-x_c)^2 + (y-y_c)^2}$. Classical wave diffraction calculations lead to the following final expression:

$$u^{(3)}(x, y, z, q^{(3)}) = \sum_{m} d_f(x-x_c, y-y_c) a_{fm,n} (x-x_c, y-y_c, z) q^{(3)}(x, y) e^{i\frac{\pi}{2}} \times e^{i\frac{\pi}{4} \operatorname{sgn}(\frac{\partial^2\varphi_{m,n}(\gamma_{m,n}, \phi_c)}{\partial\nu^2})} \pi a^2 \frac{2J(\varphi_c - \varphi_m)}{\varphi_{m,n}(\gamma_{m,n}, \phi_c) a} e^{i\varphi_{m,n}(\gamma_{m,n}, \phi_c)R,} \quad (15)$$

where the amplitude is proportional to the Fourier transform of the circular aperture.

### Rectangular-shaped transducers

A rectangular transducer of aperture $l \times L$ is considered, so that radiation integral to calculate is given by:

$$u^{(3)}(x, y, z, q^{(3)}) = \sum_{m} d_f(x-x_c, y-y_c) a_{fm,n} (x-x_c, y-y_c, z) q^{(3)}(x_c, y_c) e^{i\frac{\pi}{2}} \times e^{i\frac{\pi}{4} \operatorname{sgn}(\frac{\partial^2\varphi_{m,n}(\gamma_{m,n}, \phi_c)}{\partial\nu^2})} e^{i\varphi_{m,n}(\gamma_{m,n}, \phi_c)} \frac{\sqrt{(x-x_c)^2+(y-y_c)^2}}{\sqrt{(x-x_c)^2+(y-y_c)^2} + \sqrt{(\gamma_{m,n}-\phi_c)(\gamma_{m,n}-\phi_c)}} \times e^{i\varphi_{m,n}(\gamma_{m,n}, \phi_c)} \frac{\sqrt{(x-x_c)^2+(y-y_c)^2}}{(x-x_c)^2+(y-y_c)^2} dx^r dy^r. \quad (16)$$

Similarly, classical diffraction calculations lead to the following result:

$$u^{(3)}(x, y, z, q^{(3)}) = \sum_{m} d_f(x-x_c, y-y_c) a_{fm,n} (x-x_c, y-y_c, z) q^{(3)}(x_c, y_c) e^{i\frac{\pi}{2}} \times e^{i\frac{\pi}{4} \operatorname{sgn}(\frac{\partial^2\varphi_{m,n}(\gamma_{m,n}, \phi_c)}{\partial\nu^2})} ILSinc(\frac{\varphi_{m,n}(\gamma_{m,n}, \phi_c)(x-x_c)}{(x-x_c)^2+(y-y_c)^2}) \times \text{sinc}(\frac{\varphi_{m,n}(\gamma_{m,n}, \phi_c)(y-y_c)}{(x-x_c)^2+(y-y_c)^2}) e^{i\varphi_{m,n}(\gamma_{m,n}, \phi_c)} \sqrt{(x-x_c)^2+(y-y_c)^2}. \quad (17)$$

where the amplitude is now proportional to the Fourier transform of the rectangular aperture. If an isotropic plate is considered momentarily, above results simplify and are expressed as given by A. Raghavan and C.E.S Cesnik [1]. In this paper, the Fraunhofer approximation is shown to give very accurate results.
$$p_{fm,n}(x-x^r, y-y^r, z) = e^{i\varphi_{m,n}(\gamma_{m,n}, \phi_c)} \frac{\sqrt{(x-x_c)^2+(y-y_c)^2} e^{-\frac{i\pi}{4} \operatorname{sgn}(\frac{\partial^2\varphi_{m,n}(\gamma_{m,n}, \phi_c)}{\partial\gamma^2})}}{e^{i\varphi_{m,n}(\gamma_{m,n}, \phi_c)}(\frac{\sqrt{(x-x_c)(x-x')}}{(x-x_c)^2+(y-y_c)^2} + \frac{\sqrt{(y-y_c)(y-y')}}{(x-x_c)^2+(y-y_c)^2})} e^{\frac{i\pi}{2}} \quad (13)$$

Fields radiated by most commonly used disk or rectangular transducers are now considered.

### Disk-shaped transducers

The transducer is of radius $a$. To ease the calculation, a polar coordinate system is now used, so that the radiation integral is now given by:

$$u^{(3)}(x, y, z, q^{(3)}) = \sum_{m} d_f(x-x_c, y-y_c) a_{fm,n} (x-x_c, y-y_c, z) q^{(3)}(x, y) e^{\frac{i\pi}{2}} \\ \times e^{\frac{i\pi}{4} \operatorname{sgn}(\frac{\partial^2\varphi_{m,n}(\gamma_{m,n}, \phi_c)}{\partial\gamma^2})} \int_{-\pi}^{\pi} \int_{0}^{a} e^{i\varphi_{m,n}(\gamma_{m,n}, \phi_c)R(1-\cos(\vartheta-\phi_c))} r^2 dr d\vartheta, \quad (14)$$

where $R = \sqrt{(x-x_c)^2 + (y-y_c)^2}$. Classical wave diffraction calculations lead to the following final expression:

$$u^{(3)}(x, y, z, q^{(3)}) = \sum_{m} d_f(x-x_c, y-y_c) a_{fm,n} (x-x_c, y-y_c, z) q^{(3)}(x, y) e^{\frac{i\pi}{2}} \\ \times e^{\frac{i\pi}{4} \operatorname{sgn}(\frac{\partial^2\varphi_{m,n}(\gamma_{m,n}, \phi_c)}{\partial\gamma^2})} \pi a^2 \frac{2J(\varphi_c - \varphi_m)}{\varphi_{m,n}(\gamma_{m,n}, \phi_c) a} e^{i\varphi_{m,n}(\gamma_{m,n}, \phi_c)R,} \quad (15)$$

where the amplitude is proportional to the Fourier transform of the circular aperture.

### Rectangular-shaped transducers

A rectangular transducer of aperture $l \times L$ is considered, so that radiation integral to calculate is given by:

$$u^{(3)}(x, y, z, q^{(3)}) = \sum_{m} d_f(x-x_c, y-y_c) a_{fm,n} (x-x_c, y-y_c, z) q^{(3)}(x_c, y_c) e^{\frac{i\pi}{2}} \\ \times e^{\frac{i\pi}{4} \operatorname{sgn}(\frac{\partial^2\varphi_{m,n}(\gamma_{m,n}, \phi_c)}{\partial\gamma^2})} e^{i\varphi_{m,n}(\gamma_{m,n}, \phi_c)} \frac{\sqrt{(x-x_c)^2+(y-y_c)^2}}{\sqrt{(x-x_c)^2+(y-y_c)^2}} \\ \times \int_{x_c+L/2}^{x_c-L/2} \int_{y_c+L/2}^{y_c-L/2} e^{i\varphi_{m,n}(\gamma_{m,n}, \phi_c)\left(\frac{\sqrt{(x-x_c)(x-x')}}{(x-x_c)^2+(y-y_c)^2} + \frac{\sqrt{(y-y_c)(y-y')}}{(x-x_c)^2+(y-y_c)^2}\right)} dx^r dy^r. \quad (16)$$

Similarly, classical diffraction calculations lead to the following result:

$$u^{(3)}(x, y, z, q^{(3)}) = \sum_{m} d_f(x-x_c, y-y_c) a_{fm,n} (x-x_c, y-y_c, z) q^{(3)}(x_c, y_c) e^{\frac{i\pi}{2}} \\ \times e^{\frac{i\pi}{4} \operatorname{sgn}(\frac{\partial^2\varphi_{m,n}(\gamma_{m,n}, \phi_c)}{\partial\gamma^2})} I_1 \operatorname{sinc}\left(\frac{\varphi_{m,n}(\gamma_{m,n}, \phi_c)(x-x_c)}{(x-x_c)^2+(y-y_c)^2}\right) \\ \times \operatorname{sinc}\left(\frac{\varphi_{m,n}(\gamma_{m,n}, \phi_c)(y-y_c)}{(x-x_c)^2+(y-y_c)^2}\right) e^{i\varphi_{m,n}(\gamma_{m,n}, \phi_c)} \sqrt{(x-x_c)^2+(y-y_c)^2}. \quad (17)$$

where the amplitude is now proportional to the Fourier transform of the rectangular aperture. If an isotropic plate is considered momentarily, above results simplify and are expressed as given by A. Raghavan and C.E.S Cesnik [1]. In this paper, the Fraunhofer approximation is shown to give very accurate results.
<table><caption>**Table 1.** Plate characteristics.</caption><thead><tr><th>C11<br>(GPa)</th><th>C22=C33<br>(GPa)</th><th>C12=C13<br>(GPa)</th><th>C23<br>(GPa)</th><th>C44<br>(GPa)</th><th>C55=C66<br>(GPa)</th><th>Mass density<br>(kg/m3)</th><th>Ply thickness<br>(mm)</th></tr></thead><tbody><tr><td>123.4</td><td>11.5</td><td>5.6</td><td>6.4</td><td>2.6</td><td>4.5</td><td>1.6x10<sup>3</sup></td><td>0.25</td></tr></tbody></table>

## 3.2. *Simulations*

Simulations made using above given results are now presented for a composite plate.

### *Plate characteristics*

The plate is 1mm-thick, made of four parallel layers of fiber-reinforced polymer(CFRP) with the following symmetry [0°/90°]s(T700GC/M21 cross-ply composite). Material parameters are listed in Table 1.

Figure 3 presents the slowness curves of the three propagating modes which exist at our work frequency of 300 kHz, below the first cut-off frequency of the plate. For each of these modes the dependency of phase directions on the observation (energy) direction is plotted on Fig. 4.

<figure><img src="image_5.png" /><figcaption>**Figure 3.** Slowness curves in s.km⁻¹ in a cross-ply [0°/90°]s T700GC/M21 CFRP at a frequency of 300kHz.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Figure 4.** Phase directions function of the observation direction at a frequency of 300kHz.</figcaption></figure>

### *Limits of the Fraunhofer like approximation*

Results obtained in the previous section are now applied to the composite plate. In the slowness curves of Fig. 3, and considering the variations of phase directions with energy direction shown in Fig. 4, it appears that the Ao mode is the less anisotropic one: its slowness curve is quite circular and the direction of phase velocity is almost equal to that of the energy velocity. Therefore, the Fraunhofer approximation is expected to lead to accurate field prediction as it does in the isotropic case. In the computation, a 5-mm-diam disk transducer producing a uniform normal stress at the plate upper surface is considered. The computation is made for observation points located at 100mm of the source centre. Results computed using the approximate formula are compared to results obtained by computing the full surface convolution integral over the source active area.
<table><caption>**Table 1.** Plate characteristics.</caption><thead><tr><th>C11<br>(GPa)</th><th>C22=C33<br>(GPa)</th><th>C12=C13<br>(GPa)</th><th>C23<br>(GPa)</th><th>C44<br>(GPa)</th><th>C55=C66<br>(GPa)</th><th>Mass density<br>(kg/m3)</th><th>Ply thickness<br>(mm)</th></tr></thead><tbody><tr><td>123.4</td><td>11.5</td><td>5.6</td><td>6.4</td><td>2.6</td><td>4.5</td><td>1.6x10<sup>3</sup></td><td>0.25</td></tr></tbody></table>

## 3.2. Simulations

Simulations made using above given results are now presented for a composite plate.

### Plate characteristics

The plate is 1mm-thick, made of four parallel layers of fiber-reinforced polymer(CFRP) with the following symmetry [$0°/90°$]s(T700GC/M21 cross-ply composite). Material parameters are listed in Table 1.

Figure 3 presents the slowness curves of the three propagating modes which exist at our work frequency of 300 kHz, below the first cut-off frequency of the plate. For each of these modes the dependency of phase directions on the observation (energy) direction is plotted on Fig. 4.

<figure><img src="image_5.png" /><figcaption>**Figure 3.** Slowness curves in s.km⁻¹ in a cross-ply [$0°/90°$]s T700GC/M21 CFRP at a frequency of 300kHz.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Figure 4.** Phase directions function of the observation direction at a frequency of 300kHz.</figcaption></figure>

### Limits of the Fraunhofer like approximation

Results obtained in the previous section are now applied to the composite plate. In the slowness curves of Fig. 3, and considering the variations of phase directions with energy direction shown in Fig. 4, it appears that the Ao mode is the less anisotropic one: its slowness curve is quite circular and the direction of phase velocity is almost equal to that of the energy velocity. Therefore, the Fraunhofer approximation is expected to lead to accurate field prediction as it does in the isotropic case. In the computation, a 5-mm-diam disk transducer producing a uniform normal stress at the plate upper surface is considered. The computation is made for observation points located at 100mm of the source centre. Results computed using the approximate formula are compared to results obtained by computing the full surface convolution integral over the source active area.
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Out of plane displacement of the Ao mode at a distance of 100mm in the case of an excitation by a circular source at the frequency of 300 kHz obtained by computing the full integration (red line) and the approximate formula (blue line).</figcaption></figure>

The two kinds of results superimpose quite well but there are angular regions where approximate results are clearly wrong. Discrepancies appear in directions where small variation of the energy direction (e.g 5°) leads to large variation of phase direction (15°). The sole phase direction used in the approximation does not permit to accurately account for phase variations over the source surface. As the mode considered is the less anisotropic, such a failure is of course even worse when more anisotropic modes are considered. It is therefore obvious that Fraunhofer approximation cannot be reliably used for anisotropic plates. To overcome this difficulty, still aiming at deriving computationally efficient formulae, a new model is now proposed.

## 4. Integration along energy directions

### 4.1. Theory

We will use a new coordinate system in this part, as shown by Fig. 6. The system is centred on the calculation point and rotated so that $e_x$ is along the direction linking the calculation point to the source centre.

In the new coordinate system, the source centre is at $(R, 0)$ with $R = \frac{\sqrt{(x-x_c)^2 + (y-y_c)^2}}{m}$. The radiated field is now expressed as:

$$ u^{(3)}(e_x = 0, e_y = 0, z, q^{(3)}) = \sum_{m} \sum_{n(\phi_s)} \int_{0}^{\pi} d(-e'_{x} - e'_{y} z) a_{m,n}(-e'_{x} - e'_{y} z) \times p_{m,n}(-e'_{x} - e'_{y} z) q^{(3)}(e'_{x} - e'_{y} z) de'_{x} de'_{y} \quad (18) $$

The field expression is reformulated to explicit the angular dependence on $\vartheta$, angle formed by the line linking the calculation point to the source centre and that linking the calculation point
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Out of plane displacement of the Ao mode at a distance of 100mm in the case of an excitation by a circular source at the frequency of 300 kHz obtained by computing the full integration (red line) and the approximate formula (blue line).</figcaption></figure>

The two kinds of results superimpose quite well but there are angular regions where approximate results are clearly wrong. Discrepancies appear in directions where small variation of the energy direction (e.g 5°) leads to large variation of phase direction (15°). The sole phase direction used in the approximation does not permit to accurately account for phase variations over the source surface. As the mode considered is the less anisotropic, such a failure is of course even worse when more anisotropic modes are considered. It is therefore obvious that Fraunhofer approximation cannot be reliably used for anisotropic plates. To overcome this difficulty, still aiming at deriving computationally efficient formulae, a new model is now proposed.

## 4. Integration along energy directions

### 4.1. Theory

We will use a new coordinate system in this part, as shown by Fig. 6. The system is centred on the calculation point and rotated so that $e_x$ is along the direction linking the calculation point to the source centre.

In the new coordinate system, the source centre is at $(R, 0)$ with $R = \frac{\sqrt{(x-x_c)^2 + (y-y_c)^2}}{m}$. The radiated field is now expressed as:

$$u^{(3)}(e_x = 0, e_y = 0, z, q^{(3)}) = \sum_{m} \sum_{n(\phi_s)} \int_{0}^{\pi} d(-e'_{x} - e'_{y} z) a_{m,n}(-e'_{x} - e'_{y} z) \times p_{m,n}(-e'_{x} - e'_{y} z) q^{(3)}(e'_{x} - e'_{y} z) de'_{x} de'_{y} \quad (18)$$

The field expression is reformulated to explicit the angular dependence on $\vartheta$, angle formed by the line linking the calculation point to the source centre and that linking the calculation point
to the sourcepoint:

$$ \vartheta = \phi' - \phi_c \qquad (19) $$

with $\phi'$ the direction between the source point and the observation point and $\phi_c$ direction between the centre of the source and the observation point in the former coordinate system.

<figure><img src="image_8.png" /><figcaption>Figure 6. New coordinate system.</figcaption></figure>

In our integral $r$ is the distance between the source and calculation points. In this new formulation $d(x, y) = d(r)$:

$$ u^{(3)}(e_x = 0, e_y = 0, z, q^{(3)}) = \sum_{m} \sum_{n(\varphi)} \int_{\vartheta_{\min}}^{\vartheta_{\max}} d \frac{r_1(\vartheta) + r_2(\vartheta)}{2} \times \int_{r_2(\vartheta)}^{r_1(\vartheta)} a_{m,n}(\vartheta, z) p_{m,n}(r, \vartheta, z) q^{(3)}(r, \vartheta) r^r dr d\vartheta, \qquad (20) $$

$[\vartheta_{\min}, \vartheta_{\max}]$ and $[r_1(\vartheta), r_2(\vartheta)]$ denote the integration intervals as shown on the Fig. 7.

<figure><img src="image_9.png" /><figcaption>Figure 7. Integration limits.</figcaption></figure>

Calculation points must be outside the source. Moreover, we still consider a uniform source distribution.
to the sourcepoint:

$$ \vartheta = \phi' - \phi_c \qquad (19) $$

with $\phi'$ the direction between the source point and the observation point and $\phi_c$ direction between the centre of the source and the observation point in the former coordinate system.

<figure><img src="image_8.png" /><figcaption>Figure 6. New coordinate system.</figcaption></figure>

In our integral $r$ is the distance between the source and calculation points. In this new formulation $d(x, y) = d(r)$:

$$ u^{(3)}(e_x = 0, e_y = 0, z, q^{(3)}) = \sum_{m} \sum_{n(\phi)} \int_{\vartheta_{\min}}^{\vartheta_{\max}} d \frac{r_1(\vartheta) + r_2(\vartheta)}{2} \times \int_{r_1(\vartheta)}^{r_2(\vartheta)} a_{m,n}(\vartheta, z) p_{m,n}(r, \vartheta, z) q^{(3)}(r, \vartheta) r^r dr d\vartheta, \qquad (20) $$

$[\vartheta_{\min}, \vartheta_{\max}]$ and $[r_1(\vartheta), r_2(\vartheta)]$ denote the integration intervals as shown on the Fig. 7.

<figure><img src="image_9.png" /><figcaption>Figure 7. Integration limits.</figcaption></figure>

Calculation points must be outside the source. Moreover, we still consider a uniform source distribution.
$$ I_{m,n}(e_x = 0, e_y = 0) = \int_{\vartheta_{min}}^{\vartheta_{max}} d \frac{r_1(\vartheta) + r_2(\vartheta)}{2} \int_{r_1(\vartheta)}^{r_2(\vartheta)} a_{m,n}(\vartheta, z) p_{m,n}(r, \vartheta, z) r^r dr d\vartheta \\ = \frac{\vartheta_{max} - \vartheta_{min}}{n_{\vartheta}} d \frac{r_1(\vartheta_t) + r_2(\vartheta_t)}{2} \\ \times \int_{r_2(\vartheta_t)}^{t_1(\vartheta_t)} \int_{r_1(\vartheta_t)}^{r_2(\vartheta_t)} a_{m,n}(\vartheta_t, z) p_{m,n}(r, \vartheta_t, z) r^r dr^r, \quad (21) $$

with $\vartheta_t = \vartheta_{min} + (t - \frac{1}{2}) \frac{\vartheta_{max} - \vartheta_{min}}{n_{\vartheta}}$ and $n_{\vartheta}$ the number of discrete directions in angular integration.

For directions far from a caustic, we have:

$$ d \frac{r_1(\vartheta_t) + r_2(\vartheta_t)}{2} \int_{r_1(\vartheta_t)}^{r_2(\vartheta_t)} a_{m,n}(\vartheta_t, z) p_{m,n}(r, \vartheta_t, z) r^r dr \\ = d_f \frac{r(\vartheta) + r(\vartheta)}{2} e^{\frac{i\pi}{4} \operatorname{sgn}(\frac{\partial^2 \varphi_{m,n}(\gamma_{m,n}, \vartheta_t)}{\partial \gamma^2})} e^{\frac{\pi}{4} a_{fm,n}(\vartheta_t, z)} \\ \times e^{ir_2(\vartheta_t) \varphi_{m,n}(\gamma_{m,n}, \vartheta_t)} \frac{ir_2(\vartheta_t) \varphi_{m,n}(\gamma_{m,n}, \vartheta_t) - 1}{(i \varphi_{m,n}(\gamma_{m,n}, \vartheta_t))^2} \\ - e^{ir_1(\vartheta_t) \varphi_{m,n}(\gamma_{m,n}, \vartheta_t)} \frac{ir_1(\vartheta_t) \varphi_{m,n}(\gamma_{m,n}, \vartheta_t) - 1}{(i \varphi_{m,n}(\gamma_{m,n}, \vartheta_t))^2}. \quad (22) $$

For directions close to a caustic, we have:

$$ d \frac{r_1(\vartheta_t) + r_2(\vartheta_t)}{2} \int_{r_1(\vartheta_t)}^{r_2(\vartheta_t)} a_{m,n}(\vartheta_t, z) p_{m,n}(r, \vartheta_t, z) r^r dr \\ = d_c \frac{r_1(\vartheta_t) + r_2(\vartheta_t)}{2} \frac{e^{\frac{i\pi}{2}}}{2} A i (\kappa(\vartheta_t)) a_{cm,n}(\vartheta_t, z) \quad (23) \\ \times e^{ir_2(\vartheta_t) L(\vartheta_t)} \frac{ir_2(\vartheta_t) L(\vartheta_t) - 1}{(i L(\vartheta_t))^2} - e^{ir_1(\vartheta_t) L(\vartheta_t)} \frac{ir_1(\vartheta_t) L(\vartheta_t) - 1}{(i L(\vartheta_t))^2}, $$

$$ \text{with, } \kappa(\vartheta_t) = - \frac{r_1(\vartheta_t) + r_2(\vartheta_t)}{2} \frac{2}{3} S(\vartheta_t). $$

This general formulation involving integration along energy directions to express the field is now used for circular or rectangular sources. Integration intervals in calculation point coordinate system are easily found for both of them.

## 4.2. Results

To validate the proposed model of integration, results predicted by means of this model are compared to equivalent results computed by means of the full convolution surface integral. The cases of a disk transducer and of a rectangular one are considered. The plate considered is the same as that used in previous computations.

### Disk-shaped transducers

A disk transducer of radius 5 mm is considered. Calculation points are located at a distance of 100mm from the source centre. The amplitude displayed is the absolute value of the out-of-plane component of the displacement.
$$ I_{m,n}(e_x = 0, e_y = 0) = \int_{\vartheta_{min}}^{\vartheta_{max}} d \frac{r_1(\vartheta) + r_2(\vartheta)}{2} \int_{r_1(\vartheta)}^{r_2(\vartheta)} a_{m,n}(\vartheta, z) p_{m,n}(r, \vartheta, z) r^r dr d\vartheta \\ = \frac{\vartheta_{max} - \vartheta_{min}}{n_{\vartheta}} d \frac{r_1(\vartheta_t) + r_2(\vartheta_t)}{2} \\ \times \int_{r_2(\vartheta_t)}^{t_1(\vartheta_t)} \int_{r_1(\vartheta_t)}^{r_2(\vartheta_t)} a_{m,n}(\vartheta_t, z) p_{m,n}(r, \vartheta_t, z) r^r dr^r, \qquad (21) $$

with $\vartheta_t = \vartheta_{min} + (t - \frac{1}{2}) \frac{\vartheta_{max} - \vartheta_{min}}{n_{\vartheta}}$ and $n_{\vartheta}$ the number of discrete directions in angular integration.

For directions far from a caustic, we have:

$$ d \frac{r_1(\vartheta_t) + r_2(\vartheta_t)}{2} \int_{r_1(\vartheta_t)}^{r_2(\vartheta_t)} a_{m,n}(\vartheta_t, z) p_{m,n}(r, \vartheta_t, z) r^r dr \\ = d_f \frac{r(\vartheta) + r(\vartheta)}{2} e^{\frac{i\pi}{4} \operatorname{sgn}(\frac{\partial^2 \varphi_{m,n}(\gamma_{m,n}, \vartheta_t)}{\partial \gamma^2})} e^{\frac{\pi}{4} a_{fm,n}(\vartheta_t, z)} \\ \times e^{ir_2(\vartheta_t) \varphi_{m,n}(\gamma_{m,n}, \vartheta_t)} \frac{ir_2(\vartheta_t) \varphi_{m,n}(\gamma_{m,n}, \vartheta_t) - 1}{(i \varphi_{m,n}(\gamma_{m,n}, \vartheta_t))^2} \\ - e^{ir_1(\vartheta_t) \varphi_{m,n}(\gamma_{m,n}, \vartheta_t)} \frac{ir_1(\vartheta_t) \varphi_{m,n}(\gamma_{m,n}, \vartheta_t) - 1}{(i \varphi_{m,n}(\gamma_{m,n}, \vartheta_t))^2}. \qquad (22) $$

For directions close to a caustic, we have:

$$ d \frac{r_1(\vartheta_t) + r_2(\vartheta_t)}{2} \int_{r_1(\vartheta_t)}^{r_2(\vartheta_t)} a_{m,n}(\vartheta_t, z) p_{m,n}(r, \vartheta_t, z) r^r dr \\ = d_c \frac{r_1(\vartheta_t) + r_2(\vartheta_t)}{2} \frac{e^{\frac{i\pi}{2}}}{2} A i(\kappa(\vartheta_t)) a_{cm,n}(\vartheta_t, z) \qquad (23) \\ \times e^{ir_2(\vartheta_t) L(\vartheta_t)} \frac{ir_2(\vartheta_t) L(\vartheta_t) - 1}{(i L(\vartheta_t))^2} - e^{ir_1(\vartheta_t) L(\vartheta_t)} \frac{ir_1(\vartheta_t) L(\vartheta_t) - 1}{(i L(\vartheta_t))^2}, $$

$$ \text{with, } \kappa(\vartheta_t) = - \frac{r_1(\vartheta_t) + r_2(\vartheta_t)}{2} \frac{2}{3} S(\vartheta_t). $$

This general formulation involving integration along energy directions to express the field is now used for circular or rectangular sources. Integration intervals in calculation point coordinate system are easily found for both of them.

## 4.2. Results

To validate the proposed model of integration, results predicted by means of this model are compared to equivalent results computed by means of the full convolution surface integral. The cases of a disk transducer and of a rectangular one are considered. The plate considered is the same as that used in previous computations.

### Disk-shaped transducers

A disk transducer of radius 5 mm is considered. Calculation points are located at a distance of 100mm from the source centre. The amplitude displayed is the absolute value of the out-of-plane component of the displacement.
<figure><img src="image_10.png" /><figcaption>**Figure 8.** Out of plane displacement of the Ao mode (a), So mode (b), SHo mode (c) at a distance of 100mm in the case of an excitation by a circular source at the frequency of 300 kHz obtained by computing the full integration (solid line) and the approximate formula (+ points).</figcaption></figure>

Figure 8 (a) displays the same results as those shown on Fig. 5. In the previous results, the Fraunhofer approximation was shown to fail in some angular regions, despite the fact that the mode considered (Ao) is the less anisotropic one among the three propagative modes existing below the first cut-off frequency. Here, the proposed model leads to predictions in excellent agreement with those computed using the full integration scheme. The two other comparisons shown on Fig. 8 (b, c) now concern the two other modes of higher anisotropy So and SHo. Even for the highest anisotropy (SHo), results predicted by the proposed model perfectly agree with those computed without approximation.
<figure><img src="image_10.png" /><figcaption>**Figure 8.** Out of plane displacement of the Ao mode (a), So mode (b), SHo mode (c) at a distance of 100mm in the case of an excitation by a circular source at the frequency of 300 kHz obtained by computing the full integration (solid line) and the approximate formula (+ points).</figcaption></figure>

Figure 8 (a) displays the same results as those shown on Fig. 5. In the previous results, the Fraunhofer approximation was shown to fail in some angular regions, despite the fact that the mode considered (Ao) is the less anisotropic one among the three propagative modes existing below the first cut-off frequency. Here, the proposed model leads to predictions in excellent agreement with those computed using the full integration scheme. The two other comparisons shown on Fig. 8 (b, c) now concern the two other modes of higher anisotropy So and SHo. Even for the highest anisotropy (SHo), results predicted by the proposed model perfectly agree with those computed without approximation.
## Rectangular-shaped transducers

Similar results are shown on Fig. 9 for a rectangular transducer (actually as square source of 9-mm-side). Directivities of the various modes are different from those radiated by the disk transducer. Again, results computed using the proposed model or using the full integration superimpose.

<figure><img src="image_11.png" /><figcaption>**Figure 9.** Out of plane displacement of the Ao mode (a), So mode (b), SHo mode (c) at a distance of 100mm in the case of an excitation by a square source at the frequency of 300 kHz obtained by computing the full integration (solid line) and the approximate formula (+ points).</figcaption></figure>

We can see on the Fig. 9 that the results, as for the circular shaped transducer are similar between the convolution and the integration along energy directions.

Note that in Fig. 8 and Fig. 9, results computed by means of the approximate formulation perfectly superimpose with those computed by means of the full integration, so much so that results shown as solid lines could not be distinguished. This is the reason for presenting the
## Rectangular-shaped transducers

Similar results are shown on Fig. 9 for a rectangular transducer (actually as square source of 9-mm-side). Directivities of the various modes are different from those radiated by the disk transducer. Again, results computed using the proposed model or using the full integration superimpose.

<figure><img src="image_11.png" /><figcaption>**Figure 9.** Out of plane displacement of the Ao mode (a), So mode (b), SHo mode (c) at a distance of 100mm in the case of an excitation by a square source at the frequency of 300 kHz obtained by computing the full integration (solid line) and the approximate formula (+ points).</figcaption></figure>

We can see on the Fig. 9 that the results, as for the circular shaped transducer are similar between the convolution and the integration along energy directions.

Note that in Fig. 8 and Fig. 9, results computed by means of the approximate formulation perfectly superimpose with those computed by means of the full integration, so much so that results shown as solid lines could not be distinguished. This is the reason for presenting the
former results with “+” signs and to limit the number of angles for which they are shown.

**5. Discussion**

The proposed approximate model leads to very accurate prediction of results obtained by means of the full integral over the source aperture. The approximation is quite limited as it only concerns the amplitude term along the intersection of a ray with the source surface, the phase term being exactly taken into account and analytically integrated along this intersection. The angular integration over ray directions covers the full surface seen from the computation point. With the present model, a 1D numerical integration is involved, this leading of course to faster computation than that required for the full 2D one. Moreover, the farther the computation point, the smaller the angular range of integration and the faster the numerical evaluation of the integral.

**6. Conclusion**

Fields of guided waves radiated by finite-sized sources must be evaluated to optimize testing configurations of plate-like structures. Fraunhofer approximation is known to lead to very accurate predictions of transducer diffraction effect for a minimum computational effort, but this is only true if the plate is isotropic. As soon as the plate structure is anisotropic, all the guided modes generated by a source behave anisotropically. For a given multi-layered composite plate, the anisotropy varies a lot from one mode to another. Some results have been shown in the first part of the paper proving that, even for the less anisotropic mode, Fraunhofer approximation fails at predicting accurately the field radiated. This makes it necessary to account for all the possible energy direction. The existence of caustics is an issue that was recently solved in the literature in the form of an approximate Green's function. The model proposed in the paper uses the modal description computed by the SAFE method and expressions of the Green's function which differ if the energy direction is close to or far from a caustic [4]. A change of variables in the convolution surface integral expressing the field was proposed to optimally compute this integral (at low cost), leading to expressions easy to implement and fast computed. The model has been shown to predict very accurately the radiated field associated to various modes, even for the most anisotropic ones.

**References**

[1] Raghavan A and Cesnik C E S 2005 *Smart Mater. Struct.* **14** 1448

[2] Taupin L, Lhémery A and Inquité G 2011 *J. Phys. : Conf. Ser.* **269** 012002

[3] Velichko A and Wilcox P D 2007 *J. Acoust. Soc. Am.* **121** 60

[4] Karmazin A, Kirillova E, Seeman W and Syromyatnikov P 2013 *Ultrasonics* **53** 283
former results with “+” signs and to limit the number of angles for which they are shown.

**5. Discussion**

The proposed approximate model leads to very accurate prediction of results obtained by means of the full integral over the source aperture. The approximation is quite limited as it only concerns the amplitude term along the intersection of a ray with the source surface, the phase term being exactly taken into account and analytically integrated along this intersection. The angular integration over ray directions covers the full surface seen from the computation point. With the present model, a 1D numerical integration is involved, this leading of course to faster computation than that required for the full 2D one. Moreover, the farther the computation point, the smaller the angular range of integration and the faster the numerical evaluation of the integral.

**6. Conclusion**

Fields of guided waves radiated by finite-sized sources must be evaluated to optimize testing configurations of plate-like structures. Fraunhofer approximation is known to lead to very accurate predictions of transducer diffraction effect for a minimum computational effort, but this is only true if the plate is isotropic. As soon as the plate structure is anisotropic, all the guided modes generated by a source behave anisotropically. For a given multi-layered composite plate, the anisotropy varies a lot from one mode to another. Some results have been shown in the first part of the paper proving that, even for the less anisotropic mode, Fraunhofer approximation fails at predicting accurately the field radiated. This makes it necessary to account for all the possible energy direction. The existence of caustics is an issue that was recently solved in the literature in the form of an approximate Green's function. The model proposed in the paper uses the modal description computed by the SAFE method and expressions of the Green's function which differ if the energy direction is close to or far from a caustic [4]. A change of variables in the convolution surface integral expressing the field was proposed to optimally compute this integral (at low cost), leading to expressions easy to implement and fast computed. The model has been shown to predict very accurately the radiated field associated to various modes, even for the most anisotropic ones.

**References**

[1] Raghavan A and Cesnik C E S 2005 *Smart Mater. Struct.* **14** 1448

[2] Taupin L, Lhémery A and Inquité G 2011 *J. Phys. : Conf. Ser.* **269** 012002

[3] Velichko A and Wilcox P D 2007 *J. Acoust. Soc. Am.* **121** 60

[4] Karmazin A, Kirillova E, Seeman W and Syromyatnikov P 2013 *Ultrasonics* **53** 283