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:
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
where
and the matrices are defined as
In addition, the polynomial terms have to satisfy the following constraint to obtain a unique solution
Arranging Eqs. (2) and (5) together leads to the following set of equations
where U is the vector of nodal values and G is the generalized moment matrix. Substituting Eq. (7) into Eq. (1) leads to
Therefore the Hermite-RPIM shape function can be introduced as
Finally, the approximated function can be obtained by
The shape functions ϕ _{i} and ϕ _{i} ^{x} (i = 1, 2, ..., n) are defined by
The derivatives of shape functions can be easily expressed for the first order as
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:
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
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
where ū are the prescribed displacement on the essential boundary Γ_{u} and are tractions on the natural boundary Γ_{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
By using Eqs. (15) and (8), the discretization of Eq. (17) yields to:
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
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
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∂^{2}/∂x ^{2}{Φ}.
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
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.
2.3 Time Integration Scheme
The central difference time integration scheme is used to solve Eq. (18) as an explicit method
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})
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:
the eigen roots (wavenumber) for Eq. (23) can be yielded as
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
where each term corresponds to one of the wave modes depending on their propagating direction and properties: the positively propagating wave (A _{1}), the negatively propagating wave (A _{2}), the positively growing ringingwave (A _{3}) and the negatively exponentially decaying wave (A _{4}). 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
where the three terms in Eq. (27a) correspond to the incident wave (A _{1}), reflected propagating wave (A _{2}) and reflected exponentially decaying wave (A _{4}), respectively, while the two terms in Eq. (27b) are the damage-induced positively propagating wave (B _{1}) and growing attenuating wave (B _{2}), respectively. The wave components A _{4} and B _{2} 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})
The reflection ratio can be found as the relation between the magnitude A _{2} of the reflected wave and the magnitude A _{1} of the incident wave
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).
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})
The Hilbert transform performs a 90º phase shift to construct the so-called analytic signal Z(t) as
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
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^{3}, 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_{1} (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.
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
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
where R_{HR PIM} is the reflection ratio of HRPIM signal and R_{FE} is equal to 0.2533 from Figure (5b).
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.
α_{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.
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%.
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 (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.
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.
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} .
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.
n_{cell} | n_{g} | # G_{p} on L_{d} | RMSE | e_{r} | n_{cell} | n_{g} | # G_{p} on L_{d} | RMSE | e_{r} |
---|---|---|---|---|---|---|---|---|---|
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.
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:
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.
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.
It is better to consider α_{s} <5 for the size of support domain. Greater values make high RMSE and reflection ratio error.
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.
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.