Study on Meshfree Hermite Radial Point Interpolation Method for Flexural Wave Propagation Modeling and Damage Quantification

Abstract

This paper investigates the numerical modeling of the flexural wave propagation in Euler-Bernoulli beams using the Hermite-type radial point interpolation method (HRPIM) under the damage quantification approach. HRPIM employs radial basis functions (RBFs) and their derivatives for shape function construction as a meshfree technique. The performance of Multiquadric(MQ) RBF to the assessment of the reflection ratio was evaluated. HRPIM signals were compared with the theoretical and finite element responses. Results represent that MQ is a suitable RBF for HRPIM and wave propagation. However, the range of the proper shape parameters is notable. The number of field nodes is the main parameter for accurate wave propagation modeling using HRPIM. The size of support domain should be less thanan upper bound in order to prevent high error. With regard to the number of quadrature points, providing the minimum numbers of points are adequate for the stable solution, but the existence of more points in damage region does not leads to necessarily the accurate responses. It is concluded that the pure HRPIM, without any polynomial terms, is acceptable but considering a few terms will improve the accuracy; even though more terms make the problem unstable and inaccurate.

Keywords:
Damage quantification; Wave propagation; Meshfree method; Radial basis function; Hermite Radial point interpolation

1 INTRODUCTION

Wave propagation can be used to identify the damage of small defects. The identification methods are based on high frequency excitation, wave interaction with the damage and scattering. A variety of numerical techniques were proposed to model wave propagation in time, frequency or time-frequency domains. The frequent scheme is finite element method (FEM) that its capability in this field was illustrated in references (Moser et al. 1999, Chakraborty et al. 2002). The high frequency loads are accompanying with shorter wavelength; therefore a large number of meshes and degrees of freedom are needed for the exact consideration of the wave energy. This limitation leads to the high cost computation which makesit necessary todevelop the alternative techniques. Modified FEM (Keramat and Ahmadi 2012), boundary element method (Zhao and Rose 2003), mass-spring lattice model (MSLM) (Yim and Sohn 2000) and several spectral FEM (Mitra and Gopalakrishnan 2005, Kudela et al. 2007) have been reported as novel algorithms for wave propagation and damage identification. However mesh-based methods are not well suited to treat problems with strong inhomogeneity, mesh distortion and discontinuities that do not align with element edges especially for dynamic problems such as crack propagation (Chen et al. 2006, Nguyen et al. 2008). To overcome these drawbacks different meshfree schemes were introduced (Li and Liu 2007, Liu 2010, Li and Mulay 2013). Wen (2010) utilized the meshfree local Petrov-Galerkin (MLPG) method to solve wave propagation problem in three-dimensional poroelastic solids. Gao et al. (2007) proposed a new MLPG method to analyze stress-wave propagation in anisotropic and cracked media. Das and Kundu (2009) simulated the ultrasonic wave field models for layered half-spaces with and without an internal crack using the meshfree semi-analytical distributed point source method. Zhang and Batra (2004) applied a modified smoothed particle hydrodynamics to study the 1D wave propagation and 2D transient heat conduction problems.

Among the various meshfree techniques, radial basis functions (RBFs) or, more generally, conditionally positive definite kernels have been achieved significant attention in approximation and the solution of partial differential equations (PDEs). RBFs are powerful tools to interpolate the multivariate scattered data. They are insensitive to spatial dimension, differentiable, computationally simple and stable in the interpolation (Buhmann 2003, Liew and Chen 2004, Wendland 2010). Hardy (1971) used multi-quadric (MQ) RBF to solve the approximation of a topographic surface. Franke (1982) explored MQ and some other scattered data interpolation methods to evaluation the timing, storage, accuracy and the ease of implementation. Lazzaro et al. (2002) proposed an efficient method for the multivariate interpolation of very large scattered data sets using RBF and Shepard's method. Kansa (1990a, b) was the pioneer who used MQ for partial derivative estimates and solving of elliptic, parabolic, and hyperbolic PDEs. Kansa's method searches an approximate solution in the function space using collocation technique by a set of RBFs. His method was extended to solve various PDEs which can be referred in (Hon and Wu 2000, Fedoseyev et al. 2002, Tiago and Leitão 2006, Dehghan and Shokri 2009).

The Primary disadvantage of the traditional RBF approach is that the coefficient matrices obtained from discretization scheme are fully populated, which limit its application in large scale practical problems (Wang and Liu 2000). One effective approach is the combination of radial and polynomial basis functions which was introduced by Wang and Liu (2002) as radial point interpolation method (RPIM). Point interpolation method (Liu and Gu 2001) satisfies Kronocker delta function property of shape functions and makes the implementation of essential boundary conditions easier but the matrix may be singular. Applying mentioned combination can guarantee the stability of shape functions and the improvement of singularity for arbitrarily and irregular node distribution (Liu 2010).

The other challenge in RPIM is the choice of the proper shape parameters of RBFs that has important effect on the accuracy of results (Wang and Liu 2002). Rippa (1999) proposed an optimization algorithm for selection of a good value for shape parameters in MQ, inverse MQ and Gaussian interpolation in data fitting. Kansa (1990a) proposed a scheme that allowed the shape parameter of MQ to vary with basis functions. Wang and Liu (2002) proposed optimal shape parameters for MQ and reciprocal MQ in 2D solid mechanics. Fornberg and Wright (2004) used a technique called Contour-Padé algorithm for stable computation of RBF approximation methods. The algorithm avoids working directly with the associated ill-conditioned linear systems for all values of the shape parameters. Bozkurt et al. (2013) discussed the effects of shape parameters on the RPIM accuracy in 2D geometrically nonlinear problems. Kanber et al. (2013) studied RPIM algorithm for the solution of 2D elastoplastic problems using the multiquadric RBF. The convergence rates after yielding were investigated for perfectly plastic and hardening cases with various shape parameters. Both regular and irregular node distributions were used by considering the standard Gauss integration and nodal integration schemes.

The RPIM within the framework of weak form applied in various engineering problems: solid mechanics (Liew and Chen 2004, Liu et al. 2005, Liu et al. 2006), plate structures (Liew and Chen 2004), geometrical nonlinearity (Bozkurt et al. 2013), composite materials (Liu et al. 2008) and consolidation process (Wang et al. 2002). Ghaffarzadeh and mansouri (2013) studied the longitudinal wave propagation based on RPIM in rods. The capability of RPIM and the effect of MQ shape parameters were explored for detection of additional mass which is located in the length of the rod.

In some problems, derivatives of field functions should be considered as independent variables other than the function values. In this field, Wu (1992) introduced the Hermite RBF approach based on Hermite-Birkhoff interpolation for scattered multidimensional data. In solid mechanics, Liu et al. (2006) proposed Hermite radial point interpolation method (HRPIM) to approximate the deflection to solve the Kirchhoff plate problem using Galerkin weak form method. In their method, static analysis of different geometric shapes and boundary conditions verified the efficiency, accuracy and robustness of the method. Cui et al. (2009) developed a smoothed HRPIM using gradient smoothing operation for thin plate analysis. The approximation of deflection and rotation variables made the proposed method very effective in enforcing the essential boundary conditions even for extremely irregular background cells.

This study is conducted to illustrate the performance of Hermite RPIM weak form for the modeling of elastic wave propagation and damage quantification. The model is Euler-Bernoulli beam which has the predefined step damage as a crack in the specific position. Damage quantification algorithm is based on the analytical solution of governed differential equation that yields the relation between the reflection ratio and depth of damage. The dispersive behavior of flexural waves is considered in solution. The reflection ratios are calculated from captured signals using Hilbert transform to achieve the energy envelope of signals. The field nodes density, the size of support domain, the values of shape parameters in MQ, the number of polynomials and integration points are the most important parameters which are investigates in this study.

2 PROBLEM FORMULATIONS

2.1 Hermite Radial Point Interpolation (HRPIM)

RBFs are the functions whose values depend only on the distance between the center point and the scattered nodes. In order to avoid the singularity of the polynomial point interpolation, RBFs are used to develop the radial point interpolation method (RPIM). Adding polynomial terms up to the linear order can ensure to pass the standard patch tests and for considering derivatives of field functions as independent variables, RPIM augments with derivatives of RBF. Therefore, the function u(x) is approximated by a linear combination of RBFs, its derivatives and polynomials as following for 1D beam model which has two degrees of freedom at each node along x-axis:

(1)

where R_i (x) and R_i,x (x) are RBF and its derivative, respectively and P_j (x) is the polynomial basis functions. a_i, a_i ^x and b_j are the unknown coefficients for R_i (x), R_i,x (x) and P_j (x) respectively. n is the number of nodes in the support domain of point x = [x, y, z]^T and m is the number of monomial terms. If the derivative of field functions are equal to the derivative of approximated function, i.e. ∂u^h (x)/∂x = ∂u(x)/∂x, the unknown coefficients of Eq. (1) can be determined by enforcing the field function (nodal displacement) and its derivatives (nodal slope) at n nodes in the support domain. This leads to the 2n linear equations in 1D beam-like model as matrix form

(2)

where

(3)

and the matrices are defined as

(4)

In addition, the polynomial terms have to satisfy the following constraint to obtain a unique solution

(5)

Arranging Eqs. (2) and (5) together leads to the following set of equations

(6)

(7)

where U is the vector of nodal values and G is the generalized moment matrix. Substituting Eq. (7) into Eq. (1) leads to

(8)

Therefore the Hermite-RPIM shape function can be introduced as

(9)

Finally, the approximated function can be obtained by

(10)

The shape functions ϕ_i and ϕ_i^x (i = 1, 2, ..., n) are defined by

(11a)

(11b)

The derivatives of shape functions can be easily expressed for the first order as

(12a)

(12b)

Because of the existence of G^-1, there is no singularity in the process of shape function construction using Hermit-RPIM. The shape function and derivatives are stable and satisfy the Kronocker delta function property and partition of unity.

In this study, the MQ RBF is used to interpolate which has the following form:

(13)

where r is Euclidian distance, q and α_c are two shape parameters and d_c is the characteristic length that is usually the average nodal spacing for all the n nodes in the support domain.

2.2 Wave Propagation Analysis

The equations of motion in an elasto-dynamic problem are given in matrix form as

(14)

(15)

where L is the differential operator, σ and ε are the stress and strain, b is the body force; ϱ is the mass density, η_c is the damping parameter and D is the matrix of material constants. ü, and u are acceleration, velocity and displacement, respectively, while ∇u is the displacement gradient matrix. The natural and essential boundary conditions are

(16)

where ū are the prescribed displacement on the essential boundary Γ_u and _t. Γ = Γ_t ∪ Γ_u presents the complementary parts of the whole boundary in which n is the unit outward normal to this boundary. The Galerkin weak form of dynamic problem in Eqs. (14-16) over the global problem domain is expressed as are tractions on the natural boundary Γ

(17)

By using Eqs. (15) and (8), the discretization of Eq. (17) yields to:

(18)

where M, C, K, F and U are the global matrices of mass, damping, stiffness, time dependent excitation signal and nodal displacement, respectively. The components of matrices are derived in each support domain as

(19a)

(19b)

(19c)

(19d)

In Eq. (19), B is the strain-displacement matrix. For a special case such as the Euler-Bernoulli beam with two dimensional displacement field described by v(x) and u(x) in x and y coordinate system which are named here as axial and transverse displacement components respectively, the axial strain ε(x) can be expressed as

(20)

where, y is the transversal distance from neural axis of the cross section of the beam. Thus, for the beam model the matrix B can be constructed from the strain componentas B = -y∂²/∂x ²{Φ}.

To determine stiffness, mass and force matrices achieved by the meshfree method it is important to implement a suitable numerical scheme to carry out the integrals of Eqs (19). The well-known Gauss quadrature numerical integration approach can be utilized. From the Gauss quadrature integration the integral of a function can be evaluated in the specified Gauss pointswith their related weight coefficients. The appropriate number of Gauss points to be used depends upon the complexity of the integrand. A background cells - independent from field nodes - is necessary for numerical integrations of Eqs. (19). Using Gauss quadrature approach to perform the integrations numerically over these cells in general form leads to

(21)

where n_g is the number of Gauss points in each background cell, w_i is the Gauss weighting factor for the i ^th Gauss point at x_Qi (Figure 1) and J^D_ik is Jacobian matrix for the area integration of the background cell k, at which Gauss point x_Qi is located.

Figure 1:
A straight beam represented by a set of field nodes along its neutral axis and background cells.

2.3 Time Integration Scheme

The central difference time integration scheme is used to solve Eq. (18) as an explicit method

(22)

where t and Δt denotes time and time step of integration, respectively. This scheme is stable if Δt ≤ Δt_cr = 1/ω_max, where ω_max is the largest frequency of the system (Chen et al. 2006). This condition guarantees the accurate propagation of wave energy along the structure. The Zero initial conditions, U = 0 and = 0 at t = 0 are assumed to be implemented as initial displacement and velocity.

2.4 Damage Quantification Based on Wave Propagation for Notched Beams

This section presents how damage extent can be estimated using wave propagation by considering a uniform Euler-Bernoulli beam with step damage perpendicular to the beam axis with depth H_d . Figure 2 demonstrates the segments of the prescribed rectangular cracked beam. By neglecting the twisting effect and the assumption of the beam axis being the x-axis, the equation of motion for Euler-Bernoulli beam can be expressed in terms of the lateral displacement u(x, t) and the rotation angle ϑ(x, t) (Doyle 1997, Su et al. 2003)

(23a)

(23b)

(23c)

where, E, I, A, ρ and κ denote the elastic modulus, second moment of area, cross-section area, mass density and shear correction parameter, respectively. By considering of the dispersion relations:

(24)

the eigen roots (wavenumber) for Eq. (23) can be yielded as

(25)

Figure 2:
Damage model for quantification based on wave propagation.

Four eigen roots states that four propagating wave components are produced synchronously in the Euler-Bernoulli beam. There are four possibilities in all, allowing the complete solution to be written as

(26)

where each term corresponds to one of the wave modes depending on their propagating direction and properties: the positively propagating wave (A ₁), the negatively propagating wave (A ₂), the positively growing ringingwave (A ₃) and the negatively exponentially decaying wave (A ₄). The coefficient of each term represents the amplitude of each wave component which may be a complex problem.

Interaction between the incident wave and the damage yields to two propagating equations results from Eq. (26) before and after the damage position (referring to Figure 2), denoted with u _- and u ₊, respectively

(27a)

(27b)

where the three terms in Eq. (27a) correspond to the incident wave (A ₁), reflected propagating wave (A ₂) and reflected exponentially decaying wave (A ₄), respectively, while the two terms in Eq. (27b) are the damage-induced positively propagating wave (B ₁) and growing attenuating wave (B ₂), respectively. The wave components A ₄ and B ₂ disappear very quickly due to exponentially decaying. Applying the following continuity and equilibrium conditions yields to the amplitudes of wave components (Kim and Kim 2000)

(28)

(29)

The reflection ratio can be found as the relation between the magnitude A ₂ of the reflected wave and the magnitude A ₁ of the incident wave

(30a)

(30b)

where l and H are the width and depth of the cross-section of the beam, respectively; while H _d represents the depth of the damage. The reflection ratio is a complex number which indicates that there is a phase shifting relative to the incident wave (Graff 1991). Eq. (30) implies that the depth of damage, H _d, can be determined with a given reflection ratio. Figure 3 represents the absolute value of the reflection ratio and the damage degree (H _d/H).

Figure 3:
Theoretical relationship between reflection ratio and damage degree.

2.5 Computation of Reflection Ratio Using Hilbert Transform

The Hilbert transform (HT) is useful to the analysis of an inverse problem such as damage assessment. The application of HT to the signal provides some additional information about signal content. One of these information is the envelope function which contains important information about the energy of the signal. For an arbitrary time series X(t), the HT is defined as (Sun et al. 2009)

(31)

The Hilbert transform performs a 90º phase shift to construct the so-called analytic signal Z(t) as

(32)

whose real part is the original signal X(t) itself and imaginary part is its HT. The envelope e(t) and instantaneous phase φ(t) can be expressed by

(33a)

(33b)

Therefore, based on HT results, the reflection ratio can be obtained by driving the amplitude of the corresponding energy peak by that of the incident energy peak in the time domain.

3 NUMERICAL ANALYSIS

3.1 Numerical Model

Figure 4 represents the schematic meshfree model of a damaged Euler-Bernoulli beam according to the section 2.4 which will be studied in this paper. It is to be noted the formulation of section 2.4 were presented for one-sided step damage, however since this study investigates the 1D beam model as a straight line the mentioned formulation governs the scheme such as Figure 4 where the damage is located in two sides of the x-axis. The field nodes are distributed regularly along the x-axis. In addition, the background cells are introduced with various numbers of Gauss points independent from field nodes. Thehomogeneous beam is located on the spring support at the two ends and its length is 1000 mm. The cross-sectional dimensions are: width 25 mm and height 25 mm and Young's modulus and the mass density are equal to 70 GPa and 2700 kg/m³, respectively. Damage extent is considered as 50% reduction in the height of cross section (H_d) for more clarity of reflection ratio in captured signal. The corresponding value of reflection ratio for 50% damage extent is equal to 0.2533 based on Figure 3. This damage is located at 700 mm from the left end (P_s) and it has the length (L_d) equal to 2 mm. The excitation signal is a Hanning-windowed sinusoidal wave with the central frequency of 250 kHz, which is applied in the first degree of freedom u₁ (Figure 5). The output response signal is captured at the point S in location P_s equal to 425 mm. This point saves the incident and reflected waves from damage position. It is to be noted the damping matrix applies as C = ηM with η=0.02 in Equation (18). Figure 6 shows the benchmark signals of FE wave propagation analysis for undamaged and damaged state including signal's envelope. It is concluded in Figure 6 (b) the existence of 50% damage caused reflection ratio nearly equal to 25% corresponds to Figure 3.

Figure 4:
Schematic meshfree model of a damaged Euler-Bernoulli beam.

Figure 5:
Excitation signal in (a) time, (b) frequency domain.

Figure 6:
FE output signal in (a) undamaged, (b) damaged.

The main parameters that affects the results of wave propagation in MQ-type HRPIM includes: the number of failed nodes (n_f ), the dimensionless size of the support domain (α_s ), shape parameters of MQ (q, α_c), the number of polynomials (n_p ) and the total number of Gauss integration points (n_g ). In this study, first, Gauss points are assumed such as n_g > 3n_f for the stable solution and the prevention of singularity and other parameters are studied under this condition. In the last the effect of the arrangement of n_g is investigated independently.

3.2 Study of the Parameters of MQ (α_c, q) and Node Density

According to Eq. (13), three parameters in MQ affect the interpolation results, including q, α_c and d_c . d_c is the characteristic length as the space between two nodes, and it covers the number of field nodes (n_f ) implicitly. These parameters relate to each other, and independent evaluation of them is not rational. Therefore, Figs. 6-10 represent the contour plot of root mean square error (RMSE) and corresponding percentage error of the reflection ratio. These figures extracted for steps equal to 0.2 for q and α_c . The field node density equals to 100, 200, 300, 400 and 500. In this result two terms of polynomials (m=2) are considered for interpolation in Eq. (1) and n_g =3000 (three Gauss points on each of 1000 cells that six points are located in damage length) RMSE measures the difference between signals predicted by HRPIM and the signal which is obtained from the benchmark FE model in Figure 6 as

(34)

where f_t and f_t ^m are HRPIM and benchmark signals, respectively, and N is the number of sampling times. Also percent error (e_r ) of reflection ratio calculates as

(35)

where R_{HR PIM} is the reflection ratio of HRPIM signal and R_FE is equal to 0.2533 from Figure (5b).

Figure 7:
100 field nodes: (a-1)to(a-4): RMSE for α_s =1, 2, 4, 8. (b-1)to(b-4): percentage error of reflection ratio for α_s=1, 2, 4, 8.

Figure 8:
200 field nodes: (a-1)to(a-4): RMSE for α_s =1, 2, 4, 8. (b-1)to(b-4): percentage error of reflection ratio for α_s=1, 2, 4, 8.

Figure 9:
300 field nodes: (a-1)to(a-4): RMSE for α_s =1, 2, 4, 8. (b-1)to(b-4): percentage error of reflection ratio for α_s=1, 2, 4, 8.

Figure 10:
400 field nodes: (a-1)to(a-4): RMSE for α_s =1, 2, 4, 8. (b-1)to(b-4): percentage error of reflection ratio for α_s=1, 2, 4, 8.

Figure 11:
500 field nodes: (a-1)to(a-4): RMSE for α_s =1, 2, 4, 8. (b-1) to (b-4): percentage error of reflection ratio for α_s =1, 2, 4, 8.

The comparison of RMSE and e_r with regard to the node distribution shows the increasing of node density leads to less discrepancy and high accuracy between FE and HRPIM signals. In particular the reflection ratio is more sensitive to node density, regardless of the weight of other variables ((b-1) to (b-4) in Figures 7 through 11). Figure 12 (a) shows the variation of minimum RMSE of each node density and the corresponding e_r relative to the node densities while in Figure 12 (b) the minimum e_r and related RMSE is plotted. Both of these figures show the tendency to the least error by increasing of the node numbers. Although in Figure 12 (b) the RMSE is stable and e_r decreases. Figure 12 shows that minimum RMSE or minimum e_r are not dependant to the same values of (α_c,q) and α_s . Table 1 compares the results of (α_c,q) for different α_s in Figure 12. The important point is that the most of the values in Table 1 are in the vicinity of zones which are sensitive to any changes (q=0.2, α_c=0.2, q=5 and α_c=5) and it has relative high errors compared with amounts in the middle range of Figures 7 through 11. It is better we seek the condition that both RMSE and e_r would be acceptable.

Figure 12:
(a) Minimum RMSE and related e_r (b) Minimum e_r and related RMSE.

Table 1:
Values of (α_c,q ) for different node density in Figure 12.

αs	100 nodes		200 nodes		300 nodes		400 nodes		500 nodes
	Fig.11(a)	Fig.11(b)	Fig.11(a)	Fig.11(b)	Fig.11(a)	Fig.11(b)	Fig.11(a)	Fig.11(b)	Fig.11(a)	Fig.11(b)
1	(1.6,0.2)	(0.4,0.6)	(1.0,1.8)	(0.8,0.2)	(2.8,0.6)	(0.8,0.2)	(1.6,3.6)	(0.2,3.6)	(2.0,0.4)	(0.6,4.4)
2	(3.0,0.2)	(0.8,0.2)	(0.6,3.2)	(1.0,0.2)	(2.4,5.0)	(0.8 0.2)	(0.4,4.8)	(0.2,0.2)	(0.2,1.8)	(0.2,0.6)
4	(5.0,4.3)	(0.2,0.4)	(2.4,1.2)	(2.4 0.6)	(5.0,4.8)	(1.0,0.4)	(5.0,4.8)	(0.2,0.4)	(2.4,4.8)	(0.2,0.4)
8	(4.0,3.4)	(0.2,0.4)	(4.8,2.6)	(0.2,0.2)	(4.6,2.8)	(3.6,4.6)	(5.0,2.4)	(0.2,0.4)	(4.2,3.2)	(4.6,3.2)

In accordance with Figures 7 through 11, the integer values of q are inappropriate and the analysis would be unstable and singular. However, just q=1 leads to singularity in all of the node densitiesfor α_s =1. Figure 11 shows that n_f =500 has the least error among the other node numbers that has the acceptable errors. Figures 13 (a-1) to (d-2) shows the path of minimum RMSE and e_r along q-axis for various α_s . There is a good correlation between the minimum RMSE and related e_r (Figures (a-1) to (d-1)) compared with the other one even for α_s =1 which has the significant dispersion between results. It should be noted the results of α_s =4 has the constant difference between RMSE and e_r however it seems α_s =2 may lead to the reliable damage quantification.

Figure 13:
500 nodes: (a-1)to(d-1) minimum RMSE and related e_r relative to q for α_s =1, 2, 4, 8 (a-2) to (d-2) minimum erand related RMSE for α_s =1, 2, 4, 8 (a-3) to (d-3) distribution of minimum RMSE relative to q (× denotes the minimum value).

According to the Figs.12 (a-3) to (d-3) which are the distribution of the (α_c,q) related to the least error, the energy envelope of the captured HRPIM signal is plotted in Figure 14 (a), (b) using Equation (33a). Figure 14 (a) represents that there is a high-quality agreement among HRPIM and FE envelopes for health beam in such a way the incident and reflected energy of FE and HRPIM are matched. However, HRPIM envelopes with the same parameters lead to lower evaluation of the damage extent in Figure 14 (b). It is to be noted that incident waves are coincided, but the reflected signals by damage has lower amplitude. The computed reflection ratios are equal to 0.2025, 0.2305, 0.2345, 0.2308 for α_s =1, 2, 4, 8, respectively. Based on these reflection ratios, damage quantification by Figure 3 indicates the extent of damage as 44.62%, 47.71%, 48.12% and 23.08% respectively while the exact value is 50%.

Figure 14:
500 nodes: (a) energy envelope of health model (b) energy envelope of damaged model.

Above analysis studied the comparison of RMSE and e_r for different α_c , α_c and q which have the minimum error. Therefore the sensitivity analysis is required for investigation on the effect of each parameter on the other constant value. Figure 15 (a) shows RMSE and e_r for q while α_c =2.5 and α_s =4. The initial point of the graphs is q=0.01 and it has the reasonable error. In line with Figure 11, the integer values of q are suffering from the singularity, but even the imperceptible amount (±0.01) less or greater than the specified integer q serves the stability of numerical analysis. It is noticeable the results for two q±0.01 are close to each other. The other noticeable consequence is that the increasing of q leads to the reduction in both of RMSE and e_r, but from q=4. 99, RMSE increases until q reaches to 9.99 and two graphs are crossover and e_r surpassed RMSE. Moreover, q>15 have a propensity to high error, but there is no singularity and the computation remains stable. In this range, reflection ratio and damage quantification have an acceptable evaluation that is due to a shift in HRPIM signal relative to the benchmark FE signal without any change in amplitude. By increasing of higher than q=33.99 the analysis is singular.

Figure 15:
500 nodes, RMSE and e_r relative to (a) q (α_c =2.5) (b) α_c (q=2.5).

Figure 15 (b) shows the variation of RMSE and e_r relative to α_c when q=2.5 and α_s =4. It shows that both RMSE and e_r tends to low value, however, for α_c >7 the unstable results are achieved and quantification of damage cannot be applied.

A question that may be raised at this section is that the high density for condition q= 2.5, α_c =2.5, α_s =4 may lead to better results. Figure 16 represents the comparison of the wave propagation results for different values of n_f from 100 to 1000. by increasing of n_f from 100 to 500, the amplitude of reflected envelope wave by damage is increasing and it can be seen in Figure 16 (b), but it is obviously not only high density node distribution cannot lead to more accurate RMSE and e_r but also the enlargement of the reflected energy portion in Figure 16 (a) shows that n_f =600, 700 and 800 has the amplitude lower than n_f = 500. Although high RMSE can be seen for n_f =1000 relative to another node densities but its amplitude is close to the FE envelope. It should be noted that generally the obtained values for RMSE and e_r has the acceptable error. Since the configuration of background cells and Gauss point needs to modify for n_f >1000 and it affects the comparison with results of Figure 16 therefore it does not included.

Figure 16:
The effect of n_f (q= 2.5, α_c =2.5, α_s =4, n_g =3000) (a) envelope energy (b) RMSE and e_r .

3.3 Study of the Size of Support Domain

The number of field nodes that should be used to interpolate and shape function construction relates to the size of support domain. This parameter covers the domain around a specific Gauss points and the field nodes located within this domain contribute to interpolate. Figure 17 (b) represents as regards α_s=1 has improper RMSE equal to 34.2980 but it leads to acceptable e_r = 9.1219. It is concluded α_s=2 assess the reflection ratio a little better than α_s=4 however RMSE related to α _s=4 is the minimum quantity along α_s-axis. Moreover, the increasing of the size of support domain greater than α_s=4 does not lead to better results.

Figure 17:
The effect of α_s (q= 2.5, α_c =2.5, n_g =3000) (a) envelope energy (b) RMSE and e_r .

3.4 Study of the Number of Gauss Points

Gauss points (G_p ) have the most sensitive effect in meshfree techniques for two reasons. The first one is the determination of nodes for interpolation in the support domain. Second is the extraction of the system matrices using Equation (21). Unlike FE that a simple model of an open crack is a reduction of of a single element stiffness, but HRPIM exerts the different geometrical properties of the cracked region in the system matrices using Equation (21) for Gauss points which are placed in damaged section. As a result, the number of G_p (n_g ) is the second significant parameter after node numbers which affect the wave reflection from crack. In accordance with the influence power of n_g on the extraction of the damaged system matrices, it seems that the existence of more G_p in damaged length may lead to the exact reflection ratio and remarkable coincidence with not only FE signal as an approximated response but also to the analytical response. However, Figure 18 - as the reflected part of the energy envelope - shows that the high n_g may not necessarily have the high accuracy. In this figure the number of cells (n_cell ) is constant as n_cell =1000 and n_g increases. For set of n_g =[3,4,5,6,7,8,9,10] the number of G_p on damage length (L_d=2mm) are [6,8,10,12,14,16,18,20]. It is notable n_g =3 is the best result in comparison with n_g =10 which has the minimum RMSE and e_r .

Figure 18:
500 nodes (a) reflected energy envelope (b) RMSE and e_r nodes (q= 2.5, α_c =2.5, α_s =4, n_cell =1000).

Table 2 represents the results of different patterns of background cells and the number of G_p whereas the all number of G_p distributed in the problem domain is equivalent to 3000. The first point is that the existence of same G_p in damage length (L_d) does not guarantee the same energy envelope, RMSE and e_r . On the other hand the high G_p on L_d cannot produce the very accurate answers. For example, n_cell =10 and n_g =300 results in 28 G_p in L_d which has the acceptable RMSE equivalent to 2.5857 but it has low reflection ratio (e_r =9.5948) even compared to states that there are five G_p on L_d. It is noticeable, the patterns which has seven G_p on L_d (n_cell =375, n_g =8) and one with 5 points on L_d (n_cell =6, n_g =500) has the results better than states that there are 6 points on L_d. However, the pattern such as n_cell=8, n_g =375 produces three G_p on L_d with considerable reflection ratio error.

Table 2:
Different patterns of Gauss points for n_cell×n_g=3000 (n _f =500, q= 2.5, α_c =2.5, α_s =4).

ncell	ng	# Gp on Ld	RMSE	er	ncell	ng	# Gp on Ld	RMSE	er
3000	1	6	2.0902	7.2083	1	3000	4	2.4095	7.7646
1500	2	6	2.1530	8.3825	2	1500	4	2.3119	6.5414
1000	3	6	2.0605	8.2291	3	1000	6	2.5696	10.8060
750	4	6	2.7824	7.8160	4	750	5	1.8218	7.6374
600	5	6	3.5825	12.2180	5	600	4	1.9092	5.0506
500	6	6	1.7989	10.8070	6	500	5	2.4417	4.6012
375	8	7	2.2835	4.7081	8	375	3	3.7753	20.5920
300	10	6	2.5732	5.6499	10	300	28	2.5885	9.5948

3.5 Study of the Number of Polynomial Terms

Equation (1) is the formulation of pure HRPIM which is augmented with m polynomial basis functions. Essentially, adding polynomials is not necessary always; however it can ensure the inverse of matrix G in Equation (7), consistency, improvement of accuracy and stability of problem. Figure 19 shows the results of wave propagation relative to the number of polynomial terms (m). In this study, pure HRPIM (m=0) is stable and nonsingular. But Figure 19 (b) shows that increasing of m leads to accurate results with slight drop in RMSE and e_r . Although, for m=0, RMSE and e_r are equal to 2.7865 and 8.7022, respectively whereas m=1 evaluates RMSE=2.8071 and e_r =8.7595. The notable point is the instability of signals and envelopes for m>5. In this form, problem is not singular but detection of reflected wave is not possible as shown in Figure 19 (a) for m=6, 7.

Figure 19:
The effect of the number of polynomial terms (m) (n_f =500, q= 2.5, α_c =2.5, α_s =4, n_g =3000) (a) envelope energy (b) RMSE and e_r .

4 RESULTS

This paper investigated the numerical modeling of wave propagation for damage quantification in Euler-Bernoulli beam using Hermite RPIM by considering MQ-type RBF. The damage was considered as a reduction in the depth of the cross section in the specific position. Extraction of the relation between the reflection ratio and damage extent was utilized as a tool for the assessment of the capability of HRPIM. In addition to the reflection ratio, RMSE between the captured signals of HRPIM and FE models were determined for the sensitivity analysis of results relative to the effective parameters. Number of field nodes, shape parameters of MQ, the size of support domain, the number of Gauss points and the number of polynomial terms were the main effective variables. This study can draw the following conclusions:

1. The number of field nodes is the most significant parameter for accurate wave propagation modeling using HRPIM. It is concluded to achieve an accurate results, the consideration of a high density of field node is required.

2. MQ is a suitable RBF for HRPIM and wave propagation; however, care must be taken in the range of the proper shape parameters. Selection of the integer values for q makes the problem singular. This study showed that wave propagation and damage quantification are sensitive to high value of q and it is better q should be less than 15. There is such restriction about other shape parameter α_c and low value is preferred for this quantity.

3. It is better to consider α_s <5 for the size of support domain. Greater values make high RMSE and reflection ratio error.

4. Gauss quadrature points are important for determination of support domain and to derive the system matrices. This study investigated different patterns of Gauss points in such a way the total number of Gauss points is constant. It is concluded the existence of extra Gauss points in damage region does not guarantee the more accurate reflected wave. In general, the number of Gauss points is assumed as n_g > 3n_f for the stable solution and to prevent the singularity.

5. Pure HRPIM without any polynomial terms is acceptable and there is no singularity in solution but considering a few terms improve the accuracy; however more terms make the problem unstable.

References

• Bozkurt, O. Y., B. Kanber and M. Z. AŞIk (2013). Assessment of RPIM shape parameters for solution accuracy of 2D geometrically nonlinear problems. International Journal of Computational Methods 10(03): 1350003.
• Buhmann, M. D. (2003). Radial basis function: theory and implementation, Cambridge University Press.
• Chakraborty, A., D. R. Mahapatra and S. Gopalakrishnan (2002). Finite element analysis of free vibration and wave propagation in asymmetric composite beams with structural discontinuities. Composite Structures 55: 23-36.
• Chen, Y., J. Lee and A. Eskandarian (2006). Meshless methods in solid mechanics, Springer.
• Cui, X., G. Liu and G. Li(2009). A smoothed Hermite radial point interpolation method for thin plate analysis. Archive of Applied Mechanics 81(1): 1-18.
• Das, S. and T. Kundu(2009). Mesh-free Modeling of Ultrasonic Wave Fields in Damaged Layered Half-spaces. Structural Health Monitoring 8(5): 369-379.
• Dehghan, M. and A. Shokri (2009). Numerical solution of the nonlinear Klein-Gordon equation using radial basis functions. Journal of Computational and Applied Mathematics 230(2): 400-410.
• Doyle, J. F. (1997). Wave propagation in structures: spectral analysis using fast discrete Fourier transforms, Springer; 2nd edition.
• Fedoseyev, A. I., M. J. Friedman and E. J. Kansa (2002). Improved multiquadric method for elliptic partial differential equations via PDE collocation on the boundary. Computers & Mathematics with Applications 43(3-5): 439-455.
• Fornberg, B. and G. Wright(2004). Stable computation of multiquadric interpolants for all values of the shape parameter. Computers & Mathematics with Applications 48(5-6): 853-867.
• Franke, R. (1982). Scattered data interpolation: test of some methods. Mathematics of Computation 38: 181-200.
• Gao, L., K. Liu and Y. Liu (2007). A meshless method for stress-wave propagation in anisotropic and cracked media. International Journal of Engineering Science 45(2-8): 601-616.
• Ghaffarzadeh, H. and A. Mansouri(2013). Numerical modeling of wave propagation using RBF-based meshless method. Proc. of the Intl. Conf. on Future Trends in Structural, Civil, Environmental and Mechanical Engineering-FTSCEM, Bangkok, Thailand, IRED.
• Graff, K. F.(1991). Wave motion in elastic solids, Dover Publications.
• Hardy, R. L. (1971). Multiquadric equations of topography and other irregular surfaces. Journal of Geophysical Research 76(8): 1905-1915.
• Hon, Y. C. and Z. M. Wu(2000). A quasi-interpolation method for solving stiff ordinary differential equations. International Journal for Numerical Methods in Engineering 48: 1187-1197.
• Kanber, B., Ö. Y. Bozkurt and A. Erklig (2013). Investigation of RPIM Shape Parameter Effects on the Solution Accuracy of 2D Elastoplastic Problems. International Journal for Computational Methods in Engineering Science and Mechanics 14(4): 354-366.
• Kansa, E. J. (1990a). Multiquadrics-A scattered data approximation scheme with applications to computational fluid-dynamics-I surface approximations and partial derivative estimates. Computers & Mathematics with Applications 19(8-9): 127-145.
• Kansa, E. J. (1990b). Multiquadrics-A scattered data approximation scheme with applications to computational fluid-dynamics-II solutions to parabolic, hyperbolic and elliptic partial differential equations. Computers & Mathematics with Applications 19(8-9): 147-161.
• Keramat, A. and A. Ahmadi (2012). Axial wave propagation in viscoelastic bars using a new finite-element-based method. Journal of Engineering Mathematics 77(1): 105-117.
• Kim, Y. Y. and E. H. Kim (2000). A new damage detection method based wavelet transform. Proceedings of 18th IMAC, San Antonio, Texas, USA.
• Kudela, P., M. Krawczuk and W. Ostachowicz (2007). Wave propagation modelling in 1D structures using spectral finite elements. Journal of Sound and Vibration 300(1-2): 88-100.
• Lazzaro, D. and L. B. Montefusco (2002). Radial basis functions for the multivariate interpolation of large scattered data sets. Journal of Computational and Applied Mathematics 140: 521-536.
• Li, H. and S. S. Mulay (2013). Meshless methods and their numerical properties, CRC Press.
• Li, S. and W. K. Liu(2007). Meshfree particle methods, Springer.
• Liew, K. M. and X. L. Chen(2004). Buckling of rectangular Mindlin plates subjected to partial in-plane edge loads using the radial point interpolation method. International Journal of Solids and Structures 41(5-6): 1677-1695.
• Liew, K. M. and X. L. Chen(2004). Mesh-free radial point interpolation method for the buckling analysis of Mindlin plates subjected to in-plane point loads. International Journal for Numerical Methods in Engineering 60(11): 1861-1877.
• Liu, G. R. (2010). Mesh Free Methods Moving Beyond the Finite Element Method, CRC Press.
• Liu, G. R. and Y. T. Gu (2001). A point interpolation method for two dimensional solids. International Journal for Numerical Methods in Engineering 50(4): 937-951.
• Liu, G. R., G. Y. Zhang, Y. T. Gu and Y. Y. Wang (2005). A meshfree radial point interpolation method (RPIM) for three-dimensional solids. Computational Mechanics 36(6): 421-430.
• Liu, G. R., X. Zhao, K. Y. Dai, Z. H. Zhong, G. Y. Li and X. Han (2008). Static and free vibration analysis of laminated composite plates using the conforming radial point interpolation method. Composites Science and Technology 68(2): 354-366.
• Liu, G. R., Y. Li and Y. Dai (2006). A linearly conforming radial point interpolation method for solid mechanics problems. International Journal of Computational Methods 3(4): 401-428.
• Liu, Y., Y. C. Hon and K. M. Liew (2006). A meshfree Hermite-type radial point interpolation method for Kirchhoff plate problems. International Journal for Numerical Methods in Engineering 66(7): 1153-1178.
• Mitra, M. and S. Gopalakrishnan (2005). Spectrally formulated wavelet finite element for wave propagation and impact force identification in connected 1-D waveguides. International Journal of Solids and Structures 42(16-17): 4695-4721.
• Moser, F., L. J. Jacobs and J. Qu (1999). Modeling elastic wave propagation in waveguides with the finite element method. NDT&E International 32: 225-234.
• Nguyen, V. P., T. Rabczuk, S. Bordas and M. Duflot(2008). Meshless methods: A review and computer implementation aspects. Mathematics and Computers in Simulation 79(3): 763-813.
• Rippa, S. (1999). An algorithm for selecting a good value for the parameter c in radial basis function interpolation. Advances in Computational Mathematics 11: 193-210.
• Su, Z., L. Ye, X. Bu, X. Wang and Y. W. Mai(2003). Quantitative Assessment of Damage in a Structural Beam Based on Wave Propagation by Impact Excitation. Structural Health Monitoring . 2(1): 27-40.
• Sun, K., G. Meng, F. Li, L. Ye and Y. Lu(2009). Damage Identification in Thick Steel Beam Based on Guided Ultrasonic Waves. Journal of Intelligent Material Systems and Structures 21(3): 225-232.
• Tiago, C. M. and V. M. A. Leitão(2006). Application of radial basis functions to linear and nonlinear structural analysis problems. Computers & Mathematics with Applications 51(8): 1311-1334.
• Wang, J. G. and G. R. Liu (2000). Radial point interpolation method for elastoplastic problems. Proceedings of the First International Conference on Structural Stability and Dynamics, Taipei, Taiwan.
• Wang, J. G. and G. R. Liu(2002). A point interpolation meshless method based on radial basis functions. International Journal for Numerical Methods in Engineering 54(11): 1623-1648.
• Wang, J. G. and G. R. Liu (2002). On the optimal shape parameters of radial basis functions used for 2-D meshless methods. Computer Methods in Applied Mechanics and Engineering 191: 2611-2630.
• Wang, J. G., G. R. Liu and P. Lin(2002). Numerical analysis of Biot's consolidation process by radial point interpolation method. International Journal of Solids and Structures 39(6): 1557-1573.
• Wen, P. H. (2010). Meshless local Petrov-Galerkin (MLPG) method for wave propagation in 3D poroelastic solids. Engineering Analysis with Boundary Elements 34(4): 315-323.
• Wendland H. (2010). Scattered data approximation. Cambridge University Press.
• Wu, Z. (1992). Hermite-Birkhoff interpolation of scattered data by radial basis function. Application Theory and its Applications 8: 1-10.
• Yim, H. and Y. Sohn (2000). Numerical simulation and visualization of elastic waves using mass-spring lattice model. IEEE transactions on ultrasonics, ferroelectrics, and frequency control 47(3): 549-558.
• Zhang, G. M. and R. C. Batra (2004). Modified smoothed particle hydrodynamics method and its application to transient problems. Computational Mechanics 34(2).
• Zhao, X. G. and J. L. Rose (2003). Boundary element modeling for defect characterization potential in a wave guide. International Journal of Solids and Structures 40(11): 2645-2658.

Publication Dates

• Publication in this collection
  2016

History

• Received
  27 Feb 2016
• Reviewed
  21 July 2016
• Accepted
  25 July 2016