Semi-Analytical Finite Element Method for calculating dispersion curves of a CFRP plate

. Composite Materials are widely used thanks to their mechanical characteristics. Non destructive evaluation requires knowing the exact dispersion curves to determine the propagative waves and to resolve the phase velocity of symmetrical and antisymmetrical modes. The aim of this paper is to plot the dispersion curves of a Carbon Fibre Reinforced Polymer plate using Semi Analytical Finite Element Method algorithm. This method combines the analytical expression of the displacement and the finite element method procedure. The resultant advantage is both simplicity and rapidity. The obtained results showed that the method accuracy depends on the elements number of meshing. To ensure good precision and fast computation of the method, the elements number and the order of the interpolation functions must be optimized.


Introduction
Today, many technics of Non-Destructive testing (NDT) are used to control the quality of mechanical pieces in several fields, especially in nuclear, medical and aeronautic applications. Ultrasonic Lamb waves are among these methods. They are dispersive and sensitive to small defects, so they are used to characterize some materials or to estimate the geometry and the position of defects in isotropic and anisotropic plates. To choose an adequate frequency to excite lamb waves or to postprocess the data obtained after experiments, it is required to know precisely the dispersion curves of studied plate. These curves define the variation of phase velocity and wave number as function of frequency.
Multiple researchers have discussed the dispersion curves with their applications in the NDT. Many technics are then proposed to plot the dispersion curves such as: the dichotomy method [1], the global matrix method [2], also the numerical spectral method [3]... The Bisection algorithm permits to compute the dispersion curves considered as analytical curves; however they present in general an inaccuracy when calculating the cutoff frequencies. So, another technic giving better calculation of the dispersion curves is proposed in this study.
Among the most used methods for calculating dispersion curves in case of isotropic plates and composites structures, one finds those based on the finite element method. But, these require important capacities of storage because of the need to mesh the whole structure.
In this paper, an alternative numerical method will be used. This method combines semi-analytical expansion of the displacement component along the direction of wave propagation and finite element interpolation for the orthogonal components of displacement. The guided propagating waves and their motion in the section perpendicular to the wave direction are computed by the finite elements' technic, and then integrated through explicit analytical development performed in the propagation axe. SAFE method allows a good prediction of the propagating modes in a plate since it requires just a discretization through the plate thickness. The time of calculation is therefore shorter than with a full discretization procedure. This method is then well adapted for the case of slender structures having complex cross section geometries and allows characterizing the properties of a wave propagating along the direction of slenderness.
The Semi Analytical-FE method was studied by Waas [4] for the calculation of surface waves in a multilayered soil. Afterwards, it has been applied by many authors to carry out a study of Lamb wave propagation in isotropic and in composite plates. Ahmad et al. [5] have presented the method to plot the dispersion curves of an aluminum plate and a laminated composite having the following sequence [0º/45º/90º/-45º]. Besides, Bartoli et al. [6] have developed the SAFEM to model the Lamb waves propagation in waveguides with arbitrary section. The authors have also used SAFE to study he viscoelastic material by taking into account damping parameter. Hayashi et al. [7] have simulated by SAFE the propagation of Lamb waves in tubes. In another work, this author [8] calculated and compared the dispersion curves obtained by SAFEM with the experimental curves computed by 2D-FFT. Mukdadi et al. [9] have studied the propagation of Lamb waves in composites structures. Predoi [10] has proposed a development of the SAFEM in case of periodic structures with infinite width. Wenbo et al. [11] have developed a formulation using the SAFEM and the PML (Perfect Match Layer) method to compute the dispersion curves of immersed pipeline in a fluid. Besides, Xing et al. [12] have presented a method to localize a defect in a rail and have plotted dispersion curves by the SAFEM.
Recently, Nissabouri et al. [13] have proposed a quantitative evaluation of dispersion curves calculated by SAFE for a simple case of orthotropic plate. In this paper, SAFE numerical method is used to compute the dispersion curves of a CFRP plate. In the following, a presentation of the algorithm used is performed. Then, the obtained curves are discussed based on a comparison with DISPERSE software results. Finally, the method accuracy is evaluated.

SAFE Method Algorithm
The algorithm was proposed firstly by Nissabouri et al. [13]. It was deduced based on the theoretical model and permits to resolve the following general equation: with K is the stiffness matrix, M the mass matrix, k a wave number,  an eigenvalue circular frequency and U an eigenmode.   1 presents the flowchart algorithm structure to plot dispersion curves. In the first step, we set the material parameters of the studied plate:  , E and C . After that, we define the geometry of the plate by the thickness denoted e . In the second step, we mesh the cross section of the plate. For that, the elements number n and the order of interpolation functions N i should be fixed. In the third step, we calculate the elementary matrix e i K by Gaussian quadrature. Then, we assemble the Elementary matrices into a global matrix. In the final step, for each value of wave number k , we solve the eigenvalue problem in order to calculate the eigenvalues which corresponds to the frequencies and their eigenvectors which corresponds to the displacements.

Quantitative evaluation of SAFEM
In order to execute the algorithm presented in Fig. 1, a Matlab program was developed. This was adapted to an orthotropic plate. The properties of considered plate are listed in Table 1. To assess accuracy of the proposed numerical method, the dispersion curves computed by SAFEM are compared to those provided by DISPERSE software which is based on the Global Matrix Method.

Minimal number of elements
For the considered orthotropic plate having thickness of 1 mm , a meshing consisting of 6 elements is enough in order to calculate the first six propagating modes in the frequency range up to 5 MHz . For higher frequencies than 5 MHz , the number of elements should be increased [10]. The following condition is be satisfied: where l is the element length, 2 / T T V     is the wavelength of transverse waves propagating at the speed T V , and 4   for quadratic interpolation functions [14]. Using 2 f    and n elements through the thickness, one obtains from Eq. (2) Numerical approximation of integrals by Gaussian quadrature to compute elementary matrices.

Assemble elementary matrices into global Matrix.
Compute the eigenvalue ωi for each value of the wavenumber k.
Parameters of the considered model. . So for a specified product frequency thickness fe , the minimum number n of elements that guarantee a good precision of the method is: In the present case, 2100 / T V m s  and min 9.45 n  from Eq. (4). So, 10 n  should be used. Fig. 2 gives a comparison of the first five dispersion curves calculated by the SAFEM, for a mesh consisting of 7 elements through the plate thickness, with those provided by DISPERSE software. Fig. 3 gives these branches when using 10 n  elements.     Fig. 3 show that the precision of SAFEM depends on the elements number. The method enables to calculate the first modes with a good precision if the number of elements used is sufficient. This can be fixed according to Eq. (4). Table 2 show that the error on 0 S branch is reduced when using 10 n  .

Effect of elements number on the running time
If n the elements number is more than 10, a good precision will be reached. However, this will cause an increase of running time. Table 3 shows the changing time regarding the number of elements for 10 n  . The machine used has a processor of 2.6 GHz and RAM of 4Go . The step increment used for the wave number k was fixed to 0.5 . Table 3 shows that the variation of running time is linear with a number of elements lesser than 4. For a number of elements n greater than 4, the method will be slower. This is related to the increasing size of the elementary matrices.  Fig. 4 gives the calculated curves of phase velocity for the 1mm CFRP plate for frequencies lesser than 15 MHz . Fig. 5 and Fig. 6 give respectively the calculated SAFE group velocity curves and those of DISPERSE software for 5 f MHz  . Fig. 4. Phase velocity curves obtained by SAFE method for 1mm orthotropic CFRP plate and 10 n  .    Fig. 6 show that the calculated group velocity branches by the SAFEM and DISPERSE software are identical. This enables to asses accuracy of the proposed method.

Transversal and longitudinal displacements
The eigenvectors of the Eq. (1) are the longitudinal displacements x u and transverse displacements z u . These can be plotted as function of the plate thickness for each mode.

Conclusions
In this work we validated, by comparison with DISPERSE software the accuracy of the SAFE method that we have proposed. The results have revealed that an increase in frequency cause an increase in error of the dispersion branches. The method was found to be slow for greater number of elements. To guarantee good precision and quality of the method, we used an optimal number of elements n for meshing the cross section. The SAFEM allows to compute propagating modes with good accuracy and reduced numerical cost as it does not necessitate the meshing of the whole model.