On tracking arbitrary crack path with complex variable meshless methods

This study presents a numerical modelling framework based on complex variable meshless methods, which can accurately and efficiently track arbitrary crack paths in two-dimensional linear elastic solids. The key novelty of this work is that the proposed meshless modelling scheme enables a direct element-free approximation for the solutions of linear elastic fracture mechanics problems. The complex variable moving least-squares approximation with a group of simple complex polynomial basis is applied to implement this meshless model, with which the fracture problems with both stationary or progressive cracks are considered and studied. The effects of choosing different definitions of weighted complex variable error norm and different forms of complex polynomial basis on the computational accuracy of crack tip fields and crack paths prediction are analysed and discussed. Five benchmark numerical examples were studied to demonstrate the superiority of the present complex variable meshless framework over a standard element-free Galerkin method in tracking arbitrary crack paths in two-dimensional elastic solids.

In general, the advantages of applying meshless methods in tracking arbitrary crack propagation are mainly attributed from two aspects.On one hand, compared with the methods using predefined elements, arbitrary crack propagation paths can be simulated and tracked, easily and accurately, by a set of nodes and their local influence/support domains .For the element-based method for which naturally repel discontinuities, incorporating cracks generally involves time-consuming remeshing schemes or complicate modification to its shape functions with extended enrichments.However, in meshless methods, only the local definition of shape functions needs to be modified when approximating the discontinuous field variables in the support domain where the cracks propagate [8][9][10][11].Therefore, a lot of advanced techniques for the processing of crack discontinuity field had been developed, such as the transparency method [8,9], the diffraction method [8,9], the level set approach [12,19], the crack particles scheme [14], and the weight-enriched function [18].On the other hand, accurate computation of the crack tip stress fields is critical for the determination of crack nucleation and propagation .The node/domain-based scheme provides a much more flexible means than the element-based method for the simulation of crack-tip-singular-field.This is because either extended basis functions [8][9][10][11][12] or any other kinds of enrichment methods [10,13,26] can be straightforwardly applied by the node/domain-based meshless methods.
The complex variable meshless methods are a type of improved meshless methods that apply the complex polynomial basis functions to construct the trial functions for the unknown fields in mechanics problems [35][36][37][38][39][40].The number of terms of the complex variable basis functions that possess the same modelling degree of freedom is much less than that of the corresponding real polynomial basis functions.As such, the number of minimum nodes in each nodal support domain is largely reduced and the computational efficiency improves, substantially.This complex variable meshless modelling idea was first proposed in J.H. Li's doctoral dissertation [35], in which the complex variable moving least-squares (CVMLS) approximation and the complex variable element-free Galerkin (CVEFG) method were developed.Despite the flaw in its error functional [36,40], the CVEFG method had been successfully applied in the nonlinear analysis of solids and structures [39,41].A few years later, the error functional of the CVMLS approximation was corrected by H.P. Ren, who developed an improved CVMLS (ICVMLS) approximation [36].At the same time, another set of ICVMLS formulas based on the conjugate form of the original complex polynomial basis functions was proposed and shown a better computational accuracy for certain kinds of problems [40].Both the ICVMLS approaches with the original form and the conjugate form of complex variable polynomial basis have been successfully applied for the analysis of nonlinear materials with large deformation, such as gels [42][43][44].Recently, the complex variable meshless methods had been applied to solve three-dimensional solid problems [45,46].Moreover, the idea of using the complex polynomial basis functions had been successfully introduced into other meshless methods resulting in a series of new complex variable meshless methods, such as the complex variable boundary element-free method (CVBEFM) [47], the complex variable RKPM (CVRKPM) [48], the complex variable meshless local Petrov-Galerkin (CVMLPG) method [49], the complex variable moving Kriging interpolation (CVMKI) method [50], the complex variable interpolating moving leastsquares (CVIMLS) method [51] and the complex variable meshless manifold method (CVMMM) [52].These methods have been extensively applied in many engineering problems like the elastodynamic analysis of isotropic and functionally graded materials, the bending of thin plates on elastic foundations and the pattern transformation in hydrogels [47][48][49][50][51][52].
The complex variable methods play an important role in the early development of linear elastic fracture theory, as the analytical solutions of linear elastic fracture problems are often in the form of complex functions [53,54].However, almost all the previous numerical approaches aimed to establish the approximation for crack tip fields only using real variable functions.Obviously, real variable functions are not as accurate and efficient as complex variable functions to approximate the solutions of crack tip fields, which are originally in complex forms.Different with previous element-free methods, the complex variable meshless method applied in this work enables a direct construction of the complex form solutions for the linear elastic fracture mechanics.The CVMLS approximation method with a group of simple complex polynomial basis is applied for the meshless method to model the linear elastic fracture mechanics in a Cartesian coordinate plane.
Although the complex variable meshless methods can directly approximate the complex variable solutions at and around the crack tip region, it remains a great challenge to apply it to model arbitrary crack growth problems.From a thorough literature review, only very limited research works were found to apply the complex variable meshless methods to solve linear elastic fracture problems, e.g., predicting the stress intensity factor in elastics with a single stationary crack by Cheng and Li [37].Although the CVMLS approximation has flaws in the definition of its error norm, which could reduce the numerical accuracy [36,40], it still showed better accuracy and higher efficiency than the traditional meshless method in the computation of singular stress field at a crack tip [37].From these early encouraging results, we believe that the ICVMLS approximation with a correctly defined error functional will yield high accurate simulation results for fracture problems [36,40].Extending the complex variable meshless methods from modelling the stationary cracks problems to tracking arbitrary crack growths, is valuable and challenge.In this paper, an improved CVEFG (ICVEFG) method that combines the ICVMLS approximation and the Galerkin weak form of governing equations is developed.Moreover, the proposed ICVEFG method is, for the first time, successfully tracking arbitrary crack propagation path for linear elastic fracture problems.Compared with the conventional MLS approximation, the ICVMLS approximation applied in the ICVEFG method can substantially reduce the computational costs in tracking arbitrary crack paths.Applying an accurate definition of error norm functional ensures the ICVEFG method to yield higher accuracy than the CVEFG method in computing crack tip stress fields and tracking arbitrary crack paths.The effectiveness, robustness, and efficiency of applying this proposed ICVEFG method with both original and conjugate forms of the complex variable polynomial basis functions to track arbitrary crack paths are studied and approved with five different fracture problems in this work.
The outline of this paper is as follows.In Section 2, the ICVMLS approximation method is reviewed.In Section 3, the complex variable meshless method for two-dimensional fracture problems of elasticity is derived.In the derivation, a basis-function-enriching scheme for the singular field of crack tips, criteria for the discontinuity field at crack tips, criteria for the directions of cracks, and computational methods for related fracture parameters are introduced.In Section 4, five numerical benchmark examples are presented.

Brief review of ICVMLS
One of the key steps in the implementation of complex variable meshless methods is to construct the trial functions using complex polynomial bases [35][36][37][38][39][40].Generally, a complex variable trial function ũh (x) for a meshless method within a local support domain of point x (x 1 , x 2 , x 3 ) is expressed in terms of two arbitrary independent field variables, u k (x) and u l (x) as in the following form, where p T (x) = ( p 1 (x), p 2 (x), . . ., p m (x)) is the vector of complex polynomial basis functions, and a T (x) = (a 1 (x), a 2 (x), . . ., a m (x)) is the vector of coefficients for each term in the basis, which are also in complex forms.m is the number of terms in the basis and i is the imaginary unit.The construction meaning of Eq. ( 1) is shown in Fig. 1, in which an arbitrary unknown function f (x k , x l ) defined in a Cartesian plane is completely represented by an equivalent function in a complex plane.In the previous works [35][36][37][38][39][40][41][42][43][44][45][46], the trial functions for the complex variable meshless methods are always constructed by a pair of two arbitrary real variables in a complex form.These meshless methods only benefit from the complex variable methods with reduced number of trial functions in the modelling.As considering that the analytical solution of a crack tip field is originally in the form of complex variables, in this work, ICVMLS is directly applied to construct the complex variable trial functions to establish a complex variable element-free model for linear elastic fracture problems.In doing so, this modelling method not only has reduced number of trial functions but also will benefit with high accuracy and efficiency due to a nature complex variable approximation for linear elastic fracture problems.
The complex polynomial bases are constructed as follows [35][36][37][38][39][40]: Linear basis: Quadratic basis for 2D: There is more than one form of the quadratic complex polynomial bases for 3D problems, one of which is given as: Note, most of the complex variable trial functions for meshless methods are constructed based on the above original form of complex polynomial bases.However, applying the conjugate form of complex polynomial bases, namely, pT (x) to model mechanics problems can yield better computational accuracy [40] than using the original ones.The inconsistent numerical precision of these two types of complex polynomial bases in modelling of different problems is because that the complex polynomial basis often failures to accomplish the completeness requirement for trial functions.Thereby, the accuracy of applying ICVMLS to track arbitrary crack paths using different forms of complex variable polynomial basis needs to be further studied and examined.Compared with the original polynomial basis, whether applying a complex polynomial basis or a conjugate one, the number of terms in the admissible functions is greatly reduced.This is in particular true when high order basis function is used.To approximate high order singularity of the stress field at the tip of a crack, the enriched high order basis is generally needed, therefore applying ICVMLS approximation in tracking crack paths can reduce computational costs, substantially.
The ICVMLS approximation uses a well-defined weighted complex variable error norm, which ensures that the complex form trial function can accurately approximate the two independent field variables with real and imaginary parts, simultaneously [40].This is an essential progress of complex variable meshless methods to ensure the accuracy of their approximation, especially in the modelling of the complicated crack tip fields.The weighted complex variable error norm is defined as, where w(x − x I ) is a weight function defined within the support domain of point x, x I (I = 1, 2, . . ., n) and is used to denote the nodes within the support domain.ũh (x I , x) and ũ(x I ) are the complex form trial and target functions of two arbitrary independent field variables u k (x I ) and u l (x I ) at point x I , respectively.In this work, a weight function based on the Student's t-distribution is applied [11], as given by, where β is a parameter defined for controlling the shape of weight function, d I is the distance from a sampling point to a node, as d I = ∥x − x I ∥.d m is the size of influence domain and given by d m = d max d c , in which d max is a scaling parameter, and d c is the characteristic nodal spacing distance.
When the above complex variable polynomial basis in the original form is used, the error norm is rewritten in the following matrix form as, where Minimization of J with respect to a(x) leads to [35][36][37][38][39][40], with The trial function of ICVMLS with complex variable polynomial basis in the original form is expressed as, where Φ(x) is the shape function expressed in the following form, Therefore, two arbitrary independent field variables for the ICVMLS are defined and given by, Similarly, in case that the conjugate complex variable polynomial basis needs to be used, the shape function in Eq. ( 17) is directly derived as [40], 3. Complex variable meshless method for 2D linear elastic crack growth problems

Trial functions in complex variable meshless method for the approximation of crack tip fields
The presence of cracks generally leads to discontinuous displacement fields and singular stress fields at the crack tip, both of which will inevitably give rise to challenges to the numerical modelling.Several techniques [8,9,12,14,18,19,24] have been developed to consider discontinuities of shape functions within local support domains in meshless methods.The transparency criterion [8,9] that can be directly applied to implement the complex variable meshless fracture analysis is applied in this work.As the meshless trial functions that can be easily p-refined, introducing a set of basis functions to enrich specific analytical solutions will be able to effectively simulate the singular stress field at the crack tip.In general, a full enrichment for the basis function is given by [10], where r and θ are distance and angle of a particular point away from the crack tip, respectively, as shown in Fig. 2.
Since the singular field is confined to the vicinity around the crack tip [10], the crack tip region is divided into an enriched zone, a bridging zone, and a non-enriched zone for the purpose of reducing computational costs.The trial function of complex variable meshless method, the complex polynomial basis with and without the full enrichment are adopted in the enriched and the non-enriched zone, respectively.The shape function in the crack tip region is expressed as, where φI (x) is the shape function with fully enriched basis, φ I (x) is the shape function with non-enriched basis, R is a ramp function to bridge and ensure a compatible displacement field of the two regions.Both of the following linear and quantic polynomials are used for ramp functions.
where λ = (r − r 1 )/(r 2 − r 1 ), r 1 and r 2 are the radius of the enriched zone and bridging zone, respectively.

Discretization of Galerkin weak form in meshless method
For two-dimensional linear elastic fracture problems, the independent field variables at an arbitrary point x (x 1 , x 2 ) are horizontal displacement u 1 and vertical displacement u 2 .From Eqs. ( 18), ( 19) and (22), their ICVMLS trial function is given by, and where The linear strain is derived and expressed as, where Then, with the constitutive matrix, the stress field can be computed and given by, where where E ′ and v ′ are the equivalent Young's Modulus and Poisson's ratio with respect to a specific two-dimensional stress and strain state.For a plain stress problem, E ′ = E and ν ′ = ν.For a plain strain problem, . E and ν are the Young's modulus and the Poisson's ratio of the material.Applying the penalty method to ensure the essential boundary conditions, the Galerkin weak-form formula for a two-dimensional linear elastic fracture problem is expressed as, where b is the body force, t is the surface traction, α is the penalty factor.From left to right, the items on the left side of the equal sign in the above formula are the variations of strain energy, work done by body force, work done by surface traction, and the penalty functional corresponding to the displacement boundary conditions u = ū on Γ u .Substituting Eqs. ( 24)-(31) into Eq.( 33), the final algebraic equation is written as, where K is the global stiffness matrix, K α is the global penalty stiffness matrix, U is the global displacement vector, F is the global force vector, F α is the global penalty force vector, and they are expressed as follows, and when the boundary conditions, e.g., u i = ūi , are applied, the corresponding S i is set as 1, otherwise it will be 0.

Simulation of crack growth
For two dimensional problems, the direction of in-plane crack extension can be determined using the maximum principal stress criterion [55].This theorem states that a crack will grow along the direction, which is perpendicular to the maximum hoop stress, σ θ θ .A crack often propagates along a path with the crack tip is locally under mode I, in other words, along where the hoop stress near the crack tip σ θ θ (r, θ) arrives a minimum, as, The crack propagation angle θ is obtained from ∂σ θ θ /∂θ = 0 and ∂ 2 σ θ θ /∂θ 2 < 0, which gives, where K I and K I I are stress intensity factors of mode I and mode II cracks, respectively.
To obtain stress intensity factors with high numerical precision, accurate computation of the stress field at the crack tip is critical for the prediction of crack propagation angle and path.A correctly defined error functional for ICVMLS that can ensure high accuracy of the complex variable meshless approximation is critical for modelling the crack problems.The M integral, which is also known as interaction integral [56], is used for the extraction of mixed-mode stress intensity factors in our numerical implementation.The M integral in a local domain of a crack tip, as shown in Fig. 3, is expressed as, )   i j ∂u (2)   i ∂u (1)   i where the superscripts (1) and (2) denote the actual state corresponding with given boundary condition and an auxiliary state, respectively.W (1,2) = (σ (1)  i j ε (2)  i j + σ (2)  i j ε (1)  i j ) / 2 is the mutual strain energy density of the elastic body.The function q takes the form of where 2 ) is the coordinates of the crack tip.Stress intensity factors K I and K I I can be determined as, where E ′ is the same as in Eq. ( 32), M (1,I) and M (1,II) are M integrals obtained by choosing mode I and model II solutions as auxiliary states in Eq. ( 42), respectively.

Numerical examples
In this section, one example of an elastic body with static crack and four benchmark examples of quasi-static crack growths under different conditions are investigated using the present complex variable meshless methods with ICVMLS approximation.Both the complex polynomial basis in original and conjugate forms are examined in the ICVMLS implementation.To distinguish these two types of approaches, in the following context of this paper, we refer ICVMLS approximation meshless method with the original complex variable polynomial basis as ICV(I), and use ICV(II) to denote the ICVMLS approximation meshless method with the conjugate form of complex variable polynomial basis.Numerical results obtained from the EFG method (referred as MLS), the CVEFG method (referred as CV), and other published works are compared with the results obtained using the proposed ICV(I) and ICV(II) methods.In the implementation of these meshless methods, the penalty method is used to enforce essential boundary conditions for each problem.The rectangular support/influence domain is adopted in the meshless analysis.Four-point Gaussian integration along each direction of a single quadrature cell is used for the computation of weakform variational integrals and the M integrals.The transparency method is used for discontinuity processing in the modelling.A fixed incremental length of crack ∆a in each step is assumed for each case with evolutionary cracks.The computational parameters involved in these examples include the uniform meshes N x × N y and n x × n y for generating nodes and quadrature cells, respectively, the scaling factor d max controlling the size of local influence domain of nodes, the penalty factor α for enforcing displacement boundary conditions, the assumed crack-length increment ∆a, the size of the local domain L x × L y and uniform quadrature mesh m x × m y for evaluating M integral.Both the corresponding linear polynomial basis with and without full enrichment are examined in the numerical analysis.The quantic ramp function in Eq. ( 23) is employed when the full enrichment is adopted in the local zone around the crack tip.
To indicate the types of meshless methods and associated computational parameters that are applied in each simulation, the obtained results in each example for evolutionary cracks are labelled as "ICV(I)/ICV(II)/CV/MLS-N x × N y -n x × n y -d max -α-∆a-L x × L y -m x × m y -L 1 /(L en -r 1 × r 2 )"."ICV(I)/ICV(II)/CV/MLS" indicates which method is used in the analysis."L 1 " stands for only the corresponding linear polynomial basis is employed in the shape function."L en -r 1 × r 2 " denotes that the full enrichment is adopted around the crack tip in the enriched zone with radius r 1 and the radius of bridging zone is r 2 .

A square plate with mode I crack on boundary
The first example considered in this study is a square plate with a stationary mode I crack, which initially occurs on its boundary, as shown in Fig. 4(a).The side length of plate is L = 16 mm.The length of the crack is a = 8 mm.The Young's modulus and Poisson's ratio of the plate are E = 2×10 5 MPa and ν = 0.286, respectively.The applied displacement field is given by, The analytical solutions of the stresses around the crack tip are where SIF K I = F I σ √ πa with F I = 1.12 − 0.231( a L ) + 10.55( a L ) 2 − 21.72( a L ) 3 + 30.39( a L ) 4 [57].In this example, the applied stress is assumed to be σ = 1 GPa.
The purpose of choosing this simple example is to approve the accuracy of the proposed complex variable meshless methods, by directly comparing against well-known analytical solutions.Uniformly distributed nodes and rectangular J integral domain as shown in Fig. 4(b) are used in each meshless analysis.The global enrichment using full basis is applied in the meshless analysis to study this example.The near-tip stresses obtained using the present complex variable meshless methods are presented and compared in Figs. 5 and 6, in which the displacements and stresses at points ahead of the crack tip and their relative errors are plotted, respectively.The complex variable meshless approaches with the complex polynomial basis either in original or conjugate forms achieved more accurate results for the approximations of the near-tip displacements and stresses than that of the MLS.This is attributed to the direct approximation of the complex variable field function at the crack tip by using the ICVMLS, which obtains higher accuracy than the real variable approximation of the MLS.Although the CV approach results in relatively small errors for displacement solutions compared with the ICV(I) and ICV(II) methods, its stress solutions are much poor.The poor performance of CV is attributed to the invalid definition of error norm function in its shape function.
To provide an overall comparison for the crack-tip stress field, the contour plots of the near-tip von Mises stresses given by analytical solution, MLS and ICV (I) are presented in Fig. 7.It can be seen that the contour lines of the von Mises stress in the MLS solution is not as smooth as those in the analytical solution and the ICV(I) solution.
The stress intensity factors (SIF) K I with a series of different crack lengths a are obtained using MLS, ICV (I), and ICV (II) meshless methods and listed in Table 1.The relative errors for the results given by MLS, ICV(I) and ICV(II) methods are also listed for the purpose of comparison.As the accuracy of the results given by CV method is very poor, CV solutions are not presented in Table 1.As shown in Table 1, all SIFs calculated by ICV(I) method for various crack lengths a are better than the MLS solutions, however, the results given by ICV(II) are worse than that the MLS solutions.It concludes that different complex variable meshless methods lead to different accuracy in approximating the crack tip field due to the use of different polynomial basis functions.Therefore, an appropriate     selection of complex variables polynomial basis function is necessary to analyse these two ICV(I) and (II) meshless methods.

An edge-crack plate under mixed-mode loading
The second example in this study is an edge-cracked plate, in which the bottom edge is fixed, and the top edge is subjected to a far-field unit shear stress τ = 1 MPa, as shown in Fig. 8(a).The dimensions of the plate are L = 7 mm, H = 16 mm.The initial crack length is a = 3.5 mm.Material parameters are Young's modulus E = 3 × 10 7 MPa and Poisson's ratio ν = 0.25.The plate is under the plane stress condition.The crack path for this example had been modelled by Rao et al. [58] with a FEM-meshless coupling approach and Yang et al. [59] with a non-matching finite element-scaled boundary finite element coupled method.
In this example, uniformly distributed grid nodes and rectangular M integral domain as shown in Fig. 8(b) are adopted in each meshless analysis.Firstly, the crack paths predicted by the ICV (I), ICV (II), CV and MLS meshless methods with linear polynomial basis are presented in Fig. 9, in which the influences of applying different values of scaling factor d max on the results are illustrated.To show more details of the predicted crack paths, a 10-times (along y direction) enlarged image is shown in each diagram.From the numerical results shown in Fig. 9 (with currently chosen parameters), when the same uniform grid nodes and the same value of d max are used, the crack path predicted by ICV (I) is more in line with the results given in [58,59] than the solutions given by other methods.Moreover, the results given by ICV (I) and (II) are much better than that of CV.This is due to the use of weighted error norm given by Eq. ( 5) for ICV (I) and (II) methods, which can yield much better accuracy than the one adopted in the CV method.However, when applying the same uniform nodes, a much large value of d max , and a small size of M integral domain, are needed for obtaining accurate results for ICV (II) method in this case.On the other hand, the results obtained with the use of linear basis are all deviated from the results given in [58,59].This is mainly because that the linear basis cannot accurately catch the singular stress field of the crack tip.At the same time, it also can be observed that the crack paths predicted by ICV (I) and (II) are closer to the results given by [58,59] than the MLS solution when the same computing parameters are applied.
Secondly, the effect of applying different values of the penalty factor in these four types of meshless methods with linear polynomial basis is studied.The crack paths obtained by four different kinds of shape functions with penalty factors ranging from 10 3 × E to 10 6 × E are presented in Fig. 10.From the results, it can be seen that all the methods, except CV, are not sensitive to the change of penalty factors in this case study.When linear basis functions are applied, the crack paths predicted by ICV (I) and (II) methods are closer to the literature results than the predictions given by the other two methods.
Subsequently, the meshless analysis using different sizes of M integral domain is performed to test its influence on the predicted crack paths.Numerical results given by these four meshless methods using linear polynomial basis and different sizes of local domain for M integrals are presented in Fig. 11.Although the crack paths predicted by these four methods are varying with the size of the M integral domain, the simulation results of ICV (I) and (II) are more stable than that given by other meshless methods.Moreover, as shown in Fig. 11, when M integral domain size is appropriately selected, the ICV(II) is the most robust method that results in the closest solution compared to the published results in literature.
The last parameter considered in this study is the assumed crack-length increment ∆a.Similarly, the crack paths predicted by these four different meshless methods using linear polynomial basis and different ∆a are obtained and presented in Fig. 12. Regarding to the numerical accuracy and modelling robustness, under a same assumed crack-length increment ∆a, ICV (II) is the best, CV is the worst, and ICV (I) is slightly better than MLS.
To further improve the simulation results for the crack-tip singular stress field, the full enrichment basis function is introduced to construct the crack tip shape function.The predicted crack propagation paths are presented in Fig. 13.The results, as shown in Fig. 13, are obtained based on the consideration of the above parameter analysis, from which the most appropriate parameters are selected.Although we have examined a large number of computations with various parameter combinations, CV cannot give any meaningful results, the crack path predicted by CV is omitted in Fig. 13.The present results match with the literature results much closer than the above numerical results obtained only using linear polynomial basis functions.The numerical results obtained by both MLS and ICV(I) methods are close to the results given in [59], and the ICV(I) are relatively better than MLS.Although the results given by ICV(II) is closer to the one in [58], ICV(II) needs more nodes to achieve convergent solutions.In terms of computational efficiency, the CPU time of ICV(I), MLS and ICV(II) is 708 s, 1110 s and 3886 s, respectively.Therefore, ICV(I) method has 36% higher computational efficiency than the MLS method.As ICV(II) needs more nodes for modelling this problem than that of other methods, its computational efficiency is greatly reduced.
The effects of applying different weight functions and different basis functions in ICV (I) on the predicted crack paths are demonstrated in Figs. 14 and 15, respectively.From the numerical results shown in Fig. 14, the paths obtained by Gaussian and cubic splines are worse than those obtained by other splines.According to Fig. 15, the crack path obtained by quadratic basis is more stable than that given by linear polynomial basis.Moreover, the crack paths obtained by a rectangular domain and a circular domain for M integral are presented in Fig. 16.Both the length and the width of the rectangular domain are set to 2, and the radius of the circular domain is set to 1, other parameters are the same with those given in Fig. 13.The results show that the shape of M integral domain has insignificant impact on crack path.The deformed nodes during the crack propagation process are illustrated in Fig. 17, which clearly shows the crack path.

A double cantilever beam
The third example is to study a double cantilever beam (DCB), which is fixed at the right boundary and has a kinked crack at its left edge, as shown in Fig. 18.The beam is subjected to a pair of opposite concentrated loads F = 100 N on the left edge.The length of the beam is L = 300 mm, the height is H = 100 mm.The kinked crack consists of two segments, one segment is horizontal with 138 mm long and the other is inclined with an angle θ = 2.6 • and 12 mm long.The material properties of the plate are Young's modulus E = 2 × 10 5 MPa and Poisson's ratio ν = 0.3.This double cantilever beam is assumed under the plane stress condition.The experimental results of crack propagation path for this problem had been provided by Sumi et al. [60].Ventura et al. [12], Ooi et al. [61] and Yang et al. [59] have conducted numerical simulation to predict the crack propagation path for this  problem using an EFG approach with a new vector level set method, a polygon scaled boundary finite element method, and a non-matching finite element-scaled boundary finite element coupled method, respectively.All these published results are used to compare the present numerical solutions in this paper.Note, there is a significant difference between the crack paths from these three numerical studies and they all have certain amount of deviation from the experimental results.
For this DCB example, the crack path is predicted using MLS, ICV (I), ICV (II) and CV methods, and their numerical results are presented and compared with each other in Fig. 20.The domain discretization is based on uniformly distributed nodes, as shown in Fig. 19.Several groups of computational parameters are tested, with which appropriate parameters are chosen for obtaining accurate modelling results.The crack paths predicted by these meshless methods are shown in Fig. 20 and compared with the experimental results.To reveal more details for the crack paths, an enlarged image (7 times along x direction and 2 times along y direction), is presented in Fig. 20(b).From the parametric study of ICV(II) and CV methods for applying to solve this DCB example, it is difficult to obtain accurate results when a linear basis function with or without full enrichment is used.For the MLS method, accurate prediction of crack path can only be achieved when the linear basis function is applied.Nevertheless, for  ICV (I) method, no matter which kind of basis functions is applied, accurate solutions can always be obtained.As shown in Fig. 20, the crack paths predicted by both of the MLS method and the ICV(I) method with linear basis are in good agreement with the numerical results obtained by Ventura et al. [12] and Ooi et al. [61].Moreover, the crack path predicted by ICV(I) matched very well with the experimental results given in [60], and it is much better than the numerical results from Yang et al. [59].When using linear basis, the CPU time for the computational process of MLS and ICV(I) is 1136 s and 457 s, respectively.Compared with MLS, the computational cost of ICV(I) is reduced by 59.7%.The pattern of deformed nodes during the crack propagation process is shown in Fig. 21.

A single-edge notched shear beam
The fourth example is to study a single-edge notched beam subjected to two concentrated loads, as shown in Fig. 22.The length of the beam is L = 916 mm, the height is H = 306 mm, and the thickness is U = 152 mm.The length of the vertical crack is 82 mm.The magnitude of the load is F = 112.1 kN (the other load on top-right corner is 0.13 F).The material properties for this beam are Young's modulus E = 24.8GPa and Poisson's ratio ν = 0.18.This problem is considered under the plane stress condition.This beam crack problem has been widely accepted   as a benchmark for the validation of mixed-mode crack propagation models.Yang and Proverbs had conducted numerical simulation for this problem using a nonlinear finite element approach [62].Yang et al. had analysed this problem with a non-matching finite element-scaled boundary finite element coupled method [59].Xie et al. [63] proposed an energy-based approach for the finite-element modelling of this problem, and their result is slightly different from that of the other two approaches.The crack paths provided in all these references are used for the purpose of comparison in this section.The uniformly distributed grid nodes, as shown in Fig. 23, are used in these four meshless methods to model this problem.Fig. 24 shows (a) the predicted crack growth paths; and (b) the enlarged image (5.5 times in x direction) obtained by these four different meshless methods using linear basis function.The results given by MLS, ICV(I) and (II) are obviously better than that obtained by CV when linear basis functions and the same computational parameters are used.The results obtained by MLS and ICV(I) are the best.The predicted crack paths and the enlarged image (5.5 times in x direction) that are obtained from these four meshless methods using linear basis with full enrichment  are presented in Fig. 25(a) and (b), respectively.The crack path predicted by ICV(I) is obviously more consistent with the published results than that obtained by the MLS method when linear basis with full enrichment, when the same computational parameters are applied.The load-crack mouth sliding displacement (CMSD) curve obtained with MLS and ICV(I) are compared with the experimental results and numerical solutions in literatures, as shown in Fig. 26.It shows that the ICV(I) leads to the results that match better with the experimental results than the other methods.Fig. 27 shows the nodal deformation after the crack propagation.

A double-edge notched shear concrete beam
The fifth example is a double-edge notched beam subjected to two concentrated loads, as shown in Fig. 28.The magnitude of the first load is F (the other load at the top-left corner is F/5).The length of the beam is L = 800 mm, the height is H = 200 mm.There are two vertical cracks in the upper mid-span and lower mid-span of the beam.The positions and the lengths of the two initial cracks are shown in Fig. 25.The material parameters of this beam are Young's modulus E = 27 GPa and Poisson's ratio ν = 0.1.The plane stress condition is considered for this problem.The experimental work for this problem had been done by Bocca et al. and presented in [64].The numerical simulations using different computational methods can be found in many previous published literatures [14,[65][66][67].Clear crack paths provided in [67] for this example are used to compare and validate the results obtained from the meshless methods.Note, the crack path from the experimental results originally obtained in [64] is taken from [67].
A 39 × uniformly distribution of grid nodes for the domain discretization, as shown in Fig. 29, is applied in the meshless analysis of this problem series of parametric studies.Similar with the previous case studies, both the linear basis functions with and full enrichment are employed in these four meshless methods to predict     the crack path for this double-edge notched beam problem.Fig. 30 shows the predicted crack paths using these four meshless methods with linear basis functions and different values of d max , compared with the experimental results.Both the crack paths predicted by MLS and ICV(I) methods converge very well within the region of experimental results given by [64], and the crack path of ICV(I) gives better results than MLS.The crack paths predicted by ICV(II) are slightly away from the experimental results, albeit the results are convergent eventually.Even with use of the complex variable approximation, the incorrectly defined weight complex error norm obviously damages the accuracy of CV.A larger support domain needs to be used, that is, d max should be larger than 3.5 to obtain appropriate results for CV.Overall, the prediction of the crack path by ICV(I) is the most accurate and efficient, followed by MLS, ICV(II) and CV when linear basis function and the same computational parameters are applied.
Fig. 31 shows the crack paths predicted by these four meshless methods with enriched basis functions and different values of d max .From the numerical results, it can be seen that ICV(I) remains to be the best method in prediction of the crack paths and follows by the MLS method.The crack paths obtained by ICV(II) become more unstable compared with the simulation results given by other methods when linear basis function with full enrichment is used.CV can only generate reasonable results when d max = 4.0.The incorrectly defined error norm in CV not only makes the numerical accuracy worse, but also indicates that the original intention of using complex polynomial basis function to reduce the minimum number of nodes in the support domain to save computational time is no longer tenable.The total load vs. the vertical displacement of the right loading point is plotted in Fig. 32 and compared with other results obtained in previous works [14,64,66].It shows that the results obtained by ICV(I) have a good agreement with experiment.Fig. 33 shows a nodal deformation pattern during the crack propagation process.

Conclusions
In this paper, the complex variable meshless methods, for the first time, are applied to model two-dimensional crack propagation problems under the theoretical framework of linear elastic fracture mechanics.The main motivation and innovation of this work is to directly construct an element-free complex variable approximation of the target function, which is originally in complex form for the linear elastic fracture problems.A numerical framework for the complex variable meshless method that applies the CVMLS approximation with a group of simple complex variable polynomial basis function and accurately defined error functionals is established to model the crack tip fields and arbitrary crack path with improved accuracy and efficiency.Five benchmark examples are studied using four different meshless methods (ICV(I), ICV(II), MLS, and CV) with linear basis function (with/without enrichment) and different computational parameters.The numerical results and predicted crack paths obtained from these four meshless methods are compared with each other and compared with the experimental and numerical results published in previous research works.In general, the proposed improved complex variable meshless methods, especially the ICV(I), have been approved to provide great simulation precision and high computational efficiency in modelling the crack propagation problems.To accurately model the progressive cracks problems, employing a precise definition of complex variable error norm in shape function is critical and the form of complex polynomial basis plays an important role in obtaining stable and accurate numerical results.Although this work successfully applies the complex variable meshless method to model crack problems, the extended terms in the enriched basis function adopted by the present method are still expressed in terms of real variables.Therefore, this proposed method is not a full complex variable method.In future work, we will develop a full complex variable meshless method and apply it to modelling three-dimensional fracture problems.

Declaration of competing interest
The authors declare the following financial interests/personal relationships which may be considered as potential competing interests: D.M. Li reports financial support was provided by EU Framework Programme for Research and Innovation Marie Sklodowska-Curie Actions.

Fig. 1 .
Fig. 1.Arbitrary function on a Cartesian plane and a Complex plane.

Fig. 3 .
Fig. 3. Local domain for M integral at a crack tip.

Fig. 4 .
Fig. 4. (a) A square plate with a stationary mode I crack; (b) Uniform mesh 20 × 20 for nodes and J integral domain.

Fig. 5 .
Fig. 5. (a) Results of u for points located ahead of crack tip and (b) relative errors.

Fig. 6 .
Fig. 6.(a) Results of σ x for points located ahead of crack tip and (b) relative errors.

Fig. 8 .
Fig. 8. (a) An edge-cracked plate subjected to a unit shear traction; (b) An example for uniform mesh 13 × 31 for nodes and 2 × 2 M integral domain for the first step.

Fig. 16 .
Fig.16.Crack propagation paths by ICV (I) using rectangular and circular domain for M integral.

Fig. 17 .
Fig. 17.Deformation of nodes during the propagation of crack.

Fig. 21 .
Fig. 21.Deformation of nodes during the propagation of crack.

Fig. 25 .
Fig. 25.Crack propagation paths with different meshless methods using linear basis with full enrichment.

Fig. 27 .
Fig. 27.Deformation of nodes during the propagation of crack.

Fig. 33 .
Fig. 33.Deformation of nodes during the propagation of crack.

Table 1
Solutions of SIF K I with different crack length a.