The local arc-length method is employed to control the incremental loading procedure for phase-field brittle fracture modeling. An improved staggered algorithm with energy and damage iterative tolerance convergence criteria is developed based on the residuals of displacement and phase-field. The improved staggered solution scheme is implemented in the commercial software ABAQUS with user-defined element subroutines. The layered system of finite elements is utilized to solve the coupled elastic displacement and phase-field fracture problem. A one-element benchmark test compared with the analytical solution was conducted to validate the feasibility and accuracy of the developed method. Our study shows that the result calculated with the developed method does not depend on the selected size of loading increments. The results of several numerical experiments show that the improved staggered algorithm is efficient for solving the more complex brittle fracture problems.
An improved staggered algorithm that uses the local arc-length method is formulated and proposed.
The energy and damage tolerance convergence criterion are employed in the numerical implementation.
The proposed algorithm results are not dependent on the selected sizes of the loading increments.
Brittle fracture is one of the most feared failure modes in solids and structures. The numerical simulation method of brittle fracture is significant in engineering applications. Numerical simulation models for crack growth problems are generally divided into discrete crack models and diffuse crack models. The discrete crack model is a geometric description method based on mesh extension, which extends on mesh and mesh boundaries, including the element deletion method [
The research field of phase-field fracture models can be divided into the field of physics and the field of mechanics. The physical basis of the phase-field fracture models are based on the Ginzburg-Laudau theory, as shown by the model proposed by Hakim et al. [
In terms of algorithm implementation, the phase-field crack model can be solved with either a monolithic algorithm [
Nevertheless, the Newton-Rapson incremental loading procedure has difficulty capturing the complete failure process and the entire load-displacement curve because of the rapid failure of the structure due to the large loading step. To obtain accurate results, however, a small loading increment is necessary for the staggered algorithm. Several loading algorithms can track the whole load-displacement curve well, such as the global arc-length method [
In this work, the phase-field method based on the rate-independent variational principle of diffuse fracture is considered. The formulations of the staggered algorithm with the local arc-length method are derived and proposed to avoid the sensitivity to the local failure of the structure and the softening of materials. The energy and damage iterative tolerance convergence criterion based on the residual norm of the displacement and phase-field is proposed to save computing resources.
The method of separating the node and degree of freedom of the element is employed to implement the improved staggered algorithm in the ABAQUS package, and they were alternating between the phase-field degree of freedom and displacement degree of freedom. An energy and damage iterative tolerance convergence criterion is embedded in each incremental loading of the local arc-length method. A one element benchmark test compared with the analytical solution was conducted to validate the feasibility and accuracy of the developed method. Several numerical experiments are conducted to simulate crack initiation, intersection, coalescence, branching, and propagation.
The paper is organized as follows. The representation of the phase-field approximation of diffuse crack topology is outlined in
First, a simple model is considered to introduce the variational principle of diffuse crack topology. A cracked 1D bar with cross section
The one-dimensional description of the diffuse crack topology can be directly extended to the multidimensional case. According to the derivation process, the surface density function
Assume that we are given a sharp crack surface topology in
Francfort et al. [
Assuming that linear elasticity exists in the unbroken solid, the elastic energy density can be calculated as:
According to the above fracture variation criterion, the surface energy of cracks results from elastic strain energy, that is, elastic strain energy drives the evolution of the phase-field variables. If the elastic strain energy is not decomposed, it will cause false branching of the crack. Given this problem, this work, based on the method proposed by Miehe et al. [
Considering the crack surface energy expression in
We can substitute part of the energy into the variational principle expression, to obtain:
By defining the change in external and internal potential energy and introducing the principle of virtual displacement, the governing equations of the model are obtained:
In order to obtain the numerical solution of the partial differential
In 2D space, the displacement field and the phase-field can be discretized as:
The strain-displacement matrices are expressed as:
According to the above expression, the residual of the discrete equation corresponding to the equilibrium condition for the displacement field and the crack phase-field can be expressed as:
To obtain the solutions of
In the case of unstable crack growth, the monolithic model is prone to convergence problems. The stress fields are rearranged due to the newly degenerated stiffness matrix when the crack propagates. Therefore, the solver cannot obtain a stable equilibrium solution. To obtain a stable model, the displacement field u and phase-field
Similarly, the residual corresponding to the displacement field is then as:
The following system of the equation can be solved iteratively by employing the Newton-Raphson method,
The new formula generates the above decoupling equations for sequential phase and displacement fields. Then, a staggered scheme of phase-field model is proposed. The displacement, phase and history fields ut,
In this section, the improved staggered algorithm with the local arc-length method formulations and the corresponding numerical implementation are established. The local arc-length method [
The arc-length method adjusts the external force load
The
After linearization by the Newton-Raphson method, the residual
The unknown increment can be obtained by the following formula:
According to the Sherman–Morrison formula,
Using the Newton-Raphson backward Euler scheme means that the numerical solution is unconditionally stable; that is, to achieve convergence means to obtain an equilibrium solution. Therefore, the solution is independent of the time increment. In the classic staggered scheme, the phase and the displacement fields are solved independently. This is unlike the monolithic case, where convergence is fulfilled for both coupled fields simultaneously, retaining unconditional stability. By solving for the phase and the displacement fields independently, the staggered scheme can track the equilibrium solution in a semi-implicit method, leading to a more robust numerical implementation. However, load increment sensitivity studies should be conducted. The residual staggered algorithm is improved to reduce the influence of the loading increment step on the solution process. The flowchart of the improved staggered solution scheme with the local arc-length method is shown in
Expression
Then, following this process and repeating the iteration by assembling the residual stiffness matrix
By introducing the energy and damage convergence criterion, the improved staggered solution scheme is implemented to further improve the computational efficiency in solving phase-field brittle fracture problems. In the computation of the staggered solution, the solutions of the displacement field and phase-field both have their convergence criteria.
At the beginning of the calculation, the damage value is fixed to solve the displacement subproblem. The damage value can be solved by the local arc-length method. When the displacement field converges, the displacement field is taken as a known quantity, and the damage subproblem is solved. There are two boundary conditions for solving the damage subproblem. One is that the phase-field is bounded, which is
The process for the improved staggered algorithm based on the double convergence criterion is shown in
The core of the traditional staggered solution scheme is to fix the phase-field Add the phase-field element and displacement element. The corresponding degrees of freedom are given in the staggered solution process. Add the phase-field element. Two submodels with the same mesh division should be established, one of which is the displacement element and the other being the phase-field element. In the staggered solution scheme, submodels are called for solutions. Separate the node and the degree of freedom of the element. Set the element type without the degree of freedom in the model file, and then give the element’s corresponding node the degree of freedom in the solution process.
Regarding program expansibility and simplicity, the third method is best. The phase-field brittle fracture model can be solved directly in the finite element framework. The commercial software ABAQUS is employed to solve the iteration process and realize the improved staggered algorithm for phase-field modeling.
The layered system of the finite element implemented in ABAQUS with UEL and UMAT is shown in
Since the phase-field equation and the balanced equation are solved separately, the implicit solver ABAQUS/Standard or the explicit solver ABAQUS/Explicit can be flexibly selected according to the model when solving the displacement field. The former is implemented through the UMAT subroutine interface provided by ABAQUS, and the latter is implemented through VUMAT.
As depicted in
The solution-dependent variables used to solve differential equations in two dimensions are listed in
SDVs numbering in subroutine | Variable | |
---|---|---|
SDV1, SDV2 | Displacement |
The first layer displacement element |
SDV3, SDV4 | Axial strains |
|
SDV5 | Shear strain |
|
SDV6, SDV7 | Axial stress |
|
SDV8 | Shear stress |
|
SDV9, SDV10 | Elastic axial stress |
|
SDV11 | Elastic shear stress |
|
SDV12 | strain energy |
|
SDV13 | Elastic strain energy |
|
SDV14 | Phase-field |
|
SDV15 | Phase-field |
The second layer phase-field element |
SDV16 | History field |
A finite element model with the UEL implementation should be developed to verify the feasibility of solving brittle fracture problems using the algorithm, as demonstrated through the most straightforward case in
This case is solved using 5 different displacements Δu to reach u = 0.1 mm. The five-step displacements are Δu = 1 × 10–5, Δu = 5 × 10–5, Δu = 1 × 10–4, Δu = 5 × 10–4, and Δu = 1 × 10–3, corresponding to 10000 steps, 2000 steps, 1000 steps, 200 steps, and 100 steps.
According to
The single edge notched test is the most common test used to validate phase-field brittle fracture models. The geometry and the boundary conditions of the square specimen are depicted in
In addition to the quasi-static crack propagation of the 2D model, the 3D model of the single edge notch test is solved using the same parameter. The geometry and boundary conditions for the 3D single edge notched tensile modeling are depicted in
This case study is a tensile sample of two materials. The upper material is harder, while the lower material is 10 times softer and has an initial notch. The sample is tensioned parallel to the separation line. The geometry and boundary conditions are shown in
The purpose of this benchmark test is to show that our implementation can model the crack branching phenomenon.
Since the purpose of this benchmark test is to show crack branching, the center of the finite element mesh is dense. As the local amplification region is showm in
After the simple example, we will show a few more complex cases to demonstrate the utility of the new implementation. A symmetric and asymmetric double-notched specimen is used to investigate coalescence of cracks.
The crack patterns for symmetric and asymmetric doubled notched specimen tension tests at different displacements are shown in
The L-shaped specimen in [
In addition to the quasi-static crack propagation examples, dynamic crack propagation is also handled by means of the improved staggered algorithm. The geometry and boundary conditions of the dynamic crack branching are shown in
The crack patterns for the dynamic crack branching are shown in
In this study, an improved staggered scheme with the local arc-length method is proposed to study phase-field brittle fracture. The energy and damage iterative tolerance convergence criteria based on the residuals of displacement and phase-field is developed and implemented in ABAQUS. We adopt a one element benchmark numerical test and several numerical examples to verify the feasibility and accuracy of the phase-field model. From the presented work, the following conclusions can be drawn:
The local arc-length method applied to the staggered solution scheme of the phase-field brittle fracture can track the crack propagation process and load-displacement curve well. The improved staggered scheme for brittle fracture is not sensitive to the size of the load increment. The accuracy and efficiency are greatly improved. The results of one element benchmark test and several complex numerical examples are in good agreement with the other numerical solutions. The improved staggered scheme for the brittle fracture model is an effective method to accurately simulate crack initiation, coalescence, intersection, and branching.