Extended finite element method (XFEM) modeling of fracture in additively manufactured polymers

Journal Pre-proof Extended finite element method (XFEM) modeling of fracture in additively manufactured polymers R. Ghandriz, K. Hart, J. Li

PII:

S2214-8604(19)30518-4

DOI:

https://doi.org/10.1016/j.addma.2019.100945

Reference:

ADDMA 100945

To appear in:

Additive Manufacturing

Received Date:

7 May 2019

Revised Date:

6 November 2019

Accepted Date:

11 November 2019

Please cite this article as: Ghandriz R, Hart K, Li J, Extended finite element method (XFEM) modeling of fracture in additively manufactured polymers, Additive Manufacturing (2019), doi: https://doi.org/10.1016/j.addma.2019.100945

This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain. © 2019 Published by Elsevier.

Extended finite element method (XFEM) modeling of fracture in additively manufactured polymers R. Ghandriz1, K. Hart2, J. Li1* 1 2

Department of Mechanical Engineering, University of Massachusetts, Dartmouth, MA 02747 Milwaukee School of Engineering, Milwaukee, WI 53202

*

of

Corresponding author: [email protected]

ro

Highlights

An anisotropic cohesive zone model with XFEM is developed to capture fracture in additively manufactured polymer materials.



The XFEM is able to model crack propagations in 3D printed materials without knowing a priori the crack path.



Parametric studies demonstrate that the ratio of cohesive strengths between inter-layer failure and max principal stress failure largely affects the kinked crack growth paths in 3D printed samples.

ur na

ABSTRACT

lP

re

-p



Jo

The fracture of additively manufactured polymer materials with various layer orientations is studied using the extended finite element method (XFEM) in an anisotropic cohesive zone model (CZM). The single edge notched bending (SENB) specimens made of acrylonitrile-butadienestyrene (ABS) materials through fused filament fabrications with various crack tip/layer orientations are considered. The XFEM coupled with anisotropic CZM is employed to model the brittle fracture (fracture between layers), ductile fracture (fracture through layers), as well as kinked fracture behaviors of ABS specimens printed with vertical, horizontal, and oblique layer orientations, respectively. Both elastic and elastoplastic fracture models, coupled with linear or exponential traction-separation laws, are developed for the inter-layer and cross-layer fracture, respectively. For mixed inter-/cross- layer fracture, an anisotropic cohesive zone model is developed to predict the kinked crack propagations. Two crack initiation and evolution criteria are defined to include both crack propagation between layers (weak plane failure) and crack penetration through layers (maximum principal stress failure) that jointly determine the zig-zag crack growth paths. The anisotropic cohesive zone model with XFEM developed in this study is

1

able to capture different fracture behaviors of additively manufactured ABS samples with different layer orientations. Keywords: Fused filament fabrication; Fracture behavior; Extended finite element method; Cohesive zone model; Anisotropic damage

1. Introduction

re

-p

ro

of

Additive manufacturing (AM), or 3D printing, is being increasingly used in a wide range of areas such as aircraft and automotive components, electronic devices, civil infrastructure, and customized medical implants, etc. AM offers significant advantages over conventional methods for model prototyping. However, the reduced fracture resistance often observed in 3D printed components limits their applications to functional end-use products. AM is a general term encompassing a variety of systems used to create three-dimensional physical parts through a layerby-layer additive process. This permits the manufacturing of intricate shapes and surfaces more easily than by conventional subtractive processes. Fused Filament Fabrication (FFF) is a popular AM technique that uses thermoplastic materials to produce designed parts through the extrusion of filament materials [3, 30], which will be the focus of this study. Several recent studies reported that the tensile properties especially the ultimate strength of AM polymer parts are anisotropic and depend significantly on the layer orientation through which the specimens were built [5,10,12,15,17, 25].

ur na

lP

There are fewer studies on the fracture properties of AM polymer parts. Recently some experiments were conducted to investigate the inter-layer and intra-layer fracture toughness of FFF polymers [4,13,24,26]. In particular, contrasting brittle and ductile fracture behaviors were observed in Acrylonitrile-Butadiene-Styrene (ABS) specimens, demonstrating the inter-laminar and cross-laminar crack propagations according to the crack tip/layer orientations [15].

Jo

On the other hand, very few studies exist on modeling fracture behaviors of AM polymers possibly due to the complexity of crack propagation patterns affected by AM microstructures and interface bondings. Recently a lattice spring model was developed to predict fracture and failure of 3D printed bone-inspired materials [22]. Li et.al [18] proposed an anisotropic cohesive damage model combined with the extended finite element method (XFEM) to predict solution-dependent kinked crack paths in FFF ABS samples printed at 45/-45 raster angles. XFEM is based on the partition-of-unity method and capable of dealing with arbitrary displacement discontinuities due to cracking. XFEM extends the piecewise polynomial function space of conventional finite element methods with extra “enrichment functions” to model crack propagation. The advantage is that the enriched element permits cracking inside and no remeshing or pre-defined crack paths are needed. The crack initiation and propagation process can be modeled by a cohesive traction-separation law defined in a cohesive zone model (CZM). CZM allows decoupling the bulk constitutive behavior and the fracture characteristics of the material. The cohesive damages bridge nascent surfaces of decohesion and govern their separation through an irreversible cohesive law, which 2

also encodes the microscopic physical processes near the crack tip. The macroscopic CZM can include multiple cohesive damage criteria to model the kinked fracture in AM polymers as well as in composite materials [18, 13]. In this study, the approach is generalized to consider layer effects. In particular, both the brittle and ductile fracture behaviors due to inter- or intra-layer fracture in AM polymers will be integrated into a unified computational modeling framework.

2. Model formulation

re

2.1 Problem description and finite element model

-p

ro

of

In this paper, three different FFF ABS layered orientation samples: vertical, horizontal and oblique layer models were considered for different fracture behaviors. A 2D plane strain anisotropic damage model combined with XFEM was developed in the commercial software package ABAQUS (Dassault Systemes Simulia Inc.) through the user subroutines. The elastic material with linear softening is used to model brittle fracture. Two material models: a linear elastic with exponential softening and a nonlinear elastoplastic with linear softening were both considered to capture ductile fracture. After that, a unified model using elastic material with linear or exponential softening was developed to predict fracture in the oblique layered sample. The model was well verified with the cohesive surface method and validated with experimental data for all cases.

ur na

lP

The experiments of three-point single edge notch bending (SENB) fracture on 3D printed acrylonitrile-butadiene-styrene (ABS) samples with different layer orientations were performed in a previous study [15]. The samples were fabricated using fused filament fabrication (FFF) in a Taz 6 desktop printer (Lulzbot; Loveland, Co) with 1.75 mm ABS M30 filament (Stratasys; Eden Prairie, MN) and infill density of 100%. The fracture tests were conducted according to the ASTM E1820 and D6068 standards. Our computational models were built referring to the same model and loading configurations as in the experiments [15].

Jo

The finite element (FE) model of three-point SENB samples was developed in commercial software ABAQUS. For computational efficiency, 2D plane strain elements (CPE4) were used in the FE model with displacement controlled boundary conditions. A 10 mm pre-crack was inserted to the mid-bottom edge of the 100 mm x 20 mm x 10 mm full-size model. Three rigid pins of 6.35 mm in diameter were in contact with the specimen to apply the three-point bending loading. We note that the rigid pins were used instead of point loading to avoid excessive stress concentrations around loading points. We found this strategy effectively improves the numerical convergence in our study. The beam span length was 80 mm. The model was simulated under 1 mm displacement loading control. A sketch of the SENB sample is shown in Fig. 1.

3

of

Fig 1. Sketch of finite element model.

-p

ro

Three different layer orientation SENB fracture specimens with respect to crack tip direction were studied as shown in Fig. 2, including the vertically printed SENB specimen (inter-laminar fracture), horizontally printed SENB specimen (cross-laminar fracture) and obliquely printed SENB specimen (kinked fracture), where 𝜃 shows the crack tip/layer orientation angle.

lP

re

Elastic modulus (E) and Poisson's ratio (𝜐) are selected as 2180 MPa and 0.35 respectively, following from the Stratasys ABS-M30 material data sheet [2]. For fracture properties, the critical strain energy release rate 𝐽𝐼𝐶 for both brittle (inter-laminar) and ductile (cross-laminar) fracture were obtained from the experimental fracture tests [15] as listed in Table 1. These parameters will be used later in the calibration of cohesive zone models for material fracture.

ur na

Y-Direction

Jo

Fig 2. Illustration of the three different layer orientation samples.

Table 1. ABS fracture property from experiments [15]. Fracture Cases Critical Strain Energy Release Rate 𝐽𝐼𝐶 (J/m2) Brittle Ductile Kinked

256 ± 84 2260 256 ~ 2260

2.2 Extended finite element method and cohesive zone model

4

The extended finite element method (XFEM) technique is employed to model the crack propagation phenomena in AM polymers without knowing a priori the crack path in the context of layer orientations. The method can readily handle an evolving crack plane and solution dependent crack propagation direction. In addition, it permits a crack to be located in the element interior. This is made possible through XFEM enrichment functions added on conventional finite element function to incorporate displacement jumps across crack faces and model the singularity at the crack tip using Eq. (1): 𝑢ℎ (𝑥) = ∑𝐼∈𝑁 𝑁𝐼 (𝑥)[𝑢𝐼 + 𝐻(𝑥)𝑎𝐼 + ∑4𝜐=1 𝐹𝜐 (𝑥)𝑏𝐼𝜐 ]

(1)

-p

ro

of

where [𝐻(𝑥)𝑎𝐼 ] represents a Heaviside enrichment term and [∑4𝜐=1 𝐹𝜐 (𝑥)𝑏𝐼𝜐 ] depicts a crack tip enrichment term. H(x) is the Heaviside distribution, 𝑎𝐼 represents the nodal enriched degrees of freedom (DOF) for elements cut by a crack which is known by jump discontinuity, 𝐹𝜐 (𝑥) specifies crack tip asymptotic functions, and 𝑏𝐼𝜐 represents crack tip enriched nodal DOF [14]. Both real and phantom nodes are used in enriched elements to model crack initiation and propagation. Crack is initiated when the cohesive damage initiation criterion is met in an element. The crack evolution happens after the phantom node moves apart along with their corresponding real nodes. The magnitude of the separation follows a Traction-Separation Law (TSL) specified in the cohesive zone model.

𝑇 = {𝑇0(1 −

lP

re

The XFEM can be used in conjunction with the cohesive zone model (CZM) to determine the magnitude of the discontinuity- the displacement jump across the crack faces. The CZM governs this displacement jump by traction-separation cohesive law [11]. In this study, two different forms of TSL– bilinear and exponential, were used to model brittle and ductile fractures, respectively. The bilinear and exponential TSL follows from Eq. (2) [11, 31]: 𝐾𝑛𝑛 . 𝛿 1−exp{−𝛼(𝛿−𝛿𝑖𝑛𝑖𝑡 )/(𝛿0 −𝛿𝑖𝑛𝑖𝑡 )} 1−exp(−𝛼)

ur na

0

𝑖𝑓

𝛿 < 𝛿𝑖𝑛𝑖𝑡

) 𝑖𝑓 𝛿𝑖𝑛𝑖𝑡 ≤ 𝛿 ≤ 𝛿0 𝑖𝑓

(2)

𝛿 > 𝛿0

Jo

where T0 is the cohesive strength, 𝛿 represents the separation of the cohesive element, 𝛿𝑖𝑛𝑖𝑡 is the characteristic separation length when fracture initiates, and 𝛿0 denotes the complete material separation. Here α is an exponential law parameter. A larger value of 𝛼 shows higher exponential curvature, as shown in Fig. 3. Note that 𝛼 = 0 corresponds to the linear softening case. For simplicity, we used 𝛼 = 1 for exponential softening TSL. The initial cohesive stiffness 𝐾𝑛𝑛 can be obtained using the relation between 𝑇0 and 𝛿init : 𝐾𝑛𝑛 = 𝑇0 /𝛿init [7, 29].

5

𝑇0

𝐾𝑛𝑛 𝛿0 𝑒𝑥𝑝𝑜𝑛𝑒𝑛𝑡𝑖𝑎𝑙

𝛿0 𝑏𝑖𝑙𝑖𝑛𝑒𝑎𝑟

of

𝛿𝑖𝑛𝑖𝑡

Fig 3. Schematic of the bi-linear and exponential cohesive traction-separation laws

2.3 Surface-based cohesive method

re

-p

ro

The two independent cohesive parameters are: the maximum cohesive strength 𝑇0 for the peak stress, and the cohesive energy Γ0 equivalent to the area enclosed by the TSL curves. The stress on the cohesive element increases with the slope of initial cohesive stiffness 𝐾𝑛𝑛 , which denotes the elastic behavior until the stress reaches a critical point (𝑇0 ). After that, damage evolution occurs until the cohesive element loses its stress carrying capability when the separation reaches δ0 .

Jo

ur na

lP

In addition to XFEM, the cohesive method is also useful in modeling fracture and delamination between adhesives or bonded interfaces. The cohesive method can be element-based or surfacebased. The cohesive method is mesh dependent and needs the crack path to be predefined along with the cohesive elements or cohesive surfaces. In the cohesive element method, the model is inserted with cohesive elements while the cohesive surface method uses contact pairs between FE element interfaces. The surface-based cohesive method provides an alternative approach to modeling fracture without introducing additional cohesive elements. It is pertinent to note that the fracture property in the cohesive surface method is defined in the surface interaction property, not in the material property. In this study, the cohesive surface method is also employed to verify XFEM results both in brittle and ductile fracture models, where the crack paths can be pre-defined and inserted with the cohesive surface property. The cohesive surface property is defined through the traction-separation law with the same material parameters as in the XFEM-CZM method for numerical verification. 2.4 Model formulation and implementation The isotropic linear elastic with a linear softening TSL material model was used to model the brittle fracture behavior in the inter-laminar case. Note that in the cross-laminar case the plastic zone of a whitening area (also called crazing) was observed in the vicinity of the crack propagation region [15]. This event is proof of nonlinear elastic-plastic deformation. Observing that the plasticity is mainly localized in the fracture deformation region, two models were considered for the nonlinear ductile fracture: the linear elastic material with exponential softening TSL that considers plastic effects during fracture evolution (model 1), and the nonlinear elastic-hardening 6

plastic material with linear softening TSL (model 2). For model 2, the plastic properties are specified according to the Stratasys ABS-M30 material data sheet [2], i.e., yield point at 31MPa with 2% strain and ultimate point at 32MPa with 7% strain. The plastic hardening is assumed to be linear isotropic between yield and ultimate points. The results of the two models will be compared in Section 3. We note that model 1 is easier to be integrated into the unified model for mixed brittle/ductile fracture in oblique layer samples.

Table 2. Cohesive zone model (CZM) parameters CZM parameters

Brittle fracture

Ductile fracture

256

2260

re

Cohesive Energy Γ0 (J/m2)

-p

ro

of

Choosing the values of cohesive parameters plays a key role in achieving accurate model predictions of the fracture simulation results. It was suggested to select the cohesive strength T0 using the ultimate tensile strength and the cohesive energy Γ0 based on the critical strain energy release rates or fracture energy [21, 31]. However, in practice, they are often calibrated through a trial and error process by matching the numerical predictions to experimental results in order to improve the model accuracy. In this study, we found that the CZM parameters can be chosen the same as in the experiments [15]. The values of CZM parameters are listed in Table 2. In addition, a numerical viscosity parameter 𝜇 of 10-6 is applied in ABAQUS/Standard to circumvent numerical convergence difficulties due to the material softening behavior while with little effect on results accuracy.

26

lP

Cohesive Strength T0 (MPa)

31

Note that in oblique layer orientation samples, the crack kinks between and across the layers, leading to mixed-mode loading conditions [16, 19, 20, 23]. The traction-separation law can be extended to mixed-mode fracture using the quadratic function for damage initiation criterion [5]:

ur na

〈𝑇 〉 2

𝑇

2

𝑇

2

(𝑇 𝑛 ) + (𝑇 i−s ) + (𝑇 o−s ) = 1 0−𝐼

0−𝐼𝐼

0−𝐼𝐼𝐼

(4)

Jo

where 𝑇n , 𝑇i−s 𝑎𝑛𝑑 𝑇o−s are the cohesive tractions in the normal (Mode-I), in-plane shear (Mode II) and out-of-plane shear (Mode III) directions, respectively; 𝑇0−i (i=I, II, III) are the cohesive strengths in the three directions, respectively and 〈∙〉 is the Macaulay bracket operator meaning: 〈𝜎𝑛 〉 = {

𝜎𝑛 , 𝜎𝑛 > 0 0, 𝜎𝑛 ≤ 0

(5)

which is to suppress the effects under compressive loading. Also an effective separation 𝛿𝑒𝑓𝑓 is used in the general TSL as: 2 𝛿𝑒𝑓𝑓 = √𝛿I2 + 𝛿II2 + 𝛿III

(6)

where 𝛿i (i=I, II, III) are the cohesive separation in the normal, in-plane and out-of-plane shear directions, respectively. The material rupture happens when 𝛿𝑒𝑓𝑓 reaches the critical displacement at failure 𝛿0 as in Eq. (2). In this study, for simplicity the cohesive strengths 𝑇0−𝑖 (i=I, II, III) are assumed to be the same as 𝑇0 in all three modes directions. In addition, the cohesive parameters 7

are presumed as constants during crack propagation. These simplifications will significantly reduce the model complexity while they have been confirmed to be effective in the study of mixedmode stable tearing crack growth events in ductile fractures [7,8,9,28,32].

〈𝜎11 〉 2

(

𝑇𝐼

𝜎

2

𝜎

2

13 ) + ( 𝑇12 𝐼𝐼 ) + (𝑇 𝐼𝐼𝐼 ) = 1

-p

Weak plane failure criterion:

ro

of

The computational model is implemented in commercial software ABAQUS using its extended finite element method (XFEM) capabilities. XFEM combined with different isotropic and anisotropic damage criteria is able to model distinct fracture behaviors in 3D printing. Note that in the case of vertically printed layers (brittle fracture) and horizontally printed layers (ductile fracture), the crack propagates straightforward and therefore the isotropic cohesive zone model can be assumed. On the other hand, for obliquely printed layers (kinked fracture), the crack propagates alternatively between and across layers. The anisotropic cohesive zone model is thus developed. A material orientation is assigned to reflect the layer orientation in the local x-direction of each finite element. Anisotropic cohesive damage initiation and evolution criterion can thus be developed to capture the local inter-layer bonding directions (fracture weak planes). Specifically, two damage initiation criteria are implemented in ABAQUS through DMGINI user subroutines [22]. To account for the case of off-axis crack propagation, the maximum principal stress (MPS) criteria is included [13]. The weak plane and MPS failure criteria are defined as in Eq. (7):

𝜎𝑝 𝑇𝑝

=1

(7b)

re

Maximum principal stress failure criterion:

(7a)

ur na

lP

where 𝜎𝑖𝑗 (i,j=1,2,3) are the stress components in local material orientations; 𝑇 𝐼 , 𝑇 𝐼𝐼 𝑎𝑛𝑑 𝑇 𝐼𝐼𝐼 represents the cohesive strength in the three fracture modes on the weak plane ; 𝜎𝑝 is the maximum principal stress, and 𝑇𝑝 refers to the cohesive strength on the MPS plane. As noted earlier, for simplicity we consider all three 𝑇 𝐼 , 𝑇 𝐼𝐼 𝑎𝑛𝑑 𝑇 𝐼𝐼𝐼 the same as 𝑇𝑊 (cohesive strength on the weak plane). It is pertinent to note that the ratio of 𝑇𝑝 /𝑇𝑊 significantly affects the crack growth paths. This can be investigated through parametric studies of 𝑇𝑝 /𝑇𝑊 in kinked fracture, which will explain in Section 3.

3. Computational results and discussions

Jo

The numerical methods and material models for different fracture behaviors of ABS samples printed with different layer orientations are summarized in Fig. 4 below. We note that the cohesive surface method is used to verify XFEM-CZM results for brittle and ductile fracture cases. For kinked fracture, the crack path is not pre-defined but solution-dependent, in which only the XFEMCZM method can be applied. All numerical simulations are validated by experimental results of SENB fracture tests reported in [15].

8

Obliquely printed layers (kinked fracture)

• • • • •

• • • •

Methods: XFEM-CZM Cohesive surface Material models: Isotropic elastic + exponential softening TSL • Isotropic elasto-plastic + linear softening TSL

Methods: XFEM-CZM Material models: Isotropic elastic + anisotropic exponential or linear softening TSL

-p

ro

• Methods: • XFEM-CZM • Cohesive surface • Material models: • Isotropic elastic + linear softening TSL

Horizontally printed layers (ductile fracture)

of

Vertically printed layers (brittle fracture)

Fig 1. Schematic of numerical methods and material models for different layer orientation samples.

re

3.1 Brittle fracture results

Jo

ur na

lP

The isotropic elastic material model with linear softening TSL is used to predict the brittle fracture behavior of vertically layered 3D printed samples. To gain confidence in the accuracy of numerical predictions, the convergence study of the finite element mesh density is carried out carefully before the major simulations. Fig. 5 shows the mesh convergence study results of load-displacement curves using the XFEM-CZM method. Three mesh sizes have been considered for mesh convergence: 100x20, 400x80, and 600x120. It was found that the load-displacement curves were converged from above. Therefore, the mesh density of 400x80 is selected, which ensures both mesh convergence and computational efficiency. Both the XFEM-CZM and the cohesive surface methods are then performed independently for numerical verifications. In the XFEM-CZM approach, the triangular traction-separation law (TSL) with linear softening is employed. The cohesive surface method is based on the cohesive surface contacts defined along the crack path and between FE elements. The cohesive surface interaction property follows from the same TSL as in the XFEM-CZM model. Results of the loaddisplacement curves are verified between the two methods and further validated with experimental results as shown in Fig. 6. Overall good agreement between numerical and experimental results is found. The numerical predictions show a bit stiffer than the experimental measurements. We note that the macroscopic continuum model developed in this study does not include the porosity which is evident in 3D printed materials and can reduce strength and fracture properties [1, 17, 27]. In addition, the introduction of defects, residual stresses, and part distortions originated from the manufacturing and curing processes could all together contribute to the differences. A microstructure-based model could further improve the numerical predictions, which is out of the scope of the current study and

9

ro

of

will be pursued in future work. Nevertheless, the macroscopic model effectively captured the brittle fracture behavior in vertical layer orientation samples.

ur na

lP

re

-p

Fig. 5. XFEM mesh convergence study on vertically layered samples.

Fig. 6. Verification (XFEM vs. cohesive surface) and validation on vertically layered samples.

Jo

3.2 Ductile fracture results

To capture the ductile fracture behavior of horizontally layered 3D printed samples, two material models are considered in this study. One is to use the isotropic linear elastic material with exponential softening TSL (model 1). The other one employes the isotropic nonlinear elastichardening plastic material with linear softening TSL (model 2). Both the XFEM and cohesive surface methods for each material model were applied independently and compared for numerical verifications. The mesh convergence study was performed first to ensure sufficient mesh density in fracture simulations. Fig. 7 shows the mesh convergence of XFEM for model 1, indicating that the 400x80 mesh density is sufficient. Similarly, the mesh convergences were also achieved in other material models and numerical methods.

10

The load-displacement curves are shown in Fig. 8 for numerical verifications between the two numerical methods and the two material models. The XFEM and cohesive surface methods agree well with each other. Also, the predictions from the elastic with exponential TSL and elastic-plastic with linear TSL material models are very close. It should be noted that, for the horizontally layered samples with ductile fracture and crazing behaviors observed from experiments [15], the elasticplastic material model is more appropriate while the elastic material model with exponential softening TSL also captures well the overall response, possibly as the plasticity is not much but mainly accompanied by the fracture deformation and material softening in this case. Since model 1 can be readily integrated into a unified model of the elastic material with linear or exponential softening TSL for mixed brittle/ductile fracture, model 1 is thus selected for further study.

ur na

lP

re

-p

ro

of

Figure 9 shows the validation of XFEM and model 1 compared to experimental data [15]. Overall the agreement is good. The best agreement was found until the peak load, after which the XFEM simulation is slightly stiffer than the experimental data. This could be understood as the effect of porosity was not considered in the macroscopic continuum model. In addition, the crack becomes blunt after initiation, especially in the horizontally layered samples where the crazing around crack-tip is significant, resulting in more compliant results in experiments. Finally, it is pertinent to note that the fracture toughness and material strength in the numerical model were the same as those in experimental measurements [15], further validating the selection of the material model and their parameters.

Jo

Fig. 7. XFEM mesh convergence study using elastic material with exponential softening TSL for horizontally layered samples.

11

of

lP

re

-p

ro

Fig 8. Verification (XFEM vs. cohesive surface) of elastic with exponential TSL and elastoplastic with linear TSL models for horizontally layered samples.

ur na

Fig 9. Validation of XFEM and elastic with exponential TSL material model for horizontally layered samples.

3.3 Kinked fracture results

Jo

The isotropic elastic material model with anisotropic linear or exponential softening TSL is developed to study the kinked fracture behavior of obliquely layered 3D printed samples. The kinked fracture is the combination of ductile and brittle fracture cases. The brittle fracture occurs when the crack propagates between the 3D printed layers as weak planes. On the other hand, the ductile fracture could appear as cross-layer crack growth when the maximum principal stress failure criterion is satisfied. The XFEM with anisotropic cohesive damage initiation and evolution criteria is developed in ABAQUS through user subroutines to capture the kinked crack behaviors. To be consistent with the samples tested in experiments [15], a model with local material orientation (layer orientation) of 75° from the crack-tip direction is considered in this study. We note that the methods can be applied to any other obliquely layered samples, providing that the local material orientation is aligned with the layer direction. 12

The mesh convergence study is performed first before the major simulations as shown in Fig. 10 for load-displacement curves at different mesh densities. It is found that overall the mesh density of 400x80 is sufficient for mesh convergence, similar to the study on vertically printed samples. We note that at more refined meshes the numerical convergence difficulties arise possibly due to the crack kinking sensitive to finite element meshes. Therefore, the mesh of 400x80 is used in this study.

ro

of

Both the maximum principal stress (MPS) failure criterion for ductile fracture (with exponential TSL) and the weak plane failure criterion for brittle fracture (with linear TSL) were integrated together in the model to jointly determine the crack initiation and propagation. The experimental measurements of ductile fracture properties in [15] are used for MPS fracture parameters. On the other hand, the fracture properties of the weak plane (between layers) are more likely to vary from the 3D printing process. Parametric studies of varying the ratio between MPS and weak plane failure criteria are thus conducted: 𝑇0𝑃 /𝑇0𝑊 and Γ0𝑃 /Γ0𝑊 , where 𝑇0𝑃 and 𝑇0𝑊 are the cohesive strengths for MPS and weak plane, respectively; Γ0𝑃 and Γ0𝑊 are cohesive energy for MPS and weak planes, respectively. For simplicity, it was assumed that 𝑇0𝑃 /𝑇0𝑊 = Γ0𝑃 /Γ0𝑊 .

Jo

ur na

lP

re

-p

Figure 11 shows the maximum principal stress contours along with the crack paths for parametric studies of three cases: 𝑅 = 𝑇0𝑃 /𝑇0𝑊 = 1.05, 1.48, 1.65. The comparisons of their loaddisplacement curves and experimental results are shown in Fig. 12. The crack initiation and propagation is clearly affected by the ratio R of 𝑇0𝑃 /𝑇0𝑊 , i.e., the competing mechanisms of MPS failure and weak plane failure. When 𝑇0𝑃 /𝑇0𝑊 is small, crack is mainly driven by MPS failure shown in Fig. 11(a). On the other hand, crack travels mainly through inter-layer (weak plane failure) under large 𝑇0𝑃 /𝑇0𝑊 as in Fig. 11(c). During the median range of 𝑇0𝑃 /𝑇0𝑊 , crack transforms more frequently between the inter-layer and cross-layer shown in Fig. 11(b). Similar results were also observed in the parametric studies of 3D printed ABS samples with 45°/-45° raster angles [18]. We note that the zig-zag crack paths in Fig. 11(b) also qualitatively agrees with experimental observations [15]. The value of 𝑇0𝑃 /𝑇0𝑊 = 1.48 is chosen for the best agreement between numerical and experimental results further shown in Fig. 12.

Fig 10. XFEM mesh convergence study on obliquely layered samples. 13

(a)

(b)

(c)

ur na

lP

re

-p

ro

of

Fig 11. Max principal stress contours along with crack paths for parametric studies of obliquely layered 𝑾 samples: (a) R=1.05, (b) R=1.48, and (c) R=1.65, where R=𝑻𝑷 𝟎 /𝑻𝟎 , the ratio of cohesive strength between 𝑷 𝑾 max principal stress failure (𝑻𝟎 ) and weak plane failure (𝑻𝟎 ).

Fig 12. Validation on obliquely layered samples, where R=𝑇0𝑃 /𝑇0𝑊 , the ratio of cohesive strength between max principal stress failure (𝑇0𝑃 ) and weak plane failure (𝑇0𝑊 ).

3.4 Comparison of all three fracture behaviors

Jo

To compare the three different fracture behaviors (brittle, ductile, and kinked fracture) for the three types of 3D printed ABS samples (vertically, horizontally, and obliquely printed layers), the maximum principal stress contours along with the crack paths are illustrated in Fig. 13. The crack propagates upwards in a straight line from its original position (pre-crack) sharply between layers (𝜃 = 0° ) in vertically layered samples, showing a brittle fracture in Fig. 13(a). In this case, the fracture is mainly due to the delamination between layers and modeled by linear softening TSL. Notably, more energy is required to initiate the cross-layer fracture (𝜃 = 90° ) accompanied by more diffusive stress distributions in horizontally layered samples, leading to ductile fracture behaviors as shown in Fig. 13 (b). Accordingly, the crack penetrates through layers as determined by the maximum principal stress failure criterion and modeled by exponential softening TSL. 14

(a)

(b)

-p

ro

of

As for obliquely layered samples, both inter- and cross-layer fractures occur with zig-zag crack paths as shown in Fig. 13(c). This could be understood as the competing mechanisms between weak plane failure (delamination between 3D printed layers) and maximum principal stress failure (crack penetration through 3D printed layers). The kinked fracture behavior is determined jointly by the two failure criteria. We note that these observations of crack initiation and propagation processes also agree well with experimental images of cracked samples in [15].

(c)

re

Fig. 13. Numerical results of max principal stress contours along with crack paths for (a) vertically layered samples, (b) horizontally layered samples, (c) obliquely layered samples.

Jo

ur na

lP

Finally, the load-displacement curves from simulations and experiments for all three types of 3D printed ABS samples are compared in Fig. 14. The brittle fracture behavior is shown in vertically layered samples, depicting the linear-elastic fracture with a sharp drop after the peak load. On the other hand, the ductile fracture is observed in horizontally layered samples, represented by the non-linear elastoplastic fracture with a gradual decrease after peak load. The curves for obliquely layered samples lie between the two cases. Accordingly, the ductile fracture in horizontally layered samples takes the largest peak loads. Note that significant crazing deformation is observed in horizontally layered samples [15], which also indicates much more energy consumed in order to initiate and propagate the crack. In all cases, numerical simulations agree well with experimental results, which validates the computational models developed in this study to capture fracture of 3D printed ABS samples with different layer orientations.

15

of ro

-p

Fig 14. Load-displacement curves from simulations and experiments for all three types of 3D printed samples with vertical, horizontal, and oblique layer orientations.

4. Conclusion

ur na

lP

re

A computational model is developed to capture the experimental results of fracture behaviors in single-edge notched bending tests of 3D printed ABS samples with various layer orientations: vertical, horizontal, and oblique layer/crack tip directions. An anisotropic cohesive zone model (CZM) coupled with the extended finite element method (XFEM) has been developed. The model is well verified with the surface-based cohesive method and validated with experimental results. A key feature of the model is that both the inter-layer fracture (weak plane failure) and cross-layer fracture (max principal stress failure) are taken into account jointly to determine the crack initiation and propagation, which is able to capture the crack kinks observed in obliquely layered samples.

Jo

This study demonstrates the ability of XFEM to model fracture behaviors in 3D printed layered materials. Specifically, when the layers are sloped, the XFEM technique is able to handle an evolving crack propagation direction without knowing a priori the crack path. A 2D finite element model based on ABAQUS user subroutines was developed to implement both the weak plane and maximum principal stress (MPS) failures. With parametric studies of varying the ratios between fracture parameters for the two failure criteria, different crack propagation patterns were revealed and the best agreement between numerical and experimental results was then achieved. It is pertinent to note that the numerical predictions show a bit stiffer than experimental results. This could be understood as the porosity is not considered in the macroscopic continuum model developed in this study. Besides, the material defects, residual stresses, and part distortions introduced during the additive manufacturing processes could all contribute to the differences. A microstructure-based model can further improve the predictions in this regard, which will be the future work. Besides, 3D simulations of the model can reveal the effects of 3D printed layers coupled with structural features as a possible toughening mechanism, which is being actively pursued by the authors. Nevertheless, the macroscopic model developed in this study effectively 16

captured different fracture behaviors in 3D printed ABS samples with different layer orientations. The computational model can be also applied to fracture in a range of other layered materials (e.g., shales, bones, laminated composite materials, etc.).

Conflicts of Interest This manuscript has not been published elsewhere in part or in entirety, and is not under consideration by another journal. There are no conflicts of interest to declare.

Acknowledgment

Jo

ur na

lP

re

-p

ro

of

The authors gratefully acknowledge the financial support provided by the start-up funds and multidisciplinary seed funds at the University of Massachusetts Dartmouth.

17

References [1] Abbott, A. C., Tandon, G. P., Bradford, R. L., Koerner, H., & Baur, J. W. (2018). Process-structure-

property effects on ABS bond strength in fused filament fabrication. Additive Manufacturing, 19, 29–38. https://doi.org/10.1016/J.ADDMA.2017.11.002 [2] ABS-M30 Production-Graded Thermoplastic for FDM 3D Printers Data Sheet, ed: Stratasys Ltd,

2013. [3] Alafaghani, A., Qattawi, A., Alrawi, B., & Guzman, A. (2017). Experimental Optimization of Fused

of

Deposition Modelling Processing Parameters: A Design-for-Manufacturing Approach. Procedia Manufacturing, 10, 791–803. https://doi.org/10.1016/J.PROMFG.2017.07.079 [4] Aliheidari, N., Tripuraneni, R., Ameli, A., & Nadimpalli, S. (2017). Fracture resistance measurement

Polymer

Testing,

60,

94–101.

ro

of fused deposition modeling 3D printed polymers. https://doi.org/10.1016/J.POLYMERTESTING.2017.03.016

[5] Camanho, P. P., Davila, C. G., & de Moura, M. F. (2003). Numerical Simulation of Mixed-Mode

-p

Progressive Delamination in Composite Materials. Journal of Composite Materials, 37(16), 1415–1438. https://doi.org/10.1177/0021998303034505 [6] Cantrell, J., Rohde, S., Damiani, D., Gurnani, R., DiSandro, L., Anton, J., Ifju, P. (2017).

re

Experimental Characterization of the Mechanical Properties of 3D Printed ABS and Polycarbonate Parts (pp. 89–105). Springer, Cham. https://doi.org/10.1007/978-3-319-41600-7_11 [7] Chen, X., Deng, X., & Sutton, M. A. (2013). Simulation of stable tearing crack growth events using

lP

the cohesive zone model approach. Engineering https://doi.org/10.1016/J.ENGFRACMECH.2012.12.017

Fracture

Mechanics,

99,

223–238.

[8] Chen, X., Deng, X., Sutton, M. A., & Zavattieri, P. (2014). An inverse analysis of cohesive zone

ur na

model parameter values for ductile crack growth simulations. International Journal of Mechanical Sciences, 79, 206–215. https://doi.org/10.1016/J.IJMECSCI.2013.12.006 [9] Chen, X., Deng, X., Sutton, M. A., & Zavattieri, P. (2014). Simulation of mixed-mode I/III stable

tearing crack growth events using the cohesive zone model approach. International Journal of Fracture, 189(1), 59–75. https://doi.org/10.1007/s10704-014-9962-4 [10] Cordisco, F. A., Zavattieri, P. D., Hector, L. G., & Carlson, B. E. (2016). Mode I fracture along

Jo

adhesively bonded sinusoidal interfaces. International Journal of Solids and Structures, 83, 45– 64. https://doi.org/10.1016/J.IJSOLSTR.2015.12.028 [11] Elices, M., Guinea, G. V., Gómez, J., & Planas, J. (2002). The cohesive zone model: advantages,

limitations, and challenges. Engineering https://doi.org/10.1016/S0013-7944(01)00083-2

Fracture

Mechanics,

69(2),

137–163.

[12] Espalin, D., Alberto Ramirez, J., Medina, F., & Wicker, R. (2014). Multi-material, multi-technology

FDM: exploring build process variations. https://doi.org/10.1108/RPJ-12-2012-0112

Rapid Prototyping Journal, 20(3), 236–244.

18

[13] Feerick, E. M., Liu, X. (Cheryl), & McGarry, P. (2013). Anisotropic mode-dependent damage of

cortical bone using the extended finite element method (XFEM). Journal of the Mechanical Behavior of Biomedical Materials, 20, 77–89. https://doi.org/10.1016/J.JMBBM.2012.12.004 [14] Giner, E., Sukumar, N., Tarancón, J. E., & Fuenmayor, F. J. (2009). An Abaqus implementation of

the extended finite element method. Engineering Fracture Mechanics, 76(3), 347–368. https://doi.org/10.1016/J.ENGFRACMECH.2008.10.015 [15] Hart, K. R., & Wetzel, E. D. (2017). Fracture behavior of additively manufactured acrylonitrile

butadiene styrene (ABS) materials. Engineering https://doi.org/10.1016/j.engfracmech.2017.03.028

Fracture

Mechanics,

177(177),

1–13.

of

[16] Kamat, S. V., & Hirth, J. P. (1995). Mixed Mode Fracture Toughness of Engineering Materials.

Journal of Engineering Materials and Technology, 117(4), 391. https://doi.org/10.1115/1.2804731 [17] Koch, C., Van Hulle, L., & Rudolph, N. (2017). Investigation of mechanical anisotropy of the fused

ro

filament fabrication process via customized tool path generation. Additive Manufacturing, 16, 138–145. https://doi.org/10.1016/J.ADDMA.2017.06.003

-p

[18] Li, J., Yang, S., Li, D., & Chalivendra, V. (2018). Numerical and experimental studies of additively

re

manufactured polymers for enhanced fracture properties. Engineering Fracture Mechanics, 204, 557– 569. https://doi.org/10.1016/J.ENGFRACMECH.2018.11.001 [19] Li, S., Thouless, M. D., Waas, A. M., Schroeder, J. A., & Zavattieri, P. D. (2006). Mixed-mode

lP

cohesive-zone models for fracture of an adhesively bonded polymer–matrix composite. Engineering Fracture Mechanics, 73(1), 64–78. https://doi.org/10.1016/J.ENGFRACMECH.2005.07.004 [20] Li, S., Thouless, M. D., Waas, A. M., Schroeder, J. A., & Zavattieri, P. D. (n.d.). Competing failure

mechanisms in mixed-mode fracture of an adhesively-bonded polymer-matrix composite. Retrieved from http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.551.7563&rep=rep1&type=pdf

ur na

[21] Li, W., & Siegmund, T. (2002). An analysis of crack growth in thin-sheet metal via a cohesive zone

model. Engineering Fracture Mechanics, 69(18), 2073–2093. https://doi.org/10.1016/S00137944(02)00013-9 [22] Libonati, F., Cipriano, V., Vergani, L., Vergani, L., & Buehler, M. J. (2017). Computational

Jo

Framework to Predict Failure and Performance of Bone-Inspired Materials. ACS Biomater. Sci. Eng.20173123236-3243, https://doi.org/10.1021/acsbiomaterials.7b00606 [23] Liu, S., Chao, Y. J., & Zhu, X. (2004). Tensile-shear transition in mixed mode I/III fracture.

International Journal of Solids and https://doi.org/10.1016/J.IJSOLSTR.2004.04.044

Structures,

41(22–23),

6147–6172.

[24] Sugavaneswaran, M., & Arumaikkannu, G. (2018). Additive manufactured multi-material structure

with directional specific mechanical properties based upon classical lamination theory. Rapid Prototyping Journal, 24(7), 1212–1220. https://doi.org/10.1108/RPJ-06-2017-0118

19

[25] Malik, I., Mirkhalaf, M., & Barthelat, F. (2017). Bio-inspired “jigsaw”-like interlocking sutures:

Modeling, optimization, 3D printing and testing. Journal of the Mechanics and Physics of Solids, 102(2017)224-238. http://dx.doi.org/10.1016/j.jmps.2017.03.003 [26] McLouth, T. D., Severino, J. V., Adams, P. M., Patel, D. N., & Zaldivar, R. J. (2017). The impact

of print orientation and raster pattern on fracture toughness in additively manufactured ABS. Additive Manufacturing, 18, 103–109. https://doi.org/10.1016/J.ADDMA.2017.09.003 [27] Owolabi, G., Peterson, A., Habtour, E., Riddick, J., Coatney, M., Olasumboye, A., & Bolling, D.

(2016). Dynamic response of acrylonitrile butadiene styrene under impact loading. International Journal of Mechanical and Materials Engineering, 11(1), 3. https://doi.org/10.1186/s40712-016-0056-0

of

[28] Park, K., Choi, H., & Paulino, G. H. (2016). Assessment of cohesive traction-separation

relationships in ABAQUS: A comparative study. Mechanics Research Communications, 78, 71–78. https://doi.org/10.1016/J.MECHRESCOM.2016.09.004

ro

[29] Simulia Abaqus User’s Manual, Dassault Systémes Simulia Corp., Providence, RI, 2016.

[30] Vastola, G., Pei, Q. X., & Zhang, Y.-W. (2018). A predictive model for porosity in powder-bed

-p

fusion additive manufacturing at high beam energy regime. Additive Manufacturing, 22, 817–822. https://doi.org/10.1016/J.ADDMA.2018.05.042

re

[31] Xu, XP., & Needleman, A. (1994). Numerical simulations of fast crack growth in brittle solids.

Journal of the Mechanics and Physics of Solids, 42(9), 1397–1434. https://doi.org/10.1016/00225096(94)90003-5

lP

[32] Zavattieri, P. D. (2006). Modeling of Crack Propagation in Thin-Walled Structures Using a

Journal

of

Applied

Mechanics,

73(6),

948.

Jo

ur na

Cohesive Model for Shell Elements. https://doi.org/10.1115/1.2173286

20