Numerical Investigation of the Cycling Loading Behavior of 3D-Printed Poly-Lactic Acid (PLA) Cylindrical Lightweight Samples during Compression Testing

Abstract

The additive technologies widely used in recent years provide enormous flexibility in the production of cellular structures. Material extrusion (MEX) technology has become very popular due to the increasing availability of relatively inexpensive desktop 3D printers and the capability of fabricating parts with complex geometries. Poly-lactic acid (PLA) is a biodegradable and commonly applied thermoplastic material in additive manufacturing (AM). In this study, using a simulation method based on the user subroutine titled “user subroutine to redefine field variables at a material point” (USDFLD) in the finite element method (FEM) ABAQUS software, the elastic stiffness (ES) of a cylindrical lightweight cellular PLA sample with a 2.4 mm infill line distance (ILD), which was designed as a layered structure similar to the laboratory mode with a MEX method and was subjected to cyclic compressive loading, was investigated by considering the variation of the Young’s modulus depending on the variation of the equivalent plastic strain (PEEQ). It was observed that the PLA sample’s elastic stiffness increases during cyclic loading. This increase is high in the initial cycles and less in the subsequent cycles. It was also observed that the simulation results are in good agreement with the experimental results.

1. Introduction

Although 3D printing technology was developed in the early 1980s and has been used in the automotive and aerospace industries for over three decades, it is relatively new and rapidly growing. Three-dimensional printing, also known as additive manufacturing (AM), is a technique in which three-dimensional objects are made layer-by-layer, based on computer models designed for the fusion of materials [1,2,3,4].

Nowadays, AM is extensively used in various industrial products [5,6]. Compared to other common manufacturing processes, such as injection molding, the essential benefit of additive manufacturing is the capability to manufacture more complex shapes [7,8,9]. This feature is due to the layer-by-layer build process. AM technologies play a vital role in production [10,11,12]. Particularly, one of their most substantial strengths is the production of small series parts with innovative geometries [13,14].

Rapid prototype (RP) technologies have steadily improved over the past three decades as they have proliferated. Some of the benefits of the RP process include the capability to produce geometries with increased complexity in very short times without the need for additional costs for tooling [15,16]; the ability to create functional sets by integrating sub-sets into a single unit in the computer-aided design phase (CAD), reducing the number of parts, time displacement, storage needs, and connection problems; the ability to optimize material consumption by manufacturing parts that are complex or even impossible to manufacture using conventional manufacturing processes; and finally the ability to reduce waste and thus minimize the impact on the environment [17].

Rapid Prototyping (RP) and rapid production technologies have moved from the first commercial stereolithography hardware (SLA) in the late 1980s to today’s desktop-sized molten sediment modeling systems. The above types of systems, and several others, mostly vary in the concept of additive production. Recently, several factors have led to an exponential growth in the use of such systems; the most profitable in the consumer market is inexpensive thermoplastic extrusion [18].

Material extrusion (MEX) is an AM technology that uses geometric models generated directly in a CAD system to create three-dimensional components of the data obtained directly. This process begins with the slicing of the three-dimensional model of the part in the CAD system to obtain two-dimensional curves (XY direction), which determines where in each layer of material should be added (Z direction) [19]. MEX produces 3D objects by quickly and accurately melting filament of thermoplastic polymers while extruding them through a nozzle and placing the melted material onto a high-precision geometry plate [1].

Generally, the MEX technique is especially favorable for the manufacture of small-volume parts and for prolonging their operational periods through approaches such as the repair or self-assembly of parts [20,21,22,23].

MEX is flexible to a high degree with the geometrical and technological parameters that affect a part’s mechanical properties and qualities. Therefore, Wichniarek et al. studied mechanical properties and behavior for defined groups of materials used for the MEX-produced parts and a corresponding set of process parameters ensuring stable manufacturing conditions [24].

Anitha et al. [25] evaluated the effect of process parameters on the modality characteristics of prototypes produced using MEX. Lee et al. [26] analyzed a set of MEX parameters, such as layer thickness, mesh angle, and air distance, and investigated their effect on the elastic performance of an acrylonitrile butadiene styrene (ABS) prototype. According to their results, the MEX parameters (raster angle, layer thickness, and air gap) notably affect the elastic performance of the compliant ABS prototype. Crocolo et al. [27] studied the impacts of melted sediment modeling production parameters on the tensile strength and stiffness of components made of ABS M30 and examined this question both empirically and analytically. The most important results of this research are the feasibility of predicting the mechanical behavior of the MEX parts, once the dimensions, number of contours, and raster angle have been stated.

The most widely used thermoplastic materials are polyethylene, polypropylene, ABS, poly-lactic acid (PLA), high compression polystyrene resin (HiPS), polyethylene terephthalate glycol (PETG), polyamide (PA), commonly known as nylon, and polyether ether ketone (PEEK), among others [19,28]. The most common source materials for commercially available MEX printers are ABS and PLA.

PLA is a bio-plastic composed of a repetitive lactic acid chain. It can be recycled through applying conventional methods [29]. Additionally, PLA can be composed similarly to all of the other organic materials. It was observed to be one of the most important biopolymers [30] with its substantial application in the field of medicine, such as for medical implants, bone fixation and reconstruction, (bone) tissue engineering, drug delivery system, and suture material [31]. Åkerlund’s [1] research developed biocompatible and bioresorbable polymer blends and composite filaments for MEX printing. Åkerlund used two biodegradable polymers (PLA (transparent filament) and polycaprolactone (PCL)); he used an Ultimaker S5 printer with a nozzle size of 0.4 mm (AA 0.4) for printing the cylinders with 100% infill density for a compression test, with the dimensions of 6 × 6 × 12 mm³. Afterward, compression tests were performed on dense cylinders by using a Shimadzu AGS-X materials tester to investigate the mechanical characteristics of both materials. Six to ten specimens for each type of material were used for the compression test and the stiffness (Young’s modulus) and yield point were reported. Torres et al. [18] presented the effect of major processing parameters on the properties of materials obtained from poly-lactic acid components modeled by MEX tested in torsion. The findings showed that the result improved with heat treatment, and the infill effect on ductility. The thickness of the layer and density of infill [32] have significant impacts on optimizing the strength of the part. Recently, for environmental reasons, biodegradable polymers such as PLA have been used instead of conventional polymers. PLA was used in biomedical applications, containers, packages, auto parts, and so on. PLA is a crystalline polymer and its Young’s modulus and tensile strength are 2 to 3 GPa and 50 to 70 MPa, respectively. PLA is very brittle, and its low toughness limits the use of PLA. To improve this fragility, PLA is combined with other flexible polymers. The biocompatibility of PLA was approved for food and medicine applications because of its good mechanical properties and processing performance [21,33]. PLA is a thermoplastic aliphatic polyester polymer derived from renewable sources. It has been considered as a potential environmentally friendly alternative to its oil counterparts. It is used for food contact, packaging, and scaffolding applications and is among the most widely used polymers, especially because of its ability to crystallize stress, heat crystallize, impact-modified, fill, copolymerize and process in most excellent polymer processing.

Designing and optimizing plastic parts remains a very complex engineering challenge. Designing a new product that is more secure, for example, involves validating its mechanical behavior by numerically simulating FEM before it is built. These simulations require accurate knowledge of the mechanical properties of the material. These properties are affected by the manufacturing process as well as the mechanical effort that affects the design in the operational state.

In the case of MEX-produced parts, dimensional tolerances seem to have impacts such as affecting the layer height, which affects the topology component of its specific position in the printing press. Additionally, the height of the layer and the rate of deposition are affected by the required product quality. Similarly, manufacturing characteristics are affected by the target cost, for example, by minimizing printer size and decreasing production support.

The accurate simulation of components produced with MEX means knowing in advance how the material behaves based on the manufacturing specifications and the final product requirements. In addition, to complete the validation process after obtaining the experimental information, it is necessary to know the best way to process the experimental information in FEM software so that the simulation results can be completely in accordance with the experimental behavior of those materials. Domingo-Spin et al. [34] proposed a model for simulating MEX components that links finite element analysis (FEA) simulation to physical experiments. They determined that choosing the appropriate direction to manufacture a part is important; moreover, the nozzle diameter, slice height, and the diameter of the extruded filaments can seriously affect the anisotropy of MEX-produced parts.

To assess the possible impact of compressive loading on the ES of the MEX-produced cellular PLA specimens with two different ILDs of 1.6 mm and 2.4 mm, Pepelnjak et al. [35] carried out research on cycling loading and analyzed the impact of such a loading on the part’s mechanical properties. To this aim, in the published research, each sample was subjected to 10 cycles of distinct compressive loading to a definite amount of displacement ranging from 0.2 to 0.5 mm (Figure 1).

The results of Pepelnjak et al. [35] demonstrated an enhancement of the ES of the sample with the ILD of 1.6 mm for 7.5% after the initial 10 cycles by compressive cycling loading with a starting nominal displacement of 0.2 mm, resulting in observing the most significant variation after the second loading cycle. Furthermore, a rise of 8.4% in the ES of the sample was revealed by comparing the impact with the mentioned samples with the 2.4 mm ILD. Regarding cycling loading from a starting displacement of 0.4 mm, in which the increase in the ES of the sample was 10.2% and 9.8% for an ILD of 1.6 mm and 2.4 mm, respectively, ES changes were shown to be more pronounced. Finally, the ES of the sample for 11% had its most significant increase after 20 loading cycles from an initial nominal movement of 0.4 mm at 1.6 mm ILD. This study aims to validate the experimental results of altering the elastic stiffness of the polymer structure by using numerical simulations. These validation results open new capabilities for the optimization of polymer lightweight parts. This way of using a deliberately deformed polymer structure is a new model not yet observed in the scientific literature; it is expected to be of help in engineering structures in lightweight materials. In the present numerical simulation, a model with 2.4 mm ILD is researched.

The results from the numerical simulation of the polymer structure produced by additive manufacturing closely match the experimental results. Therefore, it is possible to alter the elastic stiffness by loading–unloading. Figure 2 summarizes the research.

2. Materials and Methods

In the numerical simulation with ABAQUS software, several user subroutines have been suggested to meet specific requirements. These subroutines are programmed in FORTRAN and used in the ABAQUS software [36,37]. The User subroutine to redefine field variables at a material point (USDFLD) was used in this research.

The USDFLD subroutine can be used to define the behavior of materials, and the output of the required results can be used by utility routines GETVRM. The state variables (SDVs) and field variables (FVs) could be obtained at each integration point of the finite element by this subroutine. Equivalent plastic strain (PEEQ) was chosen as a field variable. The dependence between Young’s modulus of FFF-printed PLA material and the PEEQ should be defined. The values of Young’s modulus correspond to different amounts of the field variable. The elastic stiffness of the PLA model changed after the compression loading and unloading. The properties of the material are altered by the dependence defined. The flow chart is shown in Figure 3.

Many biodegradable polymers are available to produce a large diversity of plastic products, each with properties relevant to the application. A wide range of mechanical properties is possible among biodegradable polymers. Some metallurgical alloys and polymers present softening behavior instantly after reaching the yield point. In a 3D elasto-plastic formulation implementing minor strains [38,39], the strain tensor is provided by an elastic and a plastic tensor:

$$\epsilon ={\epsilon }^{\mathrm{e}}+{\epsilon }^{\mathrm{p}}$$

(1)

For disparate engineering materials, a yield criterion indicating the stress level at which plastic flow begins must be postulated. A relationship between stress and strain is obliged to be developed for post-yield behavior, i.e., when the deformation comprises both elastic and plastic components. The yield surface separates the plastic region from the elastic region. To describe the stress–strain relation after plastic deformation, a plastic constitutive tensor needs to be established.

The yield function can be demonstrated by

$$F\left({\sigma }_{\mathrm{ij}},k\right)=0$$

(2)

where the yield surface, $F\left({\sigma }_{\mathrm{ij}},k\right)$, depends on the magnitude of the load applied and of a hardening parameter k [40]. For simplicity, the yield function in Equation (2) can be rephrased in terms of the three main stresses,

$$F=\left({\sigma }_1, {\sigma }_2, {\sigma }_3, k\right)=0$$

(3)

By differentiating Equation (2) and considering a plasticity flow law,

$$dF=\frac{\partial F}{\partial {\sigma }_1}d{\sigma }_1+\frac{\partial F}{\partial {\sigma }_2}d{\sigma }_2+\frac{\partial F}{\partial {\sigma }_3}d{\sigma }_3+\frac{\partial F}{\partial k}dk=0$$

(4)

Introducing parameter A, where $A=-\frac{1}{d\lambda }\frac{\partial F}{\partial k}dk$ (A is defined as a hardening parameter and $d\lambda$ is the plastic strain multiplier), Equation (4) can be rewritten to a form,

$$\frac{\partial F}{\partial {\sigma }_1}d{\sigma }_1+\frac{\partial F}{\partial {\sigma }_2}d{\sigma }_2+\frac{\partial F}{\partial {\sigma }_3}d{\sigma }_3-Ad\lambda =0$$

(5)

Decomposing the strain increment into the sum of infinitesimal elastic and plastic strain increments $d{\epsilon }^{\mathrm{e}}$ and $d{\epsilon }^{\mathrm{p}}$, respectively $d\epsilon =d{\epsilon }^{\mathrm{e}}+d{\epsilon }^{\mathrm{p}}$, the strain increments can then be derived as Equation (6).

$$d\left\{\epsilon \right\}={\left[D\right]}^{-1}d\left\{\sigma \right\}+\left\{\frac{\partial F}{\partial {\sigma }_{ij}}\right\}d\lambda$$

(6)

where [D] is the elastic constitutive matrix that linearly relates the six components of the stress ({$\sigma$}) with the six components of strain ({$\epsilon$}) [41]. The entire elastoplastic incremental stress–strain relation can be derived to be

$$d\left\{\sigma \right\}={\left[D\right]}_{\mathrm{ep}}d\left\{\epsilon \right\}$$

(7)

where ${\left[D\right]}_{\mathrm{ep}}$ is the elastoplastic constitutive stiffness matrix provided by

$${\left[D\right]}_{\mathrm{ep}}=\left[D\right]-\frac{\left[D\right]\left\{a\right\}{\left\{a\right\}}^{\mathrm{T}}\left[D\right]}{A+{\left\{a\right\}}^{\mathrm{T}}\left[D\right]\left\{a\right\}}$$

(8)

where {a} is equal to {$\frac{\partial F}{\partial {\sigma }_1},\frac{\partial F}{\partial {\sigma }_2},\frac{\partial F}{\partial {\sigma }_3}$}.

2.1. Equivalent Strain Approach

The equivalent plastic strain equation can be defined in a manner consistent with the definition of the Von Mises equation [39,42]. The equivalent plastic strain (${\epsilon }_0^{\mathrm{p}}$) can be obtained by the following equation:

$${\epsilon }_0^{\mathrm{p}}=\frac{\sqrt{2}}{3}{\left[{\left({\epsilon }_{\mathrm{x}}^{\mathrm{p}}-{\epsilon }_{\mathrm{y}}^{\mathrm{p}}\right)}^2+{\left({\epsilon }_{\mathrm{y}}^{\mathrm{p}}-{\mathsf{\epsilon }}_{\mathrm{z}}^{\mathrm{p}}\right)}^2+{\left({\epsilon }_{\mathrm{z}}^{\mathrm{p}}-{\epsilon }_{\mathrm{x}}^{\mathrm{p}}\right)}^2+6{\left({\epsilon }_{\mathrm{xy}}^{\mathrm{p}}\right)}^2+6{\left({\epsilon }_{\mathrm{xz}}^{\mathrm{p}}\right)}^2+6{\left({\epsilon }_{\mathrm{yz}}^{\mathrm{p}}\right)}^2\right]}^{\frac{1}{2}}$$

(9)

where are ${\epsilon }_{\mathrm{x}}^{\mathrm{p}}$,${\epsilon }_{\mathrm{y}}^{\mathrm{p}}$,${\epsilon }_z^p$, ${\epsilon }_{\mathrm{xy}}^{\mathrm{p}}$, ${\epsilon }_{\mathrm{xz}}^{\mathrm{p}}$, and ${\epsilon }_{\mathrm{yz}}^{\mathrm{p}}$ the plastic strain tensor components [43].

2.2. Numerical Simulation

The 3D modeling and simulation were performed using ABAQUS finite element code to develop a layered (PLA) model that is subjected to compressive loading. The experimental sample was in an open-top cylinder form of D = 11.42 mm (diameter) and H = 10 mm (height) built with MEX technology. The MEX process parameters used in the sample production are defined in [35]. It was determined that it is not possible to print a 100% infill sample with that technology (Figure 4). In order to decrease computational costs and save time, the numerical model used in this study was one-twenty-fourth of the entire open-top cylinder. Therefore, first, the cross-section of the sample was divided into twelve equal and symmetrical parts, and the height was halved (Figure 5).

The influence of model shape on numerical simulation results was examined. Therefore, a simplified geometrical model (Figure 6) and a layered model, following the shape of the printed sample and the cross-section of the real shape (Figure 7), were created.

The mechanical properties of the PLA used in this simulation are provided in

Table 1. Mechanical properties of the PLA applied in the numerical simulation.

Property	Unit	Value
Density (ρ)	kg/m3	1240 [5]
Initial Young’s modulus (E) (altered with PEEQ during the simulation)	MPa	2800 [35]
Poisson’s ratio (υ)	-	0.36
Yield strength (σy)	MPa	59

and Figure 8 [1].

Therefore, to achieve the plastic properties of the PLA used in the numerical simulation, the plastic properties of Åkerlund [1] were chosen as a reference. After analyzing the results of numerical simulation, the stress and strain values of the plastic properties used by Åkerlund [1] were shifted together to some different percentages in different steps. By examining the results of numerical simulation after each step of making changes according to Figure 8, the data of diagram 5 (simultaneous 22% decrease of the stress and strain values) in this numerical simulation were used.

It was not possible to print the sample with one hundred percent infill using MEX technology; therefore, references for the material (PLA) data were sought for use in the numerical simulation. The Z-direction fixed boundary condition (normal direction displacement = 0) was applied at the middle cross-section, and the normal direction fixed boundary condition (normal direction displacement = 0) was applied at both side cross-sections of the model in the local coordinate system, as shown in Figure 7. The numerical model of the specimen is composed of 308,734 C3D8R elements (linear 8 node hexahedra with reduced integration) (Figure 9). The model consists of 26 parts, of which 25 are PLA layers (each separate object) stacked on top of each other; the other part is the tool (Figure 5). The Poisson’s ratio was 0.36 [44]. In this model, Young’s modulus changes with the equivalent plastic strain, which was performed using the USDFLD subroutine. In this way, by changing the equivalent plastic strain in the model, a new Young’s modulus was assigned to the material. To validate with laboratory results, loading was applied in two conditions of three and five cycles. For this purpose, the model was compressed by 0.1 mm; next, the applied load was removed to return the model to its original state. This procedure was repeated in two conditions of three and five cycles.

3. Results and Discussion

3.1. Comparison of Numerical Simulation and Experimental Results

PLA-based cellular samples were fabricated using MEX 3D printing technology and their mechanical properties investigated. During the experimental work, cyclic compressive loading was applied to PLA samples, and the ES was altered through combined elastic and plastic deformation. The outcome of the experimental work showed that, under cyclic compressive loading, the sample alters its elastic stiffness. In the present study, the first five cycles of the experimental procedure were simulated via ABAQUS finite element code. To obtain the elastic stiffness results from numerical simulation, the true stress–strain diagram was plotted with the force–displacement diagram data, and the elastic stiffness was calculated from the slope of this diagram. The simulation results show that the elastic stiffness of the model changes, which is consistent with the results obtained from previous experimental research.

Table 2. Results of the ES of the sample in experiments and model in the numerical simulation, induced by cyclic compressive loading.

Number of Cycles	1	2	3	4	5
Elastic stiffness in the experiment [35] (MPa)	2800	2970	2998	3005	3015
Elastic stiffness in the numerical simulation (MPa)	2732.8	2936.8	3005.0	3036.5	3057.5

and

Table 3. Summary of the ES percentage changes in the experiments and numerical simulation, induced by cyclic compressive loading.

Cycles Comparison	2 → 1	3 → 1	4 → 1	5 → 1
Percentage of elastic stiffness changes in experiment compared to the first cycle (%)	6.1	7.1	7.3	7.7
Percentage of elastic stiffness changes in numerical simulation compared to the first cycle (%)	7.5	10	11.1	11.9
Percentage difference between the experiment and numerical simulation results (%)	1.4	2.9	3.8	4.2

summarize the results obtained from both studies (experimental and numerical simulation). It can be seen that the results obtained from the numerical simulation and experimental research are in good agreement, and the percentage of discrepancies in results is reasonably low. This outcome can be considered a validation of both parts of the research. Comparing the results of each cycle, it is clear that the first cycle’s elastic response varied significantly from the subsequent ones. In both studies, there is a significant increase in the elastic stiffness in the first cycle, but in the following cycles, the amount of alteration decreases.

The results of the experiment and numerical simulation research in each cycle are presented in Figure 10. As is evident, the numerical simulation results and experimental results are significantly close.

In addition to the above observations, the hysteresis loop shown in the force–displacement diagram of the numerical simulation (Figure 11) is also fully consistent with the experimental results. The hysteresis loop demonstrates the mechanism of cyclic hardening within the model.

3.2. Model Shape Influence on Simulation Results

Based on the findings from

Table 4. Summary of the ES of the geometrical simplified model changes and the layer model changes in the numerical simulation, induced by cyclic compressive loading.

Number of Cycles	1	2	3	4	5
Elastic stiffness in the geometrical simplified model (MPa)	2818.2	2943.9	2989.1	3018.2	3039.8
Elastic stiffness in the layer-by-layer model (MPa)	2732.8	2936.8	3005.0	3036.5	3057.5
Elastic stiffness in the experiment [35] (MPa)	2800	2970	2998	3005	3015

, elastic stiffness changed in the simplified geometrical model and the layered model, and its value increased. When comparing these results, data from both models (geometrical simplified model and layer-by-layer model) show a difference between experimental and numerical simulation results of less than 5%. In the present research, the layered model was chosen because it follows better the real shape of the experimental MEX-produced sample and the real cross-section of the specimen.

It is also seen from Figure 12 that the force–displacement diagrams have a slight difference in the layered and geometrical simplified models. Based on Table 4 and Figure 12, it can be concluded that the shape of the model (layered and geometrical simplified) does not significantly affect the results.

3.3. Simulation Results along a Specific Path

The exact position of the specific paths can be seen in Figure 13. In the research, the results alongside the two paths during the numerical simulation were investigated. These paths were chosen in the position that included all layers of the model.

Figure 14 shows the Mises stress alongside Paths 1 and 2 when each cycle is finished. As can be seen, at the end of each cycle, the stress values increase from top to bottom of the model alongside the paths on the primitive layers, and then start to decrease. The maximum stress is almost at the middle layers of Paths 1 and 2.

Figure 15 shows the amount of the elastic strain (maximum principle) at the end of each cycle, alongside Paths 1 and 2. The results indicate that the behavior of the model is approximately similar during both paths. As seen in Figure 15, alongside the paths, the maximum elastic strain value obtained is about 0.011, and the elastic strain value alongside the paths first increases and then starts to decrease.

Figure 16 shows the equivalent plastic strain in Path 1, specified in Figure 13. As is evident from Figure 16, the equivalent plastic strain values increase from top to bottom of the model alongside Path 1. The equivalent plastic strain has its highest value in the last layers and has the lowest value in the initial layers.

Figure 17 shows the equivalent plastic strain in Path 2, specified in Figure 13. As seen in Figure 17, the variation in the equivalent plastic strain values is smaller at Path 2 compared with Path 1, but both paths have almost identical ways of behaving. The equivalent plastic strain has its highest value in the last layers and has the lowest value in the initial layers.

Figure 18 shows the elastic strain obtained in each region of the model after five cycles of cyclic compressive loading. As can be seen, the maximum elastic strain occurs between the infill lines layers.

4. Conclusions

In this study, a simulation of a PLA model under cyclic compressive loading–unloading was performed, and the results were compared to the experimental results. At first, the effect of shape on the results and whether the modeling method is important were examined. In this regard, two simulations were performed; in the first one, the model was layered, which follows better the contour of the MEX-produced sample and the cross-section of the real shape; in the second one, the model was geometrically simplified. The results show that the model’s shape is also effective in the results, and the results of the layered models were different from the geometrically simplified ones. The elastic stiffness obtained in the full solid model from Cycle 1 to Cycle 5 was 2818.2 MPa, 2943.9 MPa, 2989.1 MPa, 3018.2 MPa, and 3039.8 MPa, respectively, while in the layer model, these numbers were 2732.8 MPa, 2936.8 MPa, 3005 MPa, 3036.5 MPa, and 3057.5 MPa, respectively. Therefore, a layer-by-layer model was used according to the experimental sample shape.

In contrast, by comparing the results of the layered model with the experimental results, it was observed that the percentage differences of changing elastic stiffness in experimental work and simulation during five cycles of loading–unloading were really low. The percentages of difference in elastic stiffness between the experimental results and numerical simulation results from cycle to cycle were 1.4, 2.9, 3.8, and 4.2 percent, respectively. The numerical simulation outcomes show that the elastic stiffness significantly increases during the first cycle, but in the following cycles, the amount of increase in elastic stiffness decreases, which is in agreement with the experimental results. The force–displacement diagrams are also in acceptable agreement with the experimental results. It is evident that plastic deformations in the sample increase from cycle to cycle and influences its elastic stiffness, which also increases. Plastic deformations developed in narrow transitions between layers. This is why the geometrical model of the sample, used in mechanical analysis, should also take into account the approximation of the real layer’s cross-section shape.

Finally, some of the results that were not obtained in experimental work were achieved in the numerical simulation, such as the Mises stress, elastic strain, and PEEQ diagrams alongside the two paths. It was observed that the Mises stress has the lowest value in the first layers (from the tool side), and the stress values increase from top to bottom of the model alongside the paths on the primitive layers, and then it will start to decrease. The maximum amount of stress is almost at the middle layers of Path 1 and Path 2. In addition, the full amount of the Mises stress was observed between the layers. The results of all observations on the numerical simulation are for investigation in future works.