Using micro-planes concrete damage model in discrete least squares meshless method for predicting the crack growth

– In the present study, the author’s previously developed micro-planes concrete damage model has been implemented in the Discrete Least Squares Meshless method (DLSM) for assessment of the crack performance in concrete domains. DLSM is a true meshless method needing no partition in the local approximation as well as in the assembling of the local approximations in order to obtain a global solution. Instead it uses a collection of the un-structured nodal points in place of the elements and discretizes the governing diﬀerential equations with the discrete least squares approach. This can reduce considerably the pre-processing eﬀorts of the calculations. However, DLSM like other meshless methods loses its eﬃciency when it is applied for the analysis of the solid domains dealing with the growing cracks. In recent years, some researchers have utilized some techniques to overcome this problem such as visibility, transparency and diﬀraction for precisely detecting the crack discontinuity in the case of a static non-growing crack, but these methods require some reciprocal corrective operations within each step of the calculation process extending the time of processing. In this study, however, a new, straightforward and general applicable method has been used for accommodating this issue using the advantages of the Moving Least Squares (MLS) shape functions constructed based on the Voronoi tessellation algorithm and micro-planes concrete damage constitutive model. Voronoi based MLS shape functions have been implemented for evaluating the crack eﬀects accurately as a discontinuity during the analysis processing and micro-planes concrete damage approach used for proper detecting the crack growth trajectory. In the proposed method, the domain of interest is divided into the Voronoi cells based on the classical Voronoi tessellation approach, a cell for each node, and then the support domain of each node located inside the domain as well as on or near the crack faces is determined based on a neighboring criterion. During each step of the analysis for all the integration points, damages are computed according to the micro-planes concrete damage model and the crack path for the next step will be deﬁned. After the end of the iterations in each step, both the stiﬀness tensor and internal load vector are modiﬁed based on the ﬁnal stress ﬁelds obtained on and near the crack. The accuracy and eﬃciency of the proposed method have been investigated by solving some benchmark examples and comparing the obtained results with the analytical or valid ﬁnite element analyses’ results.

In general, whenever meshless methods are employed to predict the crack growth in an engineering solid material like a concrete structure, fail in giving proper and accurate results.This difficulty comes from an inherent deficiency of these methods; they cannot precisely determine the support domain of the nodes located in the vicinity of the crack tip or crack faces.The reason is obvious; these methods are so formulated to be mesh-independent.Therefore, they cannot not utilize the advantages of the elements like to the well-known Finite Element Methods (FEMs) in which the faces of the crack must always coincide exactly with the boundaries before beginning of each step of the analysis.This means that, in FEMs, the continuum space must be well established before starting each step of the calculations.However, in contrast to FEMs, meshfree methods have no such limitation and construct the support domain of each integration point during the process of the analysis not prior to the beginning of the computation.Some techniques have been implemented by researchers in recent years such as visibility [11], diffraction or transparency [12][13][14] in the meshless methods to overcome such deficiency but each of them has own specific shortcomings and they do not possess a wide applicability, for instance, the visibility technique works based on the light transport theory.In this approach, a discontinuity is distinguished as a part of a medium which light cannot travel through it.With this method, the concave boundaries are detected but the continuity condition of the weight and obtained shape functions is violated near the crack tip.Two other mentioned techniques (diffraction and transparency) do not include the same drawback but in these methods, there exist some restrictions on the position of the nodes limiting their applicability and extending considerably the time of the calculations because of their reciprocal nature of the operations.Pirali and his co-workers described these techniques and mentioned their deficiencies [15].A more detail review of these techniques has been done in section two, helps in better describing the alternative approach suggested in this paper.
In this study, a new and straightforward approach has been implemented in the framework of DLSM for predicting the crack growth in a concrete solid medium.The scheme has been developed based on the combination of two theories; the classical Voronoi tessellation algorithm and the micro-planes concrete damage material model.The mentioned concrete material model had been developed before by the author in 2006 based on the combination of the damage and micro-planes theories [16].The latter will be explained briefly in section three.Section two discusses the Voronoi algorithm used for constructing the shape functions.In fact, it can be said that the Voronoi algorithm in this study has been proposed for detecting the crack as a discontinuity and incorporates its effects into the local approximations in the vicinity of the crack in DLSM framework and at the same time, the micro-planes concrete damage model has been applied for diagnosing the path of the progressive crack during the loading procedure.Equations govern the task of the crack growth in a concrete medium with an isotropic elastic damage material behavior and are formulated in the framework of the DLSM in section four.Section five deals with the description of the solution strategy based on the theories explained through sections two, three and four.Finally, two benchmark examples have been solved using the solution procedure of section five verifying the validity and accuracy of the proposed method.

Voronoi-based moving least square shape functions
In the proposed approach, crack detection and contribution of its effects in calculations are provided by implementing the concept of Voronoi tessellation algorithm used for dividing the domain of interest between some scattered nodal points.In this method the support domain of a node consists only of neighbor nodes.Neighbor nodes of a node are detected using Voronoi diagram algorithm.In this way, each node has a Voronoi cell.A Voronoi cell of a node j based on the Voronoi diagram algorithm is defined as the collection of the points that are located closer to the node j than any other node like node i (see Fig. 1).
In the mathematical formulation this definition can be stated as follows where V j is the Voronoi cell of the node j and d is the Euclidean distance between two points x and ξ i ∈R n defined as (2) Once the Voronoi cells are created, the nodes which have a common facet with the Voronoi cell of the node j are considered as the first layer of neighbor nodes (Fig. 1).We can define more than one layer of neighbor nodes if needed.
The second step after determining the support domain of the nodal points is to construct MLS (Moving Least Squares) shape functions defined as follows:  where, a(x) is the vector of coefficients and P (x) is the vector of basis functions; m is the order of basis functions and r is the total number of the basis function terms.In the proposed method, the number of the layers of neighbor nodes is equal to the order of the basis functions.
In other words, the value of m determines the number of neighbor layers used for constructing the support domain.For example, in Figure 2, two layers of neighbor nodes are selected for creating a second order basis function.This method guarantees the required number of the nodes for basis function calculations and there is no need to additional methods to determine the number of the nodes for the support domain.In this paper, the basis functions of the form P (x) = [1, x, y, x 2 , xy, y 2 ] (m = 2 and r = 6) have been used, so two layers of the neighbor nodal points have been employed.
Error or residual function "J" is then considered as below: u h j is the nodal value of the function to be approximated (in this study the plane displacements u and v) of the point j, and n s is the total number of nodes in the support domain of point j and W j is the weight function defined as below; where di = ||x − xj|| and d = d i /d wj in which d wj is the radius of the weight function of point j.d wj is set to be the distance of the farthest node of the support domain multiplied by a constant number which is considered in this study as 1.1 [19][20][21][22].Some considerations should be taken into account for excluding some points in the neighborhood of the crack from the support domain of the nodes like the node j near the crack which they are contributed in this study through the computing of d i .It should be recalled that by Voronoi approach used here, the nodes which are located in the support domain of the node j are not determined based on their straight distances from the node j.Instead, they are determined based on a neighborhood criterion.In other words, if the distance of a node i to the node j is smaller than the distance between the nodes j and k, the node i is not necessarily considered in the support domain of the node j.For example, for the case m = 1, node i in Figure 3 does not place in the j's support domain even if its distance from the node j is smaller than the distance of the node m from the node j.This is because that the Voronoi cell of the node i in Figure 3 has no common facet with j's Voronoi cell or in other words, the node i is not a neighbor node of the node j.Furthermore, note that the node k which is located behind the crack vertex is excluded from the support domain of node j even if it is located at a very close position to the node j, because the vertex which is located on the crack line is not considered as the common facet and as a result, the neighborhood between the nodes k and j will be violated.By this approach, for the case m = 1 (one layer of neighbor's Voronoi cells is considered), d i can be normally defined as the straight distance of node i from the node j and no additional operation is needed.However, for the case m ≥ 2 (two layers of neighbor's Voronoi cells are considered), the distance d i cannot be measured along a straight line connecting nodes i and j because this line may pass through a cell which is not a neighbor cell of node j.This problem can be fixed by imposing a constraint.This constraint is that for computing d i at first a path must be defined between nodes i and j.This path which connects the node i to the node j must always pass through the neighbor cells of the node j.By this way, for computing d i in Figure 4, the path which connects two nodes i and j passes from a middle node k.According to this figure, d i is the summation of P 1 and P 2 .
So d i is the summation of P 1 and P 2 and it can be written: It can be seen that the proposed technique in this study requires no geometric operations in comparison with classical Voronoi diagram like as the approach implemented by Aurenhammer et al. [21].
Before further continuing of the formulation, for better understanding the scheme which was described above, let us review briefly some techniques which have been implemented by researchers in recent years to incorporate the effects of the discontinuities like cracks into the analysis.These techniques may not be understood well by the readers before this stage of the study.
In visibility technique, the weight of node X I at point X is set to zero if the segment that joins them is cut by a crack.Hence, by this way, the weight function and subsequently the shape functions at the point X are modified considering the discontinuity opaque for rays of light coming from the node X I (see Fig. 5) [15].In other words, by this approach, a discontinuity is distinguished as a part of a medium which light cannot travel through it.With this method, the concave boundaries are detected but the continuity condition of the weight and obtained shape functions is violated near the crack tip (yields spurious discontinuity) because that at the crack tip, an artificial discontinuity inside the domain is constructed as shown in Figure 5 [22].
In diffraction method, if the segment that joins node X I to point X is cut by a crack, the normalized distance d which it was used in relation (6) in this study, is defined according to the length of the shortest path from node X I to point X passing by the point of the crack tip X C : X I −X−X C In mathematical form this can be written as: (see Fig. 5).The weight function of node X I is then continuous except across the crack.
By using transparency method, if the segment that joins node X I to point X is cut by a crack, the normalized distance d = di S0 in relation ( 6) is augmented by a transparency function that increases with the distance between the crack tip and the intersection of the segment with the crack.That is: f t is a strictly increasing function to choose (see Fig. 5).Transparency criterion gives sharply varying shape functions for the nodes near the crack [23].
From the above reviewing, it is revealed that these techniques are not guaranteeing the smoothness of the shape functions near the crack or crack tip and also are not straight-forward schemes.This is due to the fact that these methods do not construct support domain of the nodal points considering the effect of crack or other discontinuity from the first step of the analysis procedure and they are trying to modify the weight and shape functions after that the support domains are constructed.However, by the approach proposed in this study, the aforementioned deficiencies are removed completely even it may not seem to be simpler than those techniques.
Continuing the main discussion, the vector of coefficients is calculated by minimizing Equation (5) as follows: where u h is the vector of nodal values of displacements as follows: and A(x) and B(x) are defined as: Equation ( 9) is rewritten in the following form: where N T (x ) = P T (x ) A −1 (x ) B (x ) is transpose of the vectors of MLS shape functions.

Micro-planes concrete damage model
The approach used in this study for determining the crack growth trajectory within the solid during the analysis has been previously developed by the author based on the concept of detecting damages on particular planes called as micro-planes in each integration point [16].In this method, in each computational point, 26 microplanes are defined in the case of three dimensional problems.These planes are in fact the planes with different orientations which considered tangential to a unit sphere containing the point of interest (integration point) and the point in fact is located at the center (see Fig. 6).
The position and orientation of each micro-plane is specified with a unit outward normal n to that plane with components of n i , i = 1, 2, 3 (any subscript refers to the components in rectangular Cartesian coordinate axis x i ) (see Tab. 1).Also to extract shear components on the micro-planes, we are required to define two extra coordinate directions M, L which represent two orthogonal unit coordinate vectors m i , l i respectively.
According to the kinematic constraint approach, at first macroscopic strain tensor is projected on the planes.This projection yields three components of micro-strains along the plane's triplet local directions which one of them is normal (ε N ) and the others are tangential (ε M , ε L ) (Fig. 7).
The following relations describe the projection process in the mathematical form: Table 1.Micro-planes orientations and their weights on the unit sphere.
According to these obtained micro-strain components on each micro-plane, the following isotropic scalar damages are defined for concrete: For hydrostatic tension state of strain (17b) For uniaxial-compression state of strain (17d) For uniaxial-tension state of strain (17e) In the above relations, ω represents a scalar isotropic damage which its value ranges between zero (un-cracked) and unit (fully cracked) states in an integration point.Subscripts HC, HT, SH, C and T represent hydrostatic compression, hydrostatic tension, pure shear, uniaxial compression and uniaxial tension states of micro-strains' path respectively.Parameters a to k in relations ( 17) can be defined solely based on the standard tests on the concrete specimens e.g.uniaxial compression and uniaxial tension in laboratory.For determining these parameters, curve fitting method used based on the minimization of the summation of the squared differences between calculated stresses from the proposed micro-plane damage concrete model and experimental stresses derived from mentioned standard tests.Other damages from other stress-strain paths can be defined from obtained damages from uniaxial-compression and uniaxial-tension states of strain according to relations (17).In this study, the so-called parameters are obtained as: a = 0.00001, b = 0.0002, c = 0.00015, d = 300, e = 0.002, f = 86000, g = 303, h = 0.6, i = 0.003, j = 0.000549, k = 0.003000001.ε eq is a scalar equivalent strain computed from the micro-strains on each plane.Once the five above mentioned damages are defined, the total damage factor (TDF) on each micro-plane can be determined using the following relation; TDF(P ) = HT (P ) × ω HT (P ) + T (P ) × ω T (P ) + SH(P )ω SH (P ) + C(P ) × ω C (P ) + HC(P ) × ω HC (P ) (18) where P indicates the micro-planes number and its value changes between 1 to 13.Based on this calculated TDF, the damaged elastic modulus on each micro-plane is computed as follows: in which E (P ) is the damaged and E 0 is the initial uncracked elastic modulus respectively.According to these damaged elastic modulus, the total damaged elastic tensor, (D ijkl ), of each integration point is calculated from the E (P ) using the below equation: (20) or in the numerical form: In relations ( 21) and ( 22), D p ijkl is the deviatoric part of the damaged elastic tensor of each micro-plane.The second term in the right hand side of Equation ( 21) denotes the volumetric part of the damaged elastic tensor of the integration point.
Based on the above relations, when the total damage factor in an integration point reaches to unity, it is concluded that this point is cracked.So, the stiffness of this point must be adjusted and recomputed using the relation (21).Also the counterpart of this point in the internal load vector must be determined accordingly.For doing this, the updated stress tensor in every integration point is needed.Relation (23) helps to compute such stress: Please note that for deriving the above relation, by implementing the virtual work principle, the equilibrium conditions in the integration point are enforced in relation (23), Ω is the unit hemi-sphere surface, σ L and σ M are tangential stress components on each plane, σ V and σ D are volumetric and deviatoric parts of the normal micro-stress components, s ij is the deviatoric stress tensor on each micro-plane, w μ is the weight of the micro-plane defined in Table 1, N m is the number of the micro-planes on hemisphere which its value is equal to 13.The previous formulations were presented in a general 3D form and can be easily transformed to 2D form for application in this study.Hence, by implementing this approach, the path of the crack growth in the concrete domain can be determined and followed.

Discrete least square meshless formulation of the governing equations
The equations governing the problem of crack growth in a solid medium can be stated in general form as following: where, A () is a differential operator defined according to the physical phenomenon of the problem in hand.f is recognized as the body force vector, φ denotes the unknown vector (in this study the displacement vector) and Ω represents the domain of the problem.Accompanying to the Equation ( 24) which defines the governing equation within the domain there exists Equation ( 25) describing the governing equation on the borders of the domain problem as follows: in which, B () analogously marks a differential operator, t defines the traction force vector, Γ t denotes the natural boundary part, φ is the prescribed known values of the displacement vector and Γ u designates the essential boundary part.Let N Ω be the number of the nodal points within the domain, N u be the corresponding quantity located on the essential boundary N t be the same parameter but defined on the natural boundary part of the domain and finally N be the total number of the nodal points.Then, the unknown function (displacement function in this study) in a representative point of the domain can be approximated as follows using the well-known trial shape function approach for discretized domains: n is the number of the nodes trapped in the support domain of the x k , N i denotes the trial shape functions and φ i represents the values of the unknown function at the nodal points of the support domain.Substituting Equation (26) in Equations ( 24) and ( 25), the residual can be obtained as below: Equation ( 27) defines the residual of the response within the domain and Equations ( 28a) and (28b) represent the corresponding quantity on the natural and essential boundaries respectively.According to the proposed enhanced DLSM method in this study, a functional was defined as below: In the above equation, α t and α u are the penalty parameters for Neumann and Dirichlet boundary conditions respectively.By minimizing the above squared residual function with respect to the unknown parameters φ i , a system of descretized equations as follows is obtained: which can be rewritten as: in which: Note that in the proposed DLSM method, K ij would obtain as a symmetric and squared matrix.For 2D elastic domains; plane stress; the governing equation in displacement form is as follows: is the Laplacian operator.The Dirichlet and Neumann boundary conditions are defined as below: Rewriting Equations ( 24) and (25) using Equation (34) yields: Using Equations ( 36) to (38), the system of Equations (31) can be constructed and solved for unknowns φ i .

Solution algorithm
The displacement control scheme has been implemented for following the load scenario on the solved examples in this study in order to simulate the behavior of the concrete after its peak strain value is reached.The displacement is applied in a smoothed incremental manner to capture the nonlinearity behavior of the concrete due to crack initiation and growth.
In the first displacement step, the behavior of concrete is assumed to be un-cracked linear elastic.According to this assumption and following the formulations of Sections 2 and 4, the displacement field of the domain is computed.Based on this obtained displacement field, the strain tensor is calculated for each integration point.This strain tensor is then projected on the micro-planes and following the relations in Section 3, the damages and corresponding TDF are obtained.If TDF is less than unity in all points, the structure is considered as un-cracked and there is no need for extra calculations in this stage and subsequent displacement step will be applied.However, if TDF is reached to unity in a point, the point must be classified as a cracked point and additional process is required.Since, it is possible that this point is cracked under a smaller amount of external displacement, hence, for finding the exact instant of crack initiation, the displacement increment is divided by two and this new increment is applied.This process is followed again and again until that mentioned instant is captured.The behavior of concrete structure after crack initiation becomes nonlinear, so after cracking the stiffness tensor and internal load vector of this point must be re-computed according to the obtained cracked elastic tensor and cracked stress tensor (relations ( 21) and ( 23)).Accordingly, the total stiffness matrix and external force vector are then computed using relations (32) and (33) respectively and convergence criterion is controlled.In this study, the convergence criterion  is satisfied when the difference between the external load vector and internal computed load vector becomes less than a specific user defined error.If this criterion does not meet, the next iteration is followed in the current step of calculations.In the next iteration, according to the difference value between internal and external load, the internal displacement vector in the nodal points and the strain and stress tensors in the integration points are updated and the convergence is re-controlled.The iterations are so followed in current step until the convergence is satisfied.At the beginning of the next step, according to the fully damaged points recognized in the current step, the crack path is re-defined and then based on Voronoi algorithm explained in Section 2, that new crack line is incorporated in the calculations.The other steps are preceded until the un-stability appears in the calculations.

Numerical examples
In this section, two benchmark examples were solved following the solution procedure in section five; to verify the validity and accuracy of the presented method in predicting the crack growth path in a concrete medium.As a first example, a rectangular concrete beam with an edge pre-existing notch is considered (Fig. 8).This example representative for the fracture in mode I.

Concrete beam with an edge single notch
In this problem the Young modulus and Poisson's ratio were specified as E = 20 GPa and ϑ = 0.2 respectively.The uniaxial tensile and compressive strength of the concrete are considered as 2.4 and 23 MPa respectively.Figure 9 shows the discretization of the beam with the proposed method.
For verification of the presented model under the static non-growing crack, a 0.05 mm vertical displacement  was first imposed on the middle span of the beam on its upper edge.Figure 10 represents the developed stresses in horizontal (x-axis) direction obtained with the current model and that obtained with the finite element method.
As it can be seen, there is a good agreement between the results.
For controlling the crack mouth opening performance during the analysis, the above applied displacement was increased gradually.The initial displacement step was chosen as 0.05 mm and was increased up to 0.25 mm. Figure 11 illustrates the crack tip openings in the five mentioned steps obtained from the proposed method.Crack growth performance predicted with the current model has been depicted in Figure 12.The predicted crack tip opening and crack trajectory with the proposed method is in good agreement with the experimental evidence [17].
The obtained load-displacement curve from the current method has been compared with the experimental curve obtained from reference [17] in Figure 13.The vertical axis shows the values of the load which are scaled with respect to 100 kN, the horizontal axis shows the values of the displacement in mm.As it can be seen from Figure 13, the proposed model can predict well the behavior up to peak load.In softening branch of the curve, the current method could not follow the experimental curve.This is because that in the experimental specimen, after reaching the mode I crack from the lower parts of the beam to the crashed zone in the upper parts, due to the existence of the friction between two fractured beam parts, the beam continues to undergo further deformations.However, in the mathematical proposed model, after fracturing the beam into two parts, the continuity condition becomes violated and the proposed numerical model becomes unstable.

Concrete double edge notched (DEN) specimen
To investigate the performance of the presented model in the case of developing of mode II cracks, a concrete specimen with double edge notches has been selected for analysis.The specimen geometry and test setup have been shown in Figure 14.
The load which is in the horizontal direction in Figure 14 was kept as 10 kN during the analysis, however the vertical load was simulated by a controlled increasing displacement.The first step value of the displacement was considered as 0.002 mm.This vertical displacement was gradually increased but the horizontal load was kept constant as mentioned before.In this problem, the Young modulus and Poisson's ratio were specified as E = 30 GPa and ϑ = 0.2 respectively.The uniaxial tensile and compressive strengths of the concrete are considered as 3 and 49 MPa respectively.Figure 15 shows the discretization of the beam with the proposed method.
In the modeling of the concrete DEN, 1643 unstructured scattered nodal points were used.The algorithm which was described in Section 5 was implemented for the analysis of this problem and the below crack propagation paths have been obtained.The experimental crack growth which is extracted from reference [18] has been brought in Figure 17 for comparison.
The vertical load-displacement curve of the presented model and that obtained from reference [18] were compared in Figure 18.The vertical axis shows the values of the vertical load which are measured in kN, the horizontal axis shows the values of the displacement which are measured in mm.It can be seen that the formulated approach in this study can relatively well predict the mode II crack performance in concrete structures.Furthermore, in comparison to the first example, the present model can undergo larger deformations in the softening region in the current example.This is due to that mentioned instability condition is not occurred in this case.

Conclusions
A new scheme has been adopted within the framework of discrete least squares meshless (DLSM) method to enhance its performance, concerning with a concrete domain contain crack growth.Meshless methods generally fail in proper diagnosis for the crack discontinuity due to their inherent nature.They do not implement predefined elements like as the F.E. based methods; rather use scattered (unstructured) nodal points in the domain for discretizing the continuum governing equations.Hence, meshless methods in contrast to F.E. methods must at first distinguish the crack properly as a discontinuous medium and then start to analysis the domain.
DLSM is a true meshless method which works with the strong forms of the governing equations and needs not any background mesh either for local response approximations or for building up of system of equations setup and solving.In this study, a new technique has been implemented in this method which enables it to proper diagnosis for the crack and improves its ability for crack growth analysis of the elastic concrete domains.The presented technique is based on the combinations of the concept of Voronoi tessellation algorithm and author's previously developed concrete micro-planes damage theory.With such a technique, in fact, the Voronoi algorithm has been used for detecting the crack and incorporates its effects into the local approximations in the vicinity of the crack in DLSM framework.At the same time micro-planes concrete damage model has been applied for diagnosing the path of the progressive crack during the loading procedure.Two standard benchmark examples have been solved with the current approach to verify its performance in practical applications.Obtained results showed that this method can be improved and extended to capture crack growth phenomenon in concrete structures in the field of meshless methods.

Fig. 2 .
Fig. 2. Formation of the support domain of node j using two layers of neighbor nodes.

Fig. 4 .
Fig. 4. Support domain for node j when m = 2 in a cracked domain.

Fig. 6 .
Fig. 6.Position of the micro-planes on the unit sphere's surface.

Fig. 15 .
Fig. 15.Modeling of the concrete DEN specimen with the proposed method.

Fig. 16 .
Fig. 16.Crack performance obtained with the current method for concrete DEN.