XFEM Simulation of Tensile and Fracture Behavior of Ultraﬁne-Grained Al 6061 Alloy

: In the present work, the tensile and fracture behavior of ultra-ﬁne grained (UFG) Al 6061 alloy was simulated using extended ﬁnite element method (XFEM). UFG Al 6061 alloy processed by cryorolling (CR) and accumulative roll bonding (ARB) was investigated in this work. Numerical simulations of two-dimensional and three-dimensional models were performed in “Abaqus 6.14” software using an elastic-plastic approach, and the results obtained were validated with the experimental results. The specimens corresponding to the three-point bend test, compact tension test with center crack, and double edge cracks were analyzed using XFEM (eXtended Finite Element Method) approach. In XFEM, the partition of unity (PU) was used to model a crack in the standard ﬁnite element mesh. The tensile and fracture properties obtained from the simulation were in tandem with the experimental data. UFG Al alloy showed higher tensile strength and fracture toughness compared to their bulk solution treated counterparts. Fracture toughness was measured in terms of stress intensity factor and J integral. In CR Al alloys, with increasing thickness reduction, an increase in stress intensity factor and a decrease in the J integral was observed. This behavior is attributed to the increase in strength and decrease in ductility of CR samples with increasing thickness reduction. In ARB Al alloys, the strength and ductility have increased with an increase in number of cycles. It also revealed an increase in both the stress intensity factor and J integral in ARB processed Al alloys with increase in number of cycles, as evident from XFEM simulation results. on


Introduction
Aluminum alloys are the second most widely used engineering metal after steel [1]. These alloys are used as structural components in a various industries such as the automotive, aerospace industry [2,3]. Al 6061 alloy is a precipitation-hardenable alloy, which provides a combination of high specific strength, tensile, fracture toughness, good workability, and high corrosion resistance [4]. It is one of the most widely used Al alloys; its application ranges from aircraft fitting and structure such as fuselage and wings in the aerospace industry to brake pistons, hydraulic pistons, chassis (Audi A8) in the automotive industry, and various daily use products such as bike frames, beverage cans, camera lens mounts, etc. [4]. Severe plastic deformation (SPD) consists of a group of metal forming processes such as cryorolling (CR) [5][6][7], accumulative roll bonding (ARB) [8], swagging [9] used to produce ultrafine-grained (UFG) materials from the bulk materials. Grain size is a primary microstructural feature that affects the mechanical behavior of materials. The UFG materials possess better mechanical properties such as tensile strength, fatigue strength, and fracture toughness than conventional coarse grain counterparts [10]. SPD processes impart a large amount of strain to the material and result in grain refinement to nano or submicron levels [11]. Fine-grained microstructure in UFG material leads to higher strength as a consequence of the "Hall-Petch equation".

Processing
Solution treatment (ST) was performed by heating the sample to 530 • C for 2 h. followed by quenching in water at room temperature [28]. During solution treatment, the second phase particles dissolve back into metal matrix, and furthermore, the residual stress present in the material is also relieved. T6 Treatment was performed by heating the solution treated samples to 175 • C for 9 h. This leads to formation of fine spherical Mg 2 Si precipitates [28].
Cryorolling was performed by dipping the sample in to liquid nitrogen for 10 min before each rolling pass. The samples were stored in deep freezer at a temperature of −20 • C to avoid precipitation [10]. At cryogenic temperature, the dynamic recovery of dislocation during rolling is suppressed. This leads to higher dislocation density in the cryorolled samples. High dislocation density along with the ultrafine grain formed in the rolling are responsible for the increase in strength and fracture toughness in cryorolled samples [5].
Accumulative roll bonding (ARB) was performed on sheets of 120 mm × 75 mm × 1 mm. Sheets were stacked on top of each other and then roll bonded with thickness reduction of 50%. These roll bonded sheets were then cut into half, and this cycle is repeated with any lubricant for up to seven cycles [29]. Both CR and ARB do not require any sophisticated dies or machine and can be performed on a conventional rolling mill [5].

Testing
Tensile testing under different strain rates and temperature was performed on an electric universal tension-compression machine [30]. Three point bend test and compact tension test were performed under displacement controlled loading [10].

Simulation Procedure
In this study, the inputs required for simulations were obtained from the literature using "Origin Pro 9.1" software (OriginLab, Northampton, MA, USA), and the simulation was performed using "Abaqus 6.14" FEM software (Dassault system, Johnston, RI, USA).

XFEM Methodology
In conventional FEM, cracks or any other discontinuity are defined as an intrinsic part of the finite element mesh and are modelled by aligning the element boundary with the crack geometry. So, when the crack tip advances, the mesh must be redefined to ensure this alignment. In XFEM approach, the crack is independent of mesh, and hence, a crack can propagate through elements without remeshing. This significantly reduces the computational resources. The XFEM approach is based on three main factors: non-smooth solution properties, the partition of unity, and enrichment function.
Non-smooth solution properties are those properties that experience a rapid change over the domain. These properties are related to discontinuities, such as cracks, grain boundaries, holes, inclusion, etc. In real-world applications, these discontinuities can be frequently observed and are usually categorized into two types: weak discontinuities (the location where field quantities change direction) and strong discontinuities (the location where field quantities have a jump) [10]. The current study is centered on strong discontinuities.
In XFEM, the conventional finite element approximation is improved by using an extrinsic PU approach to introduce discontinuities within the problem. In the past two decades, XFEM has been widely used to model crack growth [19,31]. The cracks in XFEM mesh are defined by enriching the elements with an additional degree of freedom. The shape function in finite elements follows the PU property. This property states that the summation of a specific element's shape function remains unity in all locations within the specified element [19].
Heaviside function is used to determine the discontinuity in the displacement due to crack; this function takes the value of −1 below crack surface and +1 above crack surface.
where x * is closed point projection of x on the crack face, e 2 is normal to crack face at x * , and sign(.) is the signed distance function. Displacement (u h ) in XFEM is defined as: where u FE is the displacement in traditional FEA, u enr is the displacement due to enrichment function.
where U i is the displacement vector corresponding to the conventional finite element, a j is the additional degree of freedom (DOF) due to enrichment by Heaviside function, and b x K is a node enriched DOF corresponding to the asymptotic crack tip. The additional degree of freedom in the enrichment function is calculated only for the elements in the vicinity of cracks or any other discontinuities, as shown in Figure 1. The crack tip is modelled using enrichment function φ(x) in the displacement approximation, which is given as: where (r, θ) is a polar coordinate system with its origin at the crack tip, and (θ = 0) is tangent to the crack at the tip. In the LEFM approach, exponent l is 0.5, while in the EPFM approach, exponent l becomes (1/(1 + n)), where n is the hardening exponent [32,33]. In the absence of enrichment, the Equation (4) becomes the equation for displacement in classical FEA. Figure 1 shows the enrichment of nodes in the vicinity of the crack. To model the crack propagation in the Abaqus using XFEM, both damage initiation and evolution are considered [34]. In the present study, the maximum principal stress criterion is used for damage initiation. The damage will initiate when the principal stress surpasses the maximum principal stress value specified as part of the material properties. For damage evolution, Abaqus uses a scalar damage parameter D to predict the damage in the metal matrix due to nucleation, growth, and coalescence of the microvoids. It increases monotonically with increase in plastic deformation and ranges from 0 at undamaged sample to 1 at complete failure [35]. The damaged stress (σ) using damage coefficient is defined as [36]: where σ eq is the equivalent stress in the undamaged model considering the plastic behavior up to necking. Damage evolution can be defined using either displacement at failure (difference between the displacement at failure and displacement at damage initiation) or fracture energy (area of the curve under the load versus displacement curve). In the current study, displacement at failure is used as a parameter for defining the damage evolution. Damage variable is calculated using following equation [34]: where δ max m is the maximum displacement during loading history, δ f m is the displacement at complete failure and δ 0 m is displacement at damage initiation.

Tensile Simulation
The tensile specimen of Al 6061 was modelled as per the ASTM E8 standard [37], and geometry of the specimen is shown in Figure 2a. Samples for simulation were modelled as a quarter section of dog bone shaped circular specimen, as shown in Figure 2b. Material model in the elastic region was defined by Young's modulus (E) and Poisson's ratio, E was calculated from the slope of the stress-strain curve in the elastic region from the experimental data of each case, and Poisson's ratio was taken as 0.3; in the plastic region, the stress and plastic strain data were used; and for damage criteria, fracture strain, strain rate and displacement at fracture were used. These inputs required were taken from the literature and our earlier work [10,28,29]. A structural mesh containing 2736 C3D8R (8-node linear brick, reduced integration, hourglass control) elements was used. The elastoplastic approach was used to predict the tensile behavior in different conditions. Displacement controlled loading was applied with the boundary condition of keeping one side fixed, while displacement loading is applied on the other end, as shown in Figure 2b. Two additional boundary conditions were used based on the assumption that the cross-section of the specimen remains circular during the tensile test; in xz plane (shaded by yellow color in Figure 2b), displacement and rotation perpendicular to the plane were kept zero i.e., u y = ur x = ur y = 0 and similarly, in the yz plane, u x , ur y , and ur z were kept zero.

Flow Curve
The effect of variation in strain rate and temperature on the tensile behavior of the Al 6061 was predicted using the Johnson Cook (JC) material model. Geometry and the boundary condition used to simulate the flow curve are the same as described in the previous section. In the elastic region, Young's modulus and Poisson's ratio were used, and in the plastic region, the JC material model was used. The constitutive equation of the Johnson Cook model is as follows: where σ is equivalent plastic stress; A is yield stress of reference state, i.e., quasi-static test, B is hardening modulus, C is strain rate dependent coefficient, n is work hardening component, ε is the equivalent plastic strain, . ε and . ε 0 are the plastic strain rate and the reference plastic strain rate corresponding to the quasi-static test, respectively. T * is the normalized temperature, T 0 is the reference temperature (i.e., room temperature), T m is the melting temperature, and m is the thermal softening coefficient. The parameters in the Johnson Cook model was determined from the experimental data at different temperature and strain rate. A is the yield stress at reference strain rate (1 × 10 −3 s −1 ) at room temperature [38]. At reference strain rate and room temperature, the Equation (8) reduces to: Taking natural logarithmic on both side of the equation, the modified equation can be obtained as shown below: The plot between ln(σ − A) and ln ε was drawn, and first-order regression model was fitted to the plot. The slope of the regression model gives constant n, and the y-intercept of the equation is ln B.
At reference temperature, the Equation (8) reduces to: The values of A, B, and n obtained using Equation (10) were put into the Equation (11).
At reference strain rate (0.001 s −1 ), the Equation (8) reduces to: Taking natural logarithmic on both sides: Substituting the material constant A, B, and n into Equation (13) and fitting the firstorder regression model, the slope of the regression curve gives us material constant m. JC failure or damage model is given by [36].
where ∆ε is the increment of the equivalent plastic strain and ε f is the equivalent plastic strain at failure. The equivalent plastic strain is given by: where D 1 − D 5 are material parameters. The first term in the equation depends on pressure; the second term corresponds to the strain rate effect, and the third term is due to thermal effects.

Three-Point Bend Test Simulation
A three-point bend test specimen was modelled as per ASTM E1820−09e1 with a predefined crack length of 1.5 mm [5]. The specifications of the three-point bend test specimen are shown in Figure 3. Cylindrical rigid bodies were defined for applying load and providing support, i.e., displacement type loading in the y-direction is applied to the loading pin, and the supporting pin is kept fixed (i.e., u x = u y = u z = ur x = ur y = ur z = 0) as shown by red colour in the Figure 3b. The three-point bend test specimen was meshed using C3D8R elements, the same elements as used in the tensile simulation. The rigid bodies were defined using a 4-node, 3-D bilinear rigid quadrilateral (R3D4) elements; these elements do not allow deformation. Geometry was loaded using a quasi-static load in displacement-controlled loading. Tensile properties and elasto-plastic approach were used to predict the load versus displacement curve [19].

Compact Tension Simulation
CT specimen was modelled as per ASTM standard E647 with crack length and simulated according to Mode-I loading [29]. The model was meshed using CPS4R, a four-node bilinear plane stress quadrilateral, reduced integration, and hourglass control elements in 2D model, and C3D8R, an 8-node linear brick, reduced integration, hourglass control elements were used for 3D model. Since the computational resources required for the XFEM simulation mainly depend on the mesh size, mesh sensitivity analysis was conducted to optimize the mesh size. Two reference points located at the center of holes were defined in the model, and these reference points were kinetically coupled with the corresponding holes. Displacement controlled loading in y-direction was applied to the upper reference point, while the lower reference point was kept fixed as shown by red color in Figure 4a. The load versus displacement curve is predicted using the XFEM simulations.

Center Crack and Double Edge Cracked Specimen
The rectangular plate of 40 mm × 20 mm dimension was modelled with a crack length of 5 mm. Displacement loading with a strain rate of 10 −6 was applied to one side of the plate, and the other end was kept fixed, as shown in Figure 5. CPS4R, a four-node bilinear, quadrilateral, reduced integration, and hourglass control element, were used to mesh the model. Tensile properties and elastoplastic approach are used as inputs to predict the load versus displacement curves in different conditions, and fracture toughness is determined by load-displacement curves obtained from the simulations.

Tensile Test Simulation
Von Mises stress distribution during the different stages of tensile simulation is shown in Figure 6. When the model is loaded, uniform stress distribution can be seen throughout the gauge length; when necking starts, the stress becomes concentrated in the necking zone. The classical Poisson's effect can be observed in the model. Mesh sensitivity analysis was performed by varying the mesh size from 1.5 mm to 2.4 mm, but no significant change in the tensile behavior was observed as shown in Table 1.  Cryorolled (CR) and Accumulative Roll bonded (ARB) specimen shows higher tensile strength and lower ductility compared to the solution treated (ST) condition. This increase in tensile strength is due to the formation of UFG structure in Al alloy [39][40][41][42]. As the deformation increases, the microstructure becomes more refined, and dislocation density increases [43,44]; so, as a direct consequence of the Hall-Petch equation, the tensile strength increases [10]. This effect can be observed in both CR and ARB samples; for example, the tensile strength of CR25 is 21.5% higher compared to the ST samples, and as thickness reduction is increased to 50%, the tensile strength increases from 210.95 MPa to 242.5 MPa, and as the thickness reduction is further increased to 75%, the tensile strength increases to 285.18 MPa, which is 64.4% higher compared to the ST samples. This increase in the tensile strength of CR is accompanied by the decrease in ductility but in ARB samples, as number of cycles increases, both tensile strength and ductility of samples increases. Tensile strength of ARB-7 cycle sample is similar to the tensile strength of CR50 sample but ductility of ARB-7 cycle (7.8) is low compared to the CR25 sample (11.5)) as shown in Table 2. CR75 samples have 19% higher tensile strength and 9.51% higher ductility compared to the ARB-7 cycles. Tensile strength and ductility of ST, T6, CR and ARB samples predicted from the simulation closely matches the results obtained from the experiments, but a tail can be observed at the end of the stress-strain curve, as shown in Figure 7. This tail is due to the fact that load will decrease drastically when displacement is increased beyond a critical value.

Flow Curve
Tensile behavior of Al 6061 in different temperatures and strain rates is shown in Figure 8 and Table 3. Tensile strength of Al6061-T6 under different strain rates shows a sharp reduction after 423 K similar to the quasi-static case with a strain rate of 1 × 10 −3 s −1 . As temperature increases from 293 K to 423 K, the tensile strength decreases by 11.5%, but when the temperature is further increased to 573 K, there is 59% reduction from 348.1 MPa to 142.7 MPa. At lower temperatures, at 423 K and 573 K, the stress-strain curve predicted by simulation using the Johnson-Cook model is comparable to the experimental results, but at higher temperatures, the simulated curve varies greatly from the experimental curve. The high-temperature sensitivity of Al 6061 alloy, especially at high temperatures, might be the cause of this behavior [30]. JC model is an empirical model that can be used to predict the tensile behavior over a range of temperature and strain rates, but it fails in a hot working process, i.e., a process involving a working temperature higher than the recrystallization temperature [30].

Three-Point Bend Test
Fracture toughness of the material is measured in terms of several parameters such as stress intensity factor (SIF) (K)) and J integral. K is a linear elastic fracture mechanics (LEFM) parameter. When the stress intensity factor attains a critical value, the crack starts propagating. The critical SIF for three-point bend test specimen is calculated according to ASTM standard E399 by [5]: f (s) = 2.9s 0.5 − 4.6s 1.5 + 21.8s 2.5 − 37.6s 3.5 + 38.7s 4.5 (16) where P is the load, S is the span length, B is the specimen thickness, W is the specimen width, a is crack length, and f(s) is called as a compliance function. J integral comes under the domain of elastic-plastic fracture mechanics, and the EPFM approach is directly linked with the crack propagation phase in the ductile material [45]. J integral can be calculated as per ASTM standard 1820-15a using the following equation [46]: where A is the area under load versus crack opening displacement (COD) curve, B is the specimen thickness, b is the length of the unbroken ligament (b = W − a), and B is the width of the specimen. Von Mises stress distributions during different stages of crack propagation are shown in Figure 9. In this simulation, the tensile properties were used to predict fracture toughness. The load versus displacement curve obtained from the simulation varies from the curve obtained from the experiments, as shown in Figure 10. However, the fracture toughness obtained from the simulated data is in good agreement with the experimental values. Mesh sensitivity analysis was performed by refining mesh size from 1.2 mm to 0.5 mm, as shown in Table 4, and it reaffirms the mesh independent crack growth capabilities of the XFEM.   The stress intensity factor of UFG Al 6061 is higher compared to the ST. This increase is contributed to the increase in the tensile strength due to the formation of ultrafine grains and the high dislocation density of UFG samples. As deformation increases, tensile strength increases, and as a result, the critical stress intensity factor increases [10,29], but the increase in tensile strength of CR samples is accompanied by a decrease in the ductility. This decrease in the ductility is reflected in SEM fractography performed by Balakrishnan et al. [10]. In CR samples, as thickness reduction increases, the number of brittle facets also increases. The stress intensity factor is a direct indication of resistance offered by a material to the crack initiation. So, the reduction in ductility has no effect on the variation of SIF. J integral is an indicator of the crack propagation phase and depends on the material's ductility and strength. So, as strength and ductility increase, the J integral also increases. In UFG Al alloys, strength is higher, and ductility is lower, as shown in Table 2. The combination of these two factors determines the variation of J integral; for example, the J integral of CR25 is 127.9% higher compared to the ST sample; this increase is due to the higher strength of CR25 alloy, but even though the tensile strength of CR50 is higher compared to the CR25, the J integral of CR50 is 32.5% lower compared to that of CR50. A similar trend is observed as thickness reduction is further increased to 75%, the J integral decreases from 18.5 KJ/m 2 to 11.38 KJ/m 2 . This decrease in the J integral is attributed to the lower ductility of the CR50 compared to the CR25 condition as shown in Table 5. Plastic zone formed around the crack is a direct indicator of the resistance offered by the material for crack propagation. So, before the crack starts propagating, the size of the plastic zone determines the stress intensity factor. This effect can be observed in the CR samples; the size of the plastic zone formed before crack propagation increases as the percentage thickness reduction increases.

Compact Tension Test
Critical stress intensity factor in CT specimen is given by [10]: where P is the load, B is thickness, W is the width of the specimen, and a is the crack length, f (s) is called as the compliance function. The Von Mises stress distribution during crack initiation and propagation can be seen in Figure 11. The stress distribution during loading shows typical mode I behavior, i.e., the stress distribution is symmetrical with respect to the crack axis. The results obtained from the simulation match closely with the experimental results, as shown in Figure 12. Simulation performed by using the 3D model shows consistently better results compared to the 2D simulation, but the time required for 3D simulation with the same parameters like mesh size, step size, etc., is 2-2.5 times higher compared to 2D simulation. A mesh sensitivity analysis was performed by varying mesh size from 0.5 to 1.2 mm. Variation in mesh size in 0.7-1.7 mm range does not show any significant variation in the result as observed from Table 6.   The value of fracture toughness obtained from the simulation are in good agreement with the experimental results. The trend of fracture toughness in CR samples is the same as discussed in the three-point bend test. The increase in the percentage thickness reduction in CR samples is accompanied by the increase in maximum load and K (critical stress intensity factor) value, but the decrease in ductility causes a lower J integral in CR50 compared to the CR25. In ARB samples, with an increase in the number of cycles, the strength as well as ductility increases, as shown in Figure 7c. So, as the number of cycles increases from one to seven, the stress intensity factor increases from 16.1 MPa m 1/2 to 25.4 MPa m 1/2 , and J integral increases from 11.5 KJ/m 2 to 24 KJ/m 2 , as shown in Table 7. SEM fractography of ARB samples was performed by Rahmatabadi et al. [29]. They observed that as the number of cycles increases, the dimple size decreases (which indicates a reduction in the grain size), and the presence of lamellar structure and brittle facets decreases.

Center Crack and Double Edge Crack under Mode-I Loading
The stress intensity factor for the center crack specimen is given by [47]: where W is the width of the specimen, a is the crack length, and σ is the nominal stress. The stress intensity factor of the double edge specimen is calculated using the following equation [48]: where σ is the stress applied; a is the crack length. Mesh sensitivity analysis of the center cracked, and double edge cracked specimen is shown in Table 6, and the change in stress intensity factor is within 1%. Fracture toughness in the center crack and double edge crack is comparable to each other, as shown in Table 7. Von Mises stress distribution during loading, deformation, and finally rupture is shown in center cracked, and double edge cracked specimen can be seen in Figures 13 and 14. The plastic zone formed around the crack can be observed in Figures 13b,c and 14b,c. The plastic zone formed is symmetrical about the crack axis, which is a characteristic feature of the stress field under mode-I loading. Table 7 shows the fracture toughness in 25%, 50%, and 75% rolled CR samples and ARB-1 cycle and ARB-7 cycle samples. Fracture toughness of the center cracked samples, and double edge cracked specimen in CR, and ARB samples are comparable to each other but compared to the fracture toughness of CT, and three-point bend samples, the fracture toughness of center cracked and double edge cracked specimen is consistently lower.

Conclusions
Numerical simulation of tensile behavior and fracture toughness of CR and ARB Al 6061 was performed in the current study using the XFEM technique. Experimental data on Tensile and fracture toughness were taken from the literature and our earlier work were used to validate the simulation results. The following conclusions can be made based on the present work.
• Prediction of tensile behavior at different strain rates and temperatures was made using the JC material model. The properties at lower temperatures were comparable to the experimental properties, but at a higher temperature, 623K, the JC model failed to predict the tensile behavior. • Fracture toughness was calculated in different testing conditions such as the threepoint bend test, compact tension specimen, center-cracked specimen, and double edge cracked specimen. UFG alloys have higher fracture toughness compared to their bulk counterparts due to their finer grain size and higher dislocation density. In CR samples, as thickness reduction increases, the stress intensity factor increases from 17.45 MPa. M 1/2 in ST samples to 27.53 MPa. M 1/2 in CR25, 31.64 MPa. M 1/2 in CR50 and finally 34.01 MPa. M 1/2 in CR75 samples, but J integral in CR samples decreases as thickness reduction decreases. This decrease is due to a decrease in ductility with an increase in thickness reduction. In ARB samples, the strength as well as ductility increases, which leads to increase in both stress intensity factor and J integral.

•
The fracture toughness obtained from XFEM simulation was in tandem with the experimental results. In the case of CT specimen, the 3D simulation produced slightly better results on fracture toughness compared to the 2D model, but the time required for the 3D simulation is approximately 2-2.5 times higher compared to the 2D simulation. The fracture toughness in center cracked, and double edge cracked specimen is lower compared to the fracture toughness in CT and three-point bend test specimen.

Funding:
The authors acknowledge the funding support from Ministry of Education through Institute of Eminence (IOE).

Institutional Review Board Statement: Not applicable.
Informed Consent Statement: Not applicable.

Data Availability Statement:
The current results are a part of ongoing research work and associated data are confidential at current stage.

Conflicts of Interest:
The authors declare no conflict of interest.

Nomenclature UFG
Ultra fine grained XFEM Extended finite element method SPD Severe plastic deformation ARB Accumulative roll bonding CR Cryorolling CR25 Cryorolling with 25% thickness reduction CR50 Cryorolling with 50% thickness reduction CR75 Cryorolling with 75% thickness reduction JC Johnson Cook material model PU Partition of unity