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

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 timefrequency 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 Hosein Ghaffarzadeh a Majid Barghian a Ali Mansouri a,* Morteza.H Sadeghi b


Study on Meshfree Hermite Radial Point Interpolation Method for Flexural Wave Propagation Modeling and Damage Quantification 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 timefrequency 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 andGopalakrishnan 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, 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  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 . Rippa (1999) proposed an optimization algorithm for selection of a good value for shape parameters in MQ, inverse MQ and Gaussi-Latin American Journal of Solids and Structures 13 (2016) 2606-2627 an interpolation in data fitting. Kansa (1990a) proposed a scheme that allowed the shape parameter of MQ to vary with basis functions.  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 , Liu et al. 2005, plate structures , geometrical nonlinearity , composite materials (Liu et al. 2008) and consolidation process . 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,  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.

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: 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 Latin American Journal of Solids and Structures 13 (2016) 2606-2627 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 derivatives of shape functions can be easily expressed for the first order as Because of the existence of G  , 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 c d is the characteristic length that is usually the average nodal spacing for all the n nodes in the support domain.

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. u  , u  and u are acceleration, velocity and displacement, respectively, while u is the displacement gradient matrix. The natural and essential boundary conditions are where u are the prescribed displacement on the essential boundary Γu and t are tractions on the natural boundary Γt.
     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 Uare 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 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 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 wellknown 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 ng is the number of Gauss points in each background cell, i w is the Gauss weighting factor for the i th Gauss point at Qi x ( Figure 1) and D ik J is Jacobian matrix for the area integration of the background cell k,at which Gauss point Qi x is located.

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  is the largest frequency of the system (Chen et al. 2006). This condi- are assumed to be implemented as initial displacement and velocity.

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 Hd. 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: where each term corresponds to one of the wave modes depending on their propagating direction and properties: the positively propagating wave (A1), the negatively propagating wave (A2), the positively growing ringingwave (A3) and the negatively exponentially decaying wave (A4). The coefficient of each term represents the amplitude of each wave component which may be a complex problem.

A iI k I I k k iI I k k I I k k iI k R A I k I I k k I I k k
where l and H are the width and depth of the cross-section of the beam, respectively; while Hd 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, Hd, can be determined with a given reflection ratio. Figure 3 represents the absolute value of the reflection ratio and the damage degree (Hd/H).

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. 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 (Hd) 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 (Ps) and it has the length (Ld) 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 u1 ( Figure 5). The output response signal is captured at the point S in location Ps 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 (nf), the dimensionless size of the support domain (αs), shape parameters of MQ (q, αc), the number of polynomials (np) and the total number of Gauss integration points (ng). In this study, first, Gauss points are assumed such as

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 dc. dc is the characteristic length as the space between two nodes, and it covers the number of field nodes (nf) 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 ng=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 t f and m t f are HRPIM and benchmark signals, respectively, and N is the number of sampling times. Also percent error (er) of reflection ratio calculates as where HRPIM R is the reflection ratio of HRPIM signal and FE R is equal to 0.2533 from Figure (     The comparison of RMSE and er 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 par-ticular 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 er relative to the node densities while in Figure 12(b) the minimum er 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 er decreases. Figure 12 shows that minimum RMSE or minimum er 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 er would be acceptable.  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 Table 1: Values of (αc,q) for different node density in Figure 12.
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 nf =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 er along q-axis for various αs. There is a good correlation between the minimum RMSE and related er (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 er 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 Equa-tion (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 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 e r and related RMSE for α s =1, 2, 4, 8 (a-3) to (d-3) distribution of minimum RMSE relative to q (× denotes the minimum value). Above analysis studied the comparison of RMSE and er 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 er 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 er, but from q=4. 99, RMSE increases until q reaches to 9.99 and two graphs are crossover and er 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 er relative to αc when q=2.5 and αs=4. It shows that both RMSE and er 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 nf from 100 to 1000. by increasing of nf 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 er but also the enlargement of the reflected energy portion in Figure 16(a) shows that nf=600, 700 and 800 has the amplitude lower than nf= 500. Although high RMSE can be seen for nf=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 er has the acceptable error. Since the configuration of background cells and Gauss point needs to modify for nf>1000 and it affects the comparison with results of Figure

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 er= 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.

Study of the Number of Gauss Points
Gauss points (Gp) 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 Gp (ng) is the second significant parameter after node numbers which affect the wave reflection from crack. In accordance with the influence power of ng on the extraction of the damaged system matrices, it seems that the existence of more Gp 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 ng may not necessarily have the high accuracy. In this figure the number of cells (ncell) is constant as ncell=1000 and ng increases. For set of ng= [3,4,5,6,7,8,9,10] the number of Gp on damage length (Ld=2mm) are [6,8,10,12,14,16,18,20]. It is notable ng=3 is the best result in comparison with ng=10 which has the minimum RMSE and er. Table 2 represents the results of different patterns of background cells and the number of Gp whereas the all number of Gp distributed in the problem domain is equivalent to 3000. The first point is that the existence of same Gp in damage length (Ld) does not guarantee the same energy envelope, RMSE and er. On the other hand the high Gp on Ld cannot produce the very accurate answers. For example, ncell=10 and ng=300 results in 28 Gp in Ld which has the acceptable RMSE equivalent to 2.5857 but it has low reflection ratio (er=9.5948) even compared to states that there are five Gp on Ld. It is noticeable, the patterns which has seven Gp on Ld (ncell=375, ng=8) and one with 5 points on Ld (ncell=6, ng=500) has the results better than states that there are 6 points on Ld. However, the pattern such as ncell=8, ng=375 produces three Gp on Ld with considerable reflection ratio error.

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 er. Although, for m=0, RMSE and er are equal to 2.7865 and 8.7022, respectively whereas m=1 evaluates RMSE=2.8071 and er=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.

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: (a) 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. (b) 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. (c) It is better to consider αs<5 for the size of support domain. Greater values make high RMSE and reflection ratio error. (d) 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 3 g f n n  for the stable solution and to prevent the singularity. (e) 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.