Adaptive phase-field modeling of heterogeneous materials

Treball realitzat per: Roger Llorenç i Vilanova Dirigit per: Antonio Rodríguez-Ferran Grau en: Enginyeria Civil (en anglès) Barcelona, juny 2023 Departament d’Enginyeria Civil i Ambiental TREBALL FINAL DE GRAU Adaptive phase-field modeling of heterogeneous materials Abstract Adaptive phase-field modeling of heterogeneous materials Roger Lloren¸c i Vilanova Understanding fracture behavior is crucial to optimizing the design of engineering structures and components, from civil infrastructure to aircraft components. The phase-field model has become popular among these models because it can explain fracture formation and propagation as a continuous and diffuse process. The goal of this thesis is to study the behavior of crack propagation in heterogeneous materials by performing numerical simulations and analyzing them. Heterogeneity is defined as the ratio of the fracture toughness between different parts of the material. An adaptive phase-field code for fracture in homogeneous materials has been extended to heterogeneous materials, and two-dimensional experiments have been carried out to explore the mechanical response, capabilities, and limitations of the adaptive phase-field modeling of fracture in heterogeneous materials. A qualitative assessment of the crack’s behavior when interacting with heterogeneities is performed, and some behaviors are related to the force-displacement plot. It is determined that the peak force and the total displacement generally increase when the ratio of fracture toughness increases, but this relation may not necessarily be true infinitely. Increasing the ratio may result in a change in the mode of fracture (from penetration of the inclusion to the deflection of it) and decrease the peak force and/or the total displacement. It is important to carry out several tests to determine what ratio best suits the desired objective and to avoid any unexpected behaviors. The length of deflection and penetration have been related to whether the crack will penetrate or deflect a more resistant region in the material. The crack has been successfully guided through less-resistant thin paths, modifying its natural fracture path for homogeneous material. And we can state that for a more resistant, constant area in the material, in the case where the crack has no chance other than penetrating through all of it, splitting this area will result in an increase in the peak force and a decrease in the total iii displacement, and concentrating all of this area will result in the opposite, a decrease in the peak force but an increase in the total displacement. Finally, we have tested the limitations of this model by performing a test that subjects the crack to an extreme ratio of the fracture toughness, leaving no room for maneuver other than going backward. The crack has managed to go backward, but the damage field showcases chaotic and unnatural behavior. iv Contents Abstract iii Contents v List of Figures vii List of Tables xiii 1 Introduction 1 1.1 Motivation.................................... 1 1.2 Stateoftheart ................................. 2 1.2.1 Deterministic heterogeneity . . . . . . . . . . . . . . . . . . . . . . 2 1.2.2 Probabilistic heterogeneity . . . . . . . . . . . . . . . . . . . . . . . 3 1.3 Goals and layout of this thesis . . . . . . . . . . . . . . . . . . . . . . . . . 4 2 Brief overview of the model 5 2.1 Hybrid phase-field model . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 2.2 Adaptive refinement based on Nitsche’s method . . . . . . . . . . . . . . . 7 2.3 One-dimensional tension test . . . . . . . . . . . . . . . . . . . . . . . . . . 11 v 2.3.1 Homogeneous tension test . . . . . . . . . . . . . . . . . . . . . . . 12 2.3.2 Heterogeneous tension test . . . . . . . . . . . . . . . . . . . . . . . 13 3 Introducing the adaptive refinement 15 4 Two-dimensional experiments 19 4.1 Arrest, deflection and penetration of the crack . . . . . . . . . . . . . . . . 21 4.2 Length of penetration vs length of deflection . . . . . . . . . . . . . . . . . 32 4.3 Crack behavior over thin paths . . . . . . . . . . . . . . . . . . . . . . . . 34 4.3.1 Triangular wave of less resistance . . . . . . . . . . . . . . . . . . . 34 4.3.2 Void inclusions to guide the crack . . . . . . . . . . . . . . . . . . . 36 4.4 Improving the resistance by placing inclusions . . . . . . . . . . . . . . . . 38 4.5 Capability of the crack to move backward . . . . . . . . . . . . . . . . . . 47 5 Concluding remarks and future work 51 5.1 Concludingremarks............................... 51 5.2 Futurework................................... 52 Bibliography 55 vi List of Figures 2.1 Figures belonging to a shear test with d∗= 0.2 in Muix´ı (2020) pages 61 — 62. a) Contour plots of the damage and b) Force-displacement plots for different mesh resolutions and refinement factors. . . . . . . . . . . . . . . 8 2.2 Two-dimensional tension test with a cut at mid-height that goes to the middle with the four elements surrounding the crack tip prescribed to be refined for two given increments. a) The deformed mesh including the standard and refined elements and a zoom. b) The deformed mesh including the standard and refined elements for a further increment and a zoom of only the refined elements. The two next elements have been automatically refined. The displacements have been amplified by a factor of 10. ..... 9 2.3 Geometry and boundary conditions for the one-dimensional tension test. . 12 2.4 One-dimensional tension test mesh. . . . . . . . . . . . . . . . . . . . . . . 12 2.5 Damage field for the one-dimensional homogeneous tension test. Fracture takes place on the left side of the material due to rounding errors. . . . . . 12 2.6 Damage profiles for the one-dimensional homogeneous tension test. . . . . 13 2.7 Damage field for the one-dimensional heterogeneous tension test, GcRatio = 0.99. A 1% reduction in resistance on the right side is enough to pivot the fracture toward that side. . . . . . . . . . . . . . . . . . . . . . . . . . . . 13 2.8 Damage field for the one-dimensional heterogeneous tension test, GcRatio = 0.5. A 50% reduction in resistance on the right side allows for the other side to suffer considerably less damage. . . . . . . . . . . . . . . . . . . . . 14 vii 2.9 Damage profiles for the one-dimensional heterogeneous tension test, GcRatio = 0.99 ....................................... 14 2.10 Damage profiles for the one-dimensional heterogeneous tension test, GcRatio = 0.5......................................... 14 3.1 Geometry and boundary conditions for the two-dimensional tension test. . 16 3.2 Two-dimensional tension test. Damage fields for adaptive 24 ×24 (a) nonadaptive 240 ×240 (b), with a refinement factor m= 10 and a GcRatio = 10. ........................................ 17 3.3 The corresponding displacement - force plots for both cases in Figure 3.2. There is a slight difference but overall it behaves the same. Tolerance = 10−4........................................ 18 3.4 Two-dimensional tension test. Force-displacement plots merged for the uniformly refined 240 ×240 and the adaptive 24 ×24 tests in Figure 3.1. a) Tolerance = 10−2b) Tolerance = 10−3c) Tolerance = 10−4d) Tolerance = 10−8...................................... 18 4.1 Geometry and boundary conditions for the two-dimensional adaptive and homogeneous reference tests. Left: Tension test. Center: Mixed-mode test. Right:Sheartest................................. 19 4.2 Damage fields for two-dimensional adaptive and homogeneous reference tests. Left: Tension test. Center: Mixed-mode test. Right: Shear test. . . . 20 4.3 Geometry and boundary conditions. a) Two-dimensional tension test with two small cuts at mid-height and a hard circular inclusion. b) Twodimensional shear test with a cut at mid-height that goes to the middle andahardbottomlayer............................. 21 4.4 Two-dimensional tension test with a hard centered circular inclusion, GcRatio = 3.5. a) Damage field showing a slight penetration and a complete deflection. b) Force-displacement plot. . . . . . . . . . . . . . . . . . . . . . . . 22 viii 4.5 Damage fields and force-displacement plots for the tension test with a hard circular inclusion and a for a GcRatio = 3.5. Increments right before and after the sudden drops occur are shown. Note that the damage scale in panel a) does not range from 0 to 1. . . . . . . . . . . . . . . . . . . . . . . 23 4.6 Damage fields and force-displacement plots for the tension test with a hard circular inclusion and for a GcRatio = 2. Increments right before and after the sudden drops occur are shown. . . . . . . . . . . . . . . . . . . . . . . 24 4.7 Example of an unstructured mesh with a circle in the middle generated with ez4usoftware. .............................. 25 4.8 Damage fields and force-displacement plots for the tension test with a hard circular inclusion and for a GcRatio = 3.25. Increments right before and after the last sudden drop occurs are shown. . . . . . . . . . . . . . . . . . 26 4.9 Damage fields for the shear test with a hard bottom layer and for a) GcRatio = 2 and b) GcRatio =5........................ 27 4.10 Force-displacement plots for the shear test test with a hard bottom layer with the following Gcratios: homogeneous, GcRatio = 2 and a GcRatio = 5. 28 4.11 Damage fields for the shear test with a hard bottom layer and for a) GcRatio = 6 b) GcRatio = 7 c) GcRatio =8.................. 29 4.12 Force-displacement plots for the shear test with a hard bottom layer and for a GcRatio = 6, a GcRatio = 7 and a GcRatio =8............. 29 4.13 Damage fields and force-displacement plots for the shear test with a hard bottom layer and for GcRatio = 6. Increments right before and after the transition between deflection and penetration and right before and after the crack nucleation are shown. . . . . . . . . . . . . . . . . . . . . . . . . 31 4.14 Geometry and boundary conditions. Two-dimensional tension test with a small cut at mid-height and a more resistant interface with an angle of 45◦ with the vertical direction. The length of penetration, length of deflection, and ˆ Laredefined................................. 32 ix 2Introduction field approach, this model provides a comprehensive framework for simulating fracture phenomena in materials with varying material properties. Accounting for heterogeneityinduced variations in fracture toughness, the model offers a realistic representation of fracture behavior in real-world materials. Heterogeneity is a common trait in many materials of engineering relevance, whether natural or manufactured. Perfect homogeneity does not exist, so, to some degree, all materials are heterogeneous. The adaptive phase-field model is a strong computational technique for studying fracture in heterogeneous materials. As a result of its ability to iteratively update the computational mesh and account for the effects of material heterogeneities, it presents a viable technique for comprehending and predicting fracture behavior in real-world engineering applications. By expanding the understanding of fracture mechanics in heterogeneous materials, this model has the potential to stimulate innovation in the design and optimization of robust and resilient structures. 1.2 State of the art In this section, we give an overview of the different approaches when dealing with phasefield models for heterogeneous materials. 1.2.1 Deterministic heterogeneity Continuously varying elastic and/or fracture properties Vicentini et al. (2023) investigate phase-field modeling of brittle fracture on heterogeneous bars. Heterogeneity is characterized by the continuous variation of fracture properties along the axis of the bar. They show how the effect of heterogeneity is strictly tied to the non-local nature of the phase-field regularization. The interface approach There are some cases in which heterogeneous problems are approached by introducing an interface between the two regions with different properties within the material. Yuan et al. (2023) propose a smeared interface that allows for the critical energy release rate of the interface to be regularized and to predict complex cracking phenomena such as multiple 1.2. State of the art 3 crack initiation, merging, and branching. Unnikrishna Pillai et al. (2023) propose, as well, an approach for modeling brittle fracture in heterogeneous composites using a combined diffused material interface method. The diffused material interface avoids the difficulty due to the discontinuity in the stress at the interface. Fracture in bones Hug et al. (2022) perform a practical approach to humerus fractures using phase-field models. It considers the bone to be an elastic, isotropic material with heterogeneous material properties. QCT scans can obtain the ash density of bones; the Young’s modulus is a function of the ash density of bones. Heterogeneity takes place in the form of a spatially varying critical energy release rate, and the model is validated in three in vitro experiments. In this thesis, we approached heterogeneity with the same approach. They successfully demonstrate that the phase-field model is a suitable tool to predict fractures in bones. Fast Fourier transform-based numerical method Ma et al. (2023) focus on the implementation of the fast Fourier transform-based numerical method within the phase-field model. Results show that it has no effect on homogeneous materials but significantly influences heterogeneous materials where fracture properties vary along the material. The method considers one of the heterogeneous materials to be unbreakable and completely suppresses the diffusion of damage toward it. In heterogeneous materials, the local fracture patterns and the global stress-strain relations differ. 1.2.2 Probabilistic heterogeneity There are different views on probabilistic heterogeneity. Noii et al. (2022) perform a stochastic analysis over a phase-field model with random materials, geometric properties, and the distribution of stiff inclusions and voids. The Monte Carlo finite element method is employed to solve the proposed stochastic partial differential equation (PDE). This allows them to study the highly heterogeneous microstructure and monitor the uncertainty in fracture mechanics. Su et al. (2023) study the statistical effects of pore features on the mechanical properties and fracture behavior of heterogeneous random porous materials using phase-field modeling. Random porous materials using Gaussian random fields with stochastic pore size and 4Introduction shape are used to analyze the statistical attributes of pore sizes and shapes. Materials with such properties include silver nanoparticles, concrete materials, and natural bones, among others. 1.3 Goals and layout of this thesis The main goal of this thesis is to study crack propagation in heterogeneous materials. To do so, these specific goals have been determined: 1. Acquire familiarity with the existing adaptive phase-field model and Matlab code for homogeneous materials and extend it to heterogeneous materials. Understanding the basics of the finite element method, the phase-field models of fracture, and the computational fracture model developed in Muix´ı (2020). Upgrading the code for heterogeneous materials. 2. Design and carry out computational tests of fracture in heterogeneous materials. Explore the mechanical response (crack arrest, crack penetration, and crack deflection). Assess the capabilities and limitations of the adaptive phase-field modeling of fracture for heterogeneous materials. Simulations are performed in Matlab, and large computations are carried out in LaC`aN’s computer cluster. Chapter 2 Brief overview of the model 2.1 Hybrid phase-field model To study crack propagation in heterogeneous materials, we use the computational phasefield model developed by Alba Muix´ı in her thesis, Muix´ı (2020). Phase-field models are widely used in computational fracture mechanics. They are a continuous and diffuse model for brittle and quasi-brittle materials; cracks propagate through the less energy-demanding path. The phase-field damage variable dquantifies the damage to the material. Its value goes from zero for intact material to one for fully broken material. The characteristic parameter of the phase-field model for fracture is the length-scale parameter l. The width of the diffuse crack depends on l. To reproduce sharp cracks, lvalues need to be small, and fine mesh is required along the crack. The total energy of the body is expressed as: E(u,Γc) = ZΩ Ψ0(ε(u))dV+GcZΓc dΓ (2.1) The first term in Eq. (2.1) is the stored elastic strain energy, and the second term is the fracture energy. The crack is understood to follow the less energy-demanding path at every loading step; it propagates according to the energy minimization principle, Francfort and Marigo (1998). 5 6Brief overview of the model      ∇ · σ(ε(u), d) = 0 d−ℓ2∆d=2ℓ Gc (1 −d)Ψ0(ε(u)) (2.2) H(x, t) = max 0≤τ≤tΨ0(ε(x, τ)) (2.3) This minimization of the regularized total energy leads to a coupled system of two nonlinear partial differential equations, Eqs. (2.2), Bourdin et al. (2000). The system of PDEs is solved incrementally in a staggered manner and iteratively for each increment. To ensure irreversible damage, a history variable His used in the damage equation, see (2.3) Miehe et al. (2010). Ψ0(ε) = Ψ0(−ε) (2.4) Moreover, the phase-field model up to this point provides the same response in tension as it does in compression because the sign of the strains is lost in the strain energy density, see (2.4). There are different ways of splitting the strain energy, but some of them lead to non-linear sub-problems to solve. An easier way to split the strain energy is to do so as a post-process; sub-problems remain linear.      ∇ · σ(ε(u), d) = 0 d−ℓ2∆d=2ℓ Gc (1 −d)H+(2.5) With:          σ(ε(u), d) = g(d)∂Ψ0 ∂ε g(d) = ((1 −d)2where Ψ+ 0>Ψ− 0 1 elsewhere (2.6) The final phase-field model we end up working with is the hybrid phase-field model proposed by Ambati et al. (2015), see (2.5) and (2.6). 2.2. Adaptive refinement based on Nitsche’s method 7 2.2 Adaptive refinement based on Nitsche’s method The main inconvenience of solving this model is its computational cost. Very fine meshes are required to properly model the crack behavior. The proposed computational phasefield model by Muix´ı (2020) introduces different adaptive techniques to reduce the computational cost of the simulations. The one we use in this thesis is adaptive refinement based on Nitsche’s method, see Muix´ı et al. (2020). This technique distinguishes between standard elements, belonging to the original and background mesh, and refined elements. Standard elements become refined elements by two means:  Elements close to the weak part of the material are prescribed to be refined before the test takes place. In this thesis, cuts are introduced to create controlled starting points for damage. The elements surrounding the notch are prescribed to be refined.  The adaptive strategy automatically determines whether an element is to be refined or not according to a given damage threshold. The damage threshold for automatic refinement is given by the refinement value d∗. Elements with a damage dlarger than the d∗are refined. The precision of the refinement is given by the refinement factor m. A standard element that is prescribed to be refined or whose damage exceeds the d∗, will be refined into m×melements, see Figure 2.2. For instance, if mis equal to 10, a single standard element will be refined into 100 elements. In Muix´ı (2020), it is shown that fine elements are only needed along the crack and on elements that suffer the most damage dat some increment. In Figure 2.1, it is observed that using a background mesh of dimensions 24 ×24 and an mequal to 10, the damage contour converges with the one belonging to a background mesh of 240 ×240; the same occurs regarding the force-displacement plot. For a proper choice in spatial resolution, the adaptive technique has no significant effect on the final results and allows for efficient fracture simulations. 8Brief overview of the model (a) (b) Figure 2.1: Figures belonging to a shear test with d∗= 0.2 in Muix´ı (2020) pages 61 — 62. a) Contour plots of the damage and b) Force-displacement plots for different mesh resolutions and refinement factors. 2.2. Adaptive refinement based on Nitsche’s method 9 (a) (b) Figure 2.2: Two-dimensional tension test with a cut at mid-height that goes to the middle with the four elements surrounding the crack tip prescribed to be refined for two given increments. a) The deformed mesh including the standard and refined elements and a zoom. b) The deformed mesh including the standard and refined elements for a further increment and a zoom of only the refined elements. The two next elements have been automatically refined. The displacements have been amplified by a factor of 10. 10 Brief overview of the model 2.3. One-dimensional tension test 11 2.3 One-dimensional tension test The aim of this chapter is to present a simple case to introduce the phase-field model for heterogeneous materials. We define heterogeneity as the critical energy release rate Gcparameter relation between regions of a material. In a material with a number iof regions with different fracture toughnesses, we define the GcRatio as shown in equation (2.7) . GcRatio,i =Gc,i Gc,o (2.7) Material, adaptive and computational parameters Meaning Symbol Value Units Young’s modulus E210 GPa Poission ratio ν0 Length-scale parameter l0.015 mm Critical energy release rate Gc,o 2.7×10−3kN/mm Refinement value d∗0.2 Refinement factor m10 Increment of displacement △uD10−4mm One-dimensional test dimensions (1 ×L) 1 ×10 mm ×mm Two-dimensional test dimensions (L×L) 10 ×10 mm ×mm Table 2.1: Material and adaptive parameters used throughout the thesis. The proposed numerical experiment has the boundary conditions and geometric properties shown in Figure 2.3. These conditions allow for a uniaxial tension test. The parameters used to perform this simulation and the upcoming ones in this thesis are shown in Table 2.1 (in the case that one parameter is modified in a test, it is specified). The mesh dimensions are 1 ×10 and quadrilateral elements are used, see Figure 2.4. 18 Introducing the adaptive refinement Figure 3.3: The corresponding displacement - force plots for both cases in Figure 3.2. There is a slight difference but overall it behaves the same. Tolerance = 10−4. (a) (b) (c) (d) Figure 3.4: Two-dimensional tension test. Force-displacement plots merged for the uniformly refined 240 ×240 and the adaptive 24 ×24 tests in Figure 3.1. a) Tolerance = 10−2b) Tolerance = 10−3c) Tolerance = 10−4d) Tolerance = 10−8 Chapter 4 Two-dimensional experiments In this chapter, a series of numerical experiments are carried out to study the mechanical response of fracture in heterogeneous materials and to assess the capabilities and limitations of the adaptive phase-field model. To compare the upcoming two-dimensional heterogeneous numerical experiments with the homogeneous fractures’ behavior, we propose the tests in Figure 4.1. The specimen has a cut at mid-height that goes to the middle of the material, and the four elements surrounding the crack tip are prescribed to be refined. Cuts are introduced to weaken a zone and guide the initialization of the crack. The properties of the test are shown in Table 2.1. The test is performed for three different boundary conditions: tension, mixed-mode, and shear. Figure 4.1: Geometry and boundary conditions for the two-dimensional adaptive and homogeneous reference tests. Left: Tension test. Center: Mixed-mode test. Right: Shear test. 19 20 Two-dimensional experiments The resulting damage fields are shown in Figure 4.2. The natural fracture paths for a tension test and a shear test are shown. We will recall them during the thesis to observe how heterogeneities modify these natural paths. Furthermore, we can observe that by mixing the boundary conditions, we can produce a spectrum of crack paths between the tension test and the shear test. Figure 4.2: Damage fields for two-dimensional adaptive and homogeneous reference tests. Left: Tension test. Center: Mixed-mode test. Right: Shear test. 4.1. Arrest, deflection and penetration of the crack 21 4.1 Arrest, deflection and penetration of the crack The aim of these tests is to evaluate fracture behavior when encountering heterogeneity during a test. As mentioned before, we define heterogeneity as the ratio of fracture toughness Gcbetween different regions in a material. To do so, the following tests, with parameters from Table 2.1, are presented, see Figure 4.3. (a) (b) Figure 4.3: Geometry and boundary conditions. a) Two-dimensional tension test with two small cuts at mid-height and a hard circular inclusion. b) Two-dimensional shear test with a cut at mid-height that goes to the middle and a hard bottom layer. Tension test with a hard circular inclusion Test 4.3(a), is performed for three different Gcratios: GcRatio = 3.5, GcRatio = 3.25 and GcRatio = 2. The reason for choosing these three ratios is that they showcase penetration, deflection, and in-between fracture behavior. Performing the test for a GcRatio = 3.5, we can observe that there is a slight penetration of the crack inside the inclusion, but a deflection ends up occurring, see Figure 4.4(a). At some increment, the penetration of the inclusion is the easiest path energy-wise, but as the test carries on, circling the obstacle becomes less energy-demanding. In the force-displacement plot, we can observe that there are two sudden drops in force, which are expected; the phase-field model we use is suitable for brittle materials (no ductility), see Figure 4.4(b). To understand why these drops and rises occur on the forcedisplacement plot, we must have a look at Figure 4.5. The first rise corresponds to the 22 Two-dimensional experiments cracks reaching the inclusion. Once it does, the force drops and starts loading again. Elastic energy accumulates during this temporary crack arrest. During this process, the crack manages to push inside the inclusion, but damage also starts to concentrate at the bottom of it. The last drop corresponds to the crack circling the inclusion. (a) (b) Figure 4.4: Two-dimensional tension test with a hard centered circular inclusion, GcRatio = 3.5. a) Damage field showing a slight penetration and a complete deflection. b) Forcedisplacement plot. 4.1. Arrest, deflection and penetration of the crack 23 (a) (b) (c) (d) (e) (f) (g) (h) Figure 4.5: Damage fields and force-displacement plots for the tension test with a hard circular inclusion and a for a GcRatio = 3.5. Increments right before and after the sudden drops occur are shown. Note that the damage scale in panel a) does not range from 0 to 1. 24 Two-dimensional experiments Test 4.3(a) is now performed for a GcRatio = 2. Two damage fields and force-displacement plots of interest are shown in Figure 4.6. For this ratio in the fracture toughness, the cracks completely penetrate the inclusion. The inclusion being twice as resistant as the rest of the material is not enough for the crack to deflect it. The force-displacement figures show initially the same behavior as in the previous case, see Figure 4.5, the initial rise and drop correspond to the accumulation of elastic energy and the crack reaching the inclusion. From the consecutive increment, force starts loading again, in a less steep and shorter rise (the inclusion is less resistant), to completely penetrate the hard circular inclusion. (a) (b) (c) (d) Figure 4.6: Damage fields and force-displacement plots for the tension test with a hard circular inclusion and for a GcRatio = 2. Increments right before and after the sudden drops occur are shown. 4.1. Arrest, deflection and penetration of the crack 25 The reason for the crack not being completely straight is due to the background being a non-structured mesh. To introduce a circular inclusion in this test ez4u, external software provided by LaC`aN, has been used; see Figure 4.7. Some upcoming experiments’ background meshes have also been generated using this software. (a) Figure 4.7: Example of an unstructured mesh with a circle in the middle generated with ez4usoftware. Test 4.3(a) is finally performed for a ratio in fracture’s toughness that lies in between the two previous ones, GcRatio = 3.25. Two damage fields and force-displacement plots of interest are shown in Figure 4.8. It is observed how the cracks manage to penetrate inside the hard inclusion even more than in case 4.4(a). A crack arrest occurs when both cracks are surrounded by a more resistant region and are unable to keep going forward or change their path. The second drop in force corresponds to the formation of a new crack at the bottom of the circular inclusion. To finish the test, there is a gentle drop that corresponds to the bottom crack connecting to the initial ”stuck” cracks. 26 Two-dimensional experiments (a) (b) (c) (d) Figure 4.8: Damage fields and force-displacement plots for the tension test with a hard circular inclusion and for a GcRatio = 3.25. Increments right before and after the last sudden drop occurs are shown. We can observe how, for this phase-field model 2.1, cracks will go through the least energydemanding path available at every increment. Fracture may penetrate, deflect, or arrest in front of a more resistant obstacle, depending on the GcRatio and the geometry of the heterogeneity. 4.1. Arrest, deflection and penetration of the crack 27 Shear test with a hard bottom layer To continue assessing fracture’s behavior when reaching heterogeneity, the shear test 4.3(b), is performed for several ratios of fracture toughness to view the possible modes of fracture. As a reference, the homogeneous damage field for this test is shown on the right in Figure 4.2. The resulting damage fields for a GcRatio = 2 and a GcRatio = 5 are shown in Figure 4.9. A complete penetration occurs in both cases, the difference being the total length of deflection; the harder the layer is, the larger the deflection is. Increasing the GcRatio, increases the energy required for the penetration of the layer. The crack deflects until, at some increment, penetration becomes cheaper, energy-wise, than continuing with the deflection. In Figure 4.10, the force-displacement plots for the two heterogeneous cases and the homogeneous case have been merged for comparison. Similarly to what we saw in the hard circular inclusion test, when the crack meets the inclusion, a drop in force occurs. Then, for the heterogeneous tests, there is an accumulation of elastic energy that corresponds to the deflections. For the homogeneous case, once the peak force is reached, there is enough elastic energy to completely penetrate the layer. A non-smooth, spiky curve corresponding to the complete penetration of the layer follows, and, finally, a sudden drop indicates the end of the test. We can observer that the total displacement and the peak force increase as we increase the resistance of the bottom inclusion. (a) (b) Figure 4.9: Damage fields for the shear test with a hard bottom layer and for a) GcRatio = 2 and b) GcRatio = 5. In the range of fracture toughnesses used in this shear test, it seems that increasing the 34 Two-dimensional experiments 4.3 Crack behavior over thin paths 4.3.1 Triangular wave of less resistance To assess the fracture’s capability to propagate along thin paths of lower resistance, the two following tests are proposed, see Figure 4.16 . The less resistant region is defined by a triangular wave with an amplitude Aand a wave length λ. Test 4.16(a) has the following properties: A=L/24 and λ=L/6. Test 4.16(b) has the same λbut double the A. The triangular wave has a thickness of L/24 in order to define the less resistant zone. (a) (b) Figure 4.16: Geometry and boundary conditions. Two-dimensional tension tension test with a small cut at mid-height and a less resistant triangular wave path. The path has a thickness = L/24. a) A=L/24 and λ=L/6 b) A=L/6 and λ=L/12. The damage fields for the two tests and for GcRatio = 0.8 (first row), GcRatio = 0.5 (second row), and GcRatio = 0.1 (last row) are shown in Figure 4.17. In the first row, Figure 4.17(a) and Figure 4.17(b), we observe that, in both tests, the crack barely sees its natural path altered. In Figure 4.17(c), we can observe that the crack imprecisely follows the proposed path, while in Figure 4.17(d), the reduction in fracture toughness for the triangular wave path is not enough and the crack is barely guided. Finally, for a 90% reduction in the Gc, in the last row, Figure 4.17(e) and Figure 4.17(f), both cracks precisely follow the less-resistant proposed paths. 4.3. Crack behavior over thin paths 35 (a) (b) (c) (d) (e) (f) Figure 4.17: Damage fields for the two-dimensional tension tests with a less resistant triangular wave path proposed for a a) and b) GcRatio = 0.8, c) and d) GcRatio = 0.5 and e) and f) GcRatio = 0.1. 36 Two-dimensional experiments 4.3.2 Void inclusions to guide the crack In this section, we test the fracture’s capability of advancing along thin paths of voids, which can be considered parts of a material with a fracture toughness Gc= 0. Small voids are commonly incorporated into industrial product packaging to facilitate easy unpacking by providing pre-cut lines or perforations in the material (such as cardboard or plastic). These pre-cut lines serve as a guide for controlled tearing, allowing for convenient access to the contents without damaging the packaging. The first experiment we perform aims to resemble a triangular wave, see Figure 4.18. In Figure 4.18(a) a steep wave is proposed, and the test is carried out guided by the bottom row of voids but unsuccessfully following the desired path. In Figure 4.18(b), another triangular wave path is proposed, but this time the voids are vertically closer to each other, allowing for a less steeper wave. We can observe that the wave successfully adheres to the path. (a) (b) Figure 4.18: Damage fields for the two-dimensional tension tests performed with the aim to drive the crack through a triangular wave path using voids. The objective for the second experiment is to guide the crack upwards for a tension test whose natural path is to advance straight, see Figure 4.2. In Figure 4.19(a), larger voids placed very close aim to drive the fracture completely upwards. The test is carried out as a normal tension test would. We have not managed to force the crack to change its natural path so directly. In Figure 4.19(b), the crack is successfully guided. Even though the inclusions are separated a bit more than for case 4.19(a), the proposed path is more gentle (in the sense that the change in the direction of fracture does not change as much 4.3. Crack behavior over thin paths 37 as for the previous case with respect to its natural tension test direction of fracture), and this allows for the crack to stick to the proposed path. (a) (b) Figure 4.19: Damage fields for the two-dimensional tension tests performed with the aim to drive the crack upwards using voids. 38 Two-dimensional experiments 4.4 Improving the resistance by placing inclusions In this section, the effect that adding hard inclusions to a specimen has on the peak force and the total displacement during a given test is assessed. To do so, the following tension tests are proposed, see Panel 4.20. To minimize any numerical error and any provoked unsymmetrical behavior of the crack, the inclusions go from the top to the bottom, we leave the same space on both sides of the inclusions (between the notch and the first inclusion and between the last inclusion and the boundary of the material), there is the same space between inclusions (tests 4.20(e) and 4.20(g)), and we perform a tension test whose behavior in this case is to advance straight. We define the area occupied by the more resistant inclusions as the InclusionRatio in these tests: InclusionRatio =InclusionArea TotalArea =96 576 (4.1) The tests in Panel 4.20 are performed for a GcRatio = 4. The resulting damage fields are shown in the same Panel 4.20. We can observe that the fracture’s path remains almost the same, there is complete penetration of the inclusions, and the direction of the crack remains unaltered. The only difference between the homogeneous case and the heterogeneous cases is a very slight thickening of the damaged zone at the start of every inclusion due to temporary crack arrests. These occur when the fracture accumulates enough elastic energy to penetrate the inclusions. 4.4. Improving the resistance by placing inclusions 39 (a) (b) (c) (d) (e) (f) (g) (h) Figure 4.20: Geometry and boundary conditions and damage fields. Two-dimensional tension tests with a small cut at mid-height, a) is homogeneous, b), c) and d) are heterogeneous with a GcRatio = 4. The more resistant constant area is split between the inclusions. 40 Two-dimensional experiments In Figure 4.21 we can observe the force-displacement plots merged for all of the cases in Panel 4.20. The first peak is almost the same in all cases. At that increment, there is enough elastic energy for the crack to begin and, for the homogeneous test, to complete the test. Comparing only the heterogeneous tests, the highest peak force and the lowest total displacement are observed for the test in which the inclusion area is evenly distributed among four inclusions. The lowest peak force and the highest total displacement can be observed in the test in which all of the hard area is concentrated in one inclusion. In Figure 4.21, we can also observe that, for all of the heterogeneous tests, the number of sudden drops after the first initial drop corresponding to the crack reaching the first inclusion is the same as the number of inclusions. By breaking down the test with four obstacles, Figure 4.20(h), we are able to link the number of sudden drops after the initial one to the number of inclusions the crack penetrates. In Figures 4.22 and 4.23, snapshots of the damage field and the force-displacement plot of the increments right after and before these sudden drops occur are shown. Figure 4.21: Force-displacement plots for the two-dimensional tension tests for a GcRatio = 4 with more resistant constant area split into the number of inclusions. 4.4. Improving the resistance by placing inclusions 41 (a) (b) (c) (d) (e) (f) (g) (h) Figure 4.22: Damage fields and force-displacement plots for the increments right before and after sudden drops occur on the the two-dimensional heterogeneous test with four inclusions with a GcRatio = 4. 42 Two-dimensional experiments (a) (b) (c) (d) (e) (f) (g) (h) Figure 4.23: Damage fields and force-displacement plots for the increments right before and after sudden drops occur on the the two-dimensional heterogeneous test with four inclusions with a GcRatio = 4. 4.4. Improving the resistance by placing inclusions 43 For a GcRatio = 4, we observe temporary crack arrests that do not alter the regular tension test path. Performing test 4.20(c) with more resistant inclusions, GcRatio = 20, results in an undesired crack behavior: at a given increment, the accumulated elastic energy is not large enough for the crack to penetrate the inclusion, and the crack deflects the inclusion splitting into two; it bifurcates. Instead, crack nucleation occurs at the top of the material on the other side. The new crack propagates straight until it meets the inclusion, the crack then deflects the inclusion downwards for almost the entirety of it and, finally, an increment is reached in which the accumulated elastic energy is sufficiently large for the two cracks to connect across the inclusion, see Figure 4.24. Figure 4.24: Damage field for the two-dimensional heterogeneous tension test with small cut at mid-height with one inclusion with a GcRatio = 20. To control the nucleation of new cracks, we introduce another cut at the same height but on the other side. We redefine the problem statements so that we keep the same disposition used previously to minimize any numerical error and to avoid any provoked unsymmetrical behavior of the crack, see Figure 4.25. 50 Two-dimensional experiments Chapter 5 Concluding remarks and future work 5.1 Concluding remarks This dissertation explores the mechanical response, capabilities, and limitations of adaptive phase-field modeling of fracture in heterogeneous materials. Throughout the work, heterogeneity is defined as the ratio in the fracture tougnhess Gcbetween different parts of the material. To do so, an existing code developed for homogeneous materials by Alba Muix´ı in Muix´ı et al. (2020) is extended for heterogeneous materials, and computational tests of fracture in heterogeneous materials are carried out. 1. Qualitative analysis of fracture’s behavior when encountering heterogeneity. In Section 4.1 two different cases are presented that allow us to observe a crack arrest, deflection, and penetration for different fracture toughnesses. The different phenomena are also related to the force-displacement plots, and we determine that when a fracture meets a more resistant inclusion, there is a drop in force. We also determine that, in general, by increasing the GcRatio (i.e., making the inclusion more resistant), the peak force and the total displacement for the complete fracture to occur increase, but it may occur that the mode of fracture changes if the GcRatio keeps being increased, and that does not translate into the general behavior mentioned. Several tests must be carried out to determine the fracture’s path and the force-displacement plot. 2. Length of penetration vs length of deflection. In Section 4.2, we have deter51 52 Concluding remarks and future work mined that one of the factors that make the fracture penetrate or deflect a more resistant inclusion is the length of penetration and the length of deflection. For a given tension test, an obstacle with a given inclination is presented. Altering the length of penetration and deflection modifies the fracture’s path for a given angle of contact between the crack and a more resistant obstacle. 3. Guiding the crack. In Section 4.3 we are able to guide the cracks by offering a less resistant path, lowering the Gcof a part of the material, and strategically placing voids (which can be considered parts of the material with a Gc = 0). For a sufficiently low GcRatio or sufficiently gentle pattern of voids, the crack precisely follows very thin paths. 4. Improving the resistance by placing inclusions. In Section 4.4 we carry out a series of tension tests with a constant area of the material covered by a higher Gc parameter. We can state that to maximize the total displacement, this area must be concentrated in a single inclusion. If the objective is to maximize the peak force, it is better to spread the more resistant area into different inclusions. 5. Capability of the model for the crack to go backward. In Section 4.5 we carry out a series of tension tests with the goal of making the crack go backward following the minimization of energy principle. We carry out a test with a certain geometry regarding the Gcparameter distribution along the elements of the mesh, and we subject the material to a very extreme difference in the Gcparameter. The results show that, in some way, the crack seems to make it backward, but the damage field shows a chaotic path of the crack. Making the crack go backward in order to advance is unnatural, and we cannot conclude that the model comfortably allows the crack to do so. 5.2 Future work 1. Probabilistic inclusions. We have worked only with deterministic inclusions, but, as seen in the state of the art, heterogeneity can also be introduced using probabilistic methods to take into account the inherent variability and uncertainty related to material heterogeneities in real life. Creating procedures for producing stochastic representations of heterogeneous microstructures while taking into account statisti- 5.2. Future work 53 cal distributions of material characteristics — in this thesis, fracture toughness — is a possible line of work. 2. Extension to three-dimensional and experimental validation. All of the tests that are performed in this thesis are in two-dimensional space. The extension of this code to perform three-dimensional tests would allow meshing many reallife bodies and their corresponding experimental validations. This extension would need to be done very carefully, focusing on optimization, and would probably require high-performance computing due to the high computational cost of the simulations. 3. Automation and optimization. The extension of the code, especially regarding unstructured meshes, includes some ”hard-coded” functions that are not completely robust (i.e., when transferred to another mesh, they need some previous work). Improving their adaptability and reliability could be done, as could a general optimization of the code. 54 Concluding remarks and future work Bibliography Ambati, M., T. Gerasimov, and L. De Lorenzis (2015). A review on phase-field models of brittle fracture and a new fast hybrid formulation. Computational Mechanics 55(2), 383 — 405. Bourdin, B., G. Francfort, and J.-J. Marigo (2000). Numerical experiments in revisited brittle fracture. Journal of the Mechanics and Physics of Solids 48(4), 797–826. Djeumen, E., G. Moln´ar, N. Tardif, M. Coret, J. Desquines, T. Taurines, and M.-C. Baietto (2022). Modeling diffusive phase transformation and fracture in viscoplastic materials. International Journal of Solids and Structures 252, 111757. Francfort, G. and J.-J. Marigo (1998). Revisiting brittle fracture as an energy minimization problem. Journal of the Mechanics and Physics of Solids 46(8), 1319–1342. Hug, L., G. Dahan, S. Kollmannsberger, E. Rank, and Z. Yosibash (2022). Predicting fracture in the proximal humerus using phase field models. Journal of the Mechanical Behavior of Biomedical Materials 134, 105415. Kumar, P., P. Steinmann, and J. Mergheim (2022). A graded interphase enhanced phasefield approach for modeling fracture in polymer composites. Forces in Mechanics 9, 100135. Ma, X., Y. Chen, M. Shakoor, D. Vasiukov, S. V. Lomov, and C. H. Park (2023). Simplified and complete phase-field fracture formulations for heterogeneous materials and their solution using a fast fourier transform based numerical method. Engineering Fracture Mechanics 279, 109049. Miehe, C., M. Hofacker, and F. Welschinger (2010). A phase field model for rateindependent crack propagation: Robust algorithmic implementation based on operator splits. Computer Methods in Applied Mechanics and Engineering 199(45), 2765–2778. 55 56 BIBLIOGRAPHY Muix´ı, A. (2020). Locally adaptive phase-field models and transition to fracture. Ph. D. thesis, Universitat Polit`ecnica de Catalunya. Muix´ı, A., S. Fern´andez-M´endez, and A. Rodr´ıguez-Ferran (2020). Adaptive refinement for phase-field models of brittle fracture based on nitsche’s method. Computational Mechanics 66(1), 69 — 85. Noii, N., A. Khodadadian, and F. Aldakheel (2022). Probabilistic failure mechanisms via monte carlo simulations of complex microstructures. Computer Methods in Applied Mechanics and Engineering 399, 115358. Su, Y., J. Zhu, X. Long, L. Zhao, C. Chen, and C. Liu (2023). Statistical effects of pore features on mechanical properties and fracture behaviors of heterogeneous random porous materials by phase-field modeling. International Journal of Solids and Structures 264, 112098. Unnikrishna Pillai, A., A. K. Behera, and M. M. Rahaman (2023). Combined diffused material interface and hybrid phase-field model for brittle fracture in heterogeneous composites. Engineering Fracture Mechanics 277, 108957. Vicentini, F., P. Carrara, and L. De Lorenzis (2023). Phase-field modeling of brittle fracture in heterogeneous bars. European Journal of Mechanics - A/Solids 97, 104826. Yuan, J., L. Wang, and C. Chen (2023). Interfacial fracture analysis for heterogeneous materials based on phase field model. Computational Materials Science 220, 112066.