# NURBS-enriched semi-analytical finite element method (SAFE) for calculation of wave dispersion in heterogeneous waveguides

Fakhraddin Seyfaddini, Vu-Hieu Nguyen

► To cite this version:

Fakhraddin Seyfaddini, Vu-Hieu Nguyen. NURBS-enriched semi-analytical finite element method (SAFE) for calculation of wave dispersion in heterogeneous waveguides. 24e Congrès Français de Mécanique, Aug 2019, Brest, France. hal-02471811

HAL Id: hal-02471811

https://hal.science/hal-02471811

Submitted on 16 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 Attribution 4.0 International License
# NURBS-enriched semi-analytical finite element method (SAFE) for calculation of wave dispersion in heterogeneous waveguides

Fakhraddin Seyfaddini, Vu-Hieu Nguyen

Université Paris-Est, Laboratoire Modélisation et Simulation Multi Echelle, MSME UMR 8208
CNRS,94010 Créteil CEDEX, France; fakhraddin.seyfaddini@u-pec.fr; vu-hieu.nguyen@u-pec.fr

## Abstract

In this work a semi analytical finite element formulation using Non-Uniform Rational B-Splines (NURBS) basis function is presented for modeling the dispersed wave propagation in heterogeneous waveguides. Convergence studies are carried out by considering a 2D constant-thickness isotropic elastic plate. The results were compared with the ones obtained by analytical methods and by conventional SAFE method. For all of cases, the dispersion curves evaluated by using enriched-NURBS basis have a significant better precision than using conventional Lagrangian elements (for the same number of degrees of freedom).

**Keywords :** Isogeometric analysis (IGA), NURBS basis, Semi analytical finite element (SAFE), Guided waves, Dispersion, Elastic plate

## 1 Introduction

Ultrasonic guided wave (UGW) technologies are powerful nondestructive testing techniques to characterize near surface materials and evaluate integrity of materials or structures. Due to the presence of boundaries, the guided waves show a strong dispersive behavior, i.e. the phase velocity and attenuations vary with frequency-content of wave package. One of major issue in this problem is how to calculate efficiently the dispersion curves of all modes in the studied frequency range which will serve later to the inversion task. The Semi-Finite Element Method (SAFE) is one of most popular technique for computing the dispersion of guided waves in structures thanks to its effectiveness in studying functionally graded or coupled fluid/solid multilayer plates [1, 2]. However, at very high frequency, using conventional high-order Lagrangian interpolation function does not allow to improve the situation because it may lead to numerical issues when solving eigenvalue problems.

The Isogeometric Analysis (IGA) is a recently introduced concept that uses most commonly the NURBS basis functions as a powerful tool from the Computer Aided Design (CAD) to represent not only the complex geometries but also to construct shape functions for finite element analysis. In particular, the use
of NURBS basis functions for spectrum and dissipation analyses shows that this method is more accurate compared to the classical finite element analysis (FEA) for a fixed number of degrees of freedom [3, 4]. The objective of this paper is to study the effectiveness of using NURBS basis functions in the context of SAFE method for analyzing the wave propagation in homogeneous waveguides.

## 2 Governing Equation

A two dimensional solid layer of thickness $h$ having infinite extent along the horizontal direction $x_1$ represented by the domain $\Omega = [-h, 0] \times [-\infty, \infty]$ as shown in Fig. 1. The equations of motion in $\Omega$ and the constitutive relation describing elastic behavior of the solid are given by :

$$ \rho \ddot{u} - \mathbb{L}^T \sigma = 0, \qquad (1) $$

$$ \sigma = C\epsilon, \quad \epsilon = \mathbb{L}u \qquad (2) $$

where the $\rho$ is the mass density; $\sigma = \{\sigma_{11}, \sigma_{22}, \sigma_{12}\}^T$ and $\epsilon = \{\epsilon_{11}, \epsilon_{22}, 2\epsilon_{12}\}^T$ are the vectors containing the components of the stress and strain tensors, respectively; $C$ is the matrix containing the components of the elasticity tensor using Voigt's notation; the operator $\mathbb{L}$ is defined by :

$$ \mathbb{L} = \mathbf{L}_1 \partial_1 + \mathbf{L}_2 \partial_2, \quad \mathbf{L}_1 = \begin{bmatrix} 1 & 0 \\ 0 & 0 \\ 0 & 1 \end{bmatrix}, \quad \mathbf{L}_2 = \begin{bmatrix} 0 & 0 \\ 0 & 1 \\ 1 & 0 \end{bmatrix} \qquad (3) $$

Free-surfaces boundary conditions interfaces are assumed at the upper and lower surfaces of the plate, implying $\mathbf{L}_2^T \sigma = \{\sigma_{12}, \sigma_{22}\}^T = 0$.

By noting that all material properties are homogeneous in $x_1$-direction ($\rho = \rho(x_2)$, $C = C(x_2)$), we look for solution of harmonic waves propagating along the axial $x_1$ direction in the following form :

$$ \mathbf{u}(x_1, x_2, t) = \tilde{\mathbf{u}}(x_2)e^{i(k_1 x_1 - \omega t)} \qquad (4) $$

where $i^2 = -1$, $\omega$ is the angular frequency and $k_1$ is the wave number in the $x_1$ direction, $\tilde{\mathbf{u}}(x_2)$ represents the amplitudes of the displacement vector. By substituting Eq. (4) into Eq. (2), we obtain the following equation in the frequency-wavenumber domain for each values of $(\omega, k_1)$:

$$ (-\omega^2 \mathbf{A}_1 + k_1^2 \mathbf{A}_2) \tilde{\mathbf{u}} - ik_1 (\mathbf{A}_3 + \mathbf{A}_3^T) \partial_2 \tilde{\mathbf{u}} - \mathbf{A}_4 \partial_2^2 \tilde{\mathbf{u}} = 0 \qquad (5) $$

in which the matrices $\mathbf{A}_1, \mathbf{A}_2, \mathbf{A}_3, \mathbf{A}_4$ are defined by :

$$ \mathbf{A}_1 = \rho \mathbf{I}, \quad \mathbf{A}_2 = \mathbf{L}_1^T C \mathbf{L}_1, \quad \mathbf{A}_3 = \mathbf{L}_2^T C \mathbf{L}_1, \quad \mathbf{A}_4 = \mathbf{L}_2^T C \mathbf{L}_2 \qquad (6) $$

The weak form or variational form of the Eq. (5) is obtained by introducing a test function $\delta \tilde{\mathbf{u}}^*$ which is the conjugate transpose of $\delta \tilde{\mathbf{u}}$ and integrate by parts. By using the free-surfaces boundary conditions.
The resulting weak form of the problem is : Find $\tilde{u}(x_2)$ such that :

$$ \int_{0}^{h} \delta \tilde{u}^{*} (-\omega^2 \mathbf{A}_1 + k_1^2 \mathbf{A}_2 - ik_1 \mathbf{A}_3^T \partial_2) \tilde{u} dx_2 + \int_{0}^{h} \partial_2 \delta \tilde{u}^{*} (ik_1 \mathbf{A}_3 + \mathbf{A}_4 \partial_2) \tilde{u} dx_2 = 0, \quad \forall \delta \tilde{u} \in \mathcal{C}^{ad} \qquad (7) $$

### 3 Discretization by NURBS

A B-spline basis functions of order $p$ is determined in a parameter domain $\hat{\Omega}$ using a knot vector $\Xi = \{\xi_1, \xi_2, \xi_3, ..., \xi_{i+p+1}\}$, and polynomial order $p = 1, 2, 3, ...$ as :

$$ N_{i,p}(\xi) = \frac{\xi - \xi_i}{\xi_{i+p} - \xi_i} N_{i,p-1}(\xi) + \frac{\xi_{i+p+1} - \xi}{\xi_{i+p+1} - \xi_{i+1}} N_{i+1,p-1}(\xi) \qquad (8) $$

where $\xi_i \in \mathbb{R}$ is the $i^{th}$ knot, $i$ is the knot index ($i = 1, 2, ..., n + p + 1$), and $n \in \mathbb{N}$ is the number of basis functions used to construct the B-spline curve. Non-uniform rational B-spline (NURBS) basis functions are built from the B-spline functions by multiplying weighting functions $w_i$ for each basis

$$ R_{i,p}(\xi) = \frac{N_{i,p}(\xi)w_i}{\sum_{j=1}^{n} N_{j,p}(\xi)w_j} \qquad (9) $$

Given $n$ basis functions $R_{i,p}(\xi)$, and corresponding set of control points $P_i \in \mathbb{R}^d (i = 1, 2, 3, ..., n)$ the NURBS geometry in 1D is defined as a mapping from parameter domain to the physical domain as : $C(\xi) = \sum_{i=1}^{n} R_{i,p}(\xi) P_i$. Fig. 1 (right) shows an example of the discretization of plate thickness and its corresponding basis functions. According to the isogeometric concept we use the same basis functions as for geometric representation in order to approximate the solution fields.

<figure><img src="image_2.png" /><figcaption>FIGURE 1 – Geometry description of infinite plate layer and the discretization in the direction of plate thickness (left). The cubic basis functions over the corresponding knot vector (right)</figcaption></figure>

The functions defined on each element are calculated using the isoparametric NURBS basis as $\delta \tilde{u}(\xi) = R_e \delta U_e$ and $\tilde{u}(\xi) = R_e U_e$. Where $R_e$ are the NURBS basis functions, $U_e$ and $\delta U_e$ are the vector of control variables of $u_e$ and $\delta u_e$ respectively. Replacing these functions into the weak formulation (7) and assembling the elementary matrices lead to a quadratic system of eigenvalue equations :

$$ (-\omega^2 \mathbf{K}_1 + \mathbf{K}_4 + k_1^2 \mathbf{K}_2 + ik_1 \mathbf{K}_3) U = 0, \qquad (10) $$
where $U$ is the vector of amplitudes of global displacement field; the matrices $K_1, K_2, K_3$ and $K_4$ are defined by

$$K_1 = \bigcup_e \int_{\Omega_e} (R_e)^T A_1(R_e) J_a J_{\tilde{\xi}} d\tilde{\xi}, \quad K_2 = \bigcup_e \int_{\Omega_e} (R_e)^T A_2^e(R_e) J_a J_{\tilde{\xi}} d\tilde{\xi}, \quad (11)$$

$$K_3 = \bigcup_e \int_{\Omega_e} 2\{(B_e)^T A_3^e(R_e)\}_a J_a J_{\tilde{\xi}} d\tilde{\xi}, \quad K_4 = \bigcup_e \int_{\Omega_e} (B_e)^T A_4^e(B_e) J_a J_{\tilde{\xi}} d\tilde{\xi}, \quad (12)$$

where $J_a$ refers to the jacobian transformation from physical domain to parameter domain ($x, \xi$) = $\frac{\partial x}{\partial \xi} = J_a$); $J_{\tilde{\xi}}$ refers to the jacobian transformation between parent domain and parametric domain ($\frac{d\xi}{d\tilde{\xi}} = J_{\tilde{\xi}}$). In this study, Gauss quadrature rule has been used for computing the integrations over the elements. By fixing the angular frequency $\omega$ and solving (Eq. 10), one may determine the eigenvalues $k_1$ and their associated eigenvectors $V(\omega, k_1)$, which represent the wavenumber and the wave structure of guided modes.

The frequency-dependent phase velocity ($C_p$) and attenuation (att) of a guided mode are obtained from $k_1$ using the following relationships

$$c_{ph} = \frac{\omega}{\text{Re}[k_1]} (\text{m s}^{-1}), \quad \text{att} = \text{Im}[k_1] (\text{Np m}^{-1}). \quad (13)$$

where Re() and Im() are the real and imaginary parts of a complex quantity.

## 4 Numerical results

A numerical example using the proposed SAFE based IGA formulation is presented in this section to demonstrate its validity and accuracy for the analysis of guided wave in elastic free plate problem. First, an isotropic homogeneous plates is considered and the validation is performed by comparing our results to the numerical ones obtained using conventional SAFE method and to analytical ones obtained using the software DISPERSE. Then the convergence analysis is carried out.

The IGA-based SAFE formulation is validated for the case of a 4mm thickness of homogeneous aluminum plate which is assumed to have a isotropic elastic properties with the density $\rho = 2700 \text{kg/m}^3$, longitudinal wave velocity $c_L = 2344 \text{m/s}$ and shear wave velocity $c_T = 953 \text{m/s}$. For IGA-based SAFE analysis, we used NURBS basis functions over uniform knot vector with a fixed number of degree of freedom ($N_{dof} = 26$), regardless of their polynomial degree.

The dispersion of phase velocities calculated with quadratic NURBS basis functions is shown in Fig. 2(a). By comparing with the results of conventional SAFE using 6 quadratic elements ($n_{el} = 6$ and ($N_{dof} = 26$)), the IGA-based SAFE formulation using quadratic NURBS ($n_{el} = 11$ and the same degree of freedom ($N_{dof} = 26$)) shows better results, particularly for the higher modes (e.g A7, S7) at high frequencies (greater than 2.5 MHz). The results using a cubic basis functions ($p = 3$)) are presented in Fig. 2(b) and are compared with the conventional FE results using the same order of basis functions.
<figure><img src="image_3.png" /><figcaption>FIGURE 2 – Dispersion curves of homogeneous isotropic aluminum plate obtained using SAFE based IGA (red marker), conventional SAFE (blue marker) and Disperse software (green dashed line) for quadratic basis function and same degree of freedom ($N_{dof} = 26$)</figcaption></figure>

Similar statement on the effectiveness of using IGA basis may be done. Thus the IGA method allows us to increase the degree of basis functions and obtain the better agreement of the results compare with analytic solution.

**Convergence Analysis** In order to carry out a convergence analysis of the proposed method, we introduce a parameter $L_{mode}^p = \text{Re}[k_1]$ which is the corresponding value of wavenumber at a specific mode. Then the relative error of the numerical methods at a point $x$ is estimated by

$$e = \left| \frac{L_{mode}^p(x) - L^{ref}}{L^{ref}} \right| \qquad (14)$$

where $L^{ref}$ is the reference analytical value obtained from the DISPERSE. The normalized dispersion errors of the first mode ($A_0$) with respect to the polynomial order are depicted in Fig. 3 for both numerical methods. It was shown that these errors decrease faster for the NURBS order $p$.

# 5 Conclusion

Due to the strong frequency dependency of dispersion properties, determination of guided wave in plate at high frequency range may require high computational cost. The dispersion curves calculated with the IGA-based SAFE shows a better accuracy results compared to conventional SAFE method. The computational time also reduced using the proposed method. The convergence analysis shows that increasing the order of NURBS basis function leads to a much faster convergence rate in comparing with a similar procedure allying to Lagrange polynomials basis functions. For example, increasing the order of NURBS basis functions from $p = 3$ to $p = 4$ leads to an improvement of the error $e = 0.03\%$ (for $N_{dof} = 26$) whereas when using Lagrange polynomials this improvement is $e = 0.27\%$. In addition, it has been
<figure><img src="image_4.png" /><figcaption>FIGURE 3 – Convergence analysis of isotropic plate for IGA-based SAFE and conventional SAFE refering to analytical results obtained by using DISPERSE</figcaption></figure>

shown that the proposed approach gives more accuracy results for studying multilayer anisotropic plates (data not shown). This model can be developed to deal with other problems such as studying guided waves in 3D structures or in poroelastic materials [5].

Références

[1] V.-H. Nguyen and S. Naili. Ultrasonic wave propagation in viscoelastic cortical bone plate coupled with fluids : a spectral finite element study. *Computer methods in biomechanics and biomedical engineering*, Vol. **16**, pp. 963-974, 2013.

[2] V.-H. Nguyen, T. N. Tran, M. D. Sacchi, S. Naili, and L. H. Le. Computing dispersion curves of elastic/viscoelastic transversely-isotropic bone plates coupled with soft tissue and marrow using semi-analytical finite element (safe) method. *Computers in biology and medicine*, Vol. **87**, pp. 371-381, 2017.

[3] J. Cottrell, A. Reali, Y. Bazilevs, and T. Hughes. Isogeometric analysis of structural vibrations. *Computer Methods in Applied Mechanics and Engineering*, Vol. **195**, pp. 5257-5296, 2006.

[4] T. J. Hughes, J. A. Evans, and A. Reali. Finite element and nurbs approximations of eigenvalue, boundary-value, and initial-value problems. *Computer Methods in Applied Mechanics and Engineering*, Vol. **272**, pp. 290-320, 2014.

[5] V.-H. Nguyen and S. Naili. Simulation of ultrasonic wave propagation in anisotropic poroelastic bone plate using hybrid spectral/finite element method. *International journal for numerical methods in biomedical engineering*, Vol. **28**, pp. 861-876, 2012.