Abstract
This work deals with a 2D finite element simulation of nonplanar multiple cracks using fracture and crack propagation analysis. This analysis was performed by using the developed source code software written by Visual Fortran Language. This source code includes the adaptive mesh generation utilizing the advanced front method and also the mesh refinement process. In order to correctly represent the field singularity, the quarter-point singular elements are constructed around the tip of the crack. The crack growth criteria are used to predict the crack growth direction by utilizing the circumferential stress factor in calculating the yielding stress in elastic fracture assumptions. The stress intensity factor determination is one of the most critical procedures as it determines the crack initiation and propagation mechanism. Moreover, the stress intensity factor histories during the crack growth are measured with the use of equivalent domain integral methods. The crack path simulation and stress intensity factor calculations are compared with the literature and revealed that the results are in agreement with research carried in this domain.
1. Introduction
Characterisation of the mechanical behavior of the solid materials of the crack-free surface involves important influencing factors such as stress, strain, force, and loading system to which the material is subjected. The crack growth mechanism of these surfaces is analysed by the two categories of fracture mechanics, namely, linear elastic fracture mechanics and elastic-plastic fracture mechanics [1]. Stress intensity factor plays a major role in predicting and analysing the crack initiation and growth in the domain of fracture mechanics analysis. Stress intensity factors (SIFs) directly relate to the amplitude of the crack initiated in terms of singular tip. The tip of the crack has maximum probability to possess higher stress intensity factor. The stress intensity factor calculation involves consideration of variables such as stress, strain, and displacement. Several books involving the calculation methods of stress intensity factor calculation are available [2, 3] for specific geometries and loading. Due to the limitation of the analytical solution obtained in stress intensity factor determination, the vast majority of fracture problems encountered in engineering is solved with the help of numerical methods [4]. Stress intensity factors can be computed by various methods, such as the boundary element method (BEM)-finite element method (FEM), both FEM and BEM [5]. The critical problem encountered by the researcher in fracture mechanics is opting the appropriate method for accurate prediction stress intensities.
There exist two categories of commercial FEM programs: finite element programs such as ANSYS [6, 7], ABAQUS, and NASTRAN can be used to add elements manually and perform analysis on complex structures. The other type is professional FEM programs such as NASGROW and AF-GROSS, and they are comparatively very accurate in providing a high-precision calculation [8]. However, there are limitations when it comes to more complicated geometries and loading conditions. Complicated geometries require a high-density mesh and a complex element in the simulation method. Other researchers developed their own two-dimensional source code program to determine and analyse the fatigue crack growth and crack propagation under the condition, namely, static loading and also the determination of stress intensity factor using the mesh strategy [9–12]. It is difficult to opt a method for accurate calculation with minimal time and at low cost. Several techniques involving various equations in SIF calculation are available in the literature [2, 13, 14]. Recently, Fu et al. [15] proposed a crack-tip element for modelling crack propagation by taking the advantages of the symplectic analytical singular element (SASE) and the floating node method (FNM). In their method, accurate crack-tip fields (displacement and stress) can be captured, and multiple crack propagations can be modelled without remeshing. This is essentially because of the use of the crack-tip asymptotic analytical solution with higher order expanding terms. Another benefit of their method is that the stress intensity factors (SIFs) can be solved without any postprocessing. An investigation of the generalized dynamic stress intensity factors of cracked homogeneous and linear magnetoelectroelastic solids using the extended finite element method was presented by Bui and Zhang [16]. Furthermore, Bui and Zhang [17] presented a transient dynamic analysis of stationary cracks in two-dimensional, homogeneous, and linear piezoelectric solids subjected to coupled electromechanical impact loads using the extended finite element method (X-FEM). They developed a dynamic X-FEM computer code using quadrilateral elements in conjunction with the level set method to accurately describe the crack geometry. The beauty of the X-FEM lies in the fact that the discontinuity and singularity induced by the crack are effectively treated as the mesh is completely independent of the crack geometry, and more interestingly, remeshing in crack propagation is no longer required [18].
The failures due to cracks developed in mechanical bodies are due to many factors like geometrical variations as well as the differing loading conditions. There are many formulas and methods to calculate the stress intensity factor depending on the various loading conditions and body geometries. In fact, most of the fracture mechanics analysis consider the single-loading condition instead of multiple loads in analysis since the later is complex to analyse in simulated analysis [19]. Due to these reasons, the determination of SIF nonstandard geometries is required by using numerical methods like finite element methods. This research was motivated by a practical need to model the crack propagation under mixed-mode loading and fatigue life prediction in LEFM. The developed program solves two-dimensional LEFM problems using the adaptive mesh FEM. This model is provided a finite element code that produces results comparable to the currently available commercial software. As far as the acquaintanceship is concerned, using commercial software is not appropriate at least in two standpoints: first, the fundamental algorithm that lies behind it is not fully comprehended, and secondly, state of the art in the programming skill is absolutely unapprehended. Commercial software may also be employed to simulate crack propagation and fatigue life prediction, but such software is very expansive and almost cannot get the source code to make some development on it. The computational efficiency is highly dependent on many parameters such as the mesh density and number of mesh refinement, as well as the number of cracks in the geometry. The computational efficiency will significantly increase by decreasing or eliminating the remeshing process by using the extended finite element method. Huynh et al. [20] introduced a novel and effective computational approach that is based on polygonal X-FEM (named as PolyXFEM) for the analysis of 2D linear elastic fracture mechanics problems. The PolyXFEM is equipped with a new numerical integration technique that uses the concept of Cartesian transformation method over polygonal domains. This method is computationally more efficient compared to the two-level mapping integration commonly used on polygons. This study presents a two-dimensional finite element simulation of nonplanar multiple cracks using fracture and crack propagation analysis with the aid of a developed source code program written by Visual Fortran Language.
2. Developed Program Criteria and Mesh Refinements
This work relates to the advancing front method [21] which is one among the easiest mesh generation processes. The algorithm used to generate the element starts with “front” form of analysed boundary of the considered domain to create the element and then advances to the discrete regions of defined boundaries to complete the entire domain. A researcher introduced procedures in the advancing front method of the element generation algorithm [22] involving a triangulation construct using a set of predetermined points inside the domain. Another research indicates that the point generation in triangular element construction is instantaneous as part of algorithm generation [23].
The geometrical characteristics of the mesh generated can be defined using the background mesh, and graded meshes are needed to define nonuniform distribution of element sizes. Introduction of stretches in specified directions is employed to analyse directional orientation of elements. The important paths of mesh generation are given [24] as follows:(1)Node generation to define the boundary edge of the domain(2)Element and node generations within the discretised boundary(3)Element shape analysis and improvement to ensure quality of the mesh
The element shape, size, and orientation are the geometrical characteristics of the mesh and said to be the spatial functions. The front mesh generation is employed prior to triangular element creation. The creation of initial front mesh involves data of all sides of the discretised boundary edges, relating to closed loops of boundary conditions. If the predetermined region consists of multiple connected subregions of differing properties, then the front mesh generation is done for every region of differing material properties [9].
Every side of the front mesh generated is defined by two end points. The sides of analysis are arranged such that the domains to be meshed fall to the left of the generation front. In the completion stage of mesh generation, mesh smoothening is carried out to ensure the improved shape of the mesh. The interior nodes of the domains are repositioned without changing the nodal connections of the element to improvise on the shape. The Laplacian smoothening is a famous and efficient method that can be used for computational smoothening of the algorithm of the analysis [25] which indicates the repositioning of the centroid of the polygon formed by the nodes. This technique provides the effective adjustment of the shape of the element in analysis.
Adaptive subdomain remeshing techniques add each crack advance segment by modifying the current finite element mesh. This is done by replacing elements and nodes in a region ahead of the current crack tip with smaller elements along the crack faces to be added. These additional elements typically include a small structured rosette of elements at the new crack-tip location. The remaining space between is filled with additional elements. Furthermore, geometric updating of a crack path during each incremental step is carried out because cracks are not represented explicitly in an adaptive mesh. This method insists on the need for the mesh refinement among the tips of the cracks considered [5]. The maximum circumferential stress, the energy release rate, and the strain-energy criteria are used in crack path prediction [26].
In this proposed work, the circumferential stress criteria are considered for the isotropic materials subjected to multiple loading conditions [9, 11]. The following is the maximum tangential stress at which the crack propagates in the normal direction:
The aforementioned normal direction can be obtained by solving for . The following is the nontrivial solution:where the crack growth angle is given by
3. Stress Intensity Factor Method
In 1968, Rice [27] introduced the J-integral method to study nonlinear material behavior in small-scale yielding. It is a path-independent contour integral defined aswhere strain-energy density is denoted 0 by W; stresses by ; local I axis displace ui; the contour has an arc length and is expressed as s; and nj is the outward unit normal to the contour C around the crack tip (Figure 1(a)).
(a)
(b)
The integral method which utilizes the equivalent domain is more appropriate to the finite element simulation when compared to that of the finite size domain method relating to the divergence theorem by integration. For 2D problems, area integral replaces the contour integral i (Figure 1(b)). Then, equation (4) is rewritten as
In LEFM simulation, the J-integral by definition take account the translational mechanical energy balance in the tip region of the crack along x-axis. In either cases of mode I or mode II, equation (5) simplifies the SIF calculation KI or KII but fails for multiple loading conditions to calculate KI and KII separately. For such a case, the integral that can be used is as follows [28]:where k represents the index of the crack tip. The integral method is used for smaller deformations, and then the same is extended for the greater one. The following are the methods by which the SIF can be obtained. The first case uses the J-integral and SIFs. These relations are
The second methodology uses the relation between SIFs:
The natural triangle quarter-point elements [29] are used to obtain the linear elastic singularity for stresses and strain field in the vicinity of the crack tip. The schematic rosette formation of the quarter-point elements around a crack tip is shown in Figure 2.
4. Adaptive Mesh Refinement
The adaptive mesh refinement is employed as the optimization scheme. This scheme is based on a posteriori error estimator which is obtained from the solution from the previous mesh. The stress error norm is taken as the error estimator. The main idea in the h-type adaptive mesh refinement is to obtain the ratio of element norm stress error to the average norm stress error of the whole domain which is also known as relative stress norm error, and from this, ratio of the new size of the element for the refinement process can be predicted. In this procedure, the mesh size of each element is defined aswhere is the area of the triangle element. The norm stress error for each element is defined bywhile the average norm stress error for the whole domain iswhere m is the number of total elements in the whole domain and is the smoothed stress vector which consists of smoothed stresses components. In the finite element treatment, the integration with the isoparametric triangular element will be converted by the summation of quadratics following the Radau rules as follows:and similarly,where is the thickness of the element for a plane stress condition and for a plane strain condition. It is obvious that these parameters are evaluated involving the values of stresses and smoothed stresses at the sampling points only.
It is considerable then to make the relative stress norm error for each element less than some specified value (say 5% for many engineering applications [24]. Thus,and the new element relative stress error norm with permissible error of is defined as
This means any element with needs to be refined, and subsequently, the new size of mesh refinement needs to be predicted. Here, the asymptotic convergence rate criteria are used whereby it is assumed thatwhere is the polynomial order of approximation. In the present study case, since the quadratic polynomial is used for the finite element approximation. The approximate size of the new element is
The current mesh will be regarded as the new background mesh, and the advancing front method will be repeated depending on the number of mesh refinement sets by the user.
5. Results and Discussion
5.1. Plate with One Central Hole and Cracks of Interaction
In this work, nonplanar multiple crack initiation and growth under the condition of uniaxial stress and consisting a hole and three cracks are presented. The boundary conditions and load details are revealed in Figure 3. The body is subjected to cycling loading with controlled displacement of −0.005 mm. The displacement was observed at the base area of the body so that the tension prevails at the top portion. The body is made of Al-7075 for which the elastic modulus and Poisson’s ratio are 72 GPa and 0.33, respectively. The thickness of the body is fixed as 3 mm to have simulation of the plane stress condition. In Figure 4, the finite element model of the specimen is depicted in their initial states.
As shown in Figure 5(a), the source causing the shear stress is the eccentricity of the cracks so that mode II of SIF is dominating. Hence, there is change of direction which is observed in crack propagation from horizontal to vertical. The vertical distance between the second and third cracks is observed to be small leading to higher SIFs compared to that of surrounding regions. Therefore, the critical observation is that the second and third cracks grow faster than that of the first crack.
(a)
(b)
Figure 5(a) also shows the path of the crack growth predicted in the present study, and the mathematical results obtained [30] by FEM were compared with results obtained by [31], as shown in Figure 5(b). As shown in the figure, the results are matching with the literature [31].
The final crack path including the mesh as well as the final mesh with deformation is shown in Figure 6.
(a)
(b)
Figure 6 clearly shows that the crack is growing continuously, and there exists a change in propagation direction due to the presence of stress redistributions. This is clearly visible while analysing Figures 7 and 8. Here, curved crack length is the summation of all previous incremental crack growth values at a given step during propagation. As shown in Figures 7 and 8, the drops of the values of KI are the reasons for changing the crack path direction from a straight line to the curvature due to the influence of the hole on the crack direction. But, the effects of the load on crack 1 are smaller compared with the other two cracks, and there is no influence of the hole on the crack direction so that it propagates in a direction away from the existence of the hole only up to 15 mm, and the predicted SIFs are shown in Figure 9.
For multiple existence of cracks, the crack propagates one after the other according to SIF values. According to Figures 7–9, crack 3 will propagate first followed by crack 2 and finally crack 1. Figures 10(a) and 10(b) show the maximum and minimum principal stress distribution for the last step of the crack growth, respectively, as well as the magnification of the cracks paths for both stresses is depicted. As seen in the legend for both stresses, there is a slightly difference in the values (28 MPa for the maximum principal stress compared to 20 MPa for the minimum principal stress). For more visualization for the crack paths, Figure 11 shows the enlargement of the area around the cracks in the last step of crack propagation.
(a)
(b)
5.2. Nonplanar Multiple Crack Propagation in a Thin Plate
The second example deals with two eccentric through-thickness cracks that eventually propagate in a nonplanar mode. This problem was studied using numerical analysis by Price and Trevelyan in 2014 [32]. In this problem, two edge cracks in a thin plate with initial lengths of 10 mm growth under cyclic tensile loading are seen. Geometry details, load, and the conditions applied for the boundaries are shown in Figure 12. As shown in the figure, zero displacements are applied at the bottom of the plate in x, y, and z directions. The top surface of the body is subjected to 100 MPa stress. The modulus of elasticity and Poisson’s ratio utilized in the calculation are 30 GPa and 0.3, respectively.
Figure 13 shows the crack growth path predicted by using the developed program compared with the numerical results obtained by Price and Trevelyan in 2014 [32] by using the boundary element method as well as the FEM results achieved by [30]. As shown in this figure, the agreement with both results is clear. The final crack growth path as well as the mesh in the deformed form are shown in Figure 14.
(a)
(b)
(c)
(a)
(b)
The predicted SIF values for the two cracks are shown in Figures 15 and 16, respectively. As shown in both figures, the SIF of mode I factor has dominated, and the crack growth direction has a slight effect on mode II which caused the curvature path on the end. According to the values of stress intensity factors, the two cracks will propagate at the same time.
For more declaration for the stress distribution on this geometry, the maximum and minimum principal stress distribution is shown in Figure 17. The units of the stresses were in Pa.
(a)
(b)
6. Conclusions
The simulation of nonplanar multiple cracks in a two-dimensional finite element using fracture and crack propagation analysis has been performed. The strategy has been used successfully to simulate the initiation and propagation of cracks in a plate specimen with a hole and others without the hole. The existence of the hole in the plate affects the crack to change its direction towards the hole based on the hole size and position. However, for multiple cracks, the cracks propagate one after the other according to the values of the stress intensity factors, as well as the direction of the crack growth will also be affected by the existence of other cracks to attract each other depending on its position one from another.
Data Availability
The data used to support the findings of this study are included within the article.
Conflicts of Interest
The authors declare that they have no conflicts of interest.