Dynamic analysis of cracked plates based on a first-order shear deformation theory formulation by using an extended meshfree method

An extended meshfree method is employed in this paper for investigating the dynamic behaviour of cracked plates based on the first-order shear deformation theory (FSDT). The FSDT is a straightforward formulation with the assumption of first-order shear deformation as its name implies, which is appropriate for relatively thick plates. In this study, the meshfree method is chosen as an alternative to the conventional mesh-based methods to model plate structures. Among various meshfree formulations, Moving Kriging (MK) is a method that satisfies the Kronecker delta property, allowing for the easy imposition of essential boundary conditions. An extended MK formulation is proposed in this paper to model cracked plates without explicitly pre-defining the crack in the geometry domain. In the extended concept, the extrinsic enriched functions are employed to model the discontinuity due to the crack. Particularly, the Heaviside step function is employed to describe the discontinuity of the displacement fields on two sides of the crack surface. And the asymptotic enriched functions are used for stress singularity around the crack tip. In the dynamic analysis of cracked plates, one of the important factors that must be evaluated is the dynamic stress resultant intensity factor (DSRIF). In this paper, the DSRIFs are shown through many numerical examples and compared with analytical solutions and other numerical methods, showing the accuracy and efficiency of the present extended MK approach.


Introduction
Fracture analysis of thin-walled structures is an important task in many industrial fields such as mechanical, civil and aerospatial engineering since it directly affected the durability and lifetime of the structures.For that reason, more research on cracked plates is necessary, especially dynamic fracture.Dynamic cracked problems are of interest when the inertia affects the behavior of the cracked plates, such as when the plate is subjected to a sudden load for a short time.There have been plenty of references on dynamic fracture for 2D and 3D problems, but scarce papers on dynamic cracked plates have been published.Because of this, dynamic cracked analysis is still a research gap that needs to be exploited.For numerical modelling of plate structures, using a plate formulation gives better performance with high accuracy than using a 3D model.With different mechanical assumptions to simplify the model, different plate theories are introduced as follows: the classical plate theory (CPT) [1 -4] with the assumption that the plate cross-section after deformation is still perpendicular to the mid-surface, the first-order shear deformation theory (FSDT) [5 -8] assumes that after deformation, the cross-section of the plate is no longer perpendicular to the mid-surface and is determined by a first-order function of the rotation angle, and higher-order shear deformation theory (HSDT) [9 -11] assumes that after deformation, the cross-section of the plate is no longer a flat surface.For moderately thick plate, the FSDT and HSDT are appropriate to model the mechanical behavior of the plates.However, HSDT usually requires C 1 or even higher continuity, bringing sophisticated computation.Whereas FSDT requires only C 0 continuity which makes it a simple plate formulation.Despite this, when the thickness of a plate is relatively small compared to its other dimensions, a numerical issue called the shear-locking is likely to occur [12 -14].Alleviating shear-locking is important to ensure accuracy in modeling thin cracked plates, several methods such as Assumed Natural Strain (ANS) [15], Enhanced Assumed Strain (EAS) [16] and the Mixed Interpolation of Tensorial Components (MITC) [17,18] can be used.
The extended finite element method (XFEM) [19] is a powerful and versatile numerical method in fracture analysis for modeling crack discontinuity without explicitly defining it in the problem geometry.This makes it possible for XFEM to simulate crack propagation without re-meshing.To do that, the enrichment functions are added to the conventional finite element basis functions to accurately represent the crack tip singularity and the strong discontinuity in displacement fields.The extended finite element method has been used to solve a wide range of dynamic fracture problems in two and three dimensions, [20 -22], and also cracked plates using FSDT theory [23].
Besides FEM and its variants, other numerical methods are also being strongly developed, especially the group of meshfree methods.Meshfree approaches discretize the continuous domain into a set of scattered nodes, as the name indicates, and as a result, no element is required.There have been many significant contributions to meshfree methods in the literature, e.g. the Element Free Galerkin method [24,25], the Reproducing Kernel Particle Method [26,27], the Radial Point Interpolation method (RPIM) [28,29], the RPIM with strain smoothing technique [30], the Moving Kriging method (MK) [31,32], etc.Among the many meshfree methods, the MK and the RPIM methods are the few that have the desirable Kronecker delta property.This property allows for the direct imposition of boundary conditions, which is similar to conventional FEM.However, the concept of these two approaches is nonidentical, the comparisons between RPIM and MK are pointed out in [33,34].The meshfree method can also be extended to solve cracked problems like XFEM.In a similar way to formulating XFEM, the XRPIM was proposed by Nguyen et al. [35 -37] by combining RPIM and the enrichment functions.The method was first used for 2D fracture analysis.Later, Truong et al. [29] employed XRPIM for static analysis of cracked FSDT plates.Nevertheless, the number of studies on fracture analysis of cracked plates using extended meshfree methods is still limited, especially dynamic fracture analysis.Dynamic fracture analysis for bending plates is a challenging problem, a small number of publications on this topic can be mentioned as follows: dynamic fracture analysis of Kane-Mindlin plates by using the dual boundary element method (DBEM) [38], using DBEM [39] and dual reciprocity BEM (DRBEM) [40] for dynamic analysis of cracked Reissner-Mindlin plates (also known as the FSDT plate), the Peridynamic (PD) approach for solving dynamic cracked Mindlin plates is shown in [41].In the scope of this study, dynamic analysis of cracked FSDT plates is conducted by using the extended MK (XMK) for the first time.
This study evaluates the dynamic stress resultant intensity factors (DSRIFs) of the FSDT plates with through-thickness cracks which is an important feature in fracture analysis.The method used for modelling is the XMK which is formulated based on the XRPIM.The MK shape function in this study is constructed by using the Quartic correlation function.This function has the property of being independent of the correlation parameter so that the MK shape function is the same regardless of the value of the correlation parameter.The DSRIFs are calculated using the interaction integral.Many numerical examples are conducted, demonstrating the accuracy and validity of the proposed XMK method.
Here is the outline of the paper.After the Introduction, Section 2 briefly overviews the meshfree method with the Moving Kriging shape function.Then, Section 3 provides a review of the FSDT plate.After that, Section 4 presents the proposed formulation of the extended Moving Kriging method for dynamic analysis of cracked plates.The results of several numerical examples are presented in Section 5 to demonstrate the accuracy and efficiency of the proposed extended Moving Kriging method.Finally, the last Section gives some conclusions for the obtained results and outlooks for future studies.

Moving Kriging interpolation
Consider a computational domain  which is discretized into a set of scattered node i x .The function

( )
ux at an arbitrary point x can be interpolated by using the values of the neighboring nodes that are located inside the support domain  x of the interested point.The Moving Kriging interpolation (MKI) ( ) ux is described as the following expression [31] ( ) ( ) ( ) ( ) ( ) where n denotes the number of nodes inside the support domain  x , and the vector ( ) For example, the quadratic form of ( ) The vector ( ) in which ( ) Nevertheless, the correlation parameter affects the accuracy of MK interpolation, which brings up the difficulty of choosing the optimum  for each specific case.The solution of the MK meshfree method is therefore depending on the specific problem being solved [42].This drawback is addressed by the application of the Quartic function in [42,43], which is rarely influenced by the correlation parameter  .To put it another way, the stability of the computational results is not affected by changing the correlation parameter.Because of this, the Quartic function is employed in this study as the correlation function ( ) in which s l is the largest distance between any two nodes inside the support domain.The MK shape function at the th i node ( )  i x is defined from Eq. ( 1) as follows ( ) ( ) ( ) The coefficient matrices a and b in Eq. ( 1) are defined by the following equations ( ) ( ) The correlation matrix R is expressed as The MK shape function in Eq. ( 7) is now easily obtained.It is noticed that

First-order shear deformation theory
In the FSDT plate, the cross-section of the plate is assumed to remain straight after deformation but not perpendicular to the mid-surface [29].Figure 1 shows the sign convention and notation of the FSDT plate.In the Cartesian coordinate, the displacement fields of a plate is defined as in which ( ) w x x denotes the deflection of a point ( ) on the mid-surface, while 1  and 2  are, respectively, the rotation about 2 x and 1 x axes.The dynamic equilibrium equations are expressed as [40] 3 where t is the thickness of the plate and  is the mass density.M  and Q  ( )


are the components of the moment and the shear force resultants.
From Eq. ( 14), the bending strains are derived as ( ) and the shear strains are defined as For an isotropic homogeneous elastic material, the constitutive equations [23] are given as below where M and Q , in that order, denote the moment resultants and the shear force resultants.The bending stiffness matrix b D is defined as 6 ( ) and the shear stiffness matrix s D is expressed as in which E is the Young's modulus,  denotes the Poisson ratio,  is the shear modulus and the number 5/6 in Eq. ( 20) denotes the shear correction factor [23].
In dynamic problems, the inertia effect is significant and cannot be neglected.Therefore, the implicit dynamic equation without damping is given as += mu Ku F (23) The lumped mass matrix expressed in lower case m to avoid confusion with the notation of moment where Λ is the inertia tensor and given by [14]   3 3 00 0 12 0 0 0 12 The stiffness matrix is defined as [29] 3 5 dd 12 6 (26) in which the subscript letter "b" stands for bending component and "s" is for shear component, for a given node i , For a support domain with n nodes, the strain computing matrices in Eq. ( 26) are obtained by appending the matrices in Eq. ( 27) of every node within the support domain, e.g., 12 ...

XMK formulation
In the extended formulation, the enriched functions are inserted into the interpolation functions (see Eq. ( )) for deflection ( ) w and rotations ( )  as expressed below [23] ( ) where W is the set of all nodes that are in the support domain, s W is the set of nodes that are in the support domain cut by the crack, and t W is the set of scattered nodes whose support domain contains the crack tip, see Figure 2.   ( ) in which ( ) f x is the sign distance function, Figure 2 illustrate the positive and negative value of ( ) The asymptotic enrichment function ( ) Fr  for the rotations 1  and 2  contains terms proportional to 12 r [23]   ( ) where ( ) , r  defines a local polar coordinate system at the crack tip, see Figure 2.And for the deflection, the enrichment function ( ) Gr  consists of terms proportional to 32 r [23].
( ) The strain-computing matrix B for a cracked plate is now defined as The standard matrix standard B is performed in Eq. (27).And the enriched matrix enriched B for bending contribution 00 00 0 The enriched matrix for shear contribution is expressed as

Evaluation of dynamic stress resultant intensity factors
In fracture analysis of cracked plate, the stress resultant intensity factors (SRIFs) are defined as [23] ( ) ( ) ( ) where 12 , KK denote the moment intensity factors and 3 K is the shear force intensity factor.To compute the dynamic SRIFs, first, the dynamic J-integral is defined as [44] ( ) ( ) where the kinetic energy T is computed as and the strain energy W is defined as ( ) Now, insert an auxiliary state ( ) , , , ).The auxiliary states for the moment and shear force resultants, deflection and rotations are given in reference [29].The sum of two states 1 sum J is expressed as below where sum W and sum T are the sum energy of two states and expressed as Expand equation ( 42) and then rearrange it, the following relation is derived in which I is called the interaction integral and computed as In numerical computation, domain integral is more appropriate than contour integral.Hence, using an enclose contour C as shown in Figure 3  In order to obtain the interaction integral I , the term , H  need to be computed Expanding and rearranging this equation, note that the auxiliary state aux w and aux   are equal 0 so all the terms containing it are vanish, , H  is shorten as From ( 15), ( 16) and ( 53) The interaction integral I is now derived as the following It is observed that the first two integrals on the right-hand side are similar to the interaction integral proposed by Dolbow in the static case [23], the last integral indicates the effect of inertia to the interaction integral in dynamic case.
From [23] the relationship between I and 1 2 3 ,, K K K is expressed as Hence, the moment ( ) KK and shear force ( ) K intensity factors are obtained as follows In practice, the interaction integral domain used is shaped like Figure 4, which means the area 0 0 = A .The weight function is now re-defined as  In this study, all the numerical examples for transient dynamic cracked plates analysis employ the Newmark integration.More detail about the Newmark scheme can be found in reference [36].

Numerical examples
The performance of the proposed XMK method is presented in this section through various cracked plate problems.Particularly, four numerical problems with several types of loads and boundary conditions are presented as follows: -Validation of the proposed method.
-A cantilever square plate with a central crack.
-A square plate with a central crack.
-A square plate with a side crack emanating from a circular hole.The stress resultant intensity factors are evaluated in each problem, and the variation of SRIFs with respect to different crack orientations and time is also examined.All the numerical results from the present approach are compared to other numerical methods and analytical solutions.
As mentioned in Section 2, the Quartic function is not affected by the correlation parameter  .For that reason,


. The plate is discretized into a set of 50 50  scattered nodes.As mentioned in Section 4.2, the static SRIFs are computed by ignoring the third integral in Eq. (55).For comparison purposes, the analytical solutions are given as [23] ( ) ( ) where  is the angle of inclination (see Figure 5).The values of coefficients ( ) ( ) ( )  ta are determined based on the ta ratio [45].For the given dimensions of this problem, these coefficients are ( ) The variation of the intensity factors 12 , KK and 3 K with respect to the angle of inclination  is shown in Figure 6.The figure shows that 1 K decreases when the orientation angle  increases, while  It is observed that the SRIFs results obtained from the proposed XMK method fit well with the analytical solution.However, it is also noticed when  in the range of 30 o to 50 o , the shear force intensity factor 3 K has small recognizable deviations compared to the analytical one.This can be explained by the fact that 3 K is caused by only the shear force.From Eq. (20) a correction term 5/6 is added to approximate the second-order distribution of shear stress into a constant distribution, which reduces the accuracy of the FSDT plate theory.One can also see that the value of 1 K and 2 K fit well with the analytical solution and are more accurate than 3 K .Therefore, this error in 3 K is an inherent property of the FSDT plate and not due to the proposed XMK method.Firstly, the angle of inclination  is set to 0. The variation of the normalized 1 K with respect to time is shown in Figure 8.The XMK result is compared to the Peridynamic solution in reference [41], showing good agreement.In both methods, the curve has the same shape.During the examined time period in the figure, the normalized 1 K slightly increases and then dramatically decreases.The negative sign of 1 K is due to the direction in which the load is applied.9 shows the relationship between the 1 K curves when the angle  changes.It can be observed that when the angle  increases, the magnitude of 1 K decreases.In addition, the three curves have the same shape.The time at which 1 K begins to variate in all three cases is the same, although the distance from the crack tip to the force-applied edge is different in the three cases.The normalized DSIF is shown in Figure 11 for the CCCC case and Figure 12 for the SSSS case.The XMK result is compared to the Peridynamic solution in reference [41], showing good agreement.In both figures, the curves of both methods have the same shape.During the examined time period in the figures, the normalized 1 K in the CCCC case is smaller than in the SSSS case.In both cases, 1 K trivially decreases at the beginning then turns to rise sharply and peaks, after that falls sharply.At the end of the period, 1 K bounces back and tends to increase.Figure 14 shows the variation of the normalized 1 K with respect to time.The XMK result is compared to the Peridynamic solution in reference [41], showing good agreement.During the examined time period in the figure, the normalized 1 K trivially decreases at the beginning then turns to rise sharply and peaks, after that falls sharply.At the end of the period, 1 K bounces back and tends to increase.

Conclusions
In this work, the meshfree MK method has been extended to cooperate with the first-order shear deformation theory plate for dynamic analysis of cracked plate problems for the first time.One advantage of using the MK shape function is that it satisfies the Kronecker delta property, which means that boundary conditions can be easily imposed without any special treatments.Another worth noting is that the use of the Quartic correlation function helps the MK approach not to depend on the correlation coefficient.And the attractiveness of using an "extended" formulation is to model cracked plates without explicitly pre-defining the crack in the geometry domain.In the dynamic analysis of cracked plates, one of the important factors that must be evaluated is the dynamic stress resultant intensity and it can be obtained by using the interaction integral.The accuracy of the proposed method is proved in the evaluation of the dynamic stress resultant intensity factors.Numerous numerical examples demonstrate good agreement with analytical solutions and other numerical approaches.The variation of DSRIFs in different cases is also observed.

ijR
xx denotes the correlation function between two nodes i x and j x .There are many different correlation functions ( ) , ij R xx .However, the Gaussian function is the well-known form of the correlation function ji a and ki b are constants while the partial derivatives of the correlation function ( ) k r x and polynomial function ( ) j p x are uncomplicated to obtain.Therefore, it is simple to derive the first-order derivatives of the MK shape function as shown below

Figure 1 .
Figure 1.Notation and sign convention of degrees of freedom and resultants of the FSDT plate.

Figure 2 .
Figure 2. Enriched nodes and local polar coordinate systems.

Figure 3 .
Figure 3. Illustration of domain integral.Domain A is enclosed by

.
49)On the outer contour 0  the weight function has the value of zero so Now the expression of G remains only the first integral.On the inner contour  , i m takes the valueThe domain integral is derived by applying the divergence theorem

Figure 5 .
Figure 5. Square plate with an inclined central crack.

2 K and 3 K
increase and decrease symmetrically.

Figure 6 .
Figure 6.Variation of 12 , KK and 3 K with respect to the crack orientation  .

Figure 7 .
Figure 7.A cantilever square plate with a central crack.

Figure 8 .
Figure 8. Normalized 1 K of a square plate. Figure

Figure 11 .
Figure 11.Normalized 1 K of the clamped plate.

Figure 14 .
Figure 14.Normalized 1 K of the perforated plate.