Finite Fracture Mechanics from the macro- to the micro-scale. Comparison with the Phase Field model

ScienceDirect Available online at www.sciencedirect.com Procedia Structural Integrity 42 (2022) 553–560 2452-3216 © 2022 The Authors. Published by Elsevier B.V. This is an open access article under the CC BY-NC-ND license (https://creativecommons.org/licenses/by-nc-nd/4.0) Peer-review under responsibility of the scientific committee of the 23 European Conference on Fracture – ECF23 10.1016/j.prostr.2022.12.070 10.1016/j.prostr.2022.12.070 2452-3216 © 2022 The Authors. Published by Elsevier B.V. This is an open access article under the CC BY-NC-ND license ( https://creativecommons.org/licenses/by-nc-nd/4.0 ) Peer-review under responsibility of the scientific committee of the 23 European Conference on Fracture – ECF23 Available online at www.sciencedirect.com ^ĐŝĞŶĐĞŝƌĞĐƚ Structural Integrity Procedia 00 (2019) 000–000 www.elsevier.com/locate/procedia 2452-3216 © 2020 The Authors. Published by Elsevier B.V. This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/4.0/) Peer-review under responsibility of 23 European Conference on Fracture - ECF23 23 European Conference on Fracture - ECF23 Finite Fracture Mechanics from the macro- to the micro-scale. Comparison with the Phase Field model. S. Jiménez-Alfaroa, J. Reinosob, D. Leguillona, C. Maurinia aInstitut Jean Le Rond d’Alembert, Sorbonne Université – CNRS UMR 7190, 4 Place Juisseu, 75005, Paris, France bGrupo de Elasticidad y Resistencia de Materiales, Escuela Técnica Superior de Ingenieria de Sevilla, Universidad de Sevilla , Camino de los Descubrimientos, 41092, Sevilla, Spain Abstract In the Phase Field methodology there is a length parameter related to the size of the damage region, which is proportional to the Irwin length. In this work, we apply the Phase Field methodology to the analysis of fracture at the micro-scale in brittle materials, studying the role of this length parameter when the size of the specimen is of the order of the Irwin length. We compare the conclusions obtained with the ones previously presented using Finite Fracture Mechanics. Two different schemes are used, a staggered and a monolithic scheme. © 2020 The Authors. Published by Elsevier B.V. This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/4.0/) Peer-review under responsibility of 23 European Conference on Fracture - ECF23 Keywords: Micro-scale – Fracture properties – Phase Field – Brittle fracture 1. Main text One of the fundamental ingredients of the recent methodology Finite Fracture Mechanics (FFM) for the prediction of crack events in materials and structures, is that a crack is assumed to jump a given finite length at onset. This can be formulated through the invocation of the Coupled Criterion (CC) presented in Leguillon (2002), stating that this length depends on the material toughness, the tensile strength but also the geometry. Complying with a different vision, according to Pham et al. (2011) and Bourdin et al. (2014), in the Phase Field (PF) of fracture, there exist a materialrelated length that describes the size of the damaged region, called the phase field length scale. Both the PF length scale and the nucleation length obtained by the CC are proportional to the Irwin length defined from the material toughness and the tensile strength. At the macro-scale, they are small compared to any dimension of the structure, Available online at www.sciencedirect.com ^ĐŝĞŶĐĞŝƌĞĐƚ Structural Integrity Procedia 00 (2019) 000–000 www.elsevier.com/locate/procedia 2452-3216 © 2020 The Authors. Published by Elsevier B.V. This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/4.0/) Peer-review under responsibility of 23 European Conference on Fracture - ECF23 23 European Conference on Fracture - ECF23 Finite Fracture Mechanics from the macro- to the micro-scale. Comparison with the Phase Field model. S. Jiménez-Alfaroa, J. Reinosob, D. Leguillona, C. Maurinia aInstitut Jean Le Rond d’Alembert, Sorbonne Université – CNRS UMR 7190, 4 Place Juisseu, 75005, Paris, France bGrupo de Elasticidad y Resistencia de Materiales, Escuela Técnica Superior de Ingenieria de Sevilla, Universidad de Sevilla , Camino de los Descubrimientos, 41092, Sevilla, Spain Abstract In the Phase Field methodology there is a length parameter related to the size of the damage region, which is proportional to the Irwin length. In this work, we apply the Phase Field methodology to the analysis of fracture at the micro-scale in brittle materials, studying the role of this length parameter when the size of the specimen is of the order of the Irwin length. We compare the conclusions obtained with the ones previously presented using Finite Fracture Mechanics. Two different schemes are used, a staggered and a monolithic scheme. © 2020 The Authors. Published by Elsevier B.V. This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/4.0/) Peer-review under responsibility of 23 European Conference on Fracture - ECF23 Keywords: Micro-scale – Fracture properties – Phase Field – Brittle fracture 1. Main text One of the fundamental ingredients of the recent methodology Finite Fracture Mechanics (FFM) for the prediction of crack events in materials and structures, is that a crack is assumed to jump a given finite length at onset. This can be formulated through the invocation of the Coupled Criterion (CC) presented in Leguillon (2002), stating that this length depends on the material toughness, the tensile strength but also the geometry. Complying with a different vision, according to Pham et al. (2011) and Bourdin et al. (2014), in the Phase Field (PF) of fracture, there exist a materialrelated length that describes the size of the damaged region, called the phase field length scale. Both the PF length scale and the nucleation length obtained by the CC are proportional to the Irwin length defined from the material toughness and the tensile strength. At the macro-scale, they are small compared to any dimension of the structure, 554 S. Jiménez-Alfaro et al. / Procedia Structural Integrity 42 (2022) 553–560 2 Author name / Structural Integrity Procedia 00 (2019) 000–000 whereas at the micro-scale both lengths are the same order of magnitude or even larger and can interact with the dimensions of the structure. In Jiménez-Alfaro and Leguillon (2021) the authors examine the issue via the CC when descending the scales from the cm-scale to the µm-scale. In line with this previous investigation, the aim of this work is to make a comparison between the conclusions obtained using the PF model and the ones previously presented with the CC. Based on previous results, it can be argued that both the CC and the PF model provide satisfactory predictions of cracking events in solids. However, this can be a controversial issue at smaller scales of analysis due to a lack of energy because of the smallness of the specimens. This is attributed to the fact that at such scales it is seen that the corresponding results are much sensitive to the toughness but less sensitive to the tensile strength. Relying on the previous discussion, in this contribution, bending tests on notched micro-cantilever beams (Figure 1) made from a ceramic material 8Y-FSZ cubic zirconia are investigated. Particularly, we analyze how the fracture nucleation is affected considering different values of the toughness and the strength, using the PF approach as a modelling tool. A staggered scheme in the PF method is applied in the analysis. Then, the conclusions obtained using PF, in terms of the critical load and the critical displacement are compared to the ones obtained considering the CC, which were previously presented in Jiménez-Alfaro and Leguillon (2021). The influence of the phase field length scale in the crack nucleation at the micro-scale is examined, considering this parameter a property of the material because of its relation to the Irwin length, as it was already mentioned. Nomenclature Regularized energy Strain energy density 𝛾𝛾 Geometric crack function Wext Work done by external forces Displacements field Damage variable Strain tensor Degradation function Gc Critical energy release rate c Tensile strength 0+/- Strain energy density without damage for tension (+) and compression (-) Young’s modulus Poisson’s ratio Lamé parameters Phase Field length scale 2. Review of the Phase Field methodology One of the main characteristics in the Phase Field approach of fracture is the fact that the sharp crack is modelled as a diffuse discontinuity, being characterized by a continuous variable which defined the smooth transition between the completely broken and unbroken states. The main advantage of this method is its versatility since it allows predicting crack nucleation and propagation in a wide range of engineering applications. This methodology is based on the variational approach to fracture introduced in the late 90s by Francfort and Marigo (1998), and then regularized by Bourdin et al. (2000) for brittle fracture using the Ambrosio and Tortorelli’s regularization of the Mumford-Shah problem in image processing given in Ambrosio and Tortorelli (1990). This variational approach, based on Griffith (1921) vision of fracture, is described by the minimization of the following regularized energy 𝐸𝐸𝑙𝑙, presented in Eq. (1). Among others, the theoretical description about this methodology given in this work is based on the work presented in Tanné et al. (2018) and Molnar et al. (2020). S. Jiménez-Alfaro et al. / Procedia Structural Integrity 42 (2022) 553–560 555 Author name / Structural Integrity Procedia 00 (2019) 000–000 3 , (1) where denotes the strain energy density stored in the solid, is the so-called geometric crack function and is the work done by external forces. The regularized energy depends on the displacements field and the damage variable , ranging from 0 (no damage) to 1 (total damage). Therefore, the strain energy density is defined in Eq. (2). A degradation function is defined to represent the stiffness reduction when the damage is increased, in this case, by a quadratic function . A very low numerical constant is used to avoid problems in convergence during the simulation. . (2) The strain energy density without damage stored in the solid is split into and , being respectively related to traction and compression strain components. In this work the strain energy split considered is based on the volumetric-deviatoric decomposition presented in Amor et al. (2009). In Eq. (3) and Eq. (4) the formulation for and are respectively given. In these equations represents the bulk modulus and is the deviatoric part of the strain tensor. , (3) , (4) where is defined in Eq. (5). . (5) On the other hand, the geometric crack function is described in Eq. (6) by the so-called phase field length scale , which defines the size of the damage region. This function is differently defined for the AT1 and AT2 models, see Tanné et al. (2018). In the AT1 model (Pham et al. (2011)), an initial elastic threshold is considered (free of damage, it means ), whereas in the AT2 model (Bourdin et al. (2014)) no linear-elastic stage is exhibited before the peak stress. In Eq. (6) the expression is given for the AT1 model, since it is the one we are going to apply in this work. The phase field length scale is considered a property of the material, related to the critical energy release rate and the tensile strength, as it is written in Eq. (7) for plane strain and the AT1 model. , (6) . (7) Due to the complexity of this minimization problem, a staggered solution is proposed to find a solution, where the displacements field and the damage variable are separately solved. At first, the displacements are found fixing the damage variable. Then, the damage variable is solved fixing the displacements. Two methodologies at this point are proposed. On the one hand, the one proposed in Marigo et al. (2016), in which the condition of irreversibility is directly imposed. On the other hand the methodology given in Miehe et al. (2010), in which the so-called history field (Η) is defined. 3. Results In this section we show the results obtained using the PF model. At first, in section 3.1. we introduce the results with the first methodology mentioned in the previous section. In section 3.2. the second methodology is presented. In 556 S. Jiménez-Alfaro et al. / Procedia Structural Integrity 42 (2022) 553–560 4 Author name / Structural Integrity Procedia 00 (2019) 000–000 this section we set the numerical parameters of the methodology through a convergence analysis with three conditions. We define the time step of the simulation ensuring the convergence of three different conditions. Finally, in section 3.3 a comparison of the conclusions brought by the two methodologies with the ones obtained by the CC, see Jiménez- Alfaro and Leguillon (2021), is made. The problem analyzed in this work was experimentally solved by Henry et al. (2020), and then it was numerically studied in the Finite Fracture Mechanics framework by the application of the Coupled Criterion (CC) in Jiménez-Alfaro and Leguillon (2021). In all these tests a ceramic material 8Y-FSZ cubic zirconia (𝐸𝐸 = 216𝐺𝐺𝐺𝐺𝐺𝐺 and 𝜈𝜈 = 0.3) was investigated by performing bending tests on notched micro-cantilever beams. This material is used to fabricate nuclear fuels. In the experiments, 14 specimens were tested. In this work we present the results related the geometrical parameters of the first one. From numerical studies in the FFM framework, an estimation on the fracture properties was proposed, as well as an equivalent 2D model to represent the bending test (see Figure 1). In this paper we consider the same 2D model, even though the methodology applied to get the solution is the PF model. Moreover, displacement control is assumed in the whole analysis. A mesh ensuring that 5 elements define the phase field length scale in the prescribed damage region is set. 3.1. Solution with the first methodology of the PF model In the first methodology of the PF model that we have applied, the irreversibility condition in the damage variable is directly imposed in the minimization problem, see Marigo et al. (2016). The code is developed in FENICSX – DOLFINX. An example case is represented in this section, for a value of the tensile strength equal to 5000 MPa and a value of the critical energy release rate equal to 10 MPa µm. The imposed displacement range analyzed is based on the results obtained in the FFM framework, µm. In Figure 2a the damage variable is represented for µm, the critical point at which the crack is nucleated. In Figure 2b the force-displacement diagram for this case is presented. When the crack is nucleated, the force is highly reduced. Fig. 2. Analysis with the first methodology of the PF model. In (a) the force-displacement curve. In (b) the deformed shape with the damage variable, image obtained with Paraview software. These conditions are analyzed for σc=5000 MPa and Gc=10 MPa ∙ µm. 3.2. Solution with the second methodology of the PF model In this section we present the results we have obtained with the second methodology mentioned in this work, in which the so-called history variable is introduced to solve the minimization problem. We apply the existing analogy between the heat transfer equation and the staggered problem introduced in Navidtehrani et al. (2021), and developed in a UMAT subroutine for ABAQUS. Therefore, the solution is found solving a coupled displacement-temperature problem with a steady-state response, considering some particularities explained in Navidtehrani et al. (2021). Figure 1 : 2D model proposed in Jiménez-Alfaro and Leguillon (2021). 0 0.5 1 1.5 2 00.05 0.1 0.15 0.2 F [mN] U [µm] S. Jiménez-Alfaro et al. / Procedia Structural Integrity 42 (2022) 553–560 557 Author name / Structural Integrity Procedia 00 (2019) 000–000 5 Since the staggered scheme is highly dependent on the time step, ensuring the convergence in the analysis is one of the most important steps. There are three conditions that can be applied to find the critical point at which a crack nucleation occurs. The condition 1, in which we can consider that the crack nucleation is given when the maximum force in the force-displacement diagram is reached, see Nguyen et al. (2016) and Yin and Zang (2019). Such condition is applied in experiments performed in Henry et al. (2020) to find the critical force. However, with the PF model we can also assume that the failure is initiated when the damage variable reaches for the first time the value of 1, which would be the condition 2. Finally, the condition 3 arises when the stiffness starts decreasing, as it is mentioned in Marigo et al. (2016). In Figures 3a, 3b and c the conditions are analyzed for a certain case, 𝜎𝜎𝑐𝑐=9000MPa and 𝐺𝐺𝑐𝑐=10MPa ∙ µm. In each Figure the evolution of the force, the damage variable and the stiffness with respect to the imposed displacement is shown for several values of the time step 𝑑𝑑𝑑𝑑 =10−5 … 5 ∙ 10−3. The critical points corresponding to the three conditions explained above were highlighted in the graphics. Fig. 3. Representation of the three conditions to determine the damage initiation. In (a) the diagram force-displacement; (b) the diagram phase field-displacement, and (c) the diagram stiffness-displacement. These conditions are analyzed for σc=8000 MPa and Gc=10 MPa ∙ µm, for several values of the simulation time step. Finally, in Figures 4a and 4b the critical displacement and the critical force with respect to the time step for each condition is represented, in a logarithmic scale. It can be concluded that when the time step is reduced the three conditions converge to the same critical point. After making the same analysis for each possible value of the strength and the critical energy release rate we can conclude the time step dt = 10−5𝑠𝑠 is enough to ensure the convergence of the results. This is particularly important in this first study since a staggered model of the problem is considered. Fig. 4. Representation of the (a) the critical displacement and (b) the critical force according to conditions 1, 2 and 3, with respect to the simulation time step, under a logarithmic scale. a 0 1 2 3 00.08 0.16 0.24 0.32 F [mN] U [µm] b 0 0.2 0.4 0.6 0.8 1 1.2 00.08 0.16 0.24 0.32 α[-] U [µm] c 6 8 10 12 14 16 00.08 0.16 0.24 0.32 K [mN/µm] U [µm] a 0 0.05 0.1 0.15 0.2 0.25 100 1000 10000 100000 Ucrit [µm] 1/dt [1/s] Condition 1 Condition 2 Condition 3 b 1 1.5 2 2.5 3 100 1000 10000 100000 F crit [mN] 1/dt [1/s] 558 S. Jiménez-Alfaro et al. / Procedia Structural Integrity 42 (2022) 553–560 6 Author name / Structural Integrity Procedia 00 (2019) 000–000 3.3. Comparison with the Coupled Criterion In this section we compare the results obtained in Jiménez-Alfaro and Leguillon (2021) using the CC with the ones obtained applying the two methodologies of the PF model. In Figure 5 and 6 the critical displacement and force calculated for several values of the tensile strength are given for both the CC and the first methodology of the PF model. Notice that the same critical energy release rate c is fixed, in order to study the influence of the tensile strength, ranged as c . It is observed a little variation in both the critical displacement and the critical force is observed when the tensile strength is changed. For this reason, we can conclude that experiments at the micro-scale seem poorly appropriate to estimate the tensile strength, as it was explained in Jiménez- Alfaro and Leguillon (2021). Notice that when descending the tensile strength, the difference between the second methodology used in PF and the other two methods (CC and the first method in PF). When the tensile strength is lower, according to Eq. (7) the phase field length scale is higher. For the lowest value of c the phase field length scale is of the order of the dimensions of the specimen, . Therefore, the convergence could be affected in the PF model. However, the first methodology seems to show similar results to the ones predicted by the CC. Another reason that could explain this phenomenon could be the time convergence in the second method of the PF model, which is not necessary in the first one. Fig. 5. Representation of the critical displacement with respect to the tensile strength, considering both the CC and the PF model, for specimen number 1, for the same value of the critical energy release rate 𝐺𝐺𝑐𝑐=10 MPa µm. Fig. 6. Representation of the critical force with respect to the tensile strength, considering both the CC and the PF model, for specimen number 1, for the same value of the critical energy release rate 𝐺𝐺𝑐𝑐=10 MPa µm. 0 0.04 0.08 0.12 0.16 1000 2000 3000 5000 7000 9000 U [µm] 𝜎𝜎c [MPa ] CC PF-1 PF-2 0 0.5 1 1.5 2 1000 2000 3000 5000 7000 9000 F [mN] 𝜎𝜎c [MPa ] CC PF-1 PF-2 S. Jiménez-Alfaro et al. / Procedia Structural Integrity 42 (2022) 553–560 559 Author name / Structural Integrity Procedia 00 (2019) 000–000 7 On the other hand, in Figures 7 and 8 the critical displacement and the critical force for several values of the critical energy release rate are given, following the three methodologies. In this study, the range in the critical energy release rate is cFurthermore, the tensile strength is fixed here as c . It is observed that both the critical displacement and the critical force are highly dependent on the critical energy release rate, as it was concluded in Jiménez-Alfaro and Leguillon (2021). Similar results can be obtained by the three methodologies. In both analysis changing the fracture properties is equivalent to change the phase field length scale, according to the AT1 model. Therefore, for the tensile strength the phase field length scale changes from 0.01 µm to 0.8 µm, whereas in the second study, where the critical energy release rate was varied, the phase field length scale ranges from 0.01µm to 0.06µm. However, the critical force is increased between 1.8-2 times when the critical energy release rate is increased, whereas the ratio between the highest and the lowest value of the critical force when increasing the tensile strength is close to 1. Similar conclusions can be obtained for the critical displacement. Fig. 7. Representation of the critical displacement with respect to the critical energy release rate considering both the CC and the PF model, for specimen number 1. Fig. 8. Representation of the critical force with respect to the critical energy release rate considering both the CC and the PF model, for specimen number 1. 4. Conclusions To conclude, both the critical displacement and the critical force are highly dependent on the critical energy release rate, whereas they are little dependent on the tensile strength. According to the CC, the reason is that it is the energy 0 0.5 1 1.5 2 2.5 3 510 15 20 F [mN] CEER [MPa µm] CC PF-1 PF-2 0 0.04 0.08 0.12 0.16 0.2 510 15 20 U [µm] CEER [MPa µm] CC PF-1 PF-2 560 S. Jiménez-Alfaro et al. / Procedia Structural Integrity 42 (2022) 553–560 8 Author name / Structural Integrity Procedia 00 (2019) 000–000 condition the one that is governing the failure. In the PF model, a change in the fracture properties is related to a change in the phase field length scale, that is proportional to the Irwin’s length. Moreover, as a general observation the same conclusions can be obtained from the two different PF models and the FFM. Moreover, when the PF length scale is of the order of the dimensions of the specimen the difference in the results obtained with the two PF methods is higher. In this work the CC is presented as a first step of the PF model, since the load range used in the simulations are set using the results brought by the CC. However, it is not always possible to follow that order. When the prescribed crack path is not clearly defined a priori, the PF model could also be a good first step of the CC. For future studies, the authors propose to make a dimensional analysis of the problem, following the Pi- Buckingham Theorem, to study the influence of different parameters. Moreover, a comparison between the AT1 and the AT2 model could be made, to check if the same conclusions that we obtain with the AT1 model are obtained with the AT2. Finally, the influence of the notch radius on these conclusions could be studied. Acknowledgements The authors acknowledge the funding received from the European Union’s Horizon 2020 research and innovation programme under the Marie Sklodowska-Curie grant agreement No 861061 – Project NEWFRAC. The authors are also indebted to Camila Zolesi for their knowledge and their help with the first Phase Field methodology. References Ambrosio, L., Tortorelli, V.M..1990. Approximation of functionals depending on jumps by elliptic functionals via Γ-convergence. Communication on pure and applied mathematics. 43(8), 999-1036. Amor, H., Marigo, J.J., Maurini, C.. 2009. Regularized formulation of the variational brittle fracture with unilateral contact: Numerical experiments. Journal of the Mechanics and Physics of Solids. 57(8), 1209-1229. Bourdin, B., Francfort, G.A., Marigo, J.J.. 2000. Numerical experiments in revisited brittle fracture. Journal of the Mechanics and Physics of Solids. 48(2000), 797-826. Bourdin, B., Marigo, J.J., Maurini, C., Sicsic, P..2014. Morphogenesis and propagation of complex cracks induced by thermal shocks. Physical review letters. 112(1), 014301. Francfort, G.A., Marigo, J.J..1998. Revisiting brittle fracture as an energy minimization problem. Journal of the Mechanics and Physics of Solids. 46(8), 1319-1342. Griffith A.A..1921. VI. The phenomena of rupture and flow in solids. Philosophical transactions of the Royal Society of London. 221(582-593), 98-163. Henry, R., Zacharie-Aubrun, I., Blay, T., Chalal, S., Gatt, J.M., Langlois, C., Meille, S.. 2020. Fracture properties of an irradiated PWR UO2 fuel evaluated by micro-cantilever bending tests. Journal of Nuclear Materials. 538, 152209. Jiménez-Alfaro, S. Leguillon, D.. 2021. Finite Fracture Mechanics at the micro-scale. Application to bending tests of micro-cantilever beams. Engineering Fracture Mechanics, 258, 108012. Leguillon, D. 2002. Strength or toughness? A criterion for crack onset at a notch. European Journal of Mechanics-A/Solids, 21(1), 61-72. Marigo, J.J., Maurini, C., Pham, K..2016. Gradient damage models and their use in brittle fracture. Miehe, C., Hofacker, M., Welschinger, F.. 2010. A phase-field model for rate independent crack propagation: Robust algorithmic implementation based on operator splits. Computer Methods in Applied Mechanics and Engineering. 199, 2765-2778. Molnar, G., Doitrand, A., Estevez, R., Gravouil, A.. 2020. Toughness or strength? Regularization in phase-field fracture explained by the coupled criterion. Theoretical and applied fracture mechanics. 109, 102736. Navidtehrani, C., Betegón, C., Martinez-Pañeda, E.. 2021. A unified Abaqus implementation of the phase field fracture method using only a user material subroutine. Materials. 14(8), 1913. Nguyen, T.T., Yvonnet, J., Bonert, M., Chateau, C., Sab, K., Romani, R., Le Roy, R.. 2016. On the choice of parameters in the phase field method for simulating crack initiation with experimental validation. International Journal of Fracture. 197(2), 213-226. Pham, K., Amor, H., Marigo, J.J., Maurini, C.. 2011. Gradient damage models and their use to approximate brittle fracture. International Journal of Damage Mechanics. 20(4), 618-652. Tanné, E., Li, T., Bourdin, B., Marigo, J.J., Maurini, C.. 2018. Crack nucleation in variational phase-field models of brittle fracture. Journal of the Mechanics and Physics of Solids. 110, 80-99. Yin, B., Zhang, L.. 2019. Phase field method for simulating the brittle fracture of fiber reinforced composites. Engineering Fracture Mechanics. 221, 321-340. View publication stats