Code development for the computational analysis of crack propagation in structures

In this study, the main objective was the creation of a code, which gives the capability to a Finite Element Analysis Program with no builtin crack study tools, to study the propagation of a crack, in a cracked surface. For this purpose, the Finite Element Program FEMAP 11.3.2 with solver the NX NASTRAN has been used, and the proposed code was created, using the Application Program Interface (API) of the program. The Linear Elastic Fracture Mechanics (LEFM) theory has been applied to the code, and can predict, if the crack will propagate, the trajectory of the crack, as well as the number of cycle loads required for the propagation of the crack, for given boundary conditions and loads. Finally, the Stress Intensity Factors (SIF) produced by the program, were compared with results from an analytical method. Also, experimental results have been used, for the verification of the results of the trajectory of the propagation, and the cycle loads.


INTRODUCTION
t is well known that numerous thin-wall metal structures, for instance in ships and airplanes, during their lifetime, show cracks in various places. If they are not detected and treated on time, they can be fatal to the structure, with devastating consequences for the crew and passengers, and of course to the environment. Therefore, predicting, locating, and generally studying cracks is of primary importance. For these reasons, methodologies, regulations, and computational tools have been developed, that can predict the creation of a crack, as well as the period that a structural detail can withstand fatigue cycle loads. In this paper, it is assumed that the crack already exists. Therefore, this study will focus on the study of the propagation of the crack. In the numerical analysis of cracked problems, a lot of methods and procedures have already been proposed [1][2][3][4][5]. The origin of the Fracture Mechanics seems to be located in 15th-16th century. According to some authors [6,7], Leonardo Da Vinci studied the fracture strength of iron wires of different lengths, using a device described on the Codice Atlantico [8]. Nowadays, researchers use various finite element analysis (FEA) software to study fracture mechanics problems [9,10]. I The majority of these studies are conducted with FEA programs with built-in tools (for example ANSYS [11,12] and ABAQUS [13,14]), which can calculate the stress intensity factors at the tip of a crack. In this study, it is shown that with some computational effort, a person can study cracked surfaces using a FEA program with no crack study tools. For this reason, it is chosen the FEA program FEMAP 11.3.2 [15], which is a well-known Finite Element Analysis program with special features and capabilities, but with no such built-in crack study tools.
At first, in order to encounter the singularity 1/ r at the crack tip, the code is programmed to calculate the Stress Intensity Factors (KI, KII, and KIII) using the Crack Opening Displacement (COD) method [16]. Using Finite Element Method (FEM) in crack problems involves the following difficulty. If one uses regular standard elements, a very fine mesh is required at the crack tip, influencing the results. In order to deal with this singularity, the use of modified six node triangular elements (quarter-point elements) [17], has been applied (Fig. 1). Figure 1: Modified six node triangular elements [16].
Secondly, after the calculation of the SIFs, the code calculates the kinking angle, using the criteria of maximum tangential stress (MTS) for mixed-mode loading, which was suggested by Erdogan and Sih [18]. Finally, having calculated all of the above, if a cycle load is applied to the structure, the code calculates the number of cycles loads that are necessary for the crack to propagate. To verify the proposed code three applications have been made. The first one concerns the calculation of stress intensity factors in the case of a cracked plate with an edge crack. The second one is the crack propagation, using three specimens with holes, with different position and length of the crack, and the third one the verification of the fatigue lifetime estimation.

THEORY
he fracture problem in the scientific community is a problem of major importance. In this section are presented the basic theory and the crack propagation criteria applied to the code.

Stress Intensity Factors
In the general case of mixed-mode loading condition, the relative displacement of the crack faces is evaluated in relation to each other. Using a FEA program and taking into consideration the crack opening displacement, the SIF KI, KII, and K III (mode I, mode II, and mode III respectively) can be obtained from the following equations [16] (see Fig. 1): In the case of mixed-mode loading condition K I and K II , and in order to evaluate the crack propagation, an equivalent SIF ( eq K ), should be calculated, and compare it with the critical SIF (KIc, fracture toughness). This equivalent SIF can be evaluated from a number of criteria. In this paper, Tanaka's (Eqn. (4)  The KIII SIF has already been considered in various studies (for example Diaz et al. [21]). Since the verification models in our study are in mixed-mode I and II (see section 'Verification of the code'), we have considered only K I and K II SIFs. The KIII SIF may be also included in our code for out of plane sliding mode III (Fig. 12).

Crack propagation criteria.
In the bibliography several criteria can be found which can calculate the angle of the propagation. One criterion is the MTS that was suggested by Erdogan and Sih [18]. According to this criterion, the crack deflection angle φ0is found to be: For pure mode II (K I =0), the deflection angle is φ 0 =-70,5 ο Another criterion for the angle of propagation is the Richard's criterion [20]. According to this: For this criterion, for pure mode II, the deflection angle is 72.1 o . As shown in Fig.3 when KII>0 the angle φ0<0 and vice versa, while always KI>0. The comparison of the Richard's and Erdogan-Sih criteria is presented in Fig.4. It appears that Richard's criterion gives a slightly higher angle. Inside the code, the criterion of Erdogan-Sih is taken into account. It should be mentioned, that the MTS has satisfactory results if the crack growth is controlled be mode I (KI). However, in some recent experiments, large differences between crack propagation direction and the predicted kinking angle were observed, when the crack growth was driven by mode II (KII) [22]. In these cases, the maximum shear stress (MSS) criterion was proposed consequently [23]. In this study, in all of the verification models, the crack growth is controlled by KI.

Fatigue Lifetime Estimation
If an alternating load acts on a structure, a crack may propagate, even if the SIFs are below the static fracture toughness. This phenomenon is called fatigue crack growth. The alternating loads can be distinguished between constant and variable amplitudes (Fig.5). In the two left diagrams σ-t (upper diagram) and K-t (lower diagram) in Fig.5, are presented the constant amplitude alternating load, and in the right diagrams (σ-t and K-t) of the same figure are presented the variable amplitude alternating load. In this study is only considered constant amplitude. The characteristics of the cyclic load are the stress range Δσ, the mean stress m  , and the stress ratio R [16]: The range ΔΚ, which is also known as cyclic stress intensity factor is [16]: is a geometry function. With the help of stress ratio R, the range ΔΚ is written as: When the crack growth rate / d dN  is plotted in a double logarithmic scale as a function of cycle stress intensity factor ΔΚ, the crack growth curve shown in Fig.6 [16] can be obtained. In Fig.6 it is shown that the crack growth curve can be divided into three regions (A1, A2, and A3): A1: In the first region is located the . In case the ΔΚ is lower than this value, the crack will not propagate, no matter what the number of the cycle loads. Above the th K  the propagation is stable.
A2: In the second region the propagation is stable. A3: The last region is the transition to the brittle fracture. If the ΔΚ is lower than the C  then the propagation is stable.
On the other hand, if the ΔK is higher than the IC  , then the propagation is unstable.
In the bibliography, one can find various phenomenological models that can calculate the cycle loads for a given da .
Three of those are:  Paris-Erdogan equation [24]. Is the well-known Paris rule. It describes the second (A2) region (Fig.6).
The coefficient C and the exponent n are characteristics of the material.
The coefficients p and q are empirical constants describing the curvature near the threshold and the instability region of the crack growth curve, respectively [26]. In Fig.7 are shown as an example, three crack growth curves, using the three previously mentioned rules. The material properties used to create the curves in Fig.7 are shown in Tab. 1 [27]. All the above-mentioned methodologies for the calculation of the cycle loads were inserted into the proposed code. For the verification of the calculated cycle loads using experimental results, the Paris Law was used. The fatigue lifetime estimation is very important for the structure, since according to the applied criteria relative to the structure, the materials used, and the experiments, the crack propagation may be predicted.

DESCRIPTION OF THE CODE
he code was created, using the API (Application Program Interface) of FEMAP [28]. The programming language that the FEMAP uses is the WinWrap. Basically, it is a modified Visual Basic [15]. The code is a .BAS file, and it is stored in the C:\FEMAPv112\api (it is not necessary to store this file in this location, but it is customary to do that for well order matters). In Fig.8 is presented a typical Desktop of FEMAP. In Fig.9 is shown the button (CRACK) that was created and it is linked to the .BAS file. When the user presses the button (CRACK) and calls the function, the code asks at first a number of questions in order for the user to declare values-variables. Secondly, the code calculates the SIFs, prints the result, checks if the crack will propagate or not, and asks the user if he/she wishes to see the propagation direction. If the user declares his/her intention to watch the trajectory of the crack, the following procedure starts:  It opens an Excel sheet with a number of empty diagrams. In every step, these diagrams are filling with new data. In Fig.10 are presented the complete diagrams of the crack propagation of the specimen shown in Fig.17, after 99 steps with increment 0.05mm.  Erasing all elements of the structure, except the Crack Tip Elements (CTE).  Rotates the CTE in relation to a local coordinate system, according to the Erdogan-Sih criterion [18]. T  Moves the CTE according to the declared crack increment da.  Re-mesh the structure (surfaces).  Creation of the crack at the re-meshed area.  When the structure is load free, the crack is completely closed (Fig. 11,a). Because of the re-meshing, the elements that are located at both sides of the crack line have common nodes (blue nodes in Fig. 11,a). These nodes are coincident, and FEMAP joins them. Therefore, we have a crack line, but not a crack. In order to create at every step the new crack (the old one, plus a new with an increment da), the code duplicates the nodes located at the crack line. After the duplication, one node goes to the element above the crack line and the other one, in the element below the crack line. The result is shown in Fig. 11,b.  Analysis of the load case.  Calculation of the SIFs of the new crack.  Extraction of the new data in order to fill the Tab. 2 with the results, and the diagrams in Fig.10   End of step. If the SIFs are higher than the ΔKth, or the step number is lower from the maximum number of iterations, then the above procedure will be repeated from the beginning.
In Fig.12, it is presented the flowchart which describes the above-mentioned functionality.

VERIFICATION OF THE CODE
n this chapter, a comparison of the results of the calculation of the SIFs, with the results from an analytical method, is presented. In addition, is expounded the comparison of the trajectory and the fatigue lifetime estimation of the crack, created by the code and FEMAP, with results from an experiment, in order to check the accuracy and the percentage deviation.

Calculation of Stress Intensity Factors.
In order to check the results of the calculation of the SIF from the code, an analytical expression was used [29]. For the calculation of   and I   , the following equations were used: According to Eqn. (13) and Eqn. (14), the plate shown in Fig. 13 with the characteristics from Tab. 3, was used, in order to calculate analytically the SIFs. The results are shown in Tab. 4 with the crack length α equal to 5mm. Taking into consideration the analytical solution, five FEMAP models were created. All models were considered to be in plane strain condition. Every model has a different number of contours, in order to check whether the results are affected by locally changing the mesh at the crack tip. The first model has one contour and the fifth has five. The global mesh size is 0.5mm. In Fig.14      It is observed from Tab. 5, that the smallest percentage error, for K I is 0.90% and it is with one concentric circle, and for KII is 4.92% it is with two contours.

Verification of the trajectory of the crack.
For the verification of the trajectory of the crack, three specimens with holes, with different length and position of the crack were used, from an experimental procedure that took place at Cornell University [30]. These specimens were considered to be in plane strain condition, and the material was a commercial Plexiglas (cyro acrylic FF Plexiglas MC). The equivalent dimensions in the international system of units (SI), as well as the material properties (in SI also), are shown in Fig.16.   In Fig.17 is presented the path of the first specimen, that produced by the proposed code into FEMAP, and the comparison with the experimental path. For the FEA, about 5450 plane elements (6 node triangular, and 8 node quadrilateral) with global mesh size 0.2 mm have been used. It appears that the results are very close. In Fig.18 is shown the second specimen and the comparison of the results. In the second specimen, about 6560 plane elements with global mesh size 0.1mm have been used, and the trajectories are almost identical.  In the Fig.19 is displayed the last (third) specimen using about 7680 plane elements with global mesh size 0.1mm. The comparison of the trajectory of this specimen, as the second one is almost identical to the experimentally obtained trajectory. Studying all of the above comparisons, the conclusion is that the code produces reliable trajectories, due to the fact that the trajectories of all specimens (from FEMAP [28]and experiments [30]) are quite similar, and in most cases identical.

Verification of the fatigue lifetime estimation.
In order to verify the results of the lifetime estimation produced by the proposed code in FEMAP, experimental result was used [31]. In Fig.20 it is shown the specimen from the experiment that was used. The modified compact specimen was originally proposed by Miranda et al. [32]. The idea was to put a through-thickness hole above the possible crack path, in order to produce mixed-mode behavior of the crack. The specimens used in the experiments of Lu et al. [31], were similar to Miranda et al. [32], but not the same since different materials and dimensions have been used. It should be mentioned, that the initial crack length was 10.5 mm [31]. The specimen used, is considered to be in plane strain condition, and the material was Al2024-T3. In addition to the geometrical characteristics, described in Fig.20, the data shown in Tab. 6 have also been used. Furthermore, the algorithm for calculating the number of loading cycles was based on the Paris rule [24]. 1.41*10  n 3 Table 6: Additional data that were used in the experiment [31].
In Fig.21 the cracked specimen and the results of the experiment and the proposed code in FEMAP are shown.  The comparison of the trajectories is presented in Fig.21. The one has been created by the experiment (Fig. 21,a), and the other has been created by the proposed code in FEMAP (Fig. 21,b) using about 4190 elements with global mesh size 0.8mm. In addition, the comparison of the cycle loads from the experiment and the proposed code in FEMAP is shown in Fig. 22. From Fig.22, up to (about) 56000 cycles and for crack length about 17mm, the results were identical. After this point, there was a difference in cycle loads. In the experimental procedure, the propagation stops at 58000 cycles, and in FEMAP for the same crack length, the propagation stops at 61900 cycles, which means a difference of 5.99%.

CONCLUSION
code was created in order for the FEA program FEMAP 11.3.2, to be able to study fracture mechanics problems. In order to verify the results produced by the proposed code in FEMAP, analytical calculations and experimental results were used. From the analysis it is observed: At first, the comparison of the calculated SIFs from FEMAP and the analytical method shows a difference of 0.90% for KI and 4.92% for KII.
Secondly, the trajectories from all cases from FEMAP were checked with experimental results. They show almost identical behavior.
Finally, the estimation of the fatigue lifetime was also checked with experimental results. The results were identical up to a certain point. After that, a difference of 5.99% occurred. The differences of the calculations for the SIFs are small, and possibly can be reduced, replacing the COD method with the J integral method, which introduced by James Rice [33]. It is more complicated to apply this method to the FEA, but many studies in the bibliography show that using this method better results can be produced. For the crack propagation criteria already used, Diaz [34] found that SED, although is computationally more expensive, it describes the processes of the calculation of the Keq and the kinking angle, very well for ductile materials. Moreover, the problem with the selected Keq criteria is that negative K II has the same effect as a positive one. Rozumek and Macha [35], have published a comprehensive review on the calculation of the Keq. Finally, the path of the crack and the estimation of the lifetime which produced by the proposed code in FEMAP, are almost identical with the experimental results. The deviation at the end of the calculation of the lifetime is considered small. A cause of this deviation can be the difference of the calculations for the KII, as was mentioned previously in (i) or even in an internal defect of the material of the specimen which used in the experiment. According to Diaz [34], the difference between experimental and numerical results in the modified compact specimen may also be attributed to crack closure and roughness. Consequently, the proposed code gives reliable results and can be used in order to study cracked surfaces. Furthermore, apart from the simple cases of cracked surfaces, this code can be used as a base for more complicated and complex studies. Two examples of these studies can be the study of Crack Growth Behavior of Welded Stiffened Panel, and the effects due to overloads (for example [36,37]). In addition, out-of-plane sliding mode III crack problems may be also studied with the proposed code. The procedure proposed by this code is the capability of the finite element program FEMAP to confront successfully crack mechanics problems. In addition, the proposed code has the possibility with a computational effort to study cracked problems with no built-in tools. Displacement at the x-direction of the local coordinate system placed at the tip of the crack u2

NOMENCLATURE
Displacement at the y-direction of the local coordinate system placed at the tip of the crack u 3 Displacement at the z-direction of the local coordinate system placed at the tip of the crack Δσ stress amplitude ΔΚ cyclic stress intensity factor range ΔΚth cyclic stress intensity factor range threshold θ the angle between the initial direction and the direction of the new crack growth increment σm mean stress level σ max maximum stress level σmin minimum stress level v Poisson's ratio ASTM American society of testing materials COD crack opening displacement CTE crack tip elements FEA finite element analysis LEFM linear elastic fracture mechanics MSS maximum shear stress MTS maximum tangential stress SED Strain energy density criterion SIF stress intensity factor