Characterization and Simulation of Shear-Induced Damage in Selective-Laser-Sintered Polyamide 12

This paper presents the characterisation of selective-laser-sintered (SLS) samples of polyamide 12 (PA12) under shear loading. PA12 is a semi-crystalline thermoplastic and is used in various industries. Its behaviour under shear stress, which is particularly important for product reliability, has not yet been sufficiently investigated. This research focuses on understanding the material and damage behaviour of PA12 under shear-induced stress conditions. The study included quasi-static experiments and numerical simulations. Samples were prepared via SLS and tested according to ASTM standards. Digital image correlation (DIC) was used for precise deformation measurements. The Chaboche material model was used for the viscoplastic behaviour in the numerical simulations. Due to existing material discontinuities in the form of voids, the material model was coupled with the Gurson–Tvergaard–Needleman (GTN) damage model. A modified approach of the GTN model was used to account for low stress triaxiality under shear loading. These models were implemented in MATLAB and integrated into Abaqus via a User Material (UMAT) subroutine. The results of the experiments and simulations showed a high degree of accuracy. An important finding was the significant influence of the shear factor kw on the damage behaviour, especially during failure. This factor proved to be essential for the accurate prediction of material behaviour under shear-induced stress conditions. The integration of the modified GTN model with the Chaboche material model in UMAT enables an accurate prediction of the material and damage behaviour and thus makes an important contribution to the understanding of the mechanical material behaviour of SLS PA12 specimens.


Introduction
The advancement of technology necessitates the creation of innovative materials.Polymers are becoming increasingly important as a construction material with additive manufacturing (AM) methods being utilised.By using AM technology, a digital computeraided design (CAD) is converted into real objects, enabling a rapid transition from design to end product.This is achieved by building up the objects layer by layer, without the need for moulding or mechanical processing [1].During their service life, 3D-printed parts produced in this way must be able to withstand various mechanical loads.Knowing the required strength under different loads is essential [1,2].A frequently used AM process for polymers is powder bed sintering, also known as SLS [3].In this process, one or more lasers selectively fuse powder particles on the surface, layer by layer [4][5][6].PA12 is one of the most commonly used materials for the SLS process.This semi-crystalline thermoplastic is ideal for functional prototypes or end-use parts due to its mechanical properties, flexibility, and heat resistance [3,7,8].The adaptability of SLS PA12 is demonstrated by its successful integration in a wide range of industries, from the manufacture of high-performance ball valves and medical devices to customised automotive parts and beyond [3,[9][10][11][12][13].In order to utilise the full potential of SLS PA12, a thorough understanding of the material properties under a wide range of loading conditions is essential.In recent years, studies have been conducted on the material behaviour of SLS PA12 under both quasi-static and cyclic loading [14][15][16][17][18]. Due to the lack of comprehensive standards in the field of AM, the procurement and interpretation of material data remain difficult [19].One particular area to be investigated is the behaviour of SLS PA12 under shear-induced stress conditions, an area that is relatively under-researched [20][21][22].The challenge is not only to identify the details of the material behaviour, but also to determine the most suitable test geometries that can provide reliable and reproducible results, especially in shear-dominated scenarios [23][24][25][26].Even though the current state of the art in SLS is advanced, it also has its limitations.The process can lead to structural inhomogeneities in the final product in the form of voids, microcracks, or even unsintered areas [16,[27][28][29][30].These microstructural features in combination with the described complex material properties pose challenges for numerical modelling.In the work of [15,16], the material model of Chaboche [31] was coupled with the damage model of Gurson-Tvergaard-Needleman [32], and a high accuracy between experiment and simulation was achieved.The GTN model does not take shear stress into account.In this respect, attempts have been made in recent years to modify this model, following the work of Xue [33] and Nahshon and Hutchinson [34], who presented a framework that potentially provides higher accuracy under combined tensile and shear loading.
The present work has three objectives.The first one is the experimental characterisation of the material and damage behaviour of SLS PA12 under shear load.The second aim of the study is to couple Chaboche's material model with the modified GTN damage model based on Nahshon and Hutschinson's approach to ensure that it captures the details of shear-induced behaviour with optimal precision.Due to the limited data available for SLS PA12 in common FE software (ABAQUS 2022) material libraries, the third objective of this work is to implement the material and damage model in a UMAT.

Sample Preparation
The specimens were fabricated utilizing SLS technique.A PA2200 powder was deployed for the process.A sPro 230 3D printer, manufactured by 3D Systems and functioning at a power output of 70 W, was utilized in the fabrication.The printer operated at a laser scanning speed of 10 m/s for the infill and 5 m/s for the contour delineation.The nominal layer thickness was maintained in the range of 0.08 to 0.15 mm.The process and the chamber temperature during the print cycle were sustained at 200 °C and 170 °C, respectively.
Given the current lack of standardized testing procedures for AM polymers under shear loading [19], the testing geometry was chosen in line with the ASTM B831 standard [35], as illustrated in Figure 1a.Based on prior research, it has been observed that the resultant material properties demonstrate only minimal dependence on the printing direction [15].Therefore, all subsequent tests were performed on specimens layered along the z-axis, Figure 1b.

Digital Image Correlation
For the execution of DIC, the ARAMIS 3D Camera measurement system was deployed for comprehensive area-based and discrete point measurements.DIC is an optical noncontact method, allowing for full-field measurement of deformations.The technique relies on a reference image to serve as a baseline for deformation measurements.Consequently, a unique binary (black and white) pattern is applied to the specimen surface (Figure 1c).These stochastic patterns facilitate the identification of discrete image areas, with subpixel accuracy, through the analysis of image-based information, Figure 2. The methodology generates point and area measurement results which are subsequently used in material research and component testing to evaluate the static behaviour of specimens [36].The fundamental correlation function is represented as [37]: where C is a function of the coordinates x and y of the reference image.The displacement and disparity are defined as u and v, respectively.I represents the pre-deformation image and is a function of the pixel values x + i and y + j.The post-correlation, I * , signifies the image function of the pixel values after the application of displacements.The function C(x, y, u, v) represents the cumulative squared differences between the reference and the deformed images, where the function is summed over the subset size N. Using the Hegewald & Peschke Inspekt-Table 10 kN testing machine (Nossen, Germany), five shear tests were conducted at a 1 mm/min displacement rate and a temperature of 23 °C.

Chaboche Model
Previous investigations have shown that SLS PA12 behaves viscoplastically [15].This study focuses on the development of material-specific adjustments of the Chaboche material model, with limitations set on purely isothermal conditions.The Chaboche model is constituted of a system of differential equations, originating from the division of the total strain rate εtot into an elastic εel and a viscoplastic εvp component [38]: with the stress σ and Young's Modulus E. Herein, the definition of the viscoplastic strain rate εvp , composed of the viscoplastic multiplier λ and the direction ∂ f /∂σ, is decisive: The index (•) v indicates that it is an equivalent stress and (•) H hydrostatic stress.The relative stress β is defined by [39] The flow function f delineates the closure of the elastic domain and X k the non-linear kinematic hardening.Equation ( 6) defines the flow function of viscoplastic behaviour [32]: The parameters q 1 and q 2 are damage parameters and f * is the void volume fraction.The flow function is dependent on the non-linear isotropic hardening σ M [38]: with the accumulated plastic strain rate εvp , where its derivative is the modulus of plastic strain rate: εvp = σ : εvp The isotropic hardening σ M describes the extension of the flow surface and is defined by the yield stress σ y , the viscosity factor K, the viscosity exponent n, and the hardening parameters Q ∞ and b.In addition to the isotropic hardening, the flow function is influenced by the non-linear kinematic hardening X.The kinematic hardening depends on the strain hardening parameter C, the dynamic recovery factor γ, the static recovery factor R kin , the decay constant ω, and the saturation parameter φ ∞ .The first term in Equation ( 9) describes the strain hardening and the second and third ones describe the dynamic and static recovery, respectively [40]: with:

GTN Damage Model
The results of the micromechanical analysis from a previous investigation [15] were taken into account for the selection of the damage model.The virgin samples have a large portion of nearly spherical voids.Due to the quasi-static and cyclic loading, the voids are growing.When the spaces between the voids become small enough coalescence occurs.This coalescence needs to be embedded in the GTN model [32] so that it takes into account pore growth, nucleation, and coalescence.
In Equation (11), the effective void volume fraction f * is presented.This models the correlation between the reduction in the material's load-bearing ability and coalescence [41].
The initial porosity, f0 , is set at t = 0.As the void volume approaches the critical limit f c , coalescence begins.This continues until the material reaches the failure void volume fraction f F and subsequently fails.To describe the void volume fraction f , modifications were made to the original GTN model.Nahshon and Hutchinson's [34] approach was utilized, which considers the impact of shear load on damage behaviour.As a result, the shear term fShear was added [34]: so the void volume fraction f , is characterized by the sum of the initial void volume f0 , void growth fGrowth , void nucleation fnucleation , and void volume fraction due to shear fShear .The increment in the void volume fraction due to pore growth ˙f Growth is described by the effective void volume fraction f * and the viscoplastic strain rate εvp [32]: Assuming a strain-controlled load, the nucleation of pores follows a normal Gaussian distribution.This has a mean value ε N and a standard deviation S N .The new pore volume created via nucleation corresponds to the parameter f N .The nucleation is strongly dependent on the matrix material [32]: Nahshon and Hutchinson formulated an equation to predict the increase in void volume fraction at low stress triaxiality, expressed as follows [34]: with the shear factor k ω , the deviatoric stress tensor S, and the stress state function [34], with the third invariant of the deviatoric stress tensor J 3 = |S|.
All the described parameters of Equations ( 2)-( 16) were adjusted using the pattern search algorithm [42] and subsequently incorporated into Abaqus for the execution of the numerical simulation.

UMAT in Abaqus
Initially, the coupled Chaboche-GTN model was implemented in Matlab for one point.The return mapping algorithm is an iterative method used in plasticity theory [43,44].For each incremental step of the simulation, an elastic prediction is made first.Then, it is checked whether the resulting stress meets the yield condition.If the yield condition has been reached, it is assumed that a plastic deformation has occurred and a correction must be made.This correction, namely the determination of the corresponding plastic deformation increment to "return" the stress to the yield surface, is performed by the return mapping algorithm [43].
The model parameters of the coupled Chaboche-GTN model were determined using a pattern search algorithm [42] in MATLAB together with the experimental data from the shear tests.This process ensured that the simulation results accurately reflected the actual experimental results.The main goal of applying the pattern search algorithm was to fit the model as closely as possible to the actual physical phenomena in order to improve its accuracy and reliability.This step is crucial for validating the effectiveness of the model and its ability to accurately replicate real physical processes.Table 1 lists the parameters determined.The modulus of elasticity E, the yield stress σ y , the initial void volume fraction f0 , and the void volume fraction at failure f F were determined directly from the experiments.All other parameters were determined by optimisation.For the integration of the model into the finite element analysis software Abaqus 2022, the UMAT function was utilized.UMAT is a subroutine within Abaqus 2022 that allows users to define and implement custom material models in FORTRAN.This function provides the flexibility to model complex material behaviours that go beyond the models typically provided in Abaqus.A schematic overview of the process flow and the implementation of the applied model is depicted in Figure 3.The parameters derived from Matlab were used as input variables for the finite element analysis.For the finite element simulation, a reference model for the "simple shear" load case was created to cross-check the results of the UMAT in MATLAB.The C3D8 element was chosen for this purpose, a three-dimensional element with eight nodes and a linear shape function.With the support of constraints and the "Equations" function, these points were linked to the respective surfaces.Specifically, RP-1 was associated with the displayed surface III in the x-direction, and RP-2 in the y-direction for surface I.In the load case, the movement of RP-1 and its associated surface was fixed in the x-direction, while RP-2 underwent a defined displacement in the y-direction.For surface VI, the displacement in both the z and y directions was set to zero.Similarly, for surface IV, the displacement in the x and z directions was also set to zero.To simulate "simple shear", two reference points, RP-1 and RP-2, were introduced, as shown in Figure 4a.In addition to the reference model, a half-model of the original geometry was used, as shown in Figure 1.This half-model is depicted in Figure 4b.The boundary conditions are defined as follows: The displacement is constrained in the x-direction on Surface I, in the y-direction on Surface and in the z-direction on Surface V, Figure 4c.A displacement in the y-direction is prescribed on Surface III.For this model, the C3D8 element type was taken.

Experimental Results
The results of the five experiments are visualised in the shear stress-time diagram in Figure 5.The mean value x of the shear modulus and ultimate strength determined from these five shear tests and the corresponding standard deviations s are listed in detail in Table 2.The standard deviation values, 10.34 MPa for shear modulus and 2.51 MPa for ultimate strength, suggest a consistent performance of the material across different samples.Test 4 was used for the further fitting of the model parameters.

Simulation Results
Prior to the analysis, a thorough verification process was conducted to ensure the integrity of the transferred code from MATLAB to UMAT.After reviewing the model in MATLAB and integrating it into the UMAT, the FE model for simple shear (Figure 4) was analysed.In Figure 6, the deformed body and the associated shear stress σ 12 are depicted.A strong correlation was observed between the results from MATLAB (red line) and the UMAT (green dotted line) in comparison to the experimental data (black line), Figure 7.Both the material behaviour and the time of failure could be simulated exactly.As expected, an increased shear stress σ 12 was observed in the region of the webs (Figure 8).The analysis of the shear stress σ 12 at a node (black point) in the shear zone and its comparison to the experiment are illustrated in Figure 9.The material and damage behaviour can be mapped with a very high accuracy.The maximum deviation is 1 MPa and occurs at the point in time when the shear and tensile loads superimpose each other.Furthermore, the influence of the newly introduced shear factor k w (Equation ( 15)) was examined (Figure 10a).The value for k w determined through optimization is 0.12 (Figure 10a, dotted red line).To analyse its influence, this value was increased to 0.18 and decreased to 0. The results are depicted in the stress-time diagram, both throughout the entire course (Figure 10a) and in detail shortly before failure (Figure 10b).From Figure 10b, it becomes evident that k w has a significant influence on the characterization of the damage behaviour, especially shortly before failure.When the value of k w decreased, the stress is overestimated, while it is underestimated when increased.Misestimating this value can lead to inaccurate stress and failure predictions.Therefore, it is essential to calibrate k w correctly to ensure that model predictions align with the real damage behaviour.
In order to investigate whether pure shear or a superposition of shear and tension is present in the component, the calculation of the stress triaxiality η, provides a statement [45,46].The stress triaxiality indicates the ratio of hydrostatic stress σ H to the Mises equivalent stress σ v .A high value of the stress triaxiality η = 0.33 indicates that pure tensile or compressive stress predominates.Low values indicate pure shear stress.Figure 11 illustrates the stress triaxiality in the same node as the analysis of the shear stress (Figure 8).At the beginning of the loading, there is pure shear stress.As the load progresses, the stress triaxiality increases to a value of η = 0.17 and a combination of shear and tensile stresses occurs, with shear stress remaining the predominant factor.The results presented in Figure 9 demonstrate that the enhanced GTN model is capable of accurately representing both pure shear stresses and the combined loading scenario of shear and tension.The programmed material and damage model in UMAT allows for the representation of individual development variables.In Figure 12, the void volume fraction at the time t = 40 s is depicted as f * .The complete progression of the void volume fraction over time is shown in Figure 12.Until the yield point is reached, the void volume fraction retains its initial value f0 .As soon as the yield point is exceeded, the void volume fraction gradually increases.When the coalescence threshold fc is reached, the increase in the void volume fraction accelerates significantly.In Figure 13, the results from Abaqus and the DIC are compared.The focus is on the comparison of the deformation fields, especially the strains in the x and y directions.The evaluation is conducted for t = 40 s.An excellent agreement between the experiment and the simulation is evident.

Conclusions
This study investigated the material and damage behaviour of SLS PA12 under shear loading and the application of Chaboche's material model with the modified GTN model.The main conclusions are the following: 1.
This study provides fundamental insights into the behaviour of SLS PA12 under shear loading, which are relevant for a wide range of applications in materials science and 3D printing.

2.
The application of the Chaboche material model in combination with the modified GTN model shows that complex material behaviours, such as those of SLS PA12 under shear loading, can be successfully simulated, improving the prediction accuracy in similar material studies.

3.
The results of this study extend the understanding of the damage behaviour of 3Dprinted materials and provide a valuable contribution to the further development of reliable simulation models.
In summary, the precise determination of the shear factor k w enables an accurate simulation of the material and damage behaviour of SLS PA12, which opens perspectives for improved design and testing methods in additive manufacturing.

Figure 1 .
Figure 1.ASTM shear specimen: (a) geometry, (b) single printing layers are in the xy-plane, the build-up direction of the layers is the z-direction; and (c) black and white pattern for DIC.

Figure 2 .
Figure 2. Reference image of DIC in Aramis.

Figure 4 .
(a) Boundary condition simple shear model (b), mesh of half model, (c) boundary condition of the half model.

Figure 5 .
Figure 5. Experimental results of shear stress over time.

Figure 9 .
Figure 9. Shear stress over time of the node of full model.

Figure 10 .
(a) Shear stress full model of UMAT and (b) cutout of stress-time curve under the influence of k w .

Figure 11 .
Figure 11.Stress triaxiality of full model of UMAT.

Figure 12 .
(a) Void volume fraction of full model of UMAT, t = 40 s and (b) development of void volume fraction over time.

Figure 13 .
Figure 13.Comparison of the strain field of Abaqus and DIC in y-direction and x-direction.

Table 2 .
Mean values and standard deviation.