Finite element prediction of distortion during the fused filament fabrication process of PLA

Additive manufacturing (AM) has rapidly grown in popularity in recent years due to the availability of low-cost approaches such as Fused Filament Fabrication (FFF) (Tientcheu et al., 2023). It finds numerous applications in Biomedical Engineering (Devine, 2019), product design (Loy et al., 2023) and engineering education (Eslahi et al., 2020). FFF is the most popular AM technique as it allows the production of bespoke, complex parts from a variety of polymeric materials in a range of environments (Cattenone et al., 2019; Roy et al., 2019; Alzyod and Ficzere, 2023b). However, several drawbacks currently limit the use of FFF in industry due to slow build speeds, inconsistent or poor mechanical properties and poor dimensional accuracy (Syrlybayev et al., 2021; Alzyod and Ficzere, 2023b).

In FFF, polymer material is heated to a semi-molten state (Roy et al., 2019), extruded through a nozzle onto a part’s surface, bonding to the previous layer through a polymer reptation process. Parts are generally printed at room temperature (Antony Samy et al., 2021; Syrlybayev et al., 2021) or at some intermediate temperature if an enclosure is used (Cattenone et al., 2019; Roy et al., 2019; Trofimov et al., 2022). Enclosure temperatures are usually maintained close to the glass transition temperature, T_g, as thermal contractions below this point most likely lead to residual stress and distortion (Armillotta et al., 2018; Trofimov et al., 2022). Therefore, to reduce residual stress and distortion magnitudes, maintaining build temperatures around the T_g ensures uniform cooling when thermal contraction becomes more pronounced.

The printed part is constrained by its adherence to a build plate which leads to tensile stress magnitudes at the interface as the material cools and contracts. As each layer is added to the part, the previous layers are heated through conduction and then cool to ambient temperature (Armillotta et al., 2018; Cattenone et al., 2019). This heating and cooling process results in thermal expansion and contraction resulting in the development of a complex residual stress distribution thereby leading parts to distort on the print bed. This has the potential to result in failed prints through delamination and cracking, and also due to the print head colliding with the part (Aljarrah et al., 2019; Alzyod and Ficzere, 2023b; Hinchy, 2019). Researchers have used the residual stress distribution to create novel, complex shapes through four-dimensional printing strategies, where parts are heat treated after printing, resulting in parts warping in a predictable way (Huang et al., 2024). However, without models to predict the distortion behaviour, it is difficult to achieve accurate results.

FFF machine vendors advise on the optimum print settings for each material to achieve robust parts and minimise the effects of unwanted residual stress. These settings which include tool path (Antony Samy et al., 2021), bead width, layer thickness, printing speeds and nozzle, bed and enclosure temperatures are generally established through trial-and-error approaches (Hu and Mahadevan, 2017; Aljarrah et al., 2019). Research has shown that distortion is mostly affected by the geometry (largest dimension on print bed and total height), layer height (Armillotta et al., 2018; Aljarrah et al., 2019; Syrlybayev et al., 2021; Antony Samy et al., 2021; Alzyod and Ficzere, 2023b; Zhang and Chou, 2008), rate of cooling (Syrlybayev et al., 2021; Antony Samy et al., 2021; Zhu et al., 2021), raster pattern (Antony Samy et al., 2021; Chen et al., 2023; Huang et al., 2024; Yu et al., 2025) and print speed (Aljarrah et al., 2019; Alzyod and Ficzere, 2023b; Zhang and Chou, 2008). Previous models designed to predict residual stress have been shown to be mostly qualitative (Trofimov et al., 2022) when compared with experimental results and researchers have found that the models are computationally expensive (Cattenone et al., 2019; Aljarrah et al., 2019; Roy et al., 2019; Trofimov et al., 2022).

Aside from FE models, analytical models to predict distortion in simple FFF parts for Acrylonitrile Butadiene Styrene can indicate the impact of layer height, processing temperature and geometry (Armillotta et al., 2018). However, they lack the complexity to accurately predict distortion over a range of thickness values but are useful for quickly estimating likely distortion sources.

Given the propensity of three-dimensional (3D) printed FFF parts to distort, it is important to be able to predict this behaviour using a scientific approach. With a model, parameters can be varied to achieve an optimum print setup without wasting significant amounts of material. This paper aims to accurately predict distortion during FFF printing of Poly(lactic acid) (PLA) using the finite element method (FEM). PLA is chosen due to its propensity to distort in certain geometries as indicated in previous research (Trofimov et al., 2022). From this starting point, the model and approach can easily be updated for other polymers with new material property definitions and related thermal boundary conditions. This paper further aims to reduce the solution time and improve the accuracy of distortion predictions at room temperature acting as a guide for future researchers.

There is no single standard geometry used in the literature to verify distortion from FFF printing. Previous example geometries analysed include planar springs (Cattenone et al., 2019), bridge samples (Cattenone et al., 2019; Alzyod and Ficzere, 2023a; Trofimov et al., 2022), tensile test samples (Syrlybayev et al., 2021; Yu et al., 2025) and simple flat rectilinear shapes (Aljarrah et al., 2019; Antony Samy et al., 2021; Chen et al., 2023; Zhu et al., 2021; Alzyod and Ficzere, 2023b; Huang et al., 2024; Yu et al., 2025). Some of the samples include support material in the build, which is subsequently removed allowing the part to distort to an equilibrium stress state (Cattenone et al., 2019). Most models have 100% infill, but some models have an infill as low as 20% which can make direct comparison between the experiment and the model more difficult (Cattenone et al., 2019; Alzyod and Ficzere, 2023b).

From reviewing the results in the literature, the prediction models are inconsistent and do not always predict the measured distortion. Due to the relatively small volume of material required and the propensity for distortion, the bridge sample (see Figures 1 and 2) described by Trofimov et al. was used for analysis in this work (Trofimov et al., 2022).

To obtain an accurate prediction of distortion, a range of temperature dependent material properties are required. These properties include Density, Elastic modulus, Poisson’s ratio, Yield stress, Thermal expansion coefficient, Specific heat capacity and Thermal conductivity. The mechanical properties change rapidly at the glass transition temperature (T_g) with properties for yield strength dropping to less than one-third of room temperature values (Trofimov et al., 2022). It is therefore important to know the properties at this point and include this information in the model (Cattenone et al., 2019; Huang et al., 2024). PLA was chosen for this work given the popularity of the material, with the properties initially taken from the work by Trofimov et al. (2022).

Abaqus is a popular choice for predicting thermal distributions (Roy et al., 2019) and residual stress and distortion in FFF (Cattenone et al., 2019; Aljarrah et al., 2019; Huang et al., 2024; Trofimov et al., 2022) and is the software code used in this manuscript. Most models analysed use a two-step approach, first solving the thermal analysis problem followed by the mechanical problem in an uncoupled approach (Cattenone et al., 2019; Aljarrah et al., 2019; Trofimov et al., 2022). Sequential element activation is essential to simulate the arrival of new material to the model. Details of the Abaqus implementation can be found in supplementary information in the research repository created to accompany this manuscript (OSF Database, 2025).

Some models assume that the part is fixed to the build plate, so a Dirichlet boundary condition is applied to the base of the sample during the print and the print bed is not included in the model (Cattenone et al., 2019). In reality, strain exists between the sample and the build plate, and it is likely that some areas of the part, especially on the first layers at the corners (Antony Samy et al., 2021; Chen et al., 2023; Alzyod and Ficzere, 2023b), will not be bonded. Some models have assumed no adhesion between the part and the build plate but fully constrained the first set of elements to be activated (Chen et al., 2023; Yu et al., 2025). Others have included the build plate in the model and used a tied contact with the built part (Huang et al., 2024) or just included the build plate in the thermal model and assumed fixed boundary conditions for the stress-displacement model (Trofimov et al., 2022). At the end of the print, the model simulates removal from the build plate by removing the boundary conditions at the base of the sample, allowing the part to distort (Cattenone et al., 2019). Others have used a spring formulation to adhere the part to the build plate allowing the part to distort during the build process (Antony Samy et al., 2021).

In this manuscript, the build plate (measuring 172 × 16 × 1 mm³) was included in the thermal and stress-displacement analyses for the first step which modelled the build. The part was placed on the centre of the build plate, so that the centre of the base surface was at (0, 0, 0) which was 43 mm from the end of the part in the x-direction and 4 mm from the edge of the part in the z-direction. For the second step, which modelled the removal from the build plate, three nodes were constrained as defined in Figure 2. Altering these boundary conditions were found to have little effect on the overall distortion predicted.

The build plate was included in the thermal analysis and had a temperature of 52°C, while the starting temperature for the printed part was 150°C as per the work by Trofimov et al. (Trofimov et al., 2022). These values were verified during the print using a Flir One Digital camera (see Section 3.1). The heat transfer coefficient used was estimated to be: h¯_x = 12 W.m⁻².K⁻¹ based on the theory discussed in supplementary information which can be found in the corresponding research repository (OSF Database, 2025).

To obtain a precise temperature and stress prediction, elements are often chosen so that the height of the elements is equivalent to the layer height or a submultiple of the height; similarly, the width of the elements is chosen to be equal to the width of the extruded filament, or a submultiple of it (Cattenone et al., 2019; Antony Samy et al., 2021; Trofimov et al., 2022). For large parts, this can result in a lot of elements and with significant computing power required to solve them (Trofimov et al., 2022), even for just the temperature model (Roy et al., 2019) or for modelling the problem in only two dimensions (Zhu et al., 2021). For example, with a filament diameter of 0.4 mm and a layer thickness of 0.2 mm, a simple cubic shape with 100 mm edge length would require over thirty-one million elements. Studies of mesh density have shown that coarser meshes are less accurate at predicting the precise temperature but have little effect on the stress and distortion prediction (Cattenone et al., 2019). If researchers are just interested in distortion, then it has been proposed that it is acceptable to use a coarser mesh (Van-Der-Veen et al., 2022).

Given the significant computational cost associated with increasing the mesh density, this analysis focussed on four meshing strategies as follows:

1. Elements with a height of 0.2 mm and width of 0.4 mm, exactly corresponding to the size of the extrudate (number of elements in bridge = 75,400).

2. Cubic elements with an edge length of 0.4 mm (36,500 elements).

3. Cubic elements with an edge length of 0.3 mm (91,098 elements).

4. Cubic elements with an edge length of 0.6 mm (9,672 elements).

Three maximum time increments were run for these four mesh densities (0.3, 0.5 and 1 s, respectively – see Section 2.7). These time increments were used for Step 1 for both the thermal and stress analysis studies. A maximum time increment of up to 50 s was allowed for Step 2, as this was found to change the predicted displacement by < 1% when compared to a maximum time increment of 5 s.

The distortion in the y-direction was recorded using the Path function in Abaqus along a line on the top of the bridge part from (0, 6, 4) to (86, 6, 4) – see Figure 2. From this data the distortion was calculated as the average of the displacement in the y-direction at the two end points minus the distortion in the centre (43, 6, 4). In the results presented by Trofimov et al. (2022), their prediction model gave a total distortion of 1.22 mm and their experiments showed a total distortion of 1.34 mm. The results of this mesh density study are shown in Figure 3 and indicate that the 0.2/0.4 mm mesh and 0.6 mm mesh produced the most precise results. While the 0.6 mm mesh did produce precise results for the time increments analysed, the results differ quite substantially (>30%) from the other three mesh density models. Only the 0.2/0.4 mm mesh produced results that were both accurate and precise. The use of tetrahedral elements resulted in a distortion that was much larger than those predicted using the brick elements. For example, using a maximum time increment of 0.5 s and typical element size between 0.2 and 0.4 mm produced a distortion 0.73 mm for quadratic tetrahedral elements and 1.91 mm for linear tetrahedral elements. Both values would be outside the range of results plotted in Figure 3 for the quadratic elements.

To obtain realistic peak temperatures, the recommended time increment should be smaller than half the characteristic dimension of the heat source in the direction of travel, divided by the speed of advance (Van-Der-Veen et al., 2022). For a typical dimension of 0.4 mm and a speed of 30 mm.sec⁻¹, this would result in a maximum time increment of 0.007 s. For a simple model with a print time of 20 min, this would require 1.8 × 10⁵ increments, which would be computationally very expensive. Some researchers have previously analysed the effect of relatively low maximum time increments on the processing time for both the thermal and stress analysis (Cattenone et al., 2019; Trofimov et al., 2022) indicating the resultant computational overhead. Accurate prediction of temperatures may be required to indicate if the temperature between layers is above T_g to predict the likelihood of interlayer adhesion and these models may require small time steps (Chen et al., 2023; Van-Der-Veen et al., 2022; Costa et al., 2017). However, capturing the temperature history very precisely may not improve the distortion prediction significantly, so it is considered acceptable to use larger time increments to achieve a result with reasonable computation efficiency (Van-Der-Veen et al., 2022; Cattenone et al., 2019).

A second study was conducted to focus on the effect of the maximum time increment allowed in the temperature and stress analyses, with values varied from 0.1 to 5 s. The results of this study are shown in Figure 4 and show that the 0.2/0.4 mm mesh are more precise and accurate than the 0.4 mm mesh at maximum time increments between 0.1 and 5 s. The optimum maximum time increment chosen with the 0.2/0.4 mm mesh is 0.5 s which results in a Total central processing unit (CPU) time for the temperature and stress displacement model of 67,600 s. This compares very well when compared to the Trofimov et al. (2022) model which took approximately 210,000 s at the optimum values chosen. The model predicts a total distortion of 1.29 mm using the 0.5 s maximum time interval which differs by only 4% when using a maximum time interval of 0.1 s, with a 78% lower CPU time. The predicted distortion results are shown in Figure 5, along with the results from Trofimov et al. (2022), which have been adapted from the original paper. The results clearly show a very good match with this published data set with an R² value 0.97 between the results predicted by Trofimov and those predicted here validating this approach which has used different elements, lower mesh density, different heat transfer coefficients and the AM Modeler approach. The Abaqus input files can be downloaded from the research repository associated with this manuscript (OSF Database, 2025).

• Heat transfer coefficient: The value used in the model was 12 W.m⁻².K⁻¹ which was determined from the assumption that the air velocity was 0.4 m.s⁻¹. The model was run again to simulate an air velocity of 0.1, 0.8, 1.5, 3 and 6 m.s⁻¹ changing the heat transfer coefficients to 6, 17, 23, 32 and 46 W.m⁻².K⁻¹, respectively, as calculated using equations in supplementary information (OSF Database, 2025). This resulted in up to an 88% change in the total displacement in the y-direction. Lower values of the heat transfer coefficient resulted in a decrease in displacement – see Figure 6. The model was run with a sink (enclosure) temperature at 20 and 60°C – it seems clear that printing at room temperature will result in much more significant distortion than at 60°C.

• Radiation: The model assumed an emissivity of 0.92. Removing emissivity from the model resulted in a decrease in distortion of 20%. Clearly, as stated earlier, when the heat transfer coefficient is relatively small (<60 W.m⁻².K⁻¹), radiation has a significant influence on the cooling rate (Trofimov et al., 2022; Costa et al., 2017).

• Free surface: Removing the “Free Surface” option in the model (described in Section 2.3, resulted in a 77% increase in the total displacement. This was due to faster predicted cooling as the model calculated the heat transfer based on the full exposed areas of the element facets as opposed to active elements volumes at the start of each increment.

• Follow deformation: The model used the “Follow Deformation” option to ensure that elements that have not yet been activated follow the movement of the active elements. Turning this option off resulted in a lower overall predicted distortion of 1.25 mm. A comparison of the results is shown in Figure 7. When follow deformation is turned off, a smoother parabola is predicted without the “waviness” observed at the centre of each span which corresponds to the presence of the bridge supports. During the build, the material that is added to the vertical supports of the bridge structure is constrained and pulled towards them as it cools, while the material between the supports is not constrained in the same way and is free to move.

The results outlined in Section 2, clearly show that the model can predict similar values to those found in literature for PLA samples printed in an enclosure with a temperature maintained at 60°C. To further test the capability of this model, parts are printed at room temperature and measured after printing to compare distortion magnitudes, as outlined in Section 3.

An Original PRUSA I3 MK3S + printer was used to conduct test prints of the bridge sample at room temperature following the print pattern outlined by Trofimov et al. (2022). This pattern consisted of line infill for the legs that was parallel to the z-direction and line infill pattern for the top section parallel to the x-direction. This pattern was repeated for all layers. G-Code was generated using PrusaSlicer version 2.8.0. with an infill set to 100%, layer height of 0.2 mm and a standard diameter 0.4 mm Prusa nozzle installed on the hot end. The Prusa files and G-Code can be downloaded from the research repository accompanying this manuscript (OSF Database, 2025). The printer was calibrated using the Prusa built-in calibration system; the bed temperature was preheated to 60°C; the print cooling fan was turned off; a Prusa textured powder-coated steel sheet was used on the heated bed to ensure adhesion of the samples; The “elephant’s feet ^[1]” setting was left turned on, to account for thermal contraction of the first layer; and the filament used was RS PRO – PLA 1.75 mm. The total printing time was 14 min and 18 s. The airflow in the room was repeatedly measured to be < 0.03 m.s⁻¹ at the print bed using a calibrated Testo 405i NTC, Thermal Anemometer. Images recorded using a Flir One Digital camera verified that the temperature of the print bed was 52°C.

Each print had a resting time of one minute on the print bed upon completion. After this, the steel sheet, with the print attached, was removed from the heated bed and placed on a flat surface at room temperature. The print was left cool for designated time periods (1–10 min) before being removed from the steel sheet. The parts were removed with a wide stainless steel scraper along the long edge of the print. Simply pulling the samples from the print bed introduced distortion to the samples, due to relatively low yield strength in the minutes following printing. Four thermal images were recorded throughout the printing and cooling process using a Flir One Digital camera. The thermal imaging results show that after just one minute resting at room temperature, the top surface of the part remained quite hot at 49°C. After two mintutes, the temperature reduced to 40°C, three minutes it was 34°C and four minutes it reached 30°C.

Experimental distortion was measured using a Quadra-Chek 2000 shadowgraph with a resolution of 1 µm. Samples were placed on their top surfaces and leg height was measured, as these points should have been colinear from their contact with the print bed. Distortion was calculated by the maximum minus minimum height, as defined in Section 2.5. While one of the samples did exhibit distortion up to 0.198 mm in the same direction as the model predicted, some samples distorted in the opposite direction to a minimum of −0.115 mm. The average distortion measured was 0.06 mm with a 95% confidence interval for the mean calculated to be between −0.0173 and 0.1375 mm. Further measurements on samples left resting on the print bed for different periods of time revealed no impact on the measured distortion magnitude. Leaving the print cooling fan turned on, appeared to result in a smoother surface finish, but there was no consistent measurable change in overall distortion.

A representative sample was analysed using an xradia Versa XRM-500 x-ray microscope with a pixel size of 41.9387 µm. The full computed tomography (CT) scan can be found in the research repository created to accompany this manuscript (OSF Database, 2025), with some sample images shown in Figure 8 [the files can be viewed by opening them in ImageJ (Schneider et al., 2012) choosing < File – Import – Image Sequence > followed by < Image – Stacks – Orthogonal Views>].

Between the supports, there is significantly more porosity than in the centre portions. Using ImageJ to estimate the porosity indicates that between the columns on the YZ plane porosity approximately 5%, while approximately 0.5% in the centre of the columns. This porosity affects the ability of the model to accurately predict distortion, as it does not take porosity into account.

From Figure 8(a), the position of the base of the columns is correct relative to the computer-aided design (CAD) model, but the end columns have been pulled in towards the centre of the part in the x-direction. The end columns have been pulled towards the centre of the part more than the other columns as these columns have a contracting force acting towards the centre, while all other columns are pulled both towards and away from the centre. Also, as shown in Figure 8(b) and (c), the sample is higher at the sides than in the centre with a total difference of up to approximately 0.5 mm.

Using the CT data, the location of the 60-point vertices (n = 64) from Figure 8(a) were compared with the CAD model using a Root Mean Square Error calculation [RMSE – see equation (1)] where $xi$ and $yi$ are the locations of vertex $i$ from the scan data; $xi¯$ and $yi¯$ are the locations of vertex $i$ for the original CAD data. For this analysis, the image was centred on the central column (43, 0, 4). For the CT scan, the RSME was found to be 0.433 mm, which is relatively small for FFF parts indicating that the geometry is close to the original CAD part. The maximum (δ_max), mean (µ) and standard deviation (s) for the Euclidean distance calculated was 1.115 mm, 0.187 mm and 0.188 mm, respectively. Using this approach for the results predicted in Figure 5 gives RMSE = 0.766 mm; δ_max = 1.687 mm; µ = 0.587 mm; s = 0.454 mm, indicating the higher distortion predicted:

For the elastic-plastic finite element model described in Section 2, the toolpath defined by the event series in Abaqus was updated to match that of the Prusa generated G-Code file; also, the sink temperature was reduced to 22°C; the temperature of the build plate (52°) and starting temperature of the model (150°) were confirmed to be the same as the existing model. To improve the stability of the Abaqus model and improve convergence, a damping factor of 0.0002 was specified in the Step Module which resulted in a static dissipation energy (ALLSD) of up to 1.8 J, which was very small compared to the internal energy (ALLIE) of up to 350 kJ. These files can be found in the associated research repository (OSF Database, 2025). Otherwise, the model remained as described earlier with the “optimised parameters” chosen. This resulted in a total predicted distortion of 3.88 mm (RMSE = 1.155 mm; δ_max = 4.009 mm; µ = 1.334 mm; s = 1.245 mm), which is significantly larger than the experimentally observed values. Varying the cooling time (time between the end of the print and removal from the print bed) had no effect on the resulting distortion.

Before the part was removed from the build plate in the simulation, the maximum temperature (T_max) was 56°C, maximum von Mises stress (σ_vM) was predicted to be 30.2 MPa (elemental average σ_vM = 10.77 MPa), while the largest magnitude principal stresses were: σ₁ = 28.2 MPa, σ₂ = 17.5 MPa and σ₃ = −25.7 MPa. After removal from the build plate and fully cooled T_max = 22.5°C, σ_vM = 38.6 MPa (average = 6.9 MPa), σ₁ = 28.2 MPa, σ₂ = 15.9 MPa and σ₃ = −27.6 MPa. So, while the stresses have redistributed due to the removal of the constraint, the maximum von Mises stress increased due to the increase in material yield strength as it cooled, but the average decreased due to the distortion and redistribution of stress. The strain energy density (SENER in Abaqus), U = 56.82 kJ.m⁻³ before removal from the build plate and U = 29.77 kJ.m⁻³ after it had fully cooled to room temperature. This shows that removal from the build plate leads to a redistribution of the residual stress and reduces the stored energy in the system, even though the stress magnitudes increased in some locations due to the increase in mechanical properties.

Luo et al. (2022) highlights that values for elastic modulus in literature vary from 2.02 to 4.19 GPa, while yield strength is purported to vary from 15.5 to 89.1 MPa at room temperature – varying these values in the model did little to improve the predicted result.

Given the robustness of the model for predicting the distortion indicated by Trofimov et al. (2022), the most likely cause of the discrepancy in results for the room temperature print is due to the material properties chosen, especially at temperatures less than 60°C.

The FE model presented assumes that the mechanical response of the material is independent of strain rate effects. However, several recent studies have shown that the mechanical properties of PLA are strain rate dependent below 60°C, (Luo et al., 2022; Tezel et al., 2019; Ye et al., 2020; Adibeig et al., 2023; Gebrehiwot et al., 2023; Dogan, 2022; Amza et al., 2021; Kontou et al., 2024), with decreasing strain rates resulting in lower yield strengths and elastic modulus values. The simulation predicts strain rates are typically of the order of 10⁻⁴ to 10⁻⁵ s⁻¹. These thermally activated plastic flow (creep) effects were not accounted for in the model and are likely to have influenced the room temperature prediction results.

Several published data sets for the creep properties of PLA were analysed (Tezel et al., 2019; Ye et al., 2020; Adibeig et al., 2023; Dogan, 2022; Gebrehiwot et al., 2023), but the results from different publications gave inconsistent properties. The publication by Adibeig et al. (2023) was the only one that included results for both more than one temperature and applied stress magnitude. Also, the data was calculated for FFF manufactured samples, providing a good match to the model being developed and would therefore have included porosity effects. However, bead level anisotropy or inter-road bonding properties were not taken into account. A two-layer viscoplasticity (TLV) model was implemented in Abaqus alongside the existing material properties. The TLV model was chosen because it separates deformation into elastic, viscoplastic and creep-like contributions, therefore meaning it can handle both short-term and long-term responses (2023). Data was extracted at available stresses, temperatures and strain rates. The Microsoft Excel solver tool was used to fit this data to the Norton-Hoff rate law (2023) with temperature dependent data:

where $ϵV˙$ is the strain rate with respect to time in the viscoelastic network and A, m and n are the Norton-Hoff rate parameters. m is set to zero thereby assuming a steady-state creep approximation, indicating that the strain rate depends on stress only, not time (2023). It is understood that this is a limitation of the model, but creep data of PLA is limited thereby impeding its implementation. Abaqus recommends using the Power law model as the value for A can become very small. q₀ is calculated from:

Data taken from literature (Luo et al., 2022) at room temperature was used to estimate the modulus ratio (f = 0.5) that defines the relative stiffness of the elastic-viscous network:

where K_V and K_P are the elastic modulus of the viscoelastic and elastoplastic networks, respectively. Table 1 shows the values used for the TLV model in Abaqus based on the Adibeig (Adibeig et al., 2023) data. Using the data from Ye et al. (2020) would have resulted in A = 2.89E−10 N⁻ⁿ.mm²ⁿ.s^{−(1 + m)} and n = 2.55 at 23°C and the data from Tezel et al. (2019) would have resulted in A = 2.35E−10 N⁻ⁿ.mm²ⁿ.s^{−(1 + m)} and n = 2.64 at 25°C, which would have been in contradiction with the Adibeig data.

With the implementation of the TLV model, Abaqus predicted a total distortion in the y-direction of 0.118 mm (RMSE = 0.647 mm; δ_max = 0.834 mm; µ = 0.419 mm; s = 0.235 mm), which compares well to the average measured values presented in Section 3.2 of 0.06 mm and is within the 95% confidence limits calculated for the experimental data. Table 2 compares the Euclidean Distance calculations for the CT-Scan and two corresponding FE models. The TLV model is the closest to the original CAD geometry in terms of RMSE and maximum Euclidean distance values, but there is still quite a difference between the final predicted model and the CT scan data. The predicted result can be seen in Figure 9 where a comparison is made between the TLV model and the elastic-plastic model.

Before the part was removed from the build plate in the simulation, the maximum von Mises stress (σ_vM) was predicted to be 92.5 MPa (elemental average σ_vM = 15 MPa), while the largest magnitude principal stresses were: σ₁ = 91.3 MPa, σ₂ = 44.1 MPa and σ₃ = −42.9 MPa. The maximum von Mises stress occurred in the material in contact with the build plate due to the boundary conditions used and the sharp corners that exist in the model leading to a stress singularity. As per St. Venant’s principle, these local stress magnitudes are unlikely to greatly impact the overall distortion. After removal from the build plate, σ_vM = 60.1 MPa (average = 0.59 MPa), σ₁ = 39.6 MPa, σ₂ = −17.6 MPa and σ₃ = −63.1 MPa.

For quantitative prediction of distortion in AM processes, finite element analysis can give reliable results with careful consideration of element size, maximum time incrementation and heat transfer characteristics. In the case analysed here, the extrudate is 0.2 mm high by 0.4 mm wide, which results in reasonable element sizes and run times – typically the total processing time for both the heat transfer and stress models was approximately 6 h on a desktop Personal Computer with Intel® Xeon® 3.5 GHz 12 core processor and 128 GB Random Access Memory (RAM). These values are significantly lower than published data for similar modelling approaches (Cattenone et al., 2019; Roy et al., 2019; Trofimov et al., 2022). However, for a part that takes less than 15 min to print, the finite element analysis time is long (approximately 24 times longer to model than print). For metals processing the computational cost is even higher, as the element size would need to be approximately 0.06 mm high by approximately 0.08 mm wide (McConnell et al., 2024) to capture the deposited material, resulting in significantly more elements. Also, the extruder speed was 40 mm.s⁻¹ for FFF printing necessitating smaller maximum time increments resulting in longer run times. Methods that reduce computational time such as the Eigenstrain approach result in more approximate solutions (Abaqus, 2023). Attempts to use the Eigenstrain approach for this work resulted in much poorer predictions of distortion.

From the literature, FFF printed polymers tend to be strain rate sensitive even at room temperature with significant data available to indicate creep effects between room temperature and the glass transition temperature (Luo et al., 2022; Tezel et al., 2019; Ye et al., 2020; Adibeig et al., 2023; Gebrehiwot et al., 2023; Dogan, 2022; Amza et al., 2021; Kontou et al., 2024). Also, material properties can vary significantly from sample to sample depending on processing history, supplier and even exposure to a humid environment (Luo et al., 2022). Further work needs to be completed to obtain a broader data set of creep properties for PLA from room temperature to the glass transition temperature to obtain a more accurate prediction of residual stress and distortion.

In the experiment, it proved difficult to remove samples from the build plate without introducing permanent deformation as the samples lacked stiffness and strength due to their high temperature relative to the glass transition temperature. Initial samples were pulled from the plate in the build direction, but this resulted in significant distortion. Removing the samples using the scraper described in Section 3.1 helped to reduce experimentally induced distortion but some may have existed as indicated by the variability of the results presented in Section 3.2. If the scraper was not aligned parallel to the z-direction it was found that distortion could be introduced in the z-direction also (see Figure 1 for directions). Therefore, careful attention needs to be given to sample removal from the build plate where distortion is being measured, especially for thin samples such as planar spring (Cattenone et al., 2019), tensile test samples (Syrlybayev et al., 2021; Yu et al., 2025) and simple flat rectilinear shapes (Aljarrah et al., 2019; Antony Samy et al., 2021; Chen et al., 2023; Zhu et al., 2021; Alzyod and Ficzere, 2023b; Huang et al., 2024; Yu et al., 2025). The bridge sample appears to be the most appropriate to help reduce distortion when removing from the build plate as it limits the contact area with the plate, while still maintaining a large top section that is likely to distort (Cattenone et al., 2019; Alzyod and Ficzere, 2023a; Trofimov et al., 2022).

While the model predicted overall distortion close to the measured values it does not predict the dip in the centre of the sample shown in Figure 8(b) and (c), nor does it predict the distortion of the end columns shown in Figure 8(a). The data in Table 2 uses a mathematical approach to compare the predicted and CT scanned results. Comparing the TLV model to the scan data, results in an RMSE value of 0.766 mm and a maximum deviation of 1.687 mm which is quite large. While the introduction of the viscoplastic effects has reduced the overall distortion, it has not accounted for the deformation of the end columns as shown in the CT scan data. The is partly due to the lack of a robust data set for the TLV model with only data included for 30 and 40°C. it is likely that outside of this temperature range, the viscoplastic response is quite different to that calculated for this model, leading to the discrepancy in the final distortion predicted. Combined with the TLV modelling limitations, the model assumes the extrudate has a rectangular cross-section and therefore does not take porosity into account. This will clearly influence the true stiffness of printed parts. To account for this, mechanical properties would need to be determined parallel and perpendicular to the build plate and these results incorporated into future models.

As outlined in Section 3.4 the elastic-plastic model results in lower stress magnitudes, but higher overall distortion than the TLV model. These apparently counterintuitive results can be explained by the fundamental differences in how these material models handle rate-dependent deformation and stress evolution. In the elastic-plastic model, thermal stresses that develop during printing are accommodated through rate-independent plastic flow, resulting in lower residual stresses (at the end of the printing step the average von Mises stress was 10.8 MPa for the elastic-plastic model and 15 MPa for the TLV model). Also, the elastic-plastic model predicts lower strain energy density values as plastic deformation provides immediate stress relief. However, when boundary conditions are removed, all stored elastic energy must be released instantaneously through distortion, leading to higher part warping. The TLV model, in contrast, incorporates rate-dependent plastic flow where the material’s yield strength and deformation behaviour depend on strain rate and time. Given the importance of thermal history for the TLV model future studies should analyse the effect of time incrementation on the resulting predictions and focus on dissipated energy as a result of viscous flow. During the rapid thermal cycling of FFF printing, the high strain rates can significantly increase the material’s resistance to plastic flow, causing thermal stresses to accumulate to higher levels with correspondingly higher strain energy density throughout the analysis. This higher energy density reflects the increased elastic strain energy stored due to the rate-dependent strengthening effect that prevents stress relaxation through plastic deformation. Paradoxically, when constraints are removed, the TLV model predicts lower distortion because the time-dependent nature of plastic flow allows gradual stress redistribution and relaxation rather than instantaneous elastic springback. The material can continue to undergo viscoplastic deformation after constraint removal, effectively “flowing” towards a lower-stress configuration over time, which reduces the final distortion compared to the purely elastic response assumed in the elastic-plastic model. Also, of note, is that the TLV model predicts an increase in the equivalent plastic strain (PEEQ) after the constraints have been removed for the TLV model, but the PEEQ for the elastic-plastic model remains constant.

In the TLV model, the value for the modulus ratio (f = 0.5) is taken from literature at room temperature (Luo et al., 2022). However, at higher temperatures it is likely that this value will be lower as the viscoplastic layer tends to soften more rapidly than the elastic layer, especially near and above the glass transition temperature (Ngai and Plazek, 2007). The model was run with modulus ratio values of 0.1, 0.25, 0.375, 0.5 and 0.75 for the data at 40°C only, with total resulting distortion shown in Table 3. Along with this decrease in distortion there was a corresponding decrease in the average equivalent plastic strain (PEEQ) and average von Mises stress in the model. Meanwhile, the strain energy density remains almost constant as a function of the modulus ratio. However, the standard model with no viscoelastic effects showed an average von Mises stress σ_vm = 6.9 MPa, PEEQ = 9.29 × 10⁻³ and U = 29.76 kJ.m⁻³. The modulus ratio therefore has a large influence on final distortion with lower values leading to higher distortion. However, the features observed in the CT-Scan data were not replicated in the model with the change in the modulus ratio.

To investigate the strain sensitivity value (m), two further models were run – one with m = −0.01 at RT and −0.1 at 40°C; a second with m = −0.05 at RT and −0.2 at 40°C. These models both resulted in a distortion of approximately 0.09 mm, increased the value of strain energy density and von Mises stress, while equivalent plastic strain remained largely the same as the model with m = 0. The strain rate sensitivity parameter clearly increases the material’s flow stress (yield strength) as a function of strain rate leading to the higher von Mises stress and strain energy density for a similar amount of plastic strain.

The FEM can be used to accurately predict distortion in 3D printed parts, but a careful modelling approach is required to ensure accuracy. It was found that the element size should correspond with the size of the extrudate to ensure more accurate predictions of distortion. The optimum maximum time increment was found to be 0.5 s, during which the extruder would have moved approximately 20 mm, given the speed of 40 mm.s⁻¹ used in this study. This time increment was found to give a good compromise between model accuracy and computational time. The chosen maximum time increment was not varied for the TLV model which is something that should be studied in future work given the sensitivity of this model to thermal history. The use of eight-noded brick elements reduced the computational time when compared to previous models but maintained similar levels of accuracy.

The heat transfer coefficient and emissivity were found to have a large effect on distortion as predicted in literature. Slower cooling rates resulted in lower distortion, as expected.

Inconsistent removal of the sample from the build plate was found to introduce residual stress and distortion which could not be accounted for in the model. Using a flat edge scraper or similar on the longest edge of the sample (86 mm) was found to be the best approach to minimise the introduction of distortion.

The work presented shows the importance of a detailed material model and the inclusion of viscoplastic behaviour between room temperature and the glass transition temperatures to more accurately predict distortion. Further work needs to be completed to obtain reliable materials data as a function of temperature, and to investigate approaches to reduce computation cost. The model developed here, which is available to download, will help researchers to develop more advanced approaches to prediction of distortion and residual stresses in AM processes.

The supplementary material for this article can be found online.