Introduction

During the past few years, the importance of additive manufacturing (AM) has been increasing thoroughly due to the great benefits it presents. Once the improvement of the mechanical properties of structures through the modification of their geometries have reached a limit, topology optimization is an approach to improve their capabilities. More importantly, within industries like the automotive and the aerospace, the optimization of structures by generative shape design methods, has become a priority. This enables the reduction of weight [1] and material waste [2], with improved mechanical properties.

AM allows the use of a wide range of materials, from polymers to metals, passing through composites and biomaterials [3]. Depending on the AM technique, some materials are more appropriate than others. Specifically, recycled polymeric materials could be used effectively in Fused Granular Fabrication (FGF) [6], a type of Material Extrusion (MEX) [4] technique where material is introduced into the 3D printer as pellets. FGF enables the use of recycled pellets feedstocks, which are heated and extruded similarly through nozzles, but at substantially bigger sizes than filament-based technologies [5]. This characteristic makes FGF a perfect candidate for the manufacturing of large-scale components using recycled composite materials converted into pellets [6]. It has been demonstrated that the continuous heating and extrusion of polymers have a negative effect on their mechanical properties [8], limiting the recycling life to 5 cycles [7]. Therefore, to decrease the number of melt/extrude cycles of recycled plastic used in 3-D printing, one alternative is to eliminate the requirement for filament and print directly from recycled plastic pellets, flakes, regrind, or shreds [6]. As an example, Gigabot X [9] is already being used as a successful large-scale 3D printer using recycled materials as pellets.

Recently developed large-format additive manufacturing (LFAM) processes have revolutionized industrial production, offering tremendous potential for cost-effective and efficient manufacturing [11, 12]. As a consequence of larger extrusion systems and nozzle diameters, LFAM has a higher production speed compared to conventional AM, as a greater material volume could be deposited in less time [13]. Therefore, this manufacturing technique could be a solution for manufacturing larger molds, or components, with applications in the aerospace or automotive industries [11, 12], or even the naval industry [14].

Nevertheless, industrial components manufactured by LFAM - and in general AM - processes present a defect that should not be neglected: geometry distortion [15]. Residual stresses are generated in the structure due to thermal gradients between adjacent layers during the manufacturing process [13]. After the 3D printing of the part, as it cools down, energy stored is released due to these residual stress, producing a undesired deformation of the final product. This deformation becomes more significant in LFAM [13], aimed for industrial applications, where tolerances should be fully respected to achieve the correct assembly of final products. Therefore, there is an urge in the industry to understand how these deformations are generated and how could warpage be controlled or diminished. Some studies propose to improve geometrical accuracy in material extrusion methods, for example, by insulating the printing chamber to maintain a constantly high temperature to reduce thermal stresses [16], or by introducing expanded perlite fillings into a polypropylene matrix [17].

In this context, in light of the ongoing drive towards digitalization in the industrial sector through the advent of Industry 4.0, an analogous advancement should be considered to transform AM into a smart manufacturing method [10]. Accordingly, the development of Digital Twins (DT) is becoming a promising approach [18]. A DT is a very accurate virtual replica of a system. The depicted system might be of any type, including materials, manufacturing systems, biological components, etc. DTs can incorporate sensors that generate data on their performance, providing for feedback [19]. Although both, simulations and DT, employ digital models to recreate the various systems, a DT is essentially a virtual world, therefore it contains far more features and data [20, 21].

The present article suggests the integration of DT and AM technologies via an innovative multi-scale approach. Numerical modeling strategies are used to simulate the micro and macroscopic mechanical behavior. In this way, the mechanical properties are obtained and the warpage of the 3D printed component could be predicted. Previously, mean-field homogenization methods have been presented to predict the elastic modulus of foamed composite materials after a microscopic characterization of the material[22]. Moreover, Bacciocchi et al. [23] anticipated the comprehensive elastic attributes of shot-earth, a sustainable composite material, through a multiphase mean-field homogenization approach. Conversely, Kudimova et al. [24] introduced finite-element homogenization techniques to determine the effective properties of piezoelectric composites.

Consequently, the main purpose of this article is the proposal of a useful tool in the development of a DT, able to predict the warpage of industrial LFAM components, prior to their manufacturing, reducing the waste material. The primary aim of this DT is to alleviate structural distortions and reduce the number of printings before obtaining the desired product. The innovative methodology proposed in this study is based in three steps: (i) AM process, (ii) material characterization and (iii) mechanical testing. During the FGF process, a composite material, in pellet form, is employed. Firstly, this material is going to be characterized by means of mean-field (MF) homogenization techniques and a finite element (FE) homogenization strategy. Later, its mechanical properties are going to be studied at two different scales: micro and macroscale. By comparing the results obtained, the differences between strategies are going to be highlighted. This will serve as the foundation for comprehending the deformations exhibited by the final products subsequent to their printing.

Materials

The manufacturing of an external expansion aerodynamic package mold for a Formula Student (FS) car is proposed as a demonstrator of the described simulation model.

In this paper, Digimat and Ultimaker One will be used to validate the properties of LNP™ THERMOCOMP™ AM Compound AC004XXAR1 and develop an optimized 3D printing process for this material. Digimat will be used to create a precise digital model of the part, used for simulations and analysis. Ultimaker One will be used to generate a G-code for the 3D printing of the part. The results of this study will provide valuable information for the development and optimization of this material for specific applications.

Composite material: ABS with 20% CF

The composite material proposed to study is LNP™ THERMOCOMP™ AM COMPOUND AC004XXAR1 [25], a proven, high-quality material that is a solid choice for similar projects. This material is based on an ABS thermoset resin and contains 20% short CF. It exhibits outstanding resistive characteristics and an adequate working temperature of 101ºC at a pressure of 1.82 MPa. These properties make it suitable for a wide range of applications and processes, including thermoforming and vacuum forming, as well as 3D printing.

The most relevant mechanical properties of this material are summarized in Table 1.

Table 1 Main mechanical properties of the composite material LNP\(^{TM}\) THERMOCOMP\(^{TM}\) AC004, given by SABIC [25]
Full size table

The studied composite material is ideal for the manufacture of autoclave molds, as it meets the pressure resistance requirements necessary for this specific application. Its ability to withstand loads and pressures, prevents deformation or warping that could occur under these extreme conditions. In addition, this material has a thermal deflection temperature (HDT) that exceeds the temperature range used by the equipment. This means that it can withstand high temperatures without suffering significant deformation.

During this study, both ABS and CF composing the composite material being studied, will be taken into account separately in the microscopic modeling. Their main characteristics studied previously by Somireddy et al. [26, 27] are shown in the Tables 2 and 3, where \(\rho \) is the density, E is the Young’s modulus, \(\nu \) is the Poisson’s ratio, \(\sigma _Y\) is the yield stress, \(\delta \) is the hardening exponent, and H is the linear hardening modulus in the isotropic hardening law. The evolution with temperature of the main characteristics of the material is shown in Fig. 1.

Table 2 Main properties of ABS in the matrix of the composite material [26, 29]
Full size table

Large-format 3D printer: Super Discovery 3D

The LFAM machine used to conduct the study is the Super Discovery 3D Printer (by CNC Barcenas) [30], an industrial 3D printer with advanced features that uses direct extrusion of pellets. It has an unlimited supply of material thanks to its automatic feeder, which allows printing without restrictions on the quantity or weight of the part. It supports a wide range of thermoplastic materials and offers working dimensions of 1300x2500x1000 mm, which can be customized according to needs. In Fig. 2 this printer can be shown.

Table 3 Main properties of short CF inclusions in the composite material [27]
Full size table
Fig. 1
figure 1

Composite material (ABS/CF20) properties evolution with temperature

Full size image

The printing bed is formed by a closed chamber with a heated bed, capable of reaching temperatures of up to 150ºC and has an auto-leveling system to compensate for unevenness during printing. The extrusion head can heat materials up to 400ºC. The Super Discovery 3D Printer represents an evolution in 3D printing, as it uses direct extrusion of pellets instead of traditional FFF technology. This results in advantages such as increased productivity in larger formats, faster printing speeds, and significant cost savings.

Methods

Demonstrator manufacturing

As stated before, throughout this work, the demonstrator used to study the deformations during the LFAM process is a mold to manufacture the external expansion aerodynamic package of a FS car. This mold was manufactured using the LFAM machine described previously, the 3D Super Discovery by CNC Barcenas. The printing parameters used are stated in Table 4. The manufacturing process of this component is illustrated in Fig. 3, after the 3D-printing using a LFAM machine, the mold is polished to obtain an optimal surface finishing.

Experimental testing

The composite material principal mechanical properties given by the supplier, summarized in Table 1, are used to validate the micro-structural modeling of the material. Therefore, to verify these properties, a series of tensile tests are designed, as described in previous studies [13].

Fig. 2
figure 2

Super Discovery 3D printer by CNC Barcenas. Main elements are displayed

Full size image

The specimens were produced according to the ASTM D638 Type I standards [31], as shown in Fig. 4(a). They were manufactured with a CNC machine from a 3D printed plate. The tensile tests were performed with a universal mechanical testing machine (Instron 6800 series as shown in Fig. 4(b)), with a displacement set to 5 mm/min until failure, as conducted by the supplier. To ensure the reliability of the test, five specimens were tested under the same conditions. The average results were taken and they can be seen in Fig. 12 [13], summarized in Tables 6 and 7 [13], in the results and discussion section.

Microstructure characterization

The first step to carry out the characterization of the microstructure is the acquisition of data. Due to confidentiality reasons, the supplier does not provide details about the exact composition of the material, hence the microstructure arrangement should be determined. LNP™ THERMOCOMP™ AM Compound AC004XXAR1 is composed by two elements, ABS and 20% short CF, with unknown fiber dimensions and orientation, needed to understand the material mechanical properties.

Gemini Ultra 55 scanning electron microscope (SEM) is used to take high-resolution images of the microstructure of samples of the composite material. The SEM images of interest for the study are contained in Fig. 5. Using these images, together with an image processing software called Fiji, it is possible to determine various statistical data about the length, diameter and orientation of the short CF inclusions in different planes. Specifically, Fig. 5(a) and (b), are used to obtain the values of length, diameter, and orientation of the fibers in the plane parallel to the printing direction. On the other hand, from Fig. 5(c) and (d), the fiber orientation in the plane perpendicular to the printing direction is studied. Fibers and pores could be clearly identified. In Fig. 5(a), it is possible to observe a section of the sample in which the fibers are clearly visible, which will therefore facilitate their subsequent study. Figure 5(b) is a magnification of Fig. 5(a), which allows a greater measurement accuracy.

On their side, Fig. 5(c) shows a tensile fracture area, whereas, Fig. 5(d) is a magnification of this region. These images shows the porosity in the ABS matrix, as well as the holes left by fibers due to the fracture event. Taking a general view, two phases can be distinguished: matrix (ABS) and reinforcements (short CFs).

Table 4 3D printing parameters in the manufacturing of the external expansion mold
Full size table
Fig. 3
figure 3

Manufacturing steps of the external expansion aerodynamic part of a FS car. (1) 3D printing in ABS/20CF. (2) Thin sanding and polishing process of the 3D printed part to obtain the final mold. (3) Lay-up process of CFRP prepreg. (4) Assembly with the rest of parts confirming the lateral aerodynamics package. (5) Assembly into the FS car, beside the cockpit

Full size image
Fig. 4
figure 4

(a) Dimensions (in mm) of the tensile test specimen ASTM D638 Type I [31]. Red lines show the XZ printing direction [13]. (b) Fractured specimen after tensile test using a universal testing machine

Full size image

To determine the elastic properties of the material, it is crucial to establish the orientation tensor of the fibers. Measuring the spatial orientation of fibers relies on two fundamental principles: (1) Fibers with infinite length always possess elliptical cross-sections, and the ratio of the two principal axes of the ellipse depends on the angle \(\beta \) of inclination relative to the cross-sectional plane. (2) The orientation angle \(\phi \) within the cross-sectional plane corresponds to the alignment of the longer principal axis of the cross-section.

The inclination angle b of the fiber can be expressed by Eq. 1.

$$\begin{aligned} \sin \beta =\frac{a}{b} \end{aligned}$$
(1)

On Fig. 6 an scheme of the different parameters mentioned above is shown. Once the values of the different orientation angles are obtained, the values of \(a_{xx}\), \(a_{yy}\) and \(a_{zz}\) can be obtained from the Eq. 2; where b represents the length of the fiber, a specifies the value of the fiber diameter, \(a_{ij}\) represents the value of the fiber tensor of the material per unit area, N is the number of fibers per unit area, and \(\beta _i\) and \(\phi _i\) are the orientation angle of each of the fibers.

$$\begin{aligned} \begin{aligned}&a_{x x}=a_{11}=\frac{1}{N} \sum _{i=1}^N \sin ^2 \beta _i \\&a_{y y}=a_{22}=\frac{1}{N} \sum _{i=1}^N \cos ^2 \beta _i \sin ^2 \varphi _i \\&a_{z z}=a_{33}=\frac{1}{N} \sum _{i=1}^N \cos ^2 \beta _i \cos ^2 \varphi _i \end{aligned} \end{aligned}$$
(2)
Fig. 5
figure 5

Sample SEM images at two different regions of the ABS/20CF test samples

Full size image
Fig. 6
figure 6

Sketch of the dimensional parameters of each fiber [22]

Full size image

Fiji is used to post-process the images to obtain the mentioned parameters. As noticed in Fig. 5, many of the observed fibers are partially embedded in the ABS matrix, hence they can not be used for the measurements as their full length remains unknown. With different plugins and modules in Fiji, orientations in both planes are measured. After analyzing several images, the mean values of the parameters of interest are contained in Table 5.

Mean-field homogenization approach

ABS reinforced with 20 % CF is considered to be an heterogeneous composite material, as it is composed of a matrix (ABS) and multiple phases (reinforcements of short CF). To predict the interaction between the microstructure and the macroscopic properties, a micro-macro (two-scale) approach is considered to perform the micromechanical modeling of the material of interest.

Table 5 Short CF reinforcements mean dimensions and orientations. Parameters introduced into Digimat-FE regarding the inclusion phase
Full size table

In this approach, two scales should be distinguished: the microscopic (heterogeneous) and the macroscopic (locally homogeneous). Both scales are linked by a representative volume element (RVE), sufficiently large to represent the properties of heterogeneous microstructure, but small enough to represent an infinitesimal part of the macroscopic solid.

The RVE, is generated by separately defining the phases materials and dimensional properties. The ABS matrix is defined as an elastoplastic material modeled with the J2-plasticity model based on the Von-Mises equivalent stress. On the other hand the CF reinforcement is defined as a linear elastic material with an isotropic behavior, as the mechanical properties are considered to be equal in all directions. The material properties of ABS and short CF inserted into Digimat are summarized in Tables 2 and 3 respectively [27,28,29]. Then, the microstructure is modeled as in previous studies conducted [13], considering the parameters of the fibers obtained after the microstructure characterization, summarized in Table 5. Following, the formulation used by Digimat in the homogenization methods is reviewed.

The RVE contains a finite number of constituents, each represented by a constitutive model. An homogenization method is used to find the macro constitutive response of RVE, where continuum mechanics is applied at the macro scale with macro constitutive equations [32]. This two-scale approach allows the consideration of microstructural effects on the macroscopic response of the solid body, without the need to solve the problem at the micro-scale for the whole solid body, which would be computationally expensive.

However, the major difficulty in the two-scale approach is to solve the RVE problems. The classical analysis of the mechanics of continuous media is performed at the macroscopic level. In the domain \(\int _{\omega }\) and considering volume \(\textrm{V}\), at each macro point \(\textbf{X}\), the macro strain \(\mathbf {\varepsilon }(\textbf{X})\) is known and the macro stress \(\sigma (\textbf{X})\) must be calculated, or vice-versa.

The average quantity over a RVE is defined by Eq. 3, where integration is performed with respect to micro coordinates, considering \(f(\textbf{X}, \textbf{x})\) as the RVE micro-field.

$$\begin{aligned} \langle \textrm{f}(\textbf{X}, \textbf{x})\rangle =\frac{1}{\mathrm {~V}} \int _{\omega } \textrm{f}(\textbf{X}, \textbf{x}) \textrm{dV} \end{aligned}$$
(3)

Hence, the homogenization problem can be expressed in terms of continuum mechanics as follows: once the macro strain is known at each macromaterial location, the macro stress could be computed, and vice versa. The fundamental problem can be stated in a simpler form (Eq. 4), in terms of macro stiffness, based on linear elastic principles.

$$\begin{aligned} \langle \sigma \rangle =\overline{\textrm{C}}:\langle \varepsilon \rangle \end{aligned}$$
(4)

Basically, finding a comparable homogeneous material with the same effective macro stiffness as the actual heterogeneous composite under the same boundary conditions is the core challenge of homogenization in linear elasticity. In the case of the software being used in this work, Digimat 2023.3, it employs two different approaches to solve the fundamental problem: direct finite element analysis (FEA) or mean-field homogenization (MFH) models.

FEA method is very general and accurate, however, with complex geometries and non-linear problems, the CPU time and memory usage can become unaffordable. Hence, MFH models methods could become relevant due to its usage ease and high speed. Note that this approach is based on semi-analytical models, hence it will only give (accurate) approximations to the volume averages of stresses and strains.

One of the most common methods for carrying out this homogenization is the use of the Mori-Tanaka equations [33]. These set of equations represent the relationship between the effective properties of the composite material and the properties of the individual phases and the material geometry.

Therefore, to perform the homogenization, Mori-Tanaka equations divide into 3 parts. First, these equations relate the effective stresses in the composite material to the stresses in the individual phases and the material geometry. Then, the elastic and thermal strains are calculated. Finally, everything is combined into Eq. 5, that enables the calculation of the final effective stress.

$$\begin{aligned} \sigma _{i j}^*=\sigma _{i j}^m+\phi \cdot \frac{E_m}{E_f}\left( \varepsilon _{i j}^f-\varepsilon _{i j}^m\right) +\Delta \alpha \cdot \left( T-T_0\right) \cdot \delta _{i j} \end{aligned}$$
(5)

Equation 5 allows the computation of effective stresses in a given composite material, given the properties of the individual phases, the geometry, and the load and temperature conditions. The effective stress (\( \sigma _{ij}^* \)) in the composite material is composed of several terms. First, (\( \sigma _{ij}^m \)) represents the stress in the matrix of the composite material, while (\( \phi \)) is the volume fraction occupied by the reinforcing phase. The factor (\( \frac{E_m}{E_f} \)) compares the elastic moduli of the matrix (\( E_m \)) and the reinforcing phase (\( E_f \)), reflecting how the stiffness of the reinforcing phase affects the deformation in the composite material. The difference between the deformations in the reinforcing phase (\( \varepsilon _{ij}^f \)) and the matrix (\( \varepsilon _{ij}^m \)) is denoted as (\( \varepsilon _{ij}^f - \varepsilon _{ij}^m \)), which influences the stress distribution. Finally, the term (\( \Delta \alpha \cdot (T - T_0) \cdot \delta _{ij} \)) considers the thermal strain due to differences in the thermal expansion coefficients between the reinforcing phase and the matrix, adjusted for the current temperature (\( T \)), the reference temperature (\( T_0 \)), and the Kronecker delta (\( \delta _{ij} \)). In total, this equation provides a complete relationship of how the effective stresses in a composite material depend on the stresses in the individual phases, the volume fraction of the reinforcing phase, and the effects of elastic and thermal deformation.

It is important to note that these equations are simplifications, and that there are several variants and more complex MFH models depending on the specific characteristics of the composite material and the load conditions.

In the case study presented, a second-order homogenization method is applied. This involves both, the first and second statistical moments of the strain field for each phase. The second moment correlates with the variance, providing a richer statistical context than just using the mean value. As a result, it is anticipated that more accurate predictions will be made with second-order MFH compared to its first-order counterpart.

To better understand and justify the idea of second-order MFH, the paper examines a composite material composed of elastic phases (r), each characterized by a unique strain energy function \(W_r(\epsilon (\textbf{x}))\).

Around a reference strain value, uniform in each phase, the second-order Taylor expansion of \(W_r(\epsilon (\textbf{x}))\) is performed, given in Eq. 6.

$$\begin{aligned} \textrm{W}_{\textrm{r}}(\varepsilon (\textrm{x})),\approx & {} \textrm{W}_{\textrm{r}}\left( \tilde{\varepsilon }_{\textrm{r}}\right) +\frac{\partial \textrm{W}_{\textrm{r}}}{\partial \varepsilon }\left( \tilde{\varepsilon }_{\textrm{r}}\right) :\left( \varepsilon (\textrm{x})-\tilde{\varepsilon }_{\textrm{r}}\right) \nonumber \\{} & {} +\frac{1}{2} \frac{\partial ^2 \mathrm {~W}_{\textrm{r}}}{\partial \varepsilon \partial \varepsilon }\left( \tilde{\varepsilon }_{\textrm{r}}\right) ::[\left( \varepsilon (\textrm{x})-\tilde{\varepsilon }_{\textrm{r}}\right) \nonumber \\{} & {} \otimes \left( \varepsilon (\textrm{x})-\tilde{\varepsilon }_{\textrm{r}}\right) ] \end{aligned}$$
(6)

Computing the per-phase volume average of the strain energy function, then averaging over all-phases, the macro strain energy of the composite RVE is obtained by Eq. 7.

$$\begin{aligned} \overline{\textrm{W}}(\textrm{E})=\sum _{\textrm{r}} \textrm{v}_{\textrm{r}}\left\langle \textrm{W}_{\textrm{r}}(\varepsilon (\textrm{x}))\right\rangle _{\omega _{\textrm{r}}} \end{aligned}$$
(7)

The second-order theory is demonstrated to be significantly more accurate than the first-order formulation [34] when the following conditions are met within the composite material being tested: fiber reinforcement, high stiffness contrast between fibers and matrix, and little hardening exhibited by the elasto-plastic matrix. These three conditions are satisfied by the composite material being tested in this paper.

A numerical simulation is performed applying the previously explained approach, used by Digimat-MF solver. The results obtained are exposed and discussed in “Resultsand discussion”.

Fig. 7
figure 7

Single-layer microstructure RVE genertaed in Digimat

Full size image
Fig. 8
figure 8

Optic microscope image from the cross-section of a sample test of ABS/CF20

Full size image
Fig. 9
figure 9

Two-layer microstructure RVE generated in Digimat

Full size image

Finite element approach

The second approach proposed to analyse the mechanical properties of a composite material at a micro-level is throughout FEA. This approach is divided in two different cases: (1) considering the microstructure to be completely homogeneous and (2) considering the microstructure to be multi-layer. Again, the materials constituting the composite material - ABS and CF - are modeled separately, alongside the microstructure, as for the MFH method case. Figure 7 shows the single layer RVE generated for the first case representing an homogeneous microstructure. The method followed has been studied in past papers [13].

However, observing the images taken with the optic microscope, e.g. Fig. 8, it was appreciated a difference at the boundary between printing layers. This region supposed approximately 1/4 of the total layer width. Here, it is noticed how the majority of fibers are spread randomly, differently to the rest of fibers along the layer, which are mainly aligned in the printing direction. This specific feature could be modeled in Digimat by performing a multi-layer analysis.

Fig. 10
figure 10

Design of the mold for both external expansion aerodynamic packages of the FS car. Geometry introduced into the 3D printer

Full size image
Fig. 11
figure 11

Flow diagram to summarize the methodolgy proposed

Full size image

A two-layer RVE is generated, containing each layer a specific representative microstructure. In the first layer, the matrix and inclusions characteristics are exactly the same as for case (1), except for the orientation tensor. In this specific case, the fibers are set to be distributed randomly, hence the diagonal term in the orientation tensor is equal to 1/3. Regarding the second layer, the microstructure is exactly the same as for the single-layer case (1). Figure 9 shows the multi-layer RVE generated.

Finally, a voxeled mesh is generated representing the microstructure for each case. The CASI iterative solver is selected in the software as it allows calculating the solution of large systems at a lower computation cost compared to the direct or iterative solvers. The results obtained by the previously explained method are exposed and discussed in the following section of results, comparing the single and multi-layer solutions.

AM numerical simulation

Furthermore, the AM process of the demonstrator suggested could be simulated by numerical means using Digimat-AM tool. By introducing the geometry, in Fig. 10, as an STL file and the G-code - that contains the location history of the nozzle at each instant of time and the boundary conditions set prior to the AM process - the solver is able to reproduce the 3D-printing process of the component. This paper evaluates the usefulness of Digimat-AM tool to reproduce the LFAM process.

The Digimat-AM simulation approach could be basically described as multiscale workflow where a highly coupled thermomechanical analysis is performed at the microscale of the material, prior to a mechanical FE analysis at the macroscale. As a result, the computational cost is reduced as the macro part analysis uses the warpage behavior identified at the micro-level.

In this work, the temporal discretization strategy selected is layer-by-layer, where finite element layers are activated. This discretization method reduces the computational cost compared to the filament discretization, where chunks of material are added along each layer.

Regarding the element type selected to undergo the numerical simulations, fully integrated 3D quad elements - or voxels - are generated. Due to the anisotropic nature of the reinforcement CF in the studied composite material, full integration elements are more appropriate to avoid oscillations in the displacements calculated and inaccuracy in the results close to singularities in the geometry.

In this context, a fine costumed voxel mesh is generated, with a mesh size of 3 mm, corresponding to 2 times the real printing layer thickness and 1/2 the printing width. In total 134,000 voxels conform the mesh. Even though this finer mesh increases the computational time, it results in a more accurate prediction of the residual stresses generated during the 3D-printing process.

The material set for this simulation is modeled as an amorphous composite material. Once the material properties have been replicated up to a certain accuracy in the MF and FE Digimat tools, this composite with the same microstructure properties is introduced into Digimat-AM, resulting in the temperature dependant properties shown in Fig. 1.

Fig. 12
figure 12

Stress-strain curve for the composite material ABS/CF20, resulting from the experimental tensile tests and numerical simulations

Full size image

A high-fidelity warpage simulation is performed over the geometry of the mold. In processes involving filament deposition, this predominant simulation strategy in Digimat-AM is based on the inherent strain method and the execution of comprehensive transient simulations. The values of inherent strain, represented by strain tensors, depict the warpage characteristics of the material, meaning they illustrate how the material contracts and induces distortion post-deposition under specific process conditions.

When the inherent strain values are established for a certain combination of material attributes and process variables, a rudimentary mechanical simulation can be conducted, layer by layer, to emulate the manufacturing process. These inherent strain values are incorporated at each layer as it is constructed until the entire part is complete.

Beyond the standard inherent strain method offered in Digimat-AM for the low-fidelity strategy, comprehensive transient simulations are executed in the high-fidelity strategy. These simulations serve to analyze material contraction and the resultant distortion following deposition under specific process conditions, while also precisely calculating the complete temperature distribution. The transient simulations utilize time discretization, where elements are sequentially incorporated until the entire component is constructed.

The flow diagram in Fig. 11 summarizes the whole methodology proposed in this work. It shows the traditional approach for the LFAM of components and, on the other hand, the approach proposed to aid in the build of a DT that could reduce the material waste in the 3D printing process.

Results and discussion

In this section, the results of the numerical simulations at the micro and macroscale are presented and discussed, comparing the numerical approaches between them and with the results obtained experimentally.

Firstly, the findings from the tensile test are depicted in Fig. 12, indicating a behavior consistent with elastoplastic materials. The key data points are outlined in Table 6. It is crucial to acknowledge the initial toe region, a byproduct of the specimen’s alignment and takeup, which does not correspond to the material inherent properties. Therefore, this toe is disregarded, and the point of zero strain is adjusted. The results display an initial elastic phase where stress correlates directly with strain. Subsequent to the yield stress (\(\sigma _Y\)), the composite material exhibits plasticity attributed to isotropic hardening at a stable rate. The Young’s modulus (E) and hardening modulus (H) are inferred from Fig. 12. Lastly, the specimen’s average failure point is identified at the ultimate tensile stress (\(\sigma _U\)). The summarized key points and associated data - obtained graphically - are listed in Tables 6 and 7.

Table 6 Summary of main properties of the composite material ABS with 20% CF, obtained from experimental testing and numerical simulation
Full size table
Table 7 Summary of average mechanical properties at failure of the composite material ABS with 20% CF
Full size table
Fig. 13
figure 13

Warpage (mm) field due to the LFAM of the external expansion mold

Full size image

After plotting the results from the different approaches, it is noticeable the scale difference between the experimental and the numerical simulation results. Therefore, the team considered appropriate the application of a correction factor to the experimental stress results. This correction factor of 1.8, accounts for all the errors and uncertainties arising from the experimental method undertaken. During mechanical testing, there are many factors that cannot be controlled in the performance of ideal tests, for example, the atmospheric conditions under which the tests are performed or the ideal printing conditions of the specimens. One of the factors that affected the experimental behavior of the specimens was the high porosity level, observed in Fig. 5, which reduced drastically the mechanical properties of the material [35]. A possible solution would be a greater dehumidification time of the raw material, prior to its melting and extrusion.

Focusing on the experimental corrected results, highlight that the error with respect to the FE single-layer analysis is just 1% regarding the plastic region, as indicated by H value in Table 6. However, note in Fig. 12 the similarity of the elastic region with respect to the FE multi-layer analysis, showing a relatively low error of 7.8%. Therefore, the significant similarity of the curves, suggesting an elastoplastic behavior under tensile loading, and the relative low error regarding the mechanical properties - as stated in Table 7 - supports the applicability of the suggested numerical model as a tool to develop a DT to represent the mechanical behavior of the composite material studied.

Not only the corrected experimental results are closer to those obtained numerically, but they are also closer to the key values provided by the supplier in the revised technical data sheet, summarized in Table 1. The process difference during the manufacturing of the specimens could explain this error outcome. Due to equipment limitations, replicating the preparation of specimens was a challenge. A difference in the printing conditions most probably meant a poor layer adhesion, thus worse mechanical properties obtained experimentally.

Fig. 14
figure 14

Von Mises stress (MPa) field due to the LFAM of the external expansion mold

Full size image

Considering both numerical approaches, even though the CPU time taken to obtain the results with the MF second order homogenization is significantly low, in the order of 10,000 times, the FE approach is significantly more accurate. As the FE module considers more parameters and computations are made voxel by voxel, results have a greater accuracy with a high - but not restricting - computational time. Hence, the FE approach is considered to be more appropriate in this applicable case.

Moreover, reviewing the results from the high-fidelity warpage analysis at the macroscale, the contour plot for warpage reveals that the predominant deformation is centralized in the part’s midpoint, proximate to the area of maximal curvature of the mold. The thermal gradient between successive layer depositions results in warpage. Figure 13 illustrates the deformation along the layers of deposition, underscoring how the region of maximum deformation consistently aligns with the area of maximal curvature, culminating in a deformation of 8.51 mm in the finished product. This warpage resulted is barely significant considering the large dimensions of the printed mold. Actually there is an error of 6.7% in the maximum displacement point with respect to the real printed part. Given the composite material minimal thermal expansion coefficient utilized in the creation of the investigated part, alterations in dimension due to temperature variations are less probable.

The Von Mises stress contour plot in Fig. 14 displays the residual thermal stresses induced during the printing, originating from the thermal gradients developed between layers. The elevated temperature layers are prone to expansion while being constrained by the cooler ones, causing compressive forces on the warmer upper layers and tensile stresses on the restraining cooler ones, leaving inherent thermal stresses within the structure. During AM, regions with significant temperature gradients typically harbor concentrated residual stresses, often located on the external surfaces of the part, where rapid heat loss occurs due to exposure to atmospheric air temperature. Furthermore, areas undergoing changes in material thickness or geometry are also susceptible to residual stresses due to the variance in cooling rates and temperatures relative to neighboring material. Figure 14 highlights these regions, demonstrating the concentration of significant residual stresses at the vertical extremities of the part. These distribution results could be corroborated with different examples in a previous study [13].

Conclusion

In this manuscript, an innovative material/process modeling technology has been presented, aiming to streamline the development of composite materials molds. This technology facilitates the simulation of the LFAM processes through the fused deposition of pellets, utilizing specialized AM process simulation software, Digimat. However, this software operates optimally under specific boundary conditions and does not accommodate material deposition in out-of-plane trajectories. The early-stage development of a DT proves to be challenging, given the need to manage substantial data volumes, necessitating resource-intensive tools both in time and economical terms.

Throughout this study, experimental test were conducted to determine the mechanical properties of ABS with 20% short CF inclusions. Additionally, a numerical model, representing the microstructure of the composite material, was designed throughout two different approach strategies: (1) involving MFH schemes and (2) using FEA. This numerical model was constructed as an initial step in DT development and corroborated through various experimental outcomes. Prior to generating the RVE, the microstructure of the composite is characterized by semi-automatic means. This showed that fibers were mainly aligned with the material extrusion direction.

The comparison between stress-strain curves presented an identical elastoplastic behavior but with a significant discrepancy between the experimental and numerical models. Hence, a correction factor of 1.8 is applied to the experimental outcome to compensate for any errors induced during the manufacturing of specimens and the mechanical testing itself. As a result, the error of the FE approach in the elastic domain is drastically reduced and is closer to the actual data values provided by the supplier. On the other hand, the MFH approach has a relative small error in the elastic region however is unable to represent the plastic properties. Hence, the FE approach is concluded to be useful in the early-stage development of a DT.

On the other hand, results at the macroscale show how the LFAM process could be accurately reproduced by Digimat-AM, confirmed as a useful tool to develop a DT. The warpage obtained is relatively small compared to the dimensions of the mold being manufactured, with an error with respect to the actual printed part of 6.8% regarding the maximum deflection point.

Some of the limitations experienced using Digimat while conducting the study includes the particular material modeling (of ABS and CF) as the material models supported by the software cannot fully represent the mechanical behavior of these materials. In addition, the AM simulations are computationally costly, thus machine learning techniques are proposed to handle large experimental databases and real-time data that will enrich the DT. Probably, this will result in reduced CPU times and greater accuracy in the representation of the composite material mechanical properties and the LFAM process.