Ductile failure analysis of epoxy resin plates containing multiple circular arc cracks by means of the equivalent material concept

Ductile failure of polymeric samples weakened by circular arc cracks is studied theoretically and experimentally in this research. Various arrangements of cracks with different arc angles are considered in the specimens such that crack tips experienced the mixed mode I/II loading conditions. Fracture tests are conducted on the multi-cracked specimens and their fracture loads are achieved. To provide the results, the equivalent material concept (EMC) is used in conjunction of dislocation method and a brittle fracture criterion such that there is no necessity for performing complex and time-consuming elastic-plastic damage analyses. Theoretical and experimental stress intensity factors are computed and compared with each other by employing the fracture curves which demonstrate the appropriate efficiency of proposed method to predict the tests results.


Introduction
The damage analysis of defected components has been accomplished recently by several researchers. The major cause of structure failure is the presence of cracks in the material. In fact, the strength of components decreases in the presence of cracks. This destructive effect is mainly due the existence of stress concentration at the crack tips during the different kinds of loadings which may cause the failure initiation in the material. Cracks usually exist in components during their manufacturing process or service life; therefore, the presence of undesirable cracks is unavoidable in the members (Aliha et al. 2016;Gustafsson et al. 2017). Different shapes of cracks can be observed in several components and structures. The natural gas transportation systems made of polyethylene (PE), and the pinch clamping of these system, combustion chambers in which arc cracks-shape geometries, are embedded in the design to facilitate the fluid flow; canned food and canisters and the grenades are some examples of components in which cracks can be seen. Consequently, it is important to find a way for predicting the structure failure which begins by crack propagation.
Epoxy resin is one of the most widespread polymeric materials used in various engineering and industrial applications. This great interest for epoxy resins is due to its considerable mechanical and thermal properties such as high adhesion to many substrates, good heat and chemical resistance, and low shrinkage (Jin et al. 2015). Polymeric specimens are usually used to work under mechanical and thermal loadings. For predicting the load-carrying capacity of these components, their mechanical properties should be determined at first. One of the main properties used for evaluating the failure of polymeric members is their fracture toughness which should be specified accurately. Therefore, an appropriate estimation of the material properties should be provided to have a reliable operation for the members under different kinds of loadings.
Although most of the researches in the field of fracture mechanics have corresponded to the damage analysis of single cracked specimens, there are various samples in engineering and industrial members weakened by complicated arrangement of cracks. The main difficulty of multi-cracked problems lies in the interaction between the cracks. Therefore, simple methods of analysis are not applicable in such cases (Pourseifi and Faal 2015;Yan and Miao 2012). For special situations of aligned multiple cracks, exact solutions were found using classical methods (Erdogan 1962;Sih 1965); but in more complicated problems, approximate methods should be employed. Polynomial approximations and truncation of a conformal mapping function have been applied to problems with zig-zag crack configurations (Kitagawa and Yuuki 1975;Kitagawa and Yuuki 1978). In several researches, dislocation method has been employed to study the multi-cracked problems in which approximate ways have been used to solve the obtained integral equations (Fotuhi and Fariborz 2008;Fotuhi and Fariborz 2013;Hills et al. 1996;Pourseifi et al. 2017;Weertman and Weertman 1992).
For the first time, Kachanov (1985 andKachanov, 1987) proposed a simple procedure to analyze the crack interactions in multiple cracks problem in which superposition method was used to obtain a system of linear algebraic expression. Benveniste et al. (1989) also successfully used the superposition technique to find the stress field and the stress intensity factors (SIFs) for the situation of complicated multiple crack arrangement. In another research, the Bueckner's principle (Buckner 1958) was extended by Yan (2010) to analyze the interactions between different cracks in an infinite linear elastic plate. The SIFs were achieved by applying the displacement discontinuity approach proposed recently by Yan (2003). There are also other researches in the literature dealing with the multiple crack problems in linear elastic materials.
To the best of authors' knowledge, in the field of the failure analysis of multi-cracked components, all of the available papers in the literature are related to the fracture behavior of brittle materials. The mechanical behavior of polymers, specially their fracture property, significantly depends on the brittleness or ductility of them ). There are several parameters like hardener percentage, environment, curing procedure, and test conditions that affect the type of crack growth and consequently, the fracture behavior of polymeric material (Kanchanomaia et al. 2005;Osswald and Menges 2003;Xianyan et al. 2017). Noticeable plastic deformations exist in the material during the ductile failure; therefore, the classical linear elastic fracture mechanic (LEFM) becomes invalid for studying such problems. Although several methods have been developed by researchers to predict the ductile fracture, such as the Jintegral (Dawes 1976), the resistance curve (R-curve) (Anderson 1995), the crack tip opening displacement (CTOD) (Roberson and Tetelman 1973), and the crack tip opening angle (CTOA) (Anderson 1995), the complexity and timeconsuming procedures of these methods make them undesirable for designers. Recently, a novel approach has been proposed by Torabi (2012), named as the equivalent material concept (EMC), which has been used as a simple method to predict the ductile failure in several researches until now Cicero et al. 2017;Majidi et al. 2018;Torabi and Alaei 2015;Torabi et al. 2016). This concept was also successfully employed in conjunction of maximum tangential stress and mean stress criteria to estimate the ductile failure of polymeric specimens containing different kinds of notches Torabi et al. 2017;Torabi et al. 2018).
As a novel study, the ductile fracture of epoxy specimens weakened by multiple circular arc cracks is investigated experimentally and theoretically in this research. For this aim, epoxy material with considerable ductile behavior is utilized to fabricate square plates weakened by circular arc cracks with various arrangements. Different angles are considered for the circular arc cracks locating in collinear and parallel arrangements to study the effect of interactions between cracks under mechanical loading. Remote tension is applied to the specimens to obtain their load-carrying capacities (LCCs). Stress intensity factors (SIFs) of specimens are determined by combination of EMC and dislocation method and also by performing finite element analyses. The EMC is combined with dislocation method and generalized maximum tangential stress criterion to evaluate the theoretical and experimental results. It is demonstrated that the proposed approach can successfully predict the experimental results.

Experiments
By using the results of the previous characterization tests , the multi-cracked specimens are fabricated from the epoxy resin and the fracture tests are conducted on them. Details are presented in the following subsections.

Material
Araldite LY 5052 with desirable ductility is selected as a commercially available epoxy resin to prepare the test specimens. The type and concentration of the hardener have significant effects on the mechanical behavior of the ultimate polymer (Gensler et al. 2000). Hence, according to the previous papers published recently Torabi et al. 2017;Torabi et al. 2018), Aradur 5052 is selected as the curing agent for this epoxy resin and a mixture with 50 wt% hardener concentration is prepared. All the fracture tests are conducted under displacement control loading with the rate of 1 mm/min. Square plates are loaded uniaxially under remote tension and their load carrying capacities, the maximum loads that specimens can undergo, are recorded. The vacuum mixing procedure is used to prepare a homogenous mixture. The epoxy resin and hardener are poured in a vacuum mixing machine container, and the mixture is uniformly blended under vacuum pressure of 0.7 bars at ambient temperature. After preparing the suitable mixture, curing process is followed based on the resin supplier recommendation (2007). Therefore, the molded sample is kept in the ambient conditions for 24 h and then it is placed in the furnace for 4 h at the temperature of 100°C.
In the next step, the mechanical properties of the prepared polymer should be specified for using in the fracture analyses. Hence, these properties were determined by performing the tensile and the Poisson's ratio tests based on the standard test methods ASTM D638 (1994) and ASTM E132-04 (2010), respectively. The tensile properties are taken from previous research of the second author, Rahimi and his co-authors . For each standard test, 5 repetitions were considered. The obtained values are reported in Table 1.
The digital image correlation (DIC) technique was employed as a high accuracy method to achieve the tensile stress-strain curve for the epoxy material. This obtained stress-strain curve is presented in Fig. 1. Digital images of the specimen surface are used during the standard tensile test to obtain the displacement and strain of the specimen. For this aim, first a white background is painted on the surface of the specimen then black spots are sprayed randomly on this surface. By choosing two arbitrary spots and using the digital images taken every 5 s, the stress-strain curve is obtained.

Fracture experiments
To obtain the final geometries of multi-cracked samples from the prepared polymeric sheets of 3 mm thick, a high-precision 2D CNC water-jet cutting machine is employed. After that, a sharp razor blade is used to create the crack at the end of the arc tips. The thickness of the specimens is specified based on the previous papers of the second authors' works Torabi et al. 2017;Torabi et al. 2018) in which the thickness of 4 mm is determined to have plane stress condition in the notched samples made from Araldite LY 5052. Therefore, specimens with the thickness of 3 mm also have the plane stress condition. Based on this condition, the critical stress intensity factor coefficient (K c ) is used in the formulation instead of K Ic . This condition is also explained in the Torabi et al. paper (Torabi et al. 2019). Different arrangements are considered for circular arc cracks in the epoxy specimens to study the interactions between the cracks. As shown in Fig. 2, each specimen contains two circular arc cracks which are located horizontally or vertically in the specimens. The samples of prepared specimens and the high magnification image of crack tip are depicted in Fig. 3. The geometrical parameters of the tested multi-cracked specimens are depicted in Fig. 2. By considering various values for the arc angle β of circular cracks, different mode mixity ratios are obtained in tested specimens. The values of β contain 20, 40, 60, and 80 (deg.), and the radius of all Table 1 The mechanical properties of the tested epoxy resin

Material property Value
Elastic modulus (GPa), E 2.41 ± 0.02 Poisson's ratio, ν 0.33 ± 0.01 Ultimate tensile strength (MPa), σ u 71.2 ± 3.5 Yield strength (MPa), σ Y 61.1 ± 2.8 Elongation at ultimate tensile strength (%) 12 ± 0.9 Fig. 1 The tensile stress-strain curve for the epoxy material the circular arcs is 10 mm. Therefore, the length of each crack (2a) can be computed as a = R sin β. The height of specimens is 14 cm, and the widths of them are about 9 and 12 cm for parallel and collinear configurations, respectively. The values of 0.6 and 0.9 are also considered for 2a/d to have better realization of cracks interaction in the samples during the loading. As the geometries and loading conditions of the specimens are symmetric, only the tips of one crack are labeled by A and B in Fig. 2.
According to the ASTM D5045-99 (1999), reliable results are achieved by repeating each fracture test for three times. Consequently, 24 multi-cracked configurations are considered and 72 specimens are tested under remote tension with a uniaxial loading machine by applying a constant loading rate of 1 mm/min. The load-displacement curve of each tested sample is prepared by the testing set-up and the load-carrying capacity (LCC), the load at which cracks start propagating, of each specimen is extracted from these curves and reported in Tables 2, 3, and 4. The average of these LCCs P i (i=1, 2, 3) is presented by P avg in Tables 2, 3, and 4.

Theoretical estimation
Combination of the equivalent material concept (EMC) with fracture criterion The aim of this research is to prevent from complicated and time-consuming ductile analyses. Therefore, the equivalent material concept (EMC) which was proposed by Torabi (2012) is used in the current study. The EMC is a new approach that has been frequently applied to predict the ductile fracture of metallic and polymeric materials in several researches (Torabi and Alaei 2015;Cicero et al. 2017;Berto and Razavi 2018;Majidi et al. 2018;Torabi et al. 2016;Torabi et al. 2017;Torabi et al. 2018;Rahimi et al. 2018). It is clear that fracture evaluation of a brittle material is much simpler than that of a ductile material. Therefore, according to EMC, a virtual brittle medium with quite linear elastic behavior is studied instead of the real ductile material (Torabi 2012).
The same elastic modulus and fracture toughness are considered for both of the real ductile material and virtual brittle material. However, these materials have different tensile strengths which should be computed according to the following procedure. In the first step, the strain energy densities (SED) of both media are supposed to be equal. Therefore, as shown in Fig. 4, the SEDs of brittle and ductile materials are achieved by computing the areas under their true tensile stress-strain curves until the peak point. By equalizing the obtained values, the following expression is attained (Torabi and Alaei 2015): where E and σ Ã f are the elastic modulus and the tensile strength of the virtual brittle material, respectively.  Table 3 The experimental LCCs for multi-cracked specimens with vertical cracks alignment ( Fig. 2b)  By performing a simple mathematical work on Eq.
(1), the tensile strength of the virtual brittle material can be computed as follows: Then, by substituting the calculated area for the ductile material into Eq. (2), the tensile strength of the virtual brittle material is obtained. Now, to provide the theoretical prediction for the ductile failure of the multi-cracked components, the computed value σ Ã f is applied to the brittle fracture criterion which is described as follows.
The generalized maximum tangential criterion declares that the crack will propagate in the material when the maximum tangential stress along θ 0 attains a critical value. Now, by considering the singular and T-stress terms of tangential stress around the crack tip, the maximum tangential stress at the border of crack tip is obtained by using Eq. (3).
The parameter θ 0 and r c indicate the fracture initiation angle and the critical distance of material, respectively. Smith et al. (2001) derived Eq. (4) which should be satisfied simultaneously with Eq. (3) to have crack propagation under mixed mode I/II loading condition. The details about how these two equations are derived are presented by Smith et al. (2001).
The K Ic is the plane-strain fracture toughness of material.
The critical distance of material can be computed as r c Torabi et al. 2017). The parameter σ u is the ultimate tensile strength of material. Now, the tensile strength σ Ã f obtained based on the EMC should be applied to the critical distance expression instead of σ u . After that, the Eqs. (3) and (4) should be simultaneously solved to draw the theoretical fracture curve. In the next subsection, the dislocation method is described which is used to obtain the experimental critical SIFs K I and K II by using the modified experimental critical loads in the obtained formulations.

The dislocation method
The problem of an elastic isotropic plane containing multiple cracks subjected to in-plane external loads is considered in this paper. For this aim, an analytical procedure based on the distributed dislocation technique (Hills et al. 1996) is employed. Weertman and Weertman (1992) analyzed the stress distribution in a plane containing climb and glide edge dislocations with Burgers vectors, B x and B y , respectively. The stress field corresponding to the above-mentioned dislocations can be expressed as: From Eq. (5), we may observe that stress components exhibit the familiar Cauchy-type singularity at dislocation location. Classical stress fields of the dislocations contain singularity which results in singular integral equations in distributed dislocation technique. The stress components caused by the climb and glide edge dislocations located at a point with coordinates (x 0 , y 0 ), read σ x x; y ð Þ σ y x; y ð Þ σ xy x; y ð Þ 8 < : By using the distributed dislocation technique, arbitrary shapes of cracks can be modeled (Hills et al. 1996). Consider an isotropic plane weakened by N curved cracks which can be described in parametric form as For analyzing curved cracks, the moveable orthogonal coordinate (s, n) is chosen such that the origin may move on the crack while s-axis remains tangent to the crack face. Hence, the normal and tangent stresses in the (s,n) coordinates can be written as And the transformations of Burgers vectors on the surface of a crack are The angle between s-and x-axes is denoted by ðtÞ ¼ tan − 1 β 0 ðtÞ α 0 ðtÞ , where the prime denotes differentiation with respect to the argument. Suppose climb and glide edge dislocations with unknown densities b nj (t) and b sj (t), respectively, are distributed on the infinitesimal segment ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ½α 0 j ðtÞ 2 þ ½β 0 j ðtÞ 2 q dt at the surface of jth crack, where −1 ≤ t ≤ 1. Employing the principal of superposition, the components of traction vector at a point with coordinates (α j (ξ), β j (ξ)), where the parameter −1 ≤ ξ ≤ 1, on the surface of all cracks yield The kernels in Eq. (10) are defined by the following expressions: The crack opening displacement across the jth crack can be formulated by applying the definition of dislocation density function, as follows: where u n and u s indicate the normal and tangential components of displacement, respectively. On the other hand, the single valued condition of displacement field out of a crack surface leads to new closure conditions expressed in Eq. (13).
Now, the system of Cauchy singular integral containing Eqs. (10) and (13) must be solved simultaneously in order to attain the dislocation densities on a crack face. Due to the presence of square root singularity for the stress field at an embedded crack tip, the dislocation density for embedded cracks is formulated as (Pourseifi and Faal 2015) In order to obtain g nj (t) and g sj (t), the Eq. (14) should be substituted into Eqs. (10) and (13). Erdogan et al. (1973) developed the numerical technique of integral equations with Cauchy-type kernel to solve the obtained system. Fotuhi and Fariborz (2008) derived the Eq. (15) for computing the stress intensity factors (K I and K II ), based on the dislocation density functions.
After preparing the formulations of dislocation approach, the MATLAB code is used to obtain the experimental stress intensity factors. For this aim, the critical loads obtained from the fracture tests should be modified by EMC as it is described in the next section and then these corrected loads are used as the far field load in the MATLAB code.

Numerical analysis
To compare the experimental and theoretical results, the fracture toughness of the tested material should be determined. For this aim, the commercial code ABA-QUS/CAE 6.10 is used to perform finite element (FE) analysis. Based on the geometries of tested multicracked specimens, plane-stress models are created in the two-dimensional (2D) space. In order to achieve a suitable convergence of the numerical results, various analyses with different numbers of elements are conducted on the FE models. About 26,000 eight-node plane-stress quadratic elements with reduced integration are employed in the FE models with the minimum size about 0.01 mm. Figure 5a displays the typical mesh pattern applied in the whole model and Fig. 5b also shows the quite fine meshes used around the crack tips. In fact, singular elements in the form of a sweep mesh pattern are employed around the crack tips due to the existence of high stress gradient in these regions. The distributed tensile load is applied to the top line of the specimen model and the opposite line located at the bottom of the model is fixed to create the appropriate loading and boundary conditions. The FE software is utilized to analyze the multi-cracked models with β = 20, 40, and 60 (deg.) and vertical configuration.
To attain the critical SIF of the material (i.e., the K c value), the obtained experimental load-carrying capacity (LCC) for the single cracked sample corresponded to the pure mode I loading should be applied to the FE model. As the epoxy material employed in this work has a ductile behavior, its experimental LCC should be modified to use in FE analysis. To achieve this aim, first the ductile epoxy material is equated with the brittle one. In fact, the real nonlinear load-displacement curve of the ductile polymer is replaced with a linear curve. Therefore, the slope of the linear portion of the real loaddisplacement curve is obtained and the virtual linear curve is drawn with the same slope (see Fig. 6). Then, the value of the fracture load P f in the virtual linear curve can be computed by equalizing the strain energies absorbed by the real cracked specimen and the virtual brittle material (the area under the load-displacement curve until the peak point as shown in Fig. 6). Now, the obtained virtual fracture load should be applied to the created finite element model to achieve the critical SIF K c by using the linear elastic analysis. Consequently, by following the described procedure, the virtual value of the critical SIF K c for the tested polymer is attained to be equal to 1.346 MPa ffiffiffiffi m p .

Results and discussion
There are no analytical formulations to obtain the values of the critical SIFs K I and K II for multi-cracked specimens; therefore, the modified experimental loads are used in the dislocation method to calculate these parameters. The experimental critical SIFs for all the tested specimens are normalized with K IC and depicted in Fig. 7. It is clear that by increasing the circular arc angle, the contribution of mode II stress intensity factor also enhances. To compare the experimental and theoretical SIFs, the theoretical curves are also obtained based on the Eq. (3) and (4), which are presented in Fig. 7. For ductile multi-cracked specimens, the proposed criterion shows a good agreement with the experimental results. It can be seen in Fig. 7 that the mean values of experimental results located under the theoretical curves, meaning that the conservative predictions are obtained in the safe zone of fracture curve.
To prepare better comparison for the fracture resistance of considered specimens in various mode mixity, the effective fracture toughness, K eff , is employed which is defined as follows: Indeed, the effects of the mode I and mode II SIFs take into account simultaneously by applying the effective fracture toughness. As the three considered configurations have a similar trend for variation of K eff , only a sample of these charts is presented in Fig. 8. It is clear in Fig. 8 that K eff for the fractured specimen grows by increasing the contribution of mode II loading on the crack tips. This enhancement of effective fracture toughness also can be realized by considering the obtained LCCs (see Tables 2, 3, and 4) which shows that much more loads requires for the failure of specimens with larger circular arc angle. This observation is mainly due to the existence of larger plastic zone around the crack tips by increasing the value of K II . A considerable amount of energy dissipates in the plastic zone during the loading; therefore, a larger plastic zone leads to an increase in loading carrying capacity and fracture toughness of multi-cracked specimen. The importance of plastic deformation in the fracture of polymeric materials has also been demonstrated in other studies (Rizov 2017;Saboori and Ayatollahi 2017).
By comparing the obtained LCCs for cracked specimens with parallel cracks (Fig. 2b, c), it is found that for the same value of 2a d , the specimens with (c) configuration have greater fracture loads which means that these specimens can undergo higher loads than the specimens with (b) configuration. This observation is due the interactions between the cracks tips in specimens; such that the contribution of mode II loading condition in the specimens with (c) configuration is much more dominant than that for specimens with (b) configuration.
Shear yielding and crazing are the two main micromechanical mechanisms leading to the existence of plastic deformation in the vicinity of crack tips in components (Kinloch 1985;Liu et al. 1998). Micro voids are the main cause of crazing which exist in a plane perpendicular to the maximum principle stress. Shear yielding happens around the crack tips when the shear stress in this region reaches its critical value. Shear bands, which relate to this phenomenon, propagate from the crack tips along the direction in which the shear stress is maximum. However, crazing is rarely seen in thermosetting polymers such as epoxy; shear yielding phenomenon is the major cause of plastic deformation in the epoxy materials. As the shear stresses relating to the contribution of mode II increase around the crack tips, much more energy dissipates by formation of shear bands. Therefore, lower energy is available for crack propagation from the crack tips and consequently, the fracture resistance of cracked sample becomes greater. It is clear that by considering the higher order terms of Williams' series in GMTS criterion, fracture parameters can be predicted more accurately and the geometry effect becomes negligible. These observations were also reported in other researches in which other criteria with higher order terms of Williams' series were used to provide the fracture predictions of brittle material Razavi et al. 2018).
To have better realization of plastic deformation in tested specimens, size of the plastic zone at the neighborhood of crack tip is specified by performing the elastic-plastic FE analysis. For this aim, the elastic-plastic stress-strain curve of the epoxy material is used to model the real ductile specimen. In fact, several points of the nonlinear portion in real stress-strain curve of the ductile specimen are extracted and used as a data table to model the plastic material in the software. The values of elastic modulus, Poisson's ratio, and yield stress used in the simulation, are also reported in Table 1. The FE software is utilized to analyze the multi-cracked models with β = 20, 40, and 60 (deg.). The sizes of the plastic deformation are determined in Fig. 9. The considerable plastic zone around the crack tip in the multi-cracked specimens can be seen in Fig. 9. By comparing the sizes of determined plastic regions in Fig. 9, it is found that by increasing the angle of circular arc crack β, the plastic zone becomes greater. This observation is due to the enhancement of mode II contribution in the crack tip vicinity at the onset of the sample failure which leads to the existence of greater shear deformations in this region.

Conclusions
For the first time, the ductile failure of multi-cracked specimens containing circular arc cracks is investigated experimentally and theoretically in this study. For this purpose, the multi-cracked specimens are prepared from specified polymeric material in various cracks arrangements. Different angles are considered for circular arc cracks to cause mixed mode loading condition at the crack tips in tested specimens. Fracture tests are conducted on the prepared samples and their failure loads are obtained. To prevent from elastic-plastic analysis, theoretical predictions are provided by applying the equivalent material concept (EMC) in conjunction with the brittle fracture criterion. Fracture curves are prepared and the theoretical and experimental results are compared with each other. The shear yielding phenomenon, which is the main micromechanical reason for the formation of plastic region in the vicinity of crack tips of the tested specimens, is also described. Sample elastic-plastic analysis is also utilized to have better realization of plastic deformations. Finally, it is revealed that theoretical predictions can provide an appropriate consistency with the experimental results.