Predicting the Bending of 3D Printed Hyperelastic Polymer Components

The advancement of 3D printing has led to its widespread use. NinjaFlex®, a thermoplastic polyurethane (TPU) filament, is a highly durable and flexible material that has been used to create flexible parts. While this material has been available for nearly two decades, the mechanical properties of 3D printed NinjaFlex® parts are not well-understood, especially in bending. The focus of this research was predicting the behavior of small 3D printed NinjaFlex® components. Three-dimensionally printed rectangular specimens of varying lengths and aspect ratios were loaded as cantilevers. The deflection of these specimens was measured using a computer. The experimental results were compared to a modified form of the Euler–Bernoulli Beam Theorem (MEB), which was developed to account for nonlinearities associated with large deflection. Additionally, experimental results were compared to the finite element analysis (FEA). The results showed that both modeling approaches were overall accurate, with the average difference between experimental deflection data and MEB predictions ranging from 0.6% to 3.0%, while the FEA predictions ranged from 0.4% to 2.4%. In the case of the most flexible specimens, MEB underestimated the experimental results, while FEA led to higher retraction.


Introduction
Scientific and technological advancements have allowed using polymer synthesis to develop materials for many applications. An example of these new materials is the hyperelastic thermoplastic polyurethanes (HTPUs). HTPUs offer the manufacturing advantage of being processed as a thermoplastic while maintaining properties similar to rubber due to their molecular structure, which contains two specific regions: soft and hard. The soft regions allow for high flexibility, while the hard ones contribute to permanent deformation, tensile strength, and hysteresis of the material [1,2]. This combination provides HTPUs with a composite-like behavior that is hyperelastic and non-linear [2].
Simultaneously, hardware development has allowed 3D printing to speed up the process of developing prototypes and part production, especially when a relatively small numbers of components are needed. However, incorporating 3D printing in the industry has been limited by changes that the materials, especially polymers, experience during material deposition. These changes affect the mechanical behavior as compared to traditionally manufactured materials because 3D printing introduces a plethora of factors including print orientation, infill raster angle, infill density, and shell thickness. In the following, a brief overview of research in this area is presented.
Some mechanical properties of hyperelastic polymers can be readily determined using uniaxial tensile testing. A brief overview follows of relevant research in this area. By 3D printing dog-bone shaped specimens from Elastosil M4061 and translucent soft silicones, stress-strain curves as well as elastic moduli were obtained [3]. The effects of infill density on uniaxial tensile strength of 3D printed Polylactic acid (PLA) specimens were considered [4]. Increasing the infill density consistently increased the modulus of elasticity,

Materials and Preparation
The manufacturer of NinjaFlex ® has provided several key mechanical properties that were found following ASTM D638 standards [13] (Table 1). However, other researchers have provided different values for these characteristics. This variation may be due to the 3D printer used as well as the printing parameters. Additionally, some researchers have provided different stress-strain curves for NinjaFlex ® [8,[10][11][12]14]. In this research, experiments were conducted using 3D printed NinjaFlex ® dog-bone specimens with varying rectangular cross sections (Figure 1a). The length (l), width (w), and height (h) were combined to create eighteen unique specimens, Table A1. The diameters of the two cylindrical ends, d, depended upon the height of the specimen. Each of the eighteen specimens were printed three times for a total of fifty-four samples. In this research, experiments were conducted using 3D printed NinjaFlex ® dog-bone specimens with varying rectangular cross sections (Figure 1a). The length (l), width (w), and height (h) were combined to create eighteen unique specimens, Table A1. The diameters of the two cylindrical ends, d, depended upon the height of the specimen. Each of the eighteen specimens were printed three times for a total of fifty-four samples.
The CAD models of the specimens were converted to .stl files using the Cura Lulzbot ® Edition (model slicing software, Fargo Additive Manufacturing Equipment 3D, LLC, Fargo, ND, USA) [16] (Figure 1b). This software has a library of print parameters for numerous filament materials including NinjaFlex ® . Print quality was set to standard; print parameters are listed in Table 2. The number of wall lines for each specimen varied depending on the specimen height, h, as shown in Table A1. All specimens were printed on the same machine with the same black NinjaFlex ® filament using a Lulzbot Taz Mini 2 ® with an Aerostruder ® tool head (3D printer, Fargo Additive Manufacturing Equipment 3D, LLC, Fargo, ND, USA) [17], with an extruder diameter of 0.5 mm. A single layer of painters' tape was used on the heated print plate for print adhesion (Figure 1c). The CAD models of the specimens were converted to .stl files using the Cura Lulzbot ® Edition (model slicing software, Fargo Additive Manufacturing Equipment 3D, LLC, Fargo, ND, USA) [16] (Figure 1b). This software has a library of print parameters for numerous filament materials including NinjaFlex ® . Print quality was set to standard; print parameters are listed in Table 2. The number of wall lines for each specimen varied depending on the specimen height, h, as shown in Table A1. All specimens were printed on the same  [17], with an extruder diameter of 0.5 mm. A single layer of painters' tape was used on the heated print plate for print adhesion (Figure 1c).

Experiment
To apply the proper boundary conditions to the two ends of the specimens, ABSplus P430 (3D printer filament, Stratasys, Ltd., Eden Prairie, MN, USA) [18] custom brackets were designed ( Figure 2) and printed using a Stratasys ® Fortus 250mc (3D printer, Stratasys, Ltd., Eden Prairie, MN, USA) [19]. The left bracket, which was attached to a fixed base, allowed one end of the specimens to be completely fixed while the right bracket, which had a hole where weights were attached to a fishing line, bent the specimens as cantilevers. Three sets of brackets were created to tightly fit the three specimen diameters, d. Loads were specimen-dependent and chosen to induce large bending in the specimens (Table A2). Stiffer specimens were, therefore, subjected to heavier loads. Loads were applied in a quasistatic manner by incrementally increasing weight attached to the right bracket, as shown in Table A2. The whole experiment was recorded using an iPhone 11 Pro ® (cellphone, Apple Inc., Cupertino, CA, USA) at 4k resolution at a rate of 60 frames per second.

Experiment
To apply the proper boundary conditions to the two ends of the specimens, ABSplus P430 (3D printer filament, Stratasys, Ltd., Eden Prairie, MN, USA) [18] custom brackets were designed ( Figure 2) and printed using a Stratasys ® Fortus 250mc (3D printer, Stratasys, Ltd., Eden Prairie, MN, USA) [19]. The left bracket, which was attached to a fixed base, allowed one end of the specimens to be completely fixed while the right bracket, which had a hole where weights were attached to a fishing line, bent the specimens as cantilevers. Three sets of brackets were created to tightly fit the three specimen diameters, d. Loads were specimen-dependent and chosen to induce large bending in the specimens (Table A2). Stiffer specimens were, therefore, subjected to heavier loads. Loads were applied in a quasi-static manner by incrementally increasing weight attached to the right bracket, as shown in Table A2. The whole experiment was recorded using an iPhone 11 Pro ® (cellphone, Apple Inc., Cupertino, CA, USA) at 4k resolution at a rate of 60 frames per second. For each specimen, frames corresponding to a steady-state configuration after the load was applied were extracted from the videos of the experiment; corresponding loads were recorded. A custom code utilizing the MATLAB ® Image Processing Toolbox (software, The MathWorks, Natick, MA, USA) [20] was written for processing the images to extract the location of the tip point deflection of the free length, l. First, all images were cropped to reduce unnecessary noise. Next, each image was binarized to isolate the specimen (Figure 3a). The two cylindrical ends of a specimen were identified using morphological detection and removed to leave the middle section ( Figure 3b). The four edges of For each specimen, frames corresponding to a steady-state configuration after the load was applied were extracted from the videos of the experiment; corresponding loads were recorded. A custom code utilizing the MATLAB ® Image Processing Toolbox (software, The MathWorks, Natick, MA, USA) [20] was written for processing the images to extract the location of the tip point deflection of the free length, l. First, all images were cropped to reduce unnecessary noise. Next, each image was binarized to isolate the specimen ( Figure 3a). The two cylindrical ends of a specimen were identified using morphological detection and removed to leave the middle section ( Figure 3b). The four edges of this section were determined using a border detection function, bwboundaries. To identify the neutral axis specimen's middle section, an equal number of equally spaced points were created on the top and bottom edges. For each corresponding point on these two curves, the midpoint was identified to be on the neutral axis ( Figure 3c). Second-order polynomials were fitted to the top and bottom edges and the neutral axis to reduce any noise associated with the averaging (Figure 3d). The last point on the neutral axis of the rectangular section was identified by the intersection of the two cylindrical sections with the neutral axis curve ( Figure 3d). Finally, all relevant information was converted from pixels to millimeters based on camera parameters. this section were determined using a border detection function, bwboundaries. To identify the neutral axis specimen's middle section, an equal number of equally spaced points were created on the top and bottom edges. For each corresponding point on these two curves, the midpoint was identified to be on the neutral axis ( Figure 3c). Second-order polynomials were fitted to the top and bottom edges and the neutral axis to reduce any noise associated with the averaging (Figure 3d). The last point on the neutral axis of the rectangular section was identified by the intersection of the two cylindrical sections with the neutral axis curve ( Figure 3d). Finally, all relevant information was converted from pixels to millimeters based on camera parameters. Figure 3. Specimen image processing: (a) images were cropped and binarized; (b) the cylindrical ends were identified and removed, and the edges of the middle portion were determined; (c) two sets of equally spaced points were created on the top (blue) and bottom (red) borders and averaged to determine the neutral axis (yellow); (d) the top and bottom edges and the neutral axis were fitted with 2 nd order polynomials, and the intersection of the neutral axis and the right cylindrical section was used for the tip point deflection (pink).

NinjaFlex ® Material Model
The material model used in the finite element analysis was based on the uniaxial tensile testing data of dog-bone NinjaFlex ® specimens with wall lines that composed approximately 56% of the specimen width [12]. A hyperelastic 3rd order Mooney-Rivlin incompressible material model presented the best fit for these data ( Figure 4). This model is represented by the following equation [21,22]: where is the strain energy density function and 1 and 2 represent the first and second invariants of the deformation tensor, expressed by, Lastly, represents the stretch of the specimen, or = 1 + ε, with being the strain. 10 , 01 , and 11 are listed in Table 3. . Specimen image processing: (a) images were cropped and binarized; (b) the cylindrical ends were identified and removed, and the edges of the middle portion were determined; (c) two sets of equally spaced points were created on the top (blue) and bottom (red) borders and averaged to determine the neutral axis (yellow); (d) the top and bottom edges and the neutral axis were fitted with 2 nd order polynomials, and the intersection of the neutral axis and the right cylindrical section was used for the tip point deflection (pink).

NinjaFlex ® Material Model
The material model used in the finite element analysis was based on the uniaxial tensile testing data of dog-bone NinjaFlex ® specimens with wall lines that composed approximately 56% of the specimen width [12]. A hyperelastic 3rd order Mooney-Rivlin incompressible material model presented the best fit for these data (Figure 4). This model is represented by the following equation [21,22]: where ω is the strain energy density function and I 1 and I 2 represent the first and second invariants of the deformation tensor, expressed by,   Based on the stress strain curve of Figure 4, E was assigned a value of 9.45 MPa by calculating the slope of the linear range of the stress-strain curve: 0 to 0.2 m/m [12]. This value is lower than that reported by the manufacturer, but falls within the range reported by other researchers, between 5.24-12.2 MPa [8,9,11,12,15] (Table 1).

Analytical Modeling Using Modified Euler-Bernoulli Equations
Classical equations of material mechanics can be used to predict the component deformations under known loads. For example, traditional beam theory was used to model a multibody system composed of NinjaFlex ® and rigid links of acrylonitrile butadiene styrene (ABS) [23]. Similarly, Euler-Bernoulli beam theory was used to model hyperelastic pneumatic actuators made of silicone [3]. However, classical beam equations may not be able to accurately describe the deformation of parts such as those made of NinjaFlex ® due to the hyperelastic material behavior and the large deflection these parts may experience.
In this section, the use of a linear material model for the MEB method was deemed appropriate if the maximum strain is with the linear portion of Figure 4. We propose modifying the Euler-Bernoulli beam theory to account for the deformation of the beam in the axial direction due to bending. As such, the component of the normal load will be a function of the load point's rotation and the load variation (Figure 5a). This accounted for the forces and moments acting upon the experimental specimen as well as the load variation as the angle of deflection increased. Euler-Bernoulli beam equations were used as a basis, though through the comparison of the virtual work done on the system and the strain energy in the system, the deflections of the beams were derived. Therefore, the virtual work and virtual strain energy of the system were assumed to be equivalent. Lastly, λ i represents the stretch of the specimen, or λ = 1 + ε, with ε being the strain. C 10 , C 01 , and C 11 are listed in Table 3. Based on the stress strain curve of Figure 4, E was assigned a value of 9.45 MPa by calculating the slope of the linear range of the stress-strain curve: 0 to 0.2 m/m [12]. This value is lower than that reported by the manufacturer, but falls within the range reported by other researchers, between 5.24-12.2 MPa [8,9,11,12,15] (Table 1).

Analytical Modeling Using Modified Euler-Bernoulli Equations
Classical equations of material mechanics can be used to predict the component deformations under known loads. For example, traditional beam theory was used to model a multibody system composed of NinjaFlex ® and rigid links of acrylonitrile butadiene styrene (ABS) [23]. Similarly, Euler-Bernoulli beam theory was used to model hyperelastic pneumatic actuators made of silicone [3]. However, classical beam equations may not be able to accurately describe the deformation of parts such as those made of NinjaFlex ® due to the hyperelastic material behavior and the large deflection these parts may experience.
In this section, the use of a linear material model for the MEB method was deemed appropriate if the maximum strain is with the linear portion of Figure 4. We propose modifying the Euler-Bernoulli beam theory to account for the deformation of the beam in the axial direction due to bending. As such, the component of the normal load will be a function of the load point's rotation and the load variation (Figure 5a). This accounted for the forces and moments acting upon the experimental specimen as well as the load variation as the angle of deflection increased. Euler-Bernoulli beam equations were used as a basis, though through the comparison of the virtual work done on the system and the strain energy in the system, the deflections of the beams were derived. Therefore, the virtual work and virtual strain energy of the system were assumed to be equivalent. It was assumed that the cylindrical ends were fitted perfectly within the fixtures with no slippage during the experiment (Figure 5a). The free-end bracket had a significantly higher stiffness than the specimen and was assumed to be rigid. To analyze the deflection of the specimens, the variations in the virtual work and strain energies were equated to zero: where W is the virtual work, UA is the axial strain energy, and UB is the bending strain energy. The work is defined as: where is the load from the applied weight, is the load from the bracket weight, θ is the angle between the plane cross-section and the vertical axis, u is displacement in the x-direction, v is displacement in the y-direction, is the moment arm from the applied weight, and is the moment arm from the weight of the bracket. Or simplifying, = ( ( ( )) ( ) + cos( ( )) ( )) + cos( ( )) ( ) where the total applied load is = + and the total applied moment is = + . It should be noted that the slope of the tip point is related to the deflection derivative using the following equation: The strain energy equations for both the axial and bending strain are as follows: where E is the modulus of elasticity, A is the cross-sectional area, and I is the area moment of inertia of the rectangular cross-section. A and I are defined as: It was assumed that the cylindrical ends were fitted perfectly within the fixtures with no slippage during the experiment (Figure 5a). The free-end bracket had a significantly higher stiffness than the specimen and was assumed to be rigid. To analyze the deflection of the specimens, the variations in the virtual work and strain energies were equated to zero: where W is the virtual work, U A is the axial strain energy, and U B is the bending strain energy. The work is defined as: W = (P W + P B )sin(θ(l))u(l)+ (P W + P B )cos(θ(l))v(l)+ (P W a W + P B a B )cos(θ(l))θ(l) (5) where P W is the load from the applied weight, P B is the load from the bracket weight, θ is the angle between the plane cross-section and the vertical axis, u is displacement in the x-direction, v is displacement in the y-direction, a W is the moment arm from the applied weight, and a B is the moment arm from the weight of the bracket.
Or simplifying, W = P(sin(θ(l))u(l) + cos(θ(l))v(l)) + Mcos(θ(l))θ(l) where the total applied load is P = P W + P B and the total applied moment is M = P W a W + P B a B . It should be noted that the slope of the tip point is related to the deflection derivative using the following equation: The strain energy equations for both the axial and bending strain are as follows: where E is the modulus of elasticity, A is the cross-sectional area, and I is the area moment of inertia of the rectangular cross-section. A and I are defined as: Equations (8) and (9) lead to these two differential equations in terms of v and u: v xxxx = 0 (12) These differential equations are subject to the following displacement boundary conditions: Additionally, these differential equations have the following force boundary conditions: Pcos(θ(l)) = −EIv xxx (l) (15) −Psin(θ(l)) (16) or, ((Pu(l) + M)cos(θ(l)) − (Pv(l) + Mθ(l))sin(θ(l))) = EIv xx (l) (17) Psin(θ(l)) = EAu x (l) (18) The solutions to the deflections u and v are: Using the force boundary conditions, C = ((Pu(l) + M)cos(θ(l)) − (Pv(l) + Mθ(l))sin(θ(l))) + PLcos(θ(l)) It can be seen that B, C, and D are functions of the external forces, specimen dimensions, and u(x), v(x), and v x (x). These three equations can be solved simultaneously using the MATLAB ® function fsolve, which requires an initial guess. It was decided to use the classical beam theory to provide the initial guess for the first step, i.e., a beam is deforming under the weight of the bracket. In this case the initial guesses for B, C, and D were − P B EI , P B l+P B a B EI , and P B l EA respectively. Subsequently, the results of an iteration were used as the initial guesses for the next one. To ensure their convergence of the solution, each load step was divided into ten sub-steps. The deflection equations of the previous section only described bending and axial deflection. These equations, however, do not account for the retraction along the x direction that is typically associated with large deformation [24]. Based on Figure 5b, the differential retraction, du E , of the beam can be expressed as (24).
Integrating this equation, the retraction of the tip point, u E , can be obtained: Applying the derivative of Equation (20), the tip point retraction equation becomes: The total deformation of the tip point can be expressed as a combination of the retraction in addition to the extension caused by (P W + P B )sin(θ(l)):

Finite Element Analysis Modeling
Several researchers investigated the factors that affect the modeling of 3D printed components using FEA [25,26]. While these studies were all reasonably successful, showing the desired material parameters or predictive capabilities, they lacked focus on hyperelastic materials such as NinjaFlex ® as well as large deflection. By understanding the bending behavior of 3D printed hyperelastic polymers, they may be used in a wider range of applications.
ANSYS ® (simulation software, ANSYS, Inc., Canonsburg, PA, USA) was used for the FEA simulations [27]. The model of the specimen was imported as a 2D IGES file and processed as a plane stress problem The free end bracket was simplified from the one used in the experiments (Figure 2) to reduce the computational load of the simulations. To ensure that this model was comparable to the experiment, three issues were addressed. Firstly, since the bracket was significantly stiffer than the specimens, it was assigned structural steel material properties. Secondly, the center of gravity's location was maintained. Lastly, gravity was turned off, and the weight of the actual fixture was applied as a force in the center of gravity's location. Figure 6a shows the model. The values of said weights are found in Table A2. A mesh stability study was conducted, and it was found that an average element size of 0.10 mm was stable. The number of elements in simulations ranged from 5483 elements in Specimens 1-3 to 17,999 in Specimens 16-18. Figure 6b represents a typical mesh.

Results and Discussion
The maximum strain calculated by FEA for any model did not exceed 0.16 mm/mm, which is within the linear range of Figure 4. Therefore, the modulus of elasticity, E, in the modified Euler-Bernoulli model was assigned a value of 9.45 MPa, as shown in Section 2.2.
The experimental results for the three samples of each specimen were averaged, and their standard deviations were calculated in the vertical and horizontal, y and x, directions. The corresponding MEB and FEA results were also recorded. Figure 7 shows the resulting deflections of one of the stiffest and most flexible specimens, Specimen 15 and 1, respectively. The dimensions of the specimens are listed in Table A1.

Results and Discussion
The maximum strain calculated by FEA for any model did not exceed 0.16 mm/mm, which is within the linear range of Figure 4. Therefore, the modulus of elasticity, E, in the modified Euler-Bernoulli model was assigned a value of 9.45 MPa, as shown in Section 2.2.
The experimental results for the three samples of each specimen were averaged, and their standard deviations were calculated in the vertical and horizontal, y and x, directions. The corresponding MEB and FEA results were also recorded. Figure 7 shows the resulting deflections of one of the stiffest and most flexible specimens, Specimen 15 and 1, respectively. The dimensions of the specimens are listed in Table A1.
In the case of Specimen 15, which was relatively stiff, the experimental tip deflection deviated slightly from moving along a smooth curve, most noticeably in the x-direction. On the other hand, the tip point of Specimen 1 moved along a smoother curvature. All other tested specimens followed these trends. For example, the tip points of Specimens 7-9 and 13-15, which were relatively stiff, did not follow a smooth curvature. This may be the combined result of small magnitudes of motion, the camera quality, and the inevitable blur. Additionally, it was observed that the standard deviation increased as loading increased for all specimens, which may be the result of built-up issues with the experimental setup or variation in the 3D printer quality.
The MEB method reasonably predicted the deflection of the relatively stiff Specimen 15. The deflection was overestimated for the lower loads and then underestimated as the load increased. For Specimen 1, the MEB lagged behind the experimental results, although it followed the average experimental tip point curvature. The MEB method worked generally well for stiff specimens, such as Specimen 15. For all tested specimens, the curvature of the MEB deflection consistently followed the curvature of the experiments. In the case of Specimen 15, which was relatively stiff, the experimental tip deflectio deviated slightly from moving along a smooth curve, most noticeably in the x-directio On the other hand, the tip point of Specimen 1 moved along a smoother curvature. A other tested specimens followed these trends. For example, the tip points of Specimens 7 9 and 13-15, which were relatively stiff, did not follow a smooth curvature. This may b the combined result of small magnitudes of motion, the camera quality, and the inevitab blur. Additionally, it was observed that the standard deviation increased as loading in creased for all specimens, which may be the result of built-up issues with the experiment setup or variation in the 3D printer quality.
The MEB method reasonably predicted the deflection of the relatively stiff Specime 15. The deflection was overestimated for the lower loads and then underestimated as th load increased. For Specimen 1, the MEB lagged behind the experimental results, althoug it followed the average experimental tip point curvature. The MEB method worked gen erally well for stiff specimens, such as Specimen 15. For all tested specimens, the curvatur of the MEB deflection consistently followed the curvature of the experiments.
The FEA method slightly overestimated the deflection of Specimen 15. Similarly t MEB, the curvature of the FEA tip point defection matched the experiment well. Howeve the FEA resulted in an overestimation of the tip point retraction in the case of Specime 1, as shown in Figure 7b. A possible reason for this behavior may be due to the way AN SYS analyzes bodies experiencing large deflection. When a body experiences large strai the deformation triggers reorientation of the applied loads. It may be possible that th load distribution was less accurate, resulting in the excessive retraction that can be ob served in the figure.
To compare the experimental and modeling data, the deflection under the weight o The FEA method slightly overestimated the deflection of Specimen 15. Similarly to MEB, the curvature of the FEA tip point defection matched the experiment well. However, the FEA resulted in an overestimation of the tip point retraction in the case of Specimen 1, as shown in Figure 7b. A possible reason for this behavior may be due to the way ANSYS analyzes bodies experiencing large deflection. When a body experiences large strain, the deformation triggers reorientation of the applied loads. It may be possible that this load distribution was less accurate, resulting in the excessive retraction that can be observed in the figure.
To compare the experimental and modeling data, the deflection under the weight of the bracket alone, P B , was considered the datum. A measure of the difference of each loading case of a specimen was performed. These measures were averaged as follows: where m is the modeling method (MEB or FEA), I is the specimen number, j is the load number, J i is total number of loads, xe, x m , ye and y m are the experiment and model deflection in the x and y directions, respectively, and l i is the length of the specimen. All results from experiments, MEB, and FEA were combined using the following nondimensional stiffness parameter, Q [28][29][30]: It should be noted that lower Q values typically correspond to higher bending stiffness. Overall, MEB underestimated the deflection of most specimens when compared to the respective experimental results, while FEA overestimated them. The results of the average difference, Z m , were plotted with respect to the non-dimensional parameter Q (Figure 8). The data were split by length for clarity. For specimens with l = 10 mm length (Figure 8a), both MEB and FEA predictions were of the same order for values of Q below 4. The same observation was valid for all tested specimens in the case of l = 15 mm (Figure 8b). However, when l was equal to 10 mm and Q was above 5, the FEA was more accurate. However, the retraction along the x direction was significantly higher than that with MEB, where the displacement closely followed the experimental curve, as can be seen in Figure 8b. It was also observed that the FEA showed more stable Z values as Q increased.
All results from experiments, MEB, and FEA were combined using the following nondimensional stiffness parameter, Q [28][29][30]: It should be noted that lower Q values typically correspond to higher bending stiffness.
Overall, MEB underestimated the deflection of most specimens when compared to the respective experimental results, while FEA overestimated them. The results of the average difference, Ζ , were plotted with respect to the non-dimensional parameter Q (Figure 8). The data were split by length for clarity. For specimens with = 10 mm length (Figure 8a), both MEB and FEA predictions were of the same order for values of Q below 4. The same observation was valid for all tested specimens in the case of = 15 mm (Figure 8b). However, when was equal to 10 mm and Q was above 5, the FEA was more accurate. However, the retraction along the x direction was significantly higher than that with MEB, where the displacement closely followed the experimental curve, as can be seen in Figure 8b. It was also observed that the FEA showed more stable Z values as Q increased.

Conclusions
This research aimed to characterize the bending behavior of 3D printed hyperelastic polymer components and to develop a relatively simple model that can predict the deformation of these components accurately. A series of bending experiments were conducted on 18 different cantilever specimens with significantly different flexibilities. A model was developed based on fundamental beam theory: the modified Euler-Bernoulli (MEB)

Conclusions
This research aimed to characterize the bending behavior of 3D printed hyperelastic polymer components and to develop a relatively simple model that can predict the deformation of these components accurately. A series of bending experiments were conducted on 18 different cantilever specimens with significantly different flexibilities. A model was developed based on fundamental beam theory: the modified Euler-Bernoulli (MEB) model, to account for the large deformations and planar motion of the specimens. A finite element analysis was also conducted to predict the bending behavior. The results of this MEB model and FEA were compared to the experimental data, indicating that both approaches produced accurate results for most specimens. However, the MEB underestimated the overall deflection, while the FEA resulted in a significantly larger retraction of the specimens under load.
MEB can sufficiently provide an accurate prediction of deformation for a large range of hyperelastic 3D printed polymer specimens and has several advantages over FEA, including simplicity and speed. MEB exhibited consistently more accurate retraction prediction, with the deflection of the tip point aligning with the experimental results; however, the deflection in the y-direction was underestimated. The cause of this may be due to the fundamental theory used to derive the MEB method, the perfectly linear material model, or the assumption of a homogeneous isotropic specimen. FEA predicted more accurate deflection in the y-direction and maintained relatively accurate retraction under most circumstances. However, the highly flexible specimens modeled by FEA exhibited overestimation of retraction. This may be due to the print geometry within the specimen. As shown in Table 1, the mechanical characteristics of the material used exhibited variability caused by differences in the type of 3D printer, printer settings, and filament size, including characteristics such as modulus of elasticity, yield strength, ultimate strength, and hardness. For a more complete understanding, numerous polymers with a wide range of mechanical characteristics should be tested and modeled using the methods described. This way, the gradient of results would show a clearer behavior of hyperelastic polymers. Alternatively, the excessive retraction in the flexible specimens exhibited in the FEA results may be mitigated by different material models that represent the behavior better. Additionally, the mechanical properties of the core portion, which is printed using 45 • lines, are expected to differ from those of the wall lines. These variations pose a challenge to modeling components made of 3D printed polymers, especially those undergoing large deflection. An in-depth study of these parameters should be considered for future work.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author.

Conflicts of Interest:
The authors declare no conflict of interest.