Computing the bond strength of 3D printed polylactic acid scaffolds in mode I and II using experimental tests, finite element method and cohesive zone modeling

The advent of the Three-Dimensional (3D) printing technique, as an Additive Manufacturing technology, made the manufacture of complex porous scaffolds plausible in the tissue engineering field. In Fused Deposition Modeling based 3D printing, layer upon layer deposition of filaments produces voids and gaps, leading to a crack generation and loose bonding. Cohesive zone model (CZM), a fracture mechanics concept, is a promising theory to study the layers bond behavior. In this paper, a combination of experimental and computational investigations was proposed to obtain bond parameters and evaluate the effect of porosity and microstructure on these parameters. First, we considered two different designs for scaffolds beside a non-porous Bulk design. Then, we performed Double cantilever beam and Singe Lap Shear tests on the 3D printed samples for Modes I and II, respectively. Afterward, we developed the numerical simulations of these tests using the Finite element method (FEM) to obtain CZM bond parameters. Results demonstrate that the initial stiffness and cohesive strength were pretty similar for all designs in Mode I. However, the cohesive energy for the Bulk sample was approximately four times of porous samples. Furthermore, for Mode II, the initial stiffness and cohesive energy of the Bulk model were five and four times of porous designs while their cohesive strengths were almost the same. Also, using cohesive parameters was significantly enhanced the accuracy of FEM predictions in comparison with fully bonded assumption. It can be concluded that for the numerical analysis of 3D printed parts mechanical behavior, it is necessary to obtain and suppose the cohesive parameters. The present work illustrates the effectiveness of CZM and FEM combination to obtain the layer adhesive parameters of the 3D printed scaffold.


Introduction
The cutting-edge 3D printing technology fabricates objects by deposition of materials in the layer upon layer fashion [1,2]. According to ASTM F42, printing processes can be divide into seven sub-groups of extrusion-based, powder bed fusion, binder jetting, material jetting, directed energy deposition, sheet lamination, and vat photo-polymerization processes [3]. Extrusion-based FDM, as a straightforward and cost-effective method [4,5], provides users to fabricate complex three-dimensional parts quickly by deposition of the melted filaments through a nozzle on a building platform.
Nowadays, 3D printing technologies have been used in tissue engineering purposes since they allow us to create scaffolds with a complex 3D microstructure and a controllable porosity, pore shape, and size [6,7]. The porosity facilitates cells to attach and proliferate by providing easy access to microenvironment supplied nutrition and oxygen [8]. Scaffolds located in the body of a patient must simulate both mechanical and biological properties of the native tissue to bear the mechanical stress of the surrounding microenvironment and not to stimulate the immune system [9].
Various types of biocompatible and biodegradable polymers were used in manufacturing 3D printed scaffold [10], e.g., Polycarbonate (PC), Polyglycolic acid (PGA), Acrylonitrile Butadiene Styrene (ABS), Polycaprolactone (PCL), and Poly-Lactic Acid (PLA) [11]. PLA is the FDA approved hydrophobic aliphatic polyester for use in different medical applications [12]. Besides its biological properties, it has satisfying mechanical properties such as thermal stability, low viscosity, high strength, and highelastic modulus made PLA a well-suited polymer for the FDM technology [13].
The mechanical properties of additively manufactured PLA bone scaffolds have been investigated by several researchers. Tayton et al. provided a platform to compare mechanical shear properties of bone scaffolds made from PLA and PLA-Hydroxyapatite [14]. Naghieh et al. predicted the mechanical properties of porous PLA bone scaffolds fabricated by FDM [6]. This research reported the scaffold compressive elastic modulus by conducting the experimental compression test and FEM. Gremare et al. performed tensile tests on the 3D printed PLA bone scaffold to examine the effect of pore dimension on the mechanical properties [15]. It is noticeable that the numerical FEM approach was used repeatedly in the literature to predict the mechanical response of FDM parts in silico very well [1,6,16]. On the other hand, in the experimental approach, fabrication of numerous specimens and performing mechanical tests to find an appropriate mechanical property for FDM parts are costly and time-consuming.
Using 3D printing in tissue engineering field to fabricate scaffold made scientists to focus on investigation of mechanical properties of 3D printed scaffold. Shear, compressive and tensile properties of 3D printed scaffold have studied experimentally and numerically in different studies. However, the fracture properties of 3D printed scaffold has not been investigated very well yet and we aimed to demonstrate the fracture behavior of 3D printed scaffolds in Modes 1 and two of loading.
On the other hand, in FDM process, by deposition of materials in the layer upon layer fashion, some voids and gaps generate among bonded struts in the final object. Moreover, depositing these struts upon each other does not make a uniform part [17]. So, it is susceptible to have crack initiation and propagation between deposited layers [18], which leads to ultimate failure. This phenomenon can be modeled and analyzed by the fracture mechanicsbased CZM to predict the behavior of scaffold correctly. So, it prevents the overestimating of the fracture properties of scaffold and more accurate estimation of scaffold failure in the in-vitro and in-vivo studies. Moreover, in this study, the effect of microstructure on the fracture properties of the 3D printed bone scaffold was investigated by choosing two different designs for them.

Cohesive zone model (CZM)
The cohesive zone is a softening region located ahead of a crack tip. When stress at the tip reaches the cohesive strength, and enough energy is available for a new cracked area creation, an existing crack starts to propagate. Afterward, separation will occur between two adjacent surfaces across an extended crack [19,20]. Cohesive zone cracks initiation and propagation can be modeled based on crack opening and needed traction relationship using bilinear Traction-separation law (TSL). In TSL law, it assumes that the relationship between traction and separation consists of two linear regions, shown in Fig. 1. Three parameters describe these linear regions: initial stiffness, cohesive strength, and fracture energy, defined in the following section.

Initial stiffness
Initial stiffness relates to the initial crack response before reaching cohesive traction to cohesive strength [17]. For calculating the initial stiffness, Eq. 1 [21] or setting an artificially high value from 10 6 N/mm 3 [22] to 10 8 N/mm 3 [23] can be employed, and in this study, we use the latter method.
In Eq. 1, E is the relevant elastic modulus, and t is the thickness of plies adjacent to a bonding region.

Cohesive strength
Damage initiation determines by the cohesive strength, and after reaching cohesive traction to the cohesive strength, it starts to reduce due to damage accumulation [17]. The adhesive strength calculates from Eq. 2.
where P C is the critical load, and A is the adhesive area. In this research, we calculated the cohesive strength initially from Eq. 2 and then calibrated it.

Cohesive energy
The fracture energy is equal to the area under a traction-separation curve. The separation will happen when the cohesive zone dissipated energy reaches the fracture energy. The fracture energy release rate of Mode I (G I ) can be obtained from Eq. 3, presented in ASTM D5528 [24].
where P is the load, δ is the load point displacement, b is sample width, and a is debonding length. In this paper, we initially calculated the fracture energy from Eq. 3, and then we calibrated it. The identification of a cohesive law and its parameters is still highly controversial. There is no effective experimental method to directly measure the parameters like mechanical stresses near the crack tip because the fracture zone is tiny. Therefore, by correlating the experimental data and numerical simulation, inverse techniques such as FEM have been used recently to estimate the CZM parameters [25][26][27][28][29][30]. Furthermore, fracture mechanics theories have been implemented by researchers to investigate the fractural properties in additively manufactured parts. Kishore et al., by considering the effect of layer temperature on interlayer strength, made an infrared preheating system to increase the overall strength of parts fabricated by additive manufacturing [31]. Spoerk M. et al. extracted optimized printing parameters (layer thickness, layer design, and temperature) to boost inter and intra-layer strength of parts fabricated by FDM [32]. Seppala et al. developed a framework from the polymer interdiffusion viewpoint to understand the interlayer-layer strength in material extruded parts [33]. Recently, researchers have implemented the CZM, a promising theory in fracture mechanics, to describe interlayer delamination, study interlayer fracture, and predict cohesive parameters of AM parts. Liravi et al. developed a framework using experiment and FEM cohesive zone simulation to extract the delamination force of a bottom-up constrained surface projectionbased Stereolithography (SLA) [34]. Ahmadi et al. implemented CZM in the FEM package to understand the effects of various microstructural properties on the macroscopic mechanical properties of Selective laser melting (SLM) parts [35]. Spackman et al. developed a CZM based FEM model to capture the adhesion behavior of fiber-reinforced soft composite additive manufactured parts [36]. Park et al. characterized the fracture behavior of truss structure parts manufactured by the material extrusion process by proposing a modified DCB test and CZM [17]. The modified model consists of two adjacent layers, where the CZM is defined between these two layers to achieve interlayer strength. Fonseca et al. proposed an experimental procedure and CZM inverse procedure to characterize interlaminar fracture in Mode I loading of Fused filament fabricated (FFF) parts [37]. To this aim, they conducted the DCB experiment for 3D printed parts, which are manufactured by two types of material: pure and short fiber reinforced Polyamide 12 (PA12). This experimental procedure was followed with numerical simulation to extract CZM parameters.
Researchers have investigated the adhesive-bonded joint in additively manufacture parts by conducting Single lap shear (SLS) test. S. Kumar et al. evaluated the potential of 3D printed multi-material adhesive interfaces by exploring their compliance-tailored adhesive designs using the SLS test [38]. Kovan [42].
Previous researches can divide into two main groups: one group focused on the experimental and numerical investigation of scaffold mechanical properties manufactured by FDM. Another group concentrated on CZM modeling of FDM parts and conducting DCB and SLS tests. Our research goal is investigating the bonding strength of 3D printed scaffolds using the CZM approach instead of assuming fully bonded layers. We conducted this research for Modes I and II and considered the effect of porosity and microstructure on the bond behavior. To reach this goal, we develop a research framework composed of experimental and numerical processes, includes mechanical tests of PLA filament tension, DCB, and SLS along with FEM simulations, shown in Fig. 2.
First, the process of designing and manufacturing specimens for mechanical tests will be explained. Then, the procedures for each experimental test will be described. Afterward, the experimental method to calculate CZM for Modes I and II and the FEM based numerical modeling containing inverse framework to estimate the cohesive parameters from experimental data will be specified. Then, cohesive parameters of Modes I and II for all designs will be extracted from the curve fitting process. Finally, the mechanical properties of 3D printed parts with and without cohesive assumption in numerical models will be compared.

Specimen fabrication for DCB and SLS tests
The specimens were fabricated using commercially available white PLA filament with a mean diameter of 1.75 mm, a melting point of 190 °C, and a density of 1.25 g.cm −3 (ZF Company, China). We performed the printing process using the Quantum 2025 FDM printer (Persia 3DPrinter Co., Iran) equipped with a brass nozzle with an inner diameter of 0.4 mm and a bed size of 200*200 mm 2 .
We provide a Bulk design to investigate the porosity effect on the cohesive parameters by comparing the mechanical behavior of porous designs with it. These parts are depicted in Fig. 3. Two types of microstructures for scaffold were manufactured and expressed as Lattice and Shifted designs. The difference between these designs was the location of deposited struts in which there was an offset equal to half of the center to center distance of horizontal struts in Shifted design compare to Lattice design. The diameter of each strut is 400 μm, according to the extrusion nozzle diameter. The layer height was 300 μm for all designs to impose enough layer penetration and bonding region. The center to center distances of two adjacent horizontal struts for Bulk and porous designs were 400 and 800 μm, respectively. It is noticeable that the average pore size of a scaffold for osteons growing is in the range of 300-400 μm [43]. In this study, we assume the pore size of the scaffold in this range, 400 μm, equal to center to center distances of two adjacent horizontal struts.

Filament test
Uniaxial tensile test of extruded filament is challenging because of its small scale and the limitation of fixtures and testing machines. Here, a particular fixture was employed to stabilize the sample during testing (Fig. 4).
First, trapezoidal PLA tabs similar to the dog-bone specimens head were 3D printed to bind the filaments within the tabs using superglue. We performed the filament test using SANTAM SMT-5 (SANTAM ENG. DESIGN CO., Iran) instrumented with a 1kN load cell at a tension rate of 1 mm/ min. The filaments' lengths were 80 ± 1 mm, and three filaments were tested. We calculated the elastic modulus from Fig. 2 The research framework was exploited in this study. From the top left, initially, we used the FDM parameters to Computeraided design (CAD) of the 3D printing part. After the manufacturing of specimens, we performed the experimental tests on the filament, DCB, and SLS samples. Then, the FEM simulations with various cohesive parameters were applied, and finally, we extracted the CZM parameters for each design and test by comparing numerical and experimental data  Tensile test configuration for PLA filament. First, the 3d printed PLA dog-bone parts were glued to the filament, and then the tensile test was performed (n = 3) the slope of the initial linear region of the stress-strain curve extracted from the tensile test.

DCB test
The specimens for the DCB test were fabricated according to the ASTM D5528 recommended dimensions, depicted in Fig. 5a. A glue tape was placed between the 21st and 22nd layers during the printing process to create an initial crack (Fig. 5b) in all designs. Five specimens were manufactured for each design, and the piano hinges were installed to specimens for exerting load (Fig. 5c). The DCB test was performed using universal testing machine SANTAM SMT-5 (SANTAM ENG. DESIGN CO., Iran) instrumented with a 1kN load cell (Fig. 5d) at an extension rate of 2 mm/min.
Besides, in attempting to find the cohesive parameters of Mode I by analyzing the DCB fracture test in ABAQUS, a 3D model of DCB has been implemented, shown in Fig. 5e-g. In order to reduce the time of calculation, we modeled four layers of the porous scaffold.
For defining the isotropic elastic material model for PLA, Young's modulus was measured as 3180 MPa from the tensile filament test, explained in Sect. 2.2, and the Poisson's ratio was chosen 0.36 based on the related study [44]. Also, we simulated the loading condition in the test by the tensile displacement in the perpendicular direction to the Fig. 5 The DCB test specifications in experimental (a-d) and numerical (e-g) studies. a The specimen dimensions based on ASTM D5528, b Creation of the initial crack during the printing process, c Installing piano hinge, d Performing DCB test, e Tensile displace-ment to simulate loading regime, f Constraining all degree of freedom at the end of DCB specimen, g Using quadratic tetrahedral (C3D10) elements for meshing of Lattice, Shifted and Bulk designs initial crack plane for all designs (Fig. 5e) and constraining all active structural degrees of freedom at the other end (Fig. 5f). We meshed the DCB models with the quadratic tetrahedral (C3D10) element shown in Fig. 4g. Since the calculated properties of the cohesive surfaces are meshdependent, the effect of mesh size on FEM results has been investigated to eliminate this dependency. Initially, by keeping the value of all CZM parameters constant and changing the mesh size, the load-displacement curve was monitored. Since all parameters were constant, the reason for the fluctuation of the load-displacement curve was mesh size. Then the procedure of refining mesh size was performed until the value of the load-displacement curve has not changed dramatically.

SLS test
The SLS specimens were manufactured by FDM according to ASTM D3163 and based on its recommended dimensions shown in Fig. 6a. Manufactured specimens were pulled apart using the Universal Testing Machine (Zwick, GA, USA) in tension, equipped with a 25 kN load cell (Fig. 6b). The crosshead speed was 1.3 mm/min, and we tested five specimens for each design. Based on ASTM D3163, the grip must be a 1 by 1-inch square, and we tightened it sufficiently to prevent slipping during testing.
To simulate the behavior of the SLS test of three different designs and extract cohesive parameters in Mode II, we assigned the isotropic elastic material model for PLA to all models similar to the DCB test. We simulated the boundary condition of the SLS test by constraining all active structural degrees of freedom at one end (Fig. 6c) to simulate real clamping conditions in the machine grips and exerting a tensile displacement at another end of the specimen in the longitudinal direction (Fig. 6d). The SLS models meshed, with the quadratic tetrahedral (C3D10) element, are shown in Fig. 6e, and then the mesh independency analysis (discussed in Sect. 2.3) was performed.
Since the elastic properties of the cohesive surfaces are mesh-dependent [], the effect of mesh size on FEM results has been investigated to eliminate this dependency. By keeping constant the value of all CZM parameters and changing the mesh size, the load-displacement curve is monitored. Since all the CZM parameters were constant, the only reason for the fluctuation of the load-displacement curve was mesh size. The mesh size was selected which the value of the load-displacement curve has not changed dramatically.
To obtain the bonding parameters, we used the CZM for the numerical simulation of a DCB and SLS tests for Modes I and II, respectively, and implemented in the ABAQUS 2018 (Dassault Systemes Simulia Corp, France) FEM package. The surface-to-surface contact with the definition of damage initiation and propagation was assigned to adjacent layers. Cohesive behavior was characterized by assigning the value of initial stiffness (K n , K s , K t ) to the model. In this study, in Eq. 1, E and t were the elastic modulus and thickness of plies adjacent to a bonding region, which were 3180 MPa and 6.2 mm respectively and so, K was calculated as 5e12. An artificially high values of 1e6 [22] and 1e8 [23] or the calculated K were the initial guesses as input for FEM simulation. Then, they needed to be calibrated by trial and error procedure to fit the experimental curves. Initially, 1e6 was the first guess and input for FEM simulation. Then we increased the value of stiffness step by step. Finally, there was no significant difference between 5e12 and 1e12 in load-displacement curve, so the value of 1e12 was chosen. Moreover, to optimize the initial stiffness of Mode II, K n was calculated from the DCB test; K s and K t were considered Fig. 6 The SLS test specifications in experimental (a and b) and numerical (c-f) studies. a The specimen dimensions based on ASTM D3163, b Performing SLS test, c Constraining all degree of free-dom at the end of the specimen, d Tensile displacement to simulate loading regime, e Using quadratic tetrahedral (C3D10) elements for meshing of Lattice, Shifted and Bulk designs equally. Furthermore, the mechanism of damage was simulated by defining the damage initiation criterion and the damage evolution law.
The damage initiation is the start of the cohesive response degradation at a contact point. In this paper, we have employed the maximum nominal stress criterion as a damage initiation criterion (Eq. 4).
where t 0 n , t 0 s and t 0 t are the peak values of the contact stress when the debonding is either purely normal to the interface or in the first or the second shear direction, respectively. As the contact stresses (t n , t s , t t ) reach certain criteria of damage initiation, the degradation process begins. We used trial and error to optimize the initial guess calculated by the Eq. 2. In the T-d curve, the area under the curve demonstrates the cohesive energy, and the maximum point in ascending line is the cohesive strength. First, two CZM parameters of initial stiffness and cohesive energy were kept constant and then by changing cohesive strength value, this parameter was calibrated to have appropriate fitted numerical curve to the experimental data. Moreover, to optimize the cohesive strength of Mode II, we calibrated T n from the DCB test and considered T s and T t equally. After reaching the initiation criterion, the cohesive stiffness is degraded base on damage evolution law. This law can be described by fracture energy (G), which is the energy dissipation because of the damage process. We assume the viscosity coefficient of 0.0002 for damage stabilization to facilitate solution convergence. Then, for obtaining cohesive energy parameter, we used trial and error again to optimize the initial guess calculated by the Eq. 3. So, we kept the initial stiffness and cohesive strength constant and similar to procedure for estimating the cohesive strength, the value of cohesive energy was calibrated. After calibrating the cohesive energy for mode I in DCB test, value of cohesive energy for mode II was calibrated in SLS test. In conclusion, an inverse analysis procedure was started by initial guesses of CZM parameters and continued until appropriate parameters to fit the experimental data were obtained. Besides, to distinguish each cohesive parameter, the other two parameters were kept constant in FEM simulations.

Morphology characterization
The morphology of specimens was evaluated using Scanning electron microscopy (SEM) (AIS2100, SERON Technology, South Korea). First, using a sputter-coater, the samples were gold-coated. Then, by evaluating SEM images, the strut diameter, and the pore size of the fabricated scaffolds were measured.

3D printing and characterization
The samples were 3D printed with specified parameters in Fig. 7, and SEM images were captured from top and side views. According to this image, two different scaffold designs have rectangular-shaped axial and lateral interconnected pores. The average strut diameter, pore size from top view for porous designs, and layer height were measured from SEM images and presented in Fig. 7. It is noticeable that in the printed models, there were more extruded materials in comparison to their computational models, specifically in contact regions between layers. It is more evident in the Bulk model in Fig. 7 from the side view compare to its CAD model in Fig. 3. In its printed sample, there were lateral contact areas between filaments, while these areas did not exist in the designed model.

DCB and SLS tests
The experimental and numerical processes to analyze the DCB and SLS tests for obtaining the cohesive parameters in Modes I and II are depicted in Figs. 8 and 9, respectively. Also, the representative load-displacement curves of the specimens and the predicted numerical results by variation of cohesive parameters are shown in these figures. In both DCB and SLS tests, the load increment is linear relatively with displacement increment until the crack initiation. Then, the load increased at a nonlinear decreasing rate during the crack growth. When it reaches a maximum value, unstable and unexpected crack growth leads to the ultimate failure. As shown in Fig. 8, the Bulk sample needed more crosshead displacement and fracture energy to complete the DCB test compared to porous designs. However, the trend of load increasing with displacement increment was pretty similar in all designs.
As shown in Fig. 9, the Bulk sample showed more stiffness than porous designs in the SLS test. However, there was no significant difference between total displacements in all designs in contrast to the DCB test. At the maximum loading level, we simulated the high distortion of the substrates using the FEM model (Fig. 9a, b, d, e), and the debonding propagation in adhesive layers interfaces as clearly illustrated in Fig. 9e and f agreed with experimental curves.
In general, the Bulk model showed more resistance to debonding than porous designs as its failure occurred in higher fracture energies. Besides, there were no significant differences between the mechanical behavior of Lattice and Shifted designs during DCB and SLS tests from load-displacement curves and cohesive parameters. The cohesive parameters in Modes I and II of three different designs extracted from DCB and SLS tests are represented in Table 1.
As it shows in Table 1, the initial stiffness of all designs was similar in Mode I but was around five times for Bulk in comparison to porous designs in this Mode. Moreover, the cohesive strength of the Bulk model in Modes I and II is slightly lower than the porous designs. However, the fracture energy in both Modes for the Bulk model is around four times of two others.

Cohesive parameters utilization
The comparison between the assumption of cohesive contact between layers of specimens and fully bonded contact using tie constraint in the software is illustrated in Fig. 10. The cohesive assumption improved the accuracy of the computational simulation for all designs, obviously in both DCB and SLS tests. So, it is essential to suppose partially bonded layers using CZM in additive manufactured parts in Modes I and II to simulate void and gaps between layers.

Discussion
In the present study, we calculated the cohesive parameters of bonding layers in Mode I and II by comparing the load-displacement curves from numerical simulations and experiments through inverse FEM. The samples were 3D printed with acceptable geometrical accuracy in terms of strut diameter, pore size, and layer height, similar to design assumptions.
There was a reasonable consistency between bilinear law and the experimental data, and they showed similar initial slopes, maximum loads, and trends. The numerical curves were in the domain of experimental curves indicated that the CZM could simulate the crack initiation and propagation of 3D printed scaffolds. In other words, the numerical predictions were in good agreement with experimental results.
Moreover, we had assessed the susceptibility of the simulated mechanical responses to cohesive parameters, i.e., initial stiffness (K), cohesive strength (T), and cohesive energy (G). The initial stiffness affects the slope of the force-displacement curve to reach the maximum strength before damage initiation. Moreover, the initial stiffness value should be large enough to prevent layer penetration but not too large to cause the numerical simulation failure Fig. 7 3D printing parameters, SEM images of Bulk, Lattice, and Shifted samples from top and side views, and morphological parameters of strut diameter, pore size, and layer height measured from SEM images to match experimental data because it may cause fictitious compliance before crack growth.
Furthermore, by employing different cohesive strengths and energies, the initial slope of the load-displacement curve remained unchanged. However, by energy increment, the maximum load increased, and it led to a higher displacement. Moreover, the cohesive strength is associated with the peak value of the load-displacement curve. The maximum load extremely increased as the cohesive strength increased, but the cohesive strength showed a small effect on the displacement related to the maximum load, in contrast to the cohesive energy.
In Lattice and Shifted designs, because of the equal number of filaments that can carry the load and so the same cohesive surfaces, cohesive parameters that are extracted by the DCB and SLS tests are approximately equal. It can be concluded that the microstructure may not affect the cohesive parameters of the porous design in the same printing conditions.
The disagreement between the curves of two scaffold designs (Simple and Shifted) and a Bulk model was the different designs of extruded filaments, which can directly influence the load-displacement of FDM parts. In Simple and Shifted designs, because of the equal number of Fig. 8 The Double cantilever beam (DCB) test steps (a-c) and its computational simulation (d-f). The load-displacement curves with the variation of cohesive parameters were compared with experimental data, and the best fitting was selected. In the right-bottom figure for each design, the numerical results with selected cohesive parameters and 5 DCB test results are illustrated Fig. 8 (continued) filaments that can carry tension load. However, in the Bulk model, the number of filaments that can yield tension load is 50% more than those two scaffold designs, so it can yield more.
Moreover, in the bulk model, there would be an interlayer adhesion between adjacent vertical layers. As a result, there would be an interfacial area between adjacent horizontal struts, resulting in high fracture resistance. Weak fracture resistance of two scaffold designs (Simple and Shifted) was due to the lack of connecting adjacent filaments too. This limits the physical crosslinking of the hard domains across the interlayer fused region. It is worth mentioning that high fracture resistance allows specimens to experience high load and displacement.
FEM can simulate the initiation, propagation of the crack, and failure of the DCB test. The crack initiation occurs at about 3 mm displacement in Bulk specimens, while at about 1.5 mm in Lattice and Shifted specimens, indicating that there was more resistance to debonding in the Bulk sample. Also, despite the similar values of initial stiffness and cohesive strength, the fracture energy was around four times for the Bulk model compared to porous designs that can be due to the difference between FEM and printed model. The FEM model was an ideal model with simplified geometry, and the excess material deposited in printed samples, especially in contact regions of struts, didn't exist in it. So, the contact region in the printed specimen enhanced in comparison to the FEM model for Bulk design, as showed in Fig. 7   Fig. 9 The Single lap shear (SLS) test steps (a-c) and its computational simulation (d-f). The load-displacement curves with the variation of cohesive parameters were compared with experimental data, and the best fitting was selected. In the rightbottm figure for each design, the numerical results with selected cohesive parameters and 5 SLS test results are illustrated (Bulk model, side view). However, the porosity of Lattice and Shifted microstructures can be influential on the level of fracture energy between 3 designs, and the non-porous microstructure of the Bulk model can increase this energy. Moreover, the stress distribution during initial bonding, crack initiation, and progression was simulated in different steps of the DCB test using FEM, which was not attainable using experimental approaches.
The Bulk model showed more resistance to layers debonding than porous designs as its initial stiffness was higher. In other words, because of high initial stiffness, the stress and consumed energy in the sample reached to the cohesive strength and energy faster than porous designs in lower displacement ranges. On the other hand, in the DCB test of Bulk model, the number of filaments that can yield tension were 50% more than those two scaffold designs, so it could yield more, and so higher displacement was observed in comparison to porous designs. Furthermore, higher displacement in DCB test in comparison with the SLS test of Bulk sample could be because of fiber bridging development. As debonded fiber formed bridges between the crack surfaces, these bridges could bear the load and let the specimen to experience higher displacement.
As we can see in Fig. 6, the FEM model of bulk specimen does not have any filament and have a smooth surface, but four layers has been modeled in porous designs. In reality, the 3D printed bulk model include some voids, discontinuity and gap between deposited filaments which did not exist in the numerical model. Because of this simplification due to reducing the time of simulation, the deviation has Table 1 The cohesive parameters: initial stiffness (K), cohesive strength (T), and fracture energy (G) in Modes I and II of Bulk, Lattice, and Shifted designs Bulk 1.0e 12 5.0e 10 2.0e 6 2.6e 7 3.9e 2 4.2e 2 Lattice 1.0e 12 9.5e 9 3.0e 6 3.0e 7 9.5e 1 1.0e 2 Shifted 1.0e 12 9.5e 9 3.0e 6 3.2e 7 9.0e 1 1.2e 2 Fig. 10 The effect of assuming partially bonded layers by the Cohesive zone model (CZM) to simulate interlayer contact for Bulk and porous additive manufacturing parts in comparison to fully bonded (tie constraint in Abaqus software) layers assumption. The modeling of voids and gaps between layers enhanced the computational accuracy using CZM in both Modes I and II occurred between numerical and experimental data. Moreover, the FEM model can model the excessive distortion of the substrates at the peak loading level produced by the huge eccentricity of the adhesive layers (Fig. 9b) in the SLS test simulation. Moreover, in agreement with experimental findings, this method can simulate the debonding propagation in the adhesive layers interface. In other words, the load-cross head displacement response of the SLS configuration is in excellent agreement with the experimental data and from Fig. 9, the bending moment has happened both experimentally and numerically, which is the sign of sliding behavior of SLS specimen. The initial stiffness for the Bulk model in the SLS test was around five times higher than porous designs that can be due to porosity or strong lateral bonding between struts of the Bulk model because of excess material deposition as observed in Fig. 7 (Bulk model, side view). Moreover, the fracture energy of the Bulk sample was around four times of porous designs in SLS tests like the DCB test. It can be because of the higher contact region in printed samples compared to the FEM model or non-porous microstructure as stated before for the DCB test. Furthermore, there was a softening region in the initial load-displacement of the Lattice and Shifted designs SLS experimental test, which did not exist in the numerical curve that can be due to simplification in the FEM model and the lower number of porous layers in comparison to actual 3D printed samples.
The tensile test of filament was repeated three times and there was no inconsistency in stress-strain curves. So, it can be concluded that, although glue was used to bonding the filament to trapezoidal tabs, no separation and slip was existed there. Moreover, according to the literature, the elastic modulus of PLA is 3500 MPa and we have calculated the elastic modulus of PLA as 3180 MPa approximately like that. However, the only reason for the difference is existence of some impurities in the composition of extruded PLA filament in comparison to pure PLA. Using glue as a fixator was not only available, cheap, easy to perform but also have a reliable result. So, that is the reason for using glue-dog bone specimen to conduct the filament tensile test to extract the elastic modulus.
It is noticeable that we proved that the assumption of cohesive contact between layers was essential for the FEM analysis of the 3D printed parts in Modes I and II, as seen in Fig. 10. In other words, there was a significant difference between full and partial bonding predictions to fit experimental data. Furthermore, this inverse FEM method can be further employed to predict the bond formation in 3D printed products, and the CZM parameters derivation under different manufacturing conditions because the printing parameters affect the properties of scaffold [5]. We note that the developed FEM models introduced in this study did not fully represent the actual microstructure of 3D printed samples, and it can affect strength and fracture properties [3,4,43]. More realistic models with exact porous layer numbers and interlayer contact regions can further improve the predictions and can be pursued in future studies.
In reality, there is no generation of excessive molding pressure and shear rate in the FDM process opposite to the relevant extrusion and injection thermoplastic processing routes. Hence, it is impossible to develop appropriate cohesion with strong interactions between layers. Indeed, in this printing method, the filament adhesion is mainly because of interdiffusion phenomena that can control the mechanical response. Accordingly, this novel technology facilitates the manufacture of parts with complicated geometry [1] but with moderate bonding properties that can be improved by postprocessing procedures, e.g., heating [6]. Therefore, using CZM in Mode I and II is essential to predict the mechanical properties because of this partial bonding between layers of FDM parts.

Conclusion
In this study, we measured the bonding strength of 3D printed PLA scaffolds using CZM and inverse FEM in Modes I and II, and investigated the effect of porosity and microstructure. The cohesive parameters of porous designs were the same in Modes I and II which means that these parameters can be used for other porous designs, printed with the PLA in similar printing conditions. Also, the computational simulations illustrated the stress distribution in the crack propagation, and it means that the unmeasurable quantities in experiments can be estimated using numerical models. Moreover, results showed that it is essential to consider the partial bonding between layers for simulating the mechanical properties of 3D printed samples, especially in fracture Modes I and II. The work demonstrated that the framework of CZM and FEM methodology could be implemented to estimate the adhesion between layers of the 3D printed parts, which lead to more realistic and accurate mechanical properties predictions. Funding This research did not receive any specific grant from funding agencies in the public, commercial, or not-for-profit sectors.
Data availability Yes, the data and material are available.
Code availability Yes, the code is available.