Effect of tailored fiber deposition in 3D printed composites: application of an anisotropic phase field model

Theoretical and Applied Fracture Mechanics 127 (2023) 104030 Available online 3 August 2023 0167-8442/© 2023 The Author(s). Published by Elsevier Ltd. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/). Effect of tailored fiber deposition in 3D printed composites: application of an anisotropic phase field model Simone Sangaletti a , * , Anatoli Mitrou b , Israel G. García a , Albertino Arteiro b a Departamento de Mec´ anica de Medios Continuos y Teoría de Estructuras, E.T.S. Ingeniería, Universidad de Sevilla, Camino de los Descubrimientos s/n,41092, Seville, Spain b DEMec, Faculdade de Engenharia, Universidade do Porto, Rua Dr. Roberto Frias, s/n 4200-465 Porto, Portugal ARTICLE INFO Keywords: Fracture Computational mechanics Anisotropy 3D printing ABSTRACT Continuous Fiber 3D printing is a relatively new technology which can allow for tailored reinforcement of critical regions in structural components, i.e., stress concentrations, following principal stress lines. The influence the fiber deposition path has on the mechanical and failure behavior of such components is assessed using an anisotropic phase field model. A comparison with experimental results for notched unidirectional composite plates, available in the literature, demonstrates the ability of the method to produce satisfactory predictions for unidirectional reinforcement paths. The analysis is then extended to Open-Hole and Double Edge-Notched tension coupons of both unidirectional and variable stiffness reinforcement patterns. It is observed that the strength obtained for the components made with a reinforcement pattern that follows the principal stress lines is markedly higher than that for the equivalent unidirectionally reinforced ones. It is highlighted that the improvement in strength deriving from the tailored fiber deposition cannot be deduced solely by the analysis of the stress concentration factor but an analysis taking damage into account is necessary. In addition, the effect of the reinforcement strategy on the size effect was also explored, highlighting how the tailored fiber path leads to an increase in the failure load attainable by the specimens for all the dimensions analyzed. 1. Introduction Carbon Fiber (CF) reinforced laminates have for years been the object of interest among academia and different industrial sectors (e.g., aerospace, automotive) due to the high performance such materials can guarantee, with special interest being devoted to improving manufacturing methods and techniques used for their production. For many years, classical lamination processes via hand lay-up were the only way to produce composite laminates. In this procedure, preimpregnated layers (prepregs) of unidirectional (UD) fibers are sequentially stacked to form the final plate or component. This technique is highly dependent on the ability of the operator who oversees the lamination procedure. In the 60 ′ s, Automated Fiber Placement (AFP) was developed, allowing the automatic deposition of prepregs and also fibers (concerning tow-steering techniques), in unique and predefined paths to generate multidirectional laminates. These laminates can be traditionally cured or cured in-situ as the plies are placed. This technique is, until now, the main one used to produce components of large dimensions. With advanced robotic systems, AFP became a means to produce complex shapes and components much more efficiently, helping eliminate wasted material, and minimize the possibility of human error, allowing for more overall control of the process. However, this too comes with its weaknesses. Prepreg types (e.g., width, material system, thickness) could be restrictive, attainable shapes are still somewhat limited, and cost remains high for increase complexity. In recent years, the additive manufacturing technologies that were developed act as gamechangers in the field. One of the most common forms seen is fused deposition modeling (FDM) which is commonly referred to with the term 3D printing. 1 3D printing has evolved to now allow the deposition of continuous fiber bundles alongside a polymeric matrix to generate fibrous composite components. This paves the way for significantly tailored component design to be achieved, since complex curvilinear paths can be used per material layer. 3D printing of composites poses several advantages (Kabir et al. [1]) such as the possibility to generate complex shapes without the use of molds; an * Corresponding author. E-mail address: [email protected] (S. Sangaletti). 1 This term will be used throughout the document to refer to the process of FDM of continuous carbon fiber components assumed in this work. Contents lists available at ScienceDirect Theoretical and Applied Fracture Mechanics journal homepage: www.elsevier.com/locate/tafmec https://doi.org/10.1016/j.tafmec.2023.104030 Received 30 April 2023; Received in revised form 11 July 2023; Accepted 24 July 2023 Theoretical and Applied Fracture Mechanics 127 (2023) 104030 2 automatic process that does not require the intervention of human operators, leading to lower production costs; multidirectional and ad-hoc reinforcement geometries can be generated based on the loading history that the component is expected to have to endure. In particular, the latter is very pertinent to design since the designer can optimize components by simply acting on the reinforcement pattern. Strengths can be improved without incurring any penalty in weight. For components presenting stress concentrations, such as holes or notches, continuous fiber bundles can be strategically deposed around the stress concentrator and improve the mechanical response of the component. This highlights the position this technology holds on the production of new generation composite structures. Zhang et al. [2] demonstrated through experiments that specimens with a tailored fiber path versus those of a unidirectional (UD) pattern demonstrated increase in the maximum load the specimen could withstand. Khan et al. [3] quantified this increase in strength, in terms of hole size, testing continuous CF printed coupons of different hole diameters, concluding that the larger the hole, the higher the increase in strength. In an effort to explain this strength increase, Zhang et al. [4], in a computational framework, studied the reinforcement strategy for different loading cases on specimen containing circular holes, concurring that the maximum stress around the hole decreases with an appropriate tailored reinforcement geometry. The same was also studied by Malakhov and Polilov [5], who also focused on the attainable reduction in stress concentration factor around the hole. However, these works only referred to stress levels near the stress raiser and from that deduced the increase in strength. The main contribution of the present work is the application of a non-linear analysis that considers not only the stress alleviation related to the tailored fiber path, but also a fracture mechanics-based theory to assess failure, which is directly influenced by the stress state around the stress concentrator. To capture this phenomenon, the Phase Field (PF) method was used, a technique that has been extensively applied to fracture problems in varying structural material systems, ranging from metal fracture (Ambati et al. [6], Vakili et al. [7]) to composite materials. For the latter, in particular, analyses have been conducted in detail on different scales, from micro-mechanical analysis (Guill´ en-Hernandez at al. [8], Tan and Pa˜ neda [9]), to macro-scale analysis (Reinoso et al. [10], Mitrou at al. [11], Dean et al. [12], Kumar et al. [13,14]). Regarding 3D printing applications, Yvonnet et al. [15] focused on 3D printed nylon via selective laser sintering (SLS) and Sangaletti and García [16] on continuous CF. In the latter work, modelling consisted of matrix and fibers being considered as separated material sections, each one with its own elastic and fracture properties, within the finite element analysis (FEA) model. This, however, significantly increases the complexity of the model, the pre-processing times necessary, with a direct consequence on the computational expense required for the simulation. This becomes even more cumbersome if a scale-up analysis is to be performed and printing parameters are to be modified. In addition, depending on the 3D printer used, this might not also be the best representation for the way the fiber bundle is deposed during the production process. In fact, modelling the 3D printed component using two different separate materials (in Sangaletti and García [16], Onyx and CF bundle were modelled) is the correct approach if a production technique like the one described in Nanya et al. [17] is adopted. The process leaves a sizable distance between strands and material systems leading to low volume fractions of fiber and matrix. In the case though another process is considered, like the one used by the 3D printer MarkTwo by Markforged, this modelling technique is not efficient. In this case, a model that can consider a homogenized representation of the material is required. Thus, it would be suitable to adopt the PF formulation including considerations of anisotropic material behavior. Therefore, this work also aims to propose a method to improve the overall simulation procedure for such 3D printed composite components by using the PF formulation and a convenient way to represent the different material properties and fiber paths. The paper is organized as follows: in Section 2, the methodology adopted in this work is described, step by step, starting from the initial considerations to the generation of the Finite Element (FE) models. In Section 3, the anisotropic PF model is described. In Section 4, unidirectional models are analyzed to check the validity of the adopted PF model. Following, in Section 5, the procedure explained in Section 2 is applied to two case studies, Open-Hole Tension (OHT) and Double Edge- Notched Tension (DENT), checking the differences in the mechanical response resulting from the ad-hoc fiber deposition path. In the last part, conclusions and final remarks are mentioned. 2. Methodology As explained in the introduction, it is of key interest to find a way to efficiently represent continuous fiber reinforced components produced from 3D printing in a numerical setting. Depending on the approach and on the printing procedure in question, this can be successfully done by using different element sections per material type, i.e., fiber and matrix are modeled with separate sections, see Sangaletti and García [16]. This, however, is a tedious task requiring significant pre-processing effort for the FEA model set up. Some authors have also used a homogenized approach like in [5], however, their work did not account for failure. In this work, a homogenized approach able to model damage based on the anisotropic PF model, (explained in Section 3) all awhile keeping the level of complexity low, is used. Let us consider a plate containing a centrally located circular hole, like the one shown in Fig. 1 (a). As the tensile stress is applied in the vertical direction there is a stress concentration region around the hole caused by the deviation of the stress flow lines, due to the presence of the hole itself. Since the CFs have higher mechanical properties, the aim is to align them to these flow lines to attain the maximum advantage of their presence. The procedure starts with a linear elastic FEA on an isotropic model of the component of interest. The software ABAQUS is used in this case and the rest of the procedure is tailored to its capabilities. Within a Python script, the results of the linear elastic analysis are extracted from the data file. Then, the data related to the in-plane stresses are extracted and the directions of the maximum principal in-plane stresses of each element of the model are computed by means of the classical formula [18]: tan(2ξ) = 2 τ xy σ x− σ y (1) where ξ is the angle of the principal direction with respect to global coordinate frame, and σ x, σ y, and τ xy are the normal and shear stress components acting on the element. After this calculation is performed for each element of the model, the element is stored as a new section (*solid section) with a new individual coordinate system. This new coordinate system is noted by the reference frame x ′ ,y ′ ,z ′ , where the x ′ -axis of this local coordinate system is aligned with the direction of the maximum principal in-plane stress calculated previously according to the global reference frame from the linear elastic analysis (Fig. 1 (b)). This step alters the new input file. After this, a new model is obtained, with all the elements having a local coordinate system aligned with the stress flow resulting from the initial analysis (Fig. 1 (c)). The x ′ -axis of the new element material coordinate frame represents the fiber direction in that element and the highest Young’s Modulus (Ex) corresponds to this direction. For considerations of failure events, the fracture toughness will also, for this type of composites, vary per direction and can be modelled in a similar fashion. In the x ′ -axis (along the CF bundle) the fracture toughness, noted as GC(0⋅), is set equal to the inter-laminar fracture toughness of the composite, essentially the fracture toughness for a crack that progresses along the fiber direction. In the y ′ -axis the fracture toughness, noted as GC(90⋅), represents the fracture toughness for crack growth perpendicularly to the S. Sangaletti et al. Theoretical and Applied Fracture Mechanics 127 (2023) 104030 3 CF bundle, which for UD fiber composites is generally significantly higher than the inter-laminar value. It is important to note that all these mechanical aspects of the anisotropic nature of the UD material are taken into account solely by use of the PF theory (presented in Section 3) and internal Abaqus orientations. The CF bundle path is modelled by simply changing element material orientations without requiring intricate maneuvers or generation of partitions, new elements or new element sets in two separate FE models to do the subsequent analysis. In addition, it is important to note that only two-dimensional models are considered in this work, which refer to the laminate as an equivalent single layer, thus providing a homogenized representation of the printed component. Therefore, no delamination between layers is taken into account. The method does suffer from the limitation that different Fiber Volume Fraction (FVF) distribution cannot be simulated. But, as mentioned depending on the printing parameters, this can in fact be a more suitable approach to obtain accurate predictions of the mechanical response. 3. The anisotropic Phase field fracture model The PF approach to fracture, used to account for fracture in the above methodology, is based on the principle of minimization of the total potential energy of a cracked solid and stems from the variational expression of Griffith’s energy balance. Thus, it treats fracturing of solids in a globalized way allowing for an efficient basis to perform component level failure analysis. The theoretical aspects of the method are only in brief explained here and more in-depth details of the formulation can be found in the early works ([16,20]). It is noted that the mathematical convention used has matrices represented by a bold capital letter in nonitalic (i.e., A), vectors bold and italic (i.e., τ ) and scalars simple italic (i. e., S). So, consider the solid of volume V, that contains the crack surface Γ, with the external tractions τ , applied to part of its boundary S, defined by the normal vector  n, as seen in Fig. 2 (a). Its total potential energy will be generally expressed as the summation of the elastic energy Ue, the fracture energy Uf and the external work W (Pham et al. [19]): Π(u,Γ) = Ue+Uf−W=∫V\Γ ψ ( ε (u))dV +∫Γ GcdΓ−∫S τ ⋅ ndS (2) where ψ ( ε (u) ) = 1 2 ε T⋅C⋅ ε is the elastic energy per unit of volume, function of the infinitesimal strains ε (expressed using Voigt notation) and the stiffness matrix of the material C, and Gc is the fracture toughness. Such a functional is function of two unknowns: u, the displacement vector, and Γ, which represents the sharp crack topology. Due to the sharp nature of the crack, a numerical solution to the problem becomes rather difficult, and thus requires special tools. One of those that has been proposed is the PF method which relies on a variational expression of the problem. It consists of a regularization of the crack domain Γ, via a crack surface density function γ, and thus, a switch from a sharp topology (in Fig. 2 (a)) to a diffuse one (in Fig. 2 (b)). This is done with a scalar Phase Field variable φ, which ranges from 0, corresponding to intact material, to 1, corresponding to a fully damaged state [21]. In this way, the energetic quantities of Eq. (2) change to: Ue=∫V g(φ) ψ ( ε (u))dV (3) Uf=∫V gcγ(φ,∇φ,A)dV (4) where g(φ)is an energy degradation function which depends on the PF variable, representing the degradation of the stiffness during fracture, and can be described by the quadratic form g(φ) = (1−φ)2, and gc is a scaling parameter that, used with the scaling constants, allows a correct definition of the fracture toughness along the two principal directions of orthotropy (see Eqs. 6–8). The tensor A is a 2nd order structural tensor which enables the consideration of an anisotropic fracture energy, something imperative to consider when modeling anisotropic materials. This was firstly introduced in the seminal work of Clayton and Knap [22] and then has been used in a variety of applications in the same sense, see Pillai et al. [23] and Bleyer and Alessi [24]. When A is equal to the Fig. 1. Schematic representation of the procedure followed to obtain the final FE model for the analysis. Fig. 2. Schematic solid (a) in a typical form, and (b) with the introduction of the PF (diffuse crack). S. Sangaletti et al. Theoretical and Applied Fracture Mechanics 127 (2023) 104030 4 identity tensor, I, the isotropic case is recovered [25]. After the regularization, the total potential energy in Eq. (2) becomes a function of the displacement vector u and the PF φ and its minimization still leads to the equilibrium solution. Finally, γ is the crack surface density function, whose form, considering the AT2 model, is given as [9–10]: γ(φ,∇φ) = 1 2l(φ2+l2(∇φ⋅A⋅∇φ))(5) where l is the regularization parameter (or length scale as it is often referred to). As mentioned in many studies, care must be given when choosing the value of the length scale, as it directly affects the strength. In this work, it will be considered as a material parameter that will be obtained by fitting a certain set of experimental results [26]. Regarding the structural tensor for the case of an orthotropic material, it can be defined considering two principal directions of orthotropy, a1, a2, as [25]: A=1+a1a1⊗a1+a2a2⊗a2(6) Where a1 and a2 are two scaling constants for each direction a1,a2 and ⊗is the dyadic product. To define them, the approach proposed in Mitrou et al. [11] for multidirectional laminates is followed, where the scaling constants a1 and a2 are directly determined considering the values for the fracture toughness Gc(θ=0◦)and Gc(θ=90◦): a1=(Gc(90◦) gc)2 −1(7) a2=(Gc(0◦) gc)2 −1(8) where θ is the angle of the assumed crack with respect to the loading axis. As stated previously, for UD composites, which are of interest in this work, Gc(0◦)is the inter-laminar fracture toughness, noted by G0◦ C, and Gc(90◦)is equal to the fracture toughness for a crack growing in the direction perpendicular to the fibers, noted by G90◦ C. As stated, the parameter gc can be given an arbitrary value. This can be understood if one looks at the overarching idea that the toughness distribution acts as an input to the method, and that it is solely defined in terms of the scaling constants assuming known and fixed values for the fracture toughness at 0◦and 90◦. Thus, as long as those are properly defined, gc acts as a parameter that scales the scaling constants accordingly. Only in the case that the material is isotropic it is convenient for gc to be considered equal to the fracture toughness of the material, as any scaling would add unnecessary complexity [25]. The specific considerations of anisotropic fracture energy alleviate the method of requiring a posteriori crack path definitions in anisotropic failure. The fracture plane is rather a product of the toughness distribution. This is considered an advantage of the above formulation for the anisotropic PF model, as shown in Mitrou et al. [11]. Most works that refer to an application of the PF method to a UD lamina (or transversely isotropic material) do not adopt the same perspective in definition of the constituents of the anisotropic PF model. While they all use the structural tensor to include anisotropic fracture energy, in the context of transverse isotropy the fracture plane is considered to be along the symmetry axis. Thus, in those cases, the structural tensor only aids in restricting failure on a certain direction. The scaling constants are treated as numerical inputs, and gc as the toughness for that preferred direction [23,24,25]. Without imposing the anisotropic model under that prism, the fidelity in accounting for nonnegligible fracture toughness values along other directions of the transversely isotropic material is not lost, and so is not also the ability to predict a fracture plane that is not along the principal direction of transverse isotropy. The latter is considered very important, especially when assessing fracture of multidirectional laminates modelled through an equivalent single layer approach [11], or fracture of variable stiffness composites, whose fracture planes may not follow the (local) principal directions of orthotropy (or transverse isotropy). In the anisotropic PF method adopted here, the underlying assumption for crack path (or direction) determination lies on a maximum energy release rate (MERR) criterion as outlined in the work of Li and Maurini [27]. While the MERR criterion might not cater to all materials and forms of anisotropic fracture [28], in what follows it will be demonstrated that the internal adaptation of the MERR criterion of the PF models is sufficient to accurately predict the fracture plane angle for the components of interest. 4. Numerical validation In this Section, initial representative cases of structural testing are simulated to validate the aforementioned methodology. Two representative examples, an OHT test and a center-notched tension (CNT) test, are simulated based on the experimental evidence presented in Modniks et al. [29]. Both specimen types are shown in Fig. 3, where geometry, loading and boundary conditions are also depicted. The specimens are made from a quasi-UD flax fiber/epoxy fabric prepregs. The material properties of this system are reported in Table 1. Different fiber off-axis angles with respect to the load direction were tested for both cases: 30˚, 45˚, 90˚for OHT and 15˚, 45˚, 60˚, 90˚for CNT. Following the theory of Section 3, the PF approach to fracture was applied to both models. The length scale was selected as l=9.8 mm to match the experimental evidence for the 30◦tests and was then kept constant for all the other simulations, while gc=10N/mm. Lacking information for the intra-laminar fracture toughness of the material for a crack running perpendicular (transverse) to the fibers, G90◦ c, this value is assumed to be 50 times higher than the longitudinal intra-laminar fracture toughness, G0◦ c, as per Zhang et al. [30]. A minimum discretization length of 0.04 mm was used to create the mesh. The analyses were executed in Abaqus using the UMAT +UMATHT approach mentioned in Mitrou et al. [11] based on Navidtehrani et al. [31]. Both models were discretized with 2D Plane Strain elements (CPE4T) and the solution scheme used was the *type =SEPARATED option for a coupled temperature displacement problem in the implicit solver available in Abaqus. The respective results can be seen in Fig. 4 and Fig. 5. Good agreement against the experimental effective strength, σ f, with off-axis angles for both cases, is shown. It is also noted that, for each off-axis angle, a crack path that follows the fibers direction is found, as is expected. These results give confidence that the present model can be used in more complex cases of non-UD fiber paths, as it is able to accurately capture the effect of change in fiber angle with respect to the loading Fig. 3. Dimensions and configuration of the (a) OHT specimen and (b) CNT specimen. S. Sangaletti et al. Theoretical and Applied Fracture Mechanics 127 (2023) 104030 5 direction all awhile being able to properly capture the expected crack path. 5. Effect of fiber deposition path on strength After the above validation of the model for UD laminates, the effect a tailored fiber path could have on the strength of structural components is analyzed. A quantification of the difference in the maximum attainable load between a unidirectional (UD) specimen versus those that have the fibers distributed at an optimal path (for the given loading condition) is of prominent interest. It is noted here that the latter will be referred to as the case with variable stiffness reinforcement. Two specimen types are considered, Open Hole Tension (OHT) and the Double Edge-Notched Tension (DENT) specimen. These tests are of outmost importance especially in the aerospace field where they can be used to develop design allowables [32] and perform material characterization [33]. The study here focuses on the attainable improvement on strength and how this affects the observed size effect when a scale-up analyses is performed. The latter is interesting to study since the PF approach of fracture used here is a purely energetic criterion (if the length scale tends to zero) or an implicit combination of a stress and energy criterion (if the length scale value is arbitrarily chosen, which is the approach adopted in this work). Energy plays a role that leads intrinsically to a size effect which could not be appreciated by adopting a method that is reliant on a stress criterion alone, and this is demonstrated in what follows. 5.1. Open-Hole tension The method described in Section 2 is here applied in conjunction with the PF method to study the effect the disposition of fibers on an optimal path could play in the structural response of open-hole structures subjected to tension (OHT). To study the size effect, the specimens were scaled up keeping a constant ratio between the width and the hole diameter of w/d=4. The dimensions of the specimens are listed in Table 2 and the configuration is depicted in Fig. 6. The material system of interest is CF/epoxy HTA/6376 with properties listed in Table 3. It is reiterated here that, with the methodology presented in Section 2, it is possible to easily consider the local change of fiber orientation through a rotation of the local orientation of each element, leading to local changes in moduli and fracture toughness as would be observed in the 3D printed component featuring a variable stiffness reinforcement pattern. Using the anisotropic PF model mentioned in Section 3, the difference in fracture toughness between longitudinal and transverse direction is directly considered and the respective values are listed in Table 3, while gc=10 N/mm. The length scale is selected as l=0.04 mm and the model is discretized with elements of minimum length of 0.008 mm, to be sufficiently smaller than the length scale along the expected crack path. This length scale value was arbitrarily chosen given the lack of experimental evidence. This, however, does not limit the analysis which is used to provide qualitative results, and a different value of the length scale would essentially just scale the strength obtained in both cases [26]. For this reason, any results presented hereafter are normalized using the notched tensile strength of the material as occurring from the UD specimen for Size A. The analyses were once again executed in Abaqus using the UMAT +UMATHT approach mentioned in Mitrou et al. [11] based on Navidtehrani et al. [31]. The solution scheme used was the *type =SEPARATED available in the implicit solver available in Abaqus for a coupled temperature displacement problem. Plane strain, CPE4T, elements were used for the smaller specimens and, above 56 mm of width, plane stress, CPS4T, elements were used. The change in element type was applied due to a careful evaluation of the stress field of 3D equivalent elastic models. In the smaller cases, plane stress assumptions did not hold near the notch region, but after 56 mm of width size the stress field at the notch satisfied plane stress assumptions. Fig. 7 depicts the strength obtained per specimen width for both the UD and variable stiffness reinforcement (i.e., optimal fiber deposition). The size effect found for both types of specimens is similar and it is in agreement with the classical results for quasi-brittle materials: when the hole radius Table 1 Properties of the quasi-UD flax fiber/epoxy composite. Ex(GPa) Ex(GPa) Gxy(GPa) vxy G0◦ c(N/mm) 26.0 2.6 1.3 0.35 0.622 Fig. 4. Numerical vs experimental results of UD off-axis OHT. Fig. 5. Numerical vs experimental results of UD off-axis CNT. Table 2 Dimensions of the scaled OHT specimens. w(mm) d(mm) L(mm) Size A 12 3 185 Size B 24 6 185 Size C 40 10 185 Size D 56 14 185 Size E 72 18 185 S. Sangaletti et al. Theoretical and Applied Fracture Mechanics 127 (2023) 104030 6 decreases with respect to the material length, the specimen is less sensitive to the presence of the hole. In the extreme case for very small radii, the apparent strength of the holed specimen should be the same of a specimen without a hole subjected to the same net stress. It is seen that, throughout the full extent of sizes, the specimens with the fibers deposed at the optimal path demonstrate improved strength, with a significant difference between the two types of specimens. Moreover, this increase in strength is more prominent as the size of the specimens increases, as will be discussed and justified later on. 5.2. Double Edge-Notched tension A similar procedure as above for the OHT is followed here for double edge notched tension (DENT) specimens. The specimen geometry is shown in Fig. 8. Once again, scaled specimens are considered with the ratio of the crack length over width kept constant at 2w/a0=10/3. It is noted that, in all cases, the notch radius is fixed to 0.5 mm. The rest of the respective dimensions are reported in Table 4. The material system is the same one considered before and thus the length scale and discretization parameters are kept the same as well. The results for strength over specimen width are shown in Fig. 9. The typical size effect for a cracked plate is found for the two types of specimens, with a strong decrease of strength when the size increases. Note that a tailored fiber deposition does not modify this tendency. The reason is that the situation for large cracks is governed by an energetic balance which promotes this tendency, and this fact is not modified by a fiber deposition tailoring. However, from a quantitative point of view, fiber deposition tailoring does affect the notched strength with a sizable increase in strength being obtained for all the modelled sizes. It is interesting to notice, however, that in this case the magnitude of the increase has a slightly more pronounced growth rate with size than the one observed for the OHT case (Fig. 7). 5.3. Discussion To explain the results obtained above, a closer look is given into the stress fields. Fig. 10 shows the stresses along the fiber direction (S11), at a certain time step for the UD and the variable stiffness reinforcement OHT specimen. For each reinforcement type (i.e., (a) UD, (b) variable stiffness) a complementary closeup image of the stress concentrator Fig. 6. OHT specimen. Table 3 Lamina properties for CF/epoxy HTA/6376 [34]. Ex(GPa) Ey(GPa) Gxy(GPa) vxy G0o c(N/mm) G90o c(N/mm) 114.8 11.7 9.66 0.21 0.27 106.3 Fig. 7. Results for variable stiffness vs UD OHT specimen. Fig. 8. DENT specimen. Table 4 Dimensions of the scaled OHT specimens. 2w(mm) a0(mm) L(mm) Size A 10 3 185 Size B 20 6 185 Size C 40 12 185 Size D 60 18 185 Size E 80 24 185 Fig. 9. Results for variable stiffness vs UD DENT specimen. S. Sangaletti et al. Theoretical and Applied Fracture Mechanics 127 (2023) 104030 7 region and the corresponding material orientations (local axis vectors) are provided. The stresses in both cases have been normalized based on the maximum stress observed for the UD specimen. In the case of variable stiffness reinforcement, the maximum stress is lower than the one in the UD case (around 16% lower, as shown in Fig. 10). It is also interesting to observe the pattern of the stress concentration. It significantly differs between the two configurations. In the UD reinforcement case, stresses are more prominently concentrated in the vicinity of the hole, aligned with the load (and fiber) direction. On the other hand, for the variable stiffness reinforcement case, stresses are more evenly distributed along the ligament area. It is noted that this is shown from an analysis of a single size, since the same pattern can be observed for all other sizes as well. Similarly, the stresses are compared for the DENT specimens as shown in Fig. 11. Once again, the contour plot is obtained for only one size as the pattern repeats itself for the other sizes. As can be seen for the DENT specimen, the maximum stress in the case of variable stiffness reinforcement is close to 60% of the value in the UD, making the effect of the tailored fiber deposition even more prominent in this case. The reason for this decrease in maximum stress is a direct consequence of tailoring the fiber orientations, making them parallel to the load flow directions. As a consequence, the load flow afterwards follows a softer distribution around the hole, decreasing the stress concentration. In some manner, tailoring the fiber deposition modifies the load flow to mimic the stress state of a geometry much better adapted to the load. Thus, if the failure was governed by the local maximum stresses, the increase on strength expected would be given by this decrease of the maximum stress. However, this is only the maximum improvement attainable in quasi-brittle materials, since other metrics beyond the pointwise maximum stress value also play a role, as will be discussed later on. The obtained improvement in strength for both configurations reported in Section 5.1 and 5.2, in Fig. 7 and Fig. 9, is associated with these differences in the observed stress fields. The more pronounced concentration in the UD reinforcement case for both configurations leads to stresses being slightly higher at the concentrator for the same applied remote displacement, leading to an eventual lower remote failure stress. As pointed out briefly in Section 5.1 and 5.2, the increase in notched strength attained is not constant and for both DENT and OHT cases it increases more prominently as the specimen size increases. This difference could be attributed to the effect the size has on the evolution of stress concentration factor (SCF), in particular for DENT specimens where not all the geometric parameters are scaled (the notch radius is fixed). However, this tendency is also observed for OHT where all the geometric parameters are scaled, thus the SCF remains constant, as seen in Fig. 12, which presents the SCF as the width of the specimens is scaled-up. In the case, however, of the DENT specimens the SCF is constantly increasing but at a different rate for UD and variable stiffness reinforcement patterns, with the rate being higher for the UD Fig. 10. Stress field and local material orientations for the (a) UD, (b) variable stiffness, OHT specimen. Fig. 11. Stress field and local material orientations for the (a) UD, (b) variable stiffness, DENT specimen. S. Sangaletti et al. Theoretical and Applied Fracture Mechanics 127 (2023) 104030 8 reinforcement pattern. These patterns are not directly and analogously observed in the obtained percentage increase in apparent strength of the specimens. This should highlight the fact that drawing conclusions for strength improvement of components via the use of variables stiffness reinforcement patterns cannot be deduced solely using the SCF of a single size as the metric. Instead, the results show that other factors clearly play a role in the attainable strength improvement. This is coherent with the fact that failure in quasi-brittle materials is not only governed by the stress field, but the energetic balance can play a key role in many situations. In the particular case of this problem, the increase of improvement of apparent strength can be easily understood analyzing the effect of a fixed length scale (assumed by the PF method) 2 on the results. Since this length scale modulates the size of the damaged region, this region is smaller with respect to the hole size when considering larger specimens. Therefore, the region affected by the damage pattern is more and more concentrated near the hole for larger specimens. Now, concerning the stress state, as is well-known, the effect of the hole on the stress state vanishes far from the hole. Considering that the material orientation was set based on the stress state for the isotropic case, it becomes understood that the material orientation is altered more prominently in a region closer to the hole than far from it. Taking all the former into account, when the specimen is scaled up, the damaged region is restricted to a smaller portion of the overall region affected by the stress concentration. Thus, it is restricted to a smaller portion of the region for which the material orientations have been significantly altered for the case of a variable stiffness reinforcement pattern. As a consequence, for larger specimens, the reinforcement based on an optimal pattern holds a bigger effect in comparison to the UD cases, leading to the result that the improvement in strength is higher for larger specimens (Figs. 7 and 9). Additional experimental evidence would be necessary to provide further information on this. However, as of current this is not possible due to limitations that exist in manufacturing constraints. The patterns that would exactly follow the principal stress lines cannot at the moment be reproduced due to restrictions of the software used by the Markforged 3D printer at the disposal of the authors, which does not allow tailored fiber deposition paths, but only predefined ones from the printer manufacturer. At this point, it is also noted that all the above claims are made under the assumption that the failure pattern, for both OHT and DENT, is that of a fiber failure form with a crack traversing along a line perpendicular to the load direction as seen in Fig. 13. As discussed in [35], this damage pattern for an OHT specimen with UD reinforcement is dependent on the combination of material properties. For some combinations, the damage form observed would instead be split cracks following the fiber line. As stated in [35], for the material system taken into consideration in the analyses described in Section 5.1 and 5.2, a horizontal crack without split cracks is expected and, thus, direct extension of the these observations and methods to materials that might not fall in this category must be done carefully. 6. Conclusion and remarks In this work, an anisotropic PF fracture model has been applied to study the mechanical and fracture behavior of composite structures that could be realized by means of a recent continuous CF 3D printing technique. Components featuring both UD and variable stiffness reinforcement patterns were compared to quantify the strength for OHT and DENT. The observations presented are based on the predicted strength for each reinforcement pattern and configuration, which adds on to the existing works that have to date only compared stress concentration factors and deduced strength increase from there. An approach based on a homogenized representation of the UD material is followed, as described in Section 2, which significantly decreases the complexity and pre-processing time for the model by steering away from individual constituents modeling approach. In conjunction with the PF model, this renders the analysis of failure events in similar configurations much simpler than common practice. Of course, it must be understood that this method refers to cases where a high volume fraction exists in the component, thus homogenization assumptions hold and the final fracture occurs in the form a single fracture plane through the thickness of the plate. The reliability and applicability of the PF model presented here was tested on unidirectional OHT and CNT coupons against experimental evidence found in the literature. The satisfactory agreement obtained demonstrates the ability of the adopted anisotropic PF model to accurately predict the mechanical and fracture behavior of UD components with a varying angle with respect to the loading direction. Therefore, the same approach was applied to variable stiffness composites to quantify improvements in the mechanical response as a result of an ad-hoc fiber deposition path, following the principal in-plane stresses direction. Both case studies demonstrate the beneficial effect of an ad-hoc fiber deposition path with strength increase, for a range of sizes, when damage occurs as a crack developing through fiber breakage (i.e., not parallel to the fibers). Moreover, as observed for both OHT and DENT, this becomes even more prominent as the specimen is scaled up in size. It is recognized that such simulations could help promote new Fig. 12. Comparison of the SCF at the notch for the OHT and DENT specimens. Fig. 13. Failed pattern for (a) OHT specimen (b) DENT specimen. 2 This discussion is generically framed in the context of the PF length scale since it is dealing with anisotropic fracture, but can be traced back to the characteristic material length l=KIc2/ σ c, where KIc is the critical stress intensity factor (or fracture toughness) and σ c is the strength of the material. S. Sangaletti et al. Theoretical and Applied Fracture Mechanics 127 (2023) 104030 9 strategies to simulate the mechanical behavior for components realized in such way, always taking into account the boundaries imposed by the manufacturing method itself. Alongside adding considerations in the PF model able to deal with structural components under multi-axial load states, which would include compressive stress states, such approaches could lead to a faster inclusion of novel 3D printed components in industrial sectors where safety and reliability of mechanical components is of outmost importance. CRediT authorship contribution statement Simone Sangaletti: Conceptualization, Data curation, Methodology, Investigation, Formal analysis, Writing – original draft. Anatoli Mitrou: Software, Data curation, Investigation, Formal analysis, Writing – original draft. Israel G. García: Supervision, Writing – review & editing. Albertino Arteiro: Supervision, Writing – review & editing. Declaration of Competing Interest The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. Data availability Data will be made available on request. Acknowledgments The authors SS and AM would like to kindly acknowledge the funding received from the European Union’s Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement No. 861061 – Project NEWFRAC under which this research is being executed. The author IG is very grateful to the Spanish Ministry of Science and Innovation for the funding through the Project PID2020- 117001 GB-I00. The author AA wants to acknowledge the support of FCT — Fundaç˜ ao para a Ciˆ encia e a Tecnologia, I.P., in the scope of the project UIDB/50022/2020. Finally, the authors extend a thank you to Prof. Jos´ e Reinoso for providing his critical feedback and reviewing on this work. References [1] S.M.F. Kabir, K. Mathur, A.F.M. Seyam, A critical review on 3D printed continuous fiber-reinforced composites: history, mechanism, materials and properties, Compos. Struct. 232 (2020), 111476, https://doi.org/10.1016/j. compstruct.2019.111476. [2] H. Zhang, A. Li, J. Wu, B. Sun, C. Wang, D. Yang, Effectiveness of fibre placement in 3D printed open-hole composites under uniaxial tension, Compos. Sci. Technol. 220 (2022), 109269, https://doi.org/10.1016/j.compscitech.2022.109269. [3] S. Khan, K. Fayazbakhsh, Z. Fawaz, M. Arian Nik, Curvilinear variable stiffness 3D printing technology for improved open-hole tensile strength, Addit. Manuf. 24 (2018) 378–385, https://doi.org/10.1016/j.addma.2018.10.013. [4] H. Zhang, D. Yang, Y. Sheng, Performance-driven 3D printing of continuous curved carbon fibre reinforced polymer composites: a preliminary numerical study, Compos. Part B Eng. 151 (2018) 256–264, https://doi.org/10.1016/j. compositesb.2018.06.017. [5] A.V. Malakhov, A.N. Polilov, Design of composite structures reinforced curvilinear fibres using FEM, Compos. Part A Appl. Sci. Manuf. 87 (2016) 23–28, https://doi. org/10.1016/j.compositesa.2016.04.005. [6] M. Ambati, R. Kruse, L. De Lorenzis, A phase-field model for ductile fracture at finite strains and its experimental verification, Comput. Mech. 57 (2016) 149–167, https://doi.org/10.1007/s00466-015-1225-3. [7] S. Vakili, P. Shanthraj, F. Roters, J.R. Mianroodi, D. Raabe, Phase-Field Modeling of Coupled Brittle-Ductile Fracture in Aluminum Alloys (2022) 1–19. http://arxiv. org/abs/2204.13994. [8] T. Guill´ en-Hern´ andez, A. Quintana-Corominas, I.G. García, J. Reinoso, M. Paggi, A. Tur´ on, In-situ strength effects in long fibre reinforced composites: a micromechanical analysis using the phase field approach of fracture, Theor. Appl. Fract. Mech. 108 (2020) 102621. [9] W. Tan, E. Martínez-Pa˜ neda, Phase field fracture predictions of microscopic bridging behaviour of composite materials, Compos. Struct. 286 (2022), 115242, https://doi.org/10.1016/j.compstruct.2022.115242. [10] J. Reinoso, A. Arteiro, M. Paggi, P.P. Camanho, Strength prediction of notched thin ply laminates using finite fracture mechanics and the phase field approach, Compos. Sci. Technol. 150 (2017) 205–216, https://doi.org/10.1016/j. compscitech.2017.07.020. [11] A. Mitrou, A. Arteiro, J. Reinoso, P.P. Camanho, Modeling fracture of multidirectional thin-ply laminates using an anisotropic phase field formulation at the macro-scale, Int. J. Solids Struct. under Rev. 273 (2023) 112221. [12] A. Dean, P.K. Asur Vijaya Kumar, J. Reinoso, C. Gerendt, M. Paggi, E. Mahdi, R. Rolfes, A multi phase-field fracture model for long fiber reinforced composites based on the Puck theory of failure, Compos. Struct. 251 (2020) 112446. [13] P.K. Asur Vijaya Kumar, A. Dean, J. Reinoso, M. Paggi, A multi phase-field- cohesive zone model for laminated composites: application to delamination migration, Compos. Struct. 276 (2021) 114471. [14] P.K. Asur Vijaya Kumar, A. Dean, J. Reinoso, P. Lenarda, M. Paggi, Phase field modeling of fracture in functionally graded materials: Γ-convergence and mechanical insight on the effect of grading, Thin-Walled Struct. 159 (2021), https://doi.org/10.1016/j.tws.2020.107234. [15] T.T. Nguyen, J. Yvonnet, Q.Z. Zhu, M. Bornert, C. Chateau, A phase field method to simulate crack nucleation and propagation in strongly heterogeneous materials from direct imaging of their microstructure, Eng. Fract. Mech. 139 (2015) 18–39, https://doi.org/10.1016/j.engfracmech.2015.03.045. [16] S. Sangaletti, I.G. García, Fracture tailoring in 3D printed continuous fibre composite materials using the Phase field approach for fracture, Compos. Struct. 300 (2022), 116127, https://doi.org/10.1016/j.compstruct.2022.116127. [17] N. Li, Y. Li, S. Liu, Rapid prototyping of continuous carbon fiber reinforced polylactic acid composites by 3D printing, J. Mater. Process. Technol. 238 (2016) 218–225, https://doi.org/10.1016/j.jmatprotec.2016.07.025. [18] M.H. Sadd, Stress and equilibrium, Elasticity. (2021) 57–82. Doi: 10.1016/b978-0- 12-815987-3.00003-7. [19] K. Pham, H. Amor, J.J. Marigo, C. Maurini, Gradient damage models and their use to approximate brittle fracture, Int. J. Damage Mech. 20 (2011) 618–652, https:// doi.org/10.1177/1056789510386852. [20] J.J. Marigo, C. Maurini, K. Pham, An overview of the modelling of fracture by gradient damage models, Meccanica. 51 (2016) 3107–3128, https://doi.org/ 10.1007/s11012-016-0538-4. [21] P.K. Kristensen, C.F. Niordson, E. Martínez-Pa˜ neda, An assessment of phase field fracture: crack initiation and growth, Philos. Trans. R. Soc. A Math. Phys. Eng. Sci. 379 (2203) (2021) 20210021. [22] J.D. Clayton, J. Knap, Phase field modeling of directional fracture in anisotropic polycrystals, Comput. Mater. Sci. 98 (2015) 158–169, https://doi.org/10.1016/j. commatsci.2014.11.009. [23] U. Pillai, S.P. Triantafyllou, Y. Essa, F.M. de la Escalera, An anisotropic cohesive phase field model for quasi-brittle fractures in thin fibre-reinforced composites, Compos. Struct. 252 (2020), 112635, https://doi.org/10.1016/j. compstruct.2020.112635. [24] J. Bleyer, R. Alessi, Phase-field modeling of anisotropic brittle fracture including several damage mechanisms, Comput. Methods Appl. Mech. Eng. 336 (2018) 213–236, https://doi.org/10.1016/j.cma.2018.03.012. [25] S. Teichtmeister, D. Kienle, F. Aldakheel, M.A. Keip, Phase field modeling of fracture in anisotropic brittle solids, Int. J. Non. Linear. Mech. 97 (2017) 1–21, https://doi.org/10.1016/j.ijnonlinmec.2017.06.018. [26] E. Tann´ e, T. Li, B. Bourdin, J.J. Marigo, C. Maurini, Crack nucleation in variational phase-field models of brittle fracture, J. Mech. Phys. Solids. 110 (2018) 80–99, https://doi.org/10.1016/j.jmps.2017.09.006. [27] B. Li, C. Maurini, Crack kinking in a variational phase-field model of brittle fracture with strongly anisotropic surface energy, J. Mech. Phys. Solids. 125 (2019) 502–522, https://doi.org/10.1016/j.jmps.2019.01.010. [28] G. Duveau, J.F. Shao, J.P. Henry, Assessment of some failure criteria for strongly anisotropic geomaterials, Mech. Cohesive-Frictional Mater. 3 (1998) 1–26, https:// doi.org/10.1002/(SICI)1099-1484(199801)3:1<1::AID-CFM38>3.0.CO;2-7. [29] J. Modniks, E. Sprnin¸ˇ s, J. Andersons, W. Becker, Analysis of the effect of a stress raiser on the strength of a UD flax/epoxy composite in off-axis tension, J. Compos. Mater. 49 (2015) 1071–1080, https://doi.org/10.1177/0021998314528827. [30] P. Zhang, X. Hu, T.Q. Bui, W. Yao, Phase field modeling of fracture in fiber reinforced composite laminate, Int. J. Mech. Sci. 161-162 (2019) 105008. [31] Y. Navidtehrani, C. Beteg´ on, E. Martínez-Pa˜ neda, A unified abaqus implementation of the phase field fracture method using only a user material subroutine, Mater. (Basel). 14 (2021) 1–19, https://doi.org/10.3390/ma14081913. [32] T.R.C. Chuaqui, M.W.D. Nielsen, J. Colton, R. Butler, A.T. Rhead, Effects of ply angle and blocking on open-hole tensile strength of composite laminates: a design and certification perspective, Compos. Part B Eng. 207 (2021), 108582, https:// doi.org/10.1016/j.compositesb.2020.108582. [33] G. Catalanotti, A. Arteiro, M. Hayati, P.P. Camanho, Determination of the mode I crack resistance curve of polymer composites using the size-effect law, Eng. Fract. Mech. 118 (2014) 49–65, https://doi.org/10.1016/j.engfracmech.2013.10.021. [34] Z.Z. Pan, L.W. Zhang, K.M. Liew, A phase-field framework for failure modeling of variable stiffness composite laminae, Comput. Methods Appl. Mech. Eng. 388 (2022), 114192, https://doi.org/10.1016/j.cma.2021.114192. [35] H. Zhang, K. Zhang, A. Li, L. Wan, C. Robert, C.M.´ O. Br´ adaigh, D. Yang, 3D printing of continuous carbon fibre reinforced powder-based epoxy composites, Compos. Commun. 33 (2022), https://doi.org/10.1016/j.coco.2022.101239. S. Sangaletti et al.