Numerical modeling based on moving mesh method to simulate fast crack propagation

An analysis to show the capability of moving mesh strategy to predict dynamic crack growth phenomena in 2D continuum media is proposed. The numerical method is implemented in the framework of the finite element method, which is coupled with moving mesh strategy to simulate the geometry variation produced by the crack tip motion. In particular, a computational procedure based on the combination of Fracture Mechanics concepts and Arbitrary Lagrangian-Eulerian approach (ALE) is developed. This represents a generalization of previous authors’ works in a dynamic framework to propose a unified approach for predicting crack propagation in both static and dynamic frameworks. The crack speed is explicitly evaluated at each time step by using a proper crack tip speed criterion, which can be expressed as a function of energy release rate or stress intensity factor. Experimental and numerical results are proposed to validate the proposed approach. Mesh dependence problem, computational efficiency and numerical complexity are verified by comparative results.


INTRODUCTION
any researchers devoted their attention to investigate dynamic fracture mechanics by means of analytical or/and computational formulations. In the literature, most of the models are developed in the framework of the Finite Element Method (FEM), because of its robustness and reliability to deal with complex geometries. Smeared crack models are utilized with large success, since they are easy to be implemented and reproduce correctly the amount of energy dissipation during the crack growth [1,2]. However, regularization techniques such as nonlocal or gradient based models are required to prevent material instability problems and to capture the internal characteristic length and displacement discontinuity in presence of microcracks. In addition, mesh dependency problems are typically observed, since a large number of FEs should be introduced in the region in which the crack path is expected. Numerical models based on the Cohesive Zone Modelling (CZM) are frequently adopted for simulating crack evolution in both static and dynamic frameworks [3,4]. An important advantage of the CZMs is their ability to predict directly crack onset and propagation, without introducing preexisting material discontinuities [5,6,7]. However, the initial finite stiffness in the constitutive laws may produce, especially in brittle solids, an excess of compliance with spurious interface traction oscillations. Moreover, a refined mesh is required to predict accurately the fracture variables, which in many cases may lead to instability or nonconvergent phenomena. Such problems may be circumvented by the use of re-meshing techniques, but numerical complexities still remain due to the need of local mesh modification algorithm [8]. Despite high computational costs, remeshing based methods have shown rigorous results to predict crack path and related fracture variables. Alternatively, Embedded or Extended Finite Element Methods (EFEM, XFEM) are able to reproduce crack growth without the use of re-meshing procedures by the enrichment of displacement fields inside an element. In particular, material discontinuities propagate within the elements, leaving unaltered the mesh discretization [9]. However, due to lack of proper criteria for crack branching and interaction, reliable predictions of complex fracture patterns remain a key challenge [10]. Recently, several refined formulations are proposed to improve crack onset criterion and growth procedures with respect to conventional models based on CZM or Fracture Mechanics (FM). For instance, methodologies, based on the Scaled Boundary Finite Element Method (SBFEM) for simulating dynamic crack propagation by adopting polygon elements, are developed in [11]. Although, these numerical schemes present good capabilities for computing the fracture variables by means of re-meshing events, they change the global mesh, leading to large computation costs. Differently, efficient numerical models consistent to the Cracking Elements (CE) are proposed, in which strong discontinuity embedded approach, based on the use of disconnected cracking segments, is developed [12]. The main advantage of this approach is the absence of remeshing or enrichments events, whereas the main disadvantages are the inability to describe accurately crack path and SIF/ERR functions. Phase-Field Methods (PFMs) have rapidly spread in view of their capability to seamlessly deal with complex crack patterns like initiation, branching, merging and even fragmentation [13]. One of the main disadvantage of the PFM is the fact that the method is still computationally intensive. In order to avoid some of the issue previously mentioned, numerical models based on moving mesh methodology are developed in the literature. In [14], a model based on a combined approach developed in the framework of Fracture Mechanics and moving mesh methodology is able to predict dynamic delamination phenomena in layered structures. Similarly, numerical strategies based on a coupled approach between moving mesh and the cohesive zone modelling are presented in [15][16][17][18]. In this case, moving mesh methodology based on ALE approach is introduced only to represent process zone region, leaving the governing equations of the structural model, basically, unaltered [19]. The aim of the present work is to generalize the numerical approach developed in [19,20] to describe dynamic crack propagation in 2D structures. In particular, the approach combines concepts arising from structural mechanics and moving mesh methodology, which are implemented in a unified framework to predict crack growth on the basis of Fracture Mechanics variables. Moving computational nodes are modified starting from a fixed referential coordinate system on the basis of a crack growth criterion to predict directionality and displacement of the tip front. This is achieved by introducing both crack propagation in terms of tip speed and angle of propagation, appropriately. Numerical results demonstrating the effectiveness of the method to simulate crack growth in continuum media under dynamic loading conditions or impact phenomena are proposed.

FORMULATION OF THE MODEL
he proposed modeling is based on the combination of FM and Moving Mesh Method (MMM). The former predicts crack growth, by the use of the ERR/SIF concepts and an advancing criterion, whereas the latter is introduced to account geometry changes produced by the crack evolution.

M T
The crack growth is predicted by the use of a moving mesh methodology based on an arbitrary Lagrangian-Eulerian (ALE) formulation. In particular, two coordinate systems are introduced, known as Referential (R) and Moving (M) ones (Fig. 1).
A one to one relationship between the R and the M is defined by the following mapping operator   : where R X  and M X  identify the positions on the computational points in R and M configurations, respectively. More details on the derivation of the governing equations are reported in [21]. According to ALE formulation, the time and spatial derivatives of a generic physical field, in referential and material configurations can be related by the following relationships: where X   represents the relative velocity of the grid points in the material reference system. Analogously, starting from Eqn.(2), the second time derivative is evaluated recursively as: where   is the gradient operator function. The governing equations in the material configuration can be written by means of the principle of d'Alembert, taking into account virtual works of inertial, external and internal forces: where, n  is the unit normal vector,  is the mass density,   is the Cauchy stress tensor, t  is the traction forces vector on the free surface, f  is the volume forces vector dV and dA are the volume and the loaded area in the material configuration.
(2)-(3), consistently to ALE formulation, the governing equations given by Eqns. (4)-(5) should be reformulated to take into account the transformation rule between the Lagrangian and referential coordinate system ( Fig.4): where r V and r  are the volume and the loaded area in the referential configuration, C  is the elastic moduli matrix collecting stiffness coefficients, det( J ) is the determinant of a scalar metric representing the ratio of differential areas. It is worth noting that Eqn. (5) requires the definition of a proper advancing scheme to enforce the crack tip displacements and a rezoning procedure to move the current positions of the crack tip front, keeping the computational mesh undistorted during the whole calculation. In the proposed modeling, the crack tip motion is consistent to the crack growth criterion based on the instantaneous crack tip speed [22]. A fracture function depending from fracture tip variables, such as for instance ERR or SIF is required. According to ALE approach, the region enclosed into the  contour is able to describe the crack motion by introducing the following boundary conditions (Fig. 2): where 0  is the crack propagation angle, whereas F   is the incremental scalar quantity computed at the current iteration step, by adopting a proper dynamic crack growth criterion. Angle prediction and crack tip displacements are evaluated by solving the following equations: where F   is the fracture multiplier, G represents the fracture variable, i.e. total ERR or SIF, 0  is the angle of propagation depending from the crack kinking function F g , G f is the fracture function. Eqn. (6)-(7) are completed by additional boundary conditions applied to the geometry contour lines and time initial conditions (Fig.2): The elements of the remaining regions of the structures are stretched due to a Laplace based regularization method leading to a consistent transition of the mesh during the crack growth: The variational form of the governing equations of the ALE problem is derived starting from Eqn.
where X  is the vector containing the horizontal and vertical nodal point positions,   is the LMM vector. Governing equations, given by Eqn. (5), (7) and Eqn. (10), introduce a nonlinear set of equations, which are solved numerically, using a customized finite element program based on COMSOL Multiphysics [23]. In particular, proper script files are carried out to manage the steps involved in the procedure, regarding the geometry variation due to the crack propagation and the mesh enrichment in the process zone. Moreover, the resulting algebraic equations are solved by using an implicit time integration scheme based on a variable step-size backward differentiation formula (BDF), in which a coupled approach, in which no splitting operators in the solving procedure for plane stress and ALE formulation are considered. More details on the implementation procedure are reported in [20].

EVALUATION OF THE ERR AND CRACK GROWTH CRITERION
n the present section, the main formulas regarding the computation of the ERR and crack growth criterion are summarized. It is worth noting that the present formulation is quite general to be utilized in conjunction with other existing formulas and procedures typically developed in Fracture Mechanics to predict crack growth. In order to compute the ERR components, a path independent J integral formulation developed in [24] has been considered into the proposed numerical scheme by means of the following expression: where   is a contour enclosing the crack tip, W and K are the strain and the kinetic energy densities and k n is the outward normal, , and i i i u u u   are the displacement velocity and acceleration of the material point, respectively,  is the material density,  is an arbitrary contour going around the crack tip. From the numerical point of view, case it is convenient to consider the following expression taking the limit of 0 V   , i.e 0    : The ERR components, evaluated with reference to the global coordinates system, can be projected on the local tip coordinates by using the following coordinate transformation rule (Fig.3):  The computation of the SIFs is developed by using the component separation method, where I K and II K expressions are related to the dynamic J integral, crack velocity parameters and evolution functions [24]. In order to compute the crack propagation direction a proper crack kinking criterion has to be considered. The proposed model incorporates Maximum Energy release rate criterion which is simply defined as follows: in which 0  is measured with respect the horizontal global axis (Fig. 3).

RESULTS
n this section, results are developed with the purpose to verify the consistency and the reliability proposed methodology, by means of comparisons with numerical and experimental data. Moreover, a sensitivity analysis is carried out to verify the consistency of the proposed model in term mesh dependence and accuracy of crack tip motion. In particular, the computational performance of the proposed model is verified by means the investigation of two loading configurations based on a rectangular sample made of Araldite-B. At first, the simulations are conducted on the loading configurations experimentally tested in [25], which involve a pure mode I loading condition. The loading scheme, boundary conditions and geometric configuration are represented in Fig. 4, whereas the mechanical parameters are summarized in Table. 1. Fig. 4b reports the mesh discretization adopted in the numerical simulations, in which a relatively refined mesh is considered only around the crack tip region, whereas the remaining part of the structure presents a transition mesh with a lower element size. The mesh discretization is featured by a ratio between the characteristic length element and the moving region radius equal to in the tip region and transition mesh in the remaining part of the plate with maximum length equal to (M2), involving in a total number of DOFs equal to 7516. Numerical and experimental data are collected by means of a two-step analysis. The first one is needed to load the structure by means a static analysis, in which a vertical displacement at the pin is applied until the mode I SIF reaches its critical value [25]. Once the SIF reaches 15/1 D R I the prescribed value, a restart procedure is achieved and the crack propagation is enforced by the crack growth criterion. In particular, the curve experimentally determined in [25] is implemented to compute the crack tip velocity at each time increment of the simulation (Fig. 5). It is worth noting that, in the present loading scheme, the crack proceeds along the horizontal direction, leading the crack angle always to be equal to zero.   The dynamic behavior of the structure is analyzed by means three different loading levels, which are featured by the initial values of SIF reached during the static loading step:   [11,26] and experimental [25] data; (b) Variation of dynamic stress intensity factor vs time: comparisons with experimental data [25].
In Figs. 6a-b and 7a, results in terms of crack velocity normalized on the Rayleigh wave speed r c a function of the crack tip position normalized on the specimen length (L) obtained by the proposed model are compared with existing experimental [25] and numerical [11,26] data. The curves show high values of the crack tip speed especially during the initiation phase. Once the crack tip moves, an oscillatory behavior is observed until the crack arrest phenomenon is achieved. In Fig. 7b, the dynamic SIF time histories, computed by means of the proposed model and for the three loading configurations, are compared with experimental data [25]. The numerical results are in good agreement with the data arising from the experimental [25] and numerical [11,26] data taken from the literature. All the computations are performed on a Xeon processor running on a Windows 10 system. The governing equations are solved by using an implicit time integration scheme based on a variable step-size backward differentiation formula (BDF). The total computational time used by the proposed numerical scheme is about 700 s, while it is about 360000 s for the remeshing technique developed by [26]. Thus, this formulation allows to save much of the total computational time. In order to verify the computational efficiency of the proposed model, the influence of the mesh discretization as a function of both the mesh dependency and computational costs are discussed. To this end the following mesh discretization are considered:  Fig. 8 shows the mesh discretization around the crack tip region for M1 M2 and M3 configurations. In this case the numerical simulations are employed by considering the loading condition with   In Fig. 10a, comparisons in terms of crack tip speed as a function of the crack tip position in terms of mesh configuration, are reported. The results show how the proposed model is able to describe dynamic crack propagation phenomena by using a coarse mesh discretization. The results, reported in Fig. 10b, present a good agreement in the prediction of the dynamic SIF time history. Moreover, the solution is quite independent from the mesh utilized in the numeric model, since also using a coarse discretization (i.e. M1) no instability phenomena are observed in the prediction of the fracture variables. Fig. 11 shows a synoptic representation of the mesh motion during the crack propagation for the M1 configuration. The mesh movements of the crack tip region are enforced rigidly, ensuring accuracy in the prediction of the fracture variables.
In the remaining parts in which a lower accuracy is required the mesh is affected by distortion phenomena produced by the mesh motion. However, once the distortion tolerance is reached, a re-meshing procedure with an update of the ALE and structural variables is achieved.  Additionally, a numerical case study similar to the previous one but subjected to mixed mode crack propagation is investigated. The geometry and the loading configurations are reported in Fig. 12a, whereas material characteristics are summarized in Tab . The motion of the crack tip region is driven by the crack growth criterion, while the crack is assumed to propagate along the direction computed by using the crack propagation angle defined on the basis of the Maximum Energy release rate criterion [19]. In Fig. 12b, the initial mesh configuration used for the simulation is reported. As shown, a refined mesh discretization is considered just around the crack tip region, whereas the remaining part of the structure requires a relatively coarse mesh. Fig.13a represents the dynamic stress intensity factor time history detected by the proposed formulation. Fig. 13b shows the crack tip speed time history components along X (blue) and Y (red), whereas crack tip speed (black) is computed as 2 2 X Y c c c   .
In Fig. 14, Von Mises stress maps, coincident at several times selected at magenta points of Fig. 13 (from A to F), are reported. The analysis shows that the crack initiation angle is equal to 33° and it remains substantially unaltered during the whole propagation process. It is worth noting that the proposed procedure is strictly dependent from internal parameters, such as the limit value of angle tolerance, which may influence the solution in terms of both computation efforts and accuracy. A detailed discussion about this aspect is reported in [19]. It is worth noting that such parameter controls the maximum variation angle allowed during the crack growth and in this case is set equal to 2°. Once that the predicted crack propagation angle is larger than the tolerance value, a re-meshing procedure is required. This procedure, extensively described in a previous authors' work developed in static framework [27], is made by proper script files implemented in a MATLAB® environment, which manages the numerical procedure.

CONCLUSIONS
ynamic crack propagation phenomena are successfully simulated by using a new numerical methodology, which combines concepts arising from solid mechanics, fracture mechanics and moving mesh methodology. The numerical formulation is able to describe crack propagation just by adopting a proper crack function and kinking D criterion. The numerical model is able to describe fast crack curving phenomena. Two loading configurations involving in both pure mode I and mixed mode are successfully simulated. The validity of the proposed model was proved by means of a sensitivity analysis in terms of mesh dependence and time required for the computation. The proposed model shows capabilities to simulate fast crack phenomena although further studies are required to simulate crack branching phenomena.