An Improved Damage Mechanics Method for Solving Three Dimensional Crack Problems of Concrete

In this paper, a damage mechanical model is improved to simulate the three-dimensional (3D) mixed mode cracking process of quasi-brittle material such as concrete. The cracking and crack propagating in an element is reflected in the model by controlling the degraded stiffness matrix and the dissipated energy. The local tracking algorithm (LTA) is employed to preselect the element that will belong to the crack path to be damaged. The local tracking algorithm is extended to 3D by analyzing and overcoming defects in the previous 3D LTAs. The present LTA has very low computational cost and this fact is proved by comparison with well-known global tracking algorithm. The numerical results from the novel damage model added 3D local tracking algorithm show good agreement with the results of the previous researches, demonstrating the superior performance of the approach in predicting crack propagation and fracture strength.


Introduction
Numerical analysis of crack fracture of concrete has attracted much attention in civil engineering field [1][2][3][4]. Continuum damage mechanics (CDM) models based on smeared band theory have been widely used to simulate the material fracture behavior [5][6][7][8]. Recently, Cevera [9] has found that the CDM combined with the crack tracking algorithm (TA) can predict the crack path well in the simulation of the mixed mode fracture process of quasi-brittle materials. This method is more practical because it can successfully solve the discrete cracking problem and is simpler and faster than the embedded finite element method (E-FEM) or the extended finite element method (X-FEM).
Slobbe [10] suggested a delayed local tracking algorithm (LTA) having the postponement of crack path fix moment based on the distinction between crack paths and damage paths and combined it with the 2D isotropic damage model. Saloustros [11] extended the use of a combination of isotropic damage model and LTA to the simulation of internal cracks such as the shear cracks for structural fracture problem analysis. Yun [12] developed a damage model with crack direction parameter, which added a LTA with self-correction of the crack tip element and propagation direction to solve U-turns of the crack path. Burnett [13] proposed a mesh objective method that uses the smeared crack approach to prevent the crack direction from being dependent on the direction or size of the mesh.
Most of these studies employ the isotropic damage model and deal with the two-dimensional crack problems. This is related to the convergence difficulties of the Newton-Raphson method for solving the anisotropic model and the lack of a 3D LTA. The LTAs, which are extensively used with CDM in 2D, have the discrepancy between the crack planes of the element when they are directly extended from 2D to 3D.
In this paper, a 3D anisotropic damage model is developed to solve the discrete crack problems in concrete. The developed model will reflect the damage evolution and crack propagation in the continuous element due to tension and shear stresses based on the energy equivalence. The issues for expanding the local tracking algorithm to three dimension are also addressed in this paper. ANSYS Parameter Design Language (APDL) is used for the present damage model and tracking algorithm. The study shows that the present 3D damage model coupled with 3D LTA resolves the problem of discrete crack path prediction.

Mechanics Model
In this section, the governing equations based on 3D anisotropic damage model are detailed described. The present 3D damage model has higher crack tracking ability than the isotropic damage model.
The constitutive equation and effective stress is expressed as the following form.
, d  σ C ε σ Cε (1) where   In general, expressions of the second order damage tensor and 4th order damaged stiffness tensors are used in damage models of isotropic material, but they are lack of physical intuitive. In the present study, the expression of the second order damaged stiffness tensor, which is widely employed for composites damage, is used for isotropic material. Converting the damaged stiffness matrix from the local coordinate system (CS) to global Cartesian CS is performed using the transformation matrix of stress σ T .
T dd   σ σ C T C T (2) The damaged stiffness tensor in the local CS is expressed as follows, allowing only the reduction of material properties in the normal vector N of the crack surface. The local CS should be given in the considered element to get the transformation matrix in Eq. (2). The following three local coordinate axes are set in each damaged element.
where 0 T is the crack direction vector at the previous crack tip element adjacent to the present crack tip element. The crack propagation direction is also defined in the 3D elements, as in the 2D elements, by introducing the crack direction vector T .
The transformation matrix of stress σ T in Eq. (2) is expressed by the components of three vectors in global CS.

T S S N T N T T S T S S N S N N N T T S S N T N T T S T S S N S N N N T T S S N T N T T S T S S N S N
The exponential damage evolution law is used.
where r is the elastic domain threshold and S H is the softening parameter. The softening parameter is defined as: where the length parameters e l and mat l are given as: where c G is the energy release rate in fracture mechanics and V G is the total energy dissipated per unit volume during the damage evolution. It is noted that the characteristic length of element e l helps to reduce the mesh dependence of the simulation results.
The elastic domain threshold r in Eq. (7) is expressed as: where   is the value of loading function  at the given time step  .
The loading function is expressed as follows: The damage activation function is defined as: Since the direction of the principal stress calculated by the local stress tensors causes the oscillations of crack surfaces, the non-local stress is applied to determine the crack direction. It is noted that the non-local tensor is used only to obtain the crack direction and local CS, and not to other equations such as damage evolution equation. Commands of ANSYS for force and displacement controls are used to prevent the occurrence of convergence difficulties in the solution system.

3D Local Tracking
In most discrete crack methods, crack tracking algorithms are needed to accurately predict the crack path [10]. It is well known that X-FEM or E-FEM without TA also leads to the mesh-induced crack path biases or discontinuous path in numerical fracture analysis [2,14,15]. TAs added to damage mechanics models prevent many elements from being damaged at the same time due to large stress at the crack tip, and required that elements belonging to the crack path are damaged. The combination of TA and the damage model can also eliminate the effect of the mesh on the simulated crack path. In this section, we present an efficient 3D LTA to solve topological difficulties and reduce computational costs.
The topological difficulty of 3D LTAs is that the crack planes generated individually in each element of the crack path do not have coincident edges [16]. Jä ger [2] checked the performance of Belytschko's LTA [17] and Holzapfel's two-step TA [18], and reported that their stability and robustness are poor due to topological difficulties. Although Gasser [19] proposed LTA averaging the crack planes of each element of the crack path, continuity of the averaged crack planes of all damaged elements was not fully guaranteed.
In the present LTA, the trace points propagate along the crack propagation direction T in Eq. Track triangles with three trace points (the present trace point, the previous trace point of the same flow line, and the tip trace point of the adjacent flow line) are drawn to form the crack surface. The creation of one trace point produces two trace triangles. And then the crack surface is generated as the sum of these trace triangles, which contain all trace points and lines. The crack surface is divided into trace points and triangles, so this method is referred as crack surface discretization. The elements that have intersections with the crack surface are added to the crack path and the local CSs are set for damage evolution of these elements.
The crack surface is generated independently with the solid mesh by eliminating the constraint that one trace triangle should be placed inside an element. Trace points should only be placed on element faces.
The correction of the crack tip element and the maximum crack inclined angle criterion are applied to avoid sharp turns often encountered in the discrete crack method [10,[20][21][22] and improve the robustness of the crack path. The correction methods are described in detail in Ref. [12] for 2D problems and are almost same as 2D for 3D. Once the damage variable D of the element is greater than 0.999, the killing element technique is applied to this element, taking into account the fact that the

Numerical Examples
Several numerical experiments are performed and the results are compared with literature to validate the proposed approach in this section.
The criterion of maximum principal stress is used to predict the crack surface and the standard finite elements (4-noded linear tetrahedron elements) with continuous displacement are used for all numerical simulations.

3D Envelope Tracking
The problem that tracking the 3D envelopes of a normal vector field   Fig. 1 for different meshes. As demonstrated by Oliver [23], the scalar field method can accurately obtain these envelopes. However, it should solve the global heat-conduction like problem for tracking the envelope at each loading step and so requires a lot of computational expense. x y z    respectively. As shown in Fig. 2, our LTA tracks envelopes well in the space given the normal vector field.
The computation times consumed by two different tracking methods (our LTA and scalar field method) are compared for different element sizes. A standard PC with a 3.2 GHz quad-core Intel processor and 16 GB of Ram is used to perform tracking computations. When using the element size of 0.1, the scalar field method and our LTA require computation time of 8s and 4s respectively. And computation time of 73s and 32s is consumed by the scalar field method and our LTA respectively when using the element size of 0.05.
The present LTA tracks the total envelopes of crack surface once, but the scalar field method should get the envelopes each iteration step. If the number of iterations is 100, the entire time consumed by the scalar field method is 800s and 7300s for two meshes respectively. From this comparison, it has been confirmed that our LTA can save much more computational cost than the scalar field method.

Three-Point Bending Beam with Offset Notch
In this example, a beam with offset notch in the bottom plane is subjected to three-point bending [12]. The beam with an initial crack at the bottom of the beam offset 0.8mm from the center of the bottom plane, is subjected to a load applied by vertical displacement at the center of the top plane.  The geometry of beam specimen under three-point bending is presented in Fig. 3 Since mixed mode I -II crack propagation is studied in this problem, the curved crack path should be obtained in simulation. The capability of the present approach to track the curved crack path is tested by solving this example. At the same time, this problem is numerically analyzed using the isotropic damage model and the comparison with our method is performed in this section. Numerical results by isotropic damage model are shown in Fig. 4. It can be seen that the crack path obtained by the isotropic damage model represents a straight crack path with little curvature, which indicates that the isotropic model has no ability to predict the curve path.   Fig. 5, it is proved that the present model has a good crack path prediction ability. It indicates that the present approach has far better capability than the isotropic damage model in crack propagation prediction. The load-displacement curve by the present approach is compared to ones by isotropic model and XFEM in Fig. 6. The load-displacement curves by our approach are very similar with the one by X-FEM in prediction of magnitude of the maximum load, although the unloading parts have a slight difference. Because of inaccurate crack path caused by the isotropic model, the load-displacement curve differs from the other two methods in predicting the maximum load.

L-Shaped Concrete Panel
The third example considers a L-shaped concrete panel specimen experimentally and numerically investigated by many researchers to study the discrete crack propagation and fracture behavior of concrete [2,21,24]. The sample geometry and the load conditions are shown in Fig. 7 The crack path predicted in this problem should initially be curved upward and gradually straightened as depicted in earlier researches. 12 initial trace points are used in this simulation. The deformed shape at an applied vertical displacement of 0.8 u mm  is shown in Fig. 8 and the predicted crack surface with trace triangles is shown in Fig. 9. As expected, the crack curves upward from the lower edge of the specimen, gradually straightens and approaches the right edge as the load increases. The crack surface obtained is in good agreement with that tested by the experiment [24]. Although Jä ger [2] mentioned that LTA cannot solve the problem of L-shaped panel with curved crack surface, simulation results demonstrate that the present anisotropic model combined with 3D LTA perfectly eliminates the drawbacks of traditional LTAs.   10 shows a comparison of load-displacement curves by the present approach and experiment. The magnitude of the maximum load and tendency of the curves are consistent with the experiment, even if there are differences in the unloading state. It is noted that these differences of the numerically simulated load-displacement curves in the unloading state are also reported in other references [21]. In this example, it is fully proved the accuracy and crack tracking ability of our approach.

Conclusion
One methodology for simulating the propagation of discrete cracks and fracture processes in concrete using the finite element method has been proposed in this paper. This method can accurately simulate the crack propagation process using the damage mechanics method, while using conventional elements without nodal enrichment or elemental enrichment in the finite element method. For tracking the 3D crack path, the local 3D tracking algorithm with very low computational cost is improved by solving topological difficulties. APDL is used to simulate some benchmark problems and validate our damage model and tracking algorithm. Crack tracking ability of our approach has been fully verified in terms of calculation cost and accuracy.