Transient modeling of material extrusion by system identification

doi:10.21203/rs.3.rs-3443933/v1

Download PDF

Research Article

Transient modeling of material extrusion by system identification

https://doi.org/10.21203/rs.3.rs-3443933/v1

This work is licensed under a CC BY 4.0 License

Version 1

posted

You are reading this latest preprint version

Material extrusion is popular for its low barriers to entry and the flexibility it gives designers relative to traditional manufacturing techniques. Material extrusion is a transient process with a high frequency of starts, stops, and accelerations. This work presents transient data collected by an instrumented printhead and models the data by way of system identification. First-order and second-order control system models are proposed. The work also includes principal component analysis to determine which model coefficients correlate with the main effect, models the first-order model coefficients as a function of the experimental factors by regression, and predicts the apparent viscosity using a fitted static gain and known parameters. Flow rate, hot end temperature, nozzle diameter, and acceleration are the factors selected for the experiment. Each of these factors influences the steady state pressure, except for acceleration. The system identification models predict the melt pressure’s transient behavior well, with standard errors less than 4% of the mean melt pressure. Statistical analysis of the first-order model coefficients verifies that the static gain and time constant are statistically significant responses of the factors. The modeled apparent viscosity follows rheological expectations, showing the trends typically seen for viscosity as a function of shear rate.

Fused deposition modeling

fused filament fabrication

instrumentation

rheology

Transient behavior of the motor torque, infeed load, nozzle melt pressure, and nozzle melt temperature.
System identification of the nozzle melt pressure as a function of the set flow rate, temperature, and nozzle orifice diameter.
Regression and principal component analysis of the system identification model coefficients.
Estimation of the apparent viscosity in the nozzle bore using the system identification model’s static gain.

In material extrusion additive manufacturing of thermoplastics, a solid filament is driven into the hot end at a commanded volumetric flow rate via a gear mounted on the shaft of a stepper motor. This mechanism leads to a dependency chain where the motor torque drives the infeed force with which the filament is driven into the hot end, and the infeed force drives the melt pressure within the hot end. The melt pressure within the hot end will be the primary focus of this work, since it is the most reflective of the state in the hot end based on the authors’ prior work [1–3].

Material extrusion is of interest because of its low barrier-to-entry for manufacturers and the resulting popularity in the additive manufacturing community [4]. This process is readily accessible for manufacturers because of the simplicity and variety of the feedstock [5–9] and the wide availability of desktop material extrusion systems [10], providing the users with competitive design flexibility for their products. Thermoplastic material extrusion finds use in biomedical [11, 12], automotive [13], and consumer applications [14].

With material extrusion’s increasing use in manufacturing [15], it is critical to improve observability and modeling efforts to maximize its performance for end users. Based on current understanding, the thermoplastic material extrusion’s system performance is strongly affected by the melting capacity of the hot ends available [16–18], the slicing parameters used to generate the G-code for printing [19], a lack of observability and control [20], and the transient behavior of the thermoplastic in the hot end and its deposition [21]. It is critical that researchers investigate these areas, since they affect the properties of the final part [22, 23], and their prediction [23], and the accuracy of the part’s dimensions to the user’s design [24].

Modeling the thermal behavior within the hot end is an ongoing area of investigation [25–27]. Most of the models are at steady state [28], although there are some transient models reported in the literature [21, 29]. One predicts nozzle clogging based on the current to the stepper motor [29]. Other researchers have modeled the initial melting stages of the thermoplastic in the hot end using finite element and computational fluid dynamics simulations [21, 30]. To validate these models, an instrumented printhead is usually implemented [31–33]. Several instrumented systems can also be found within the literature, which employ an infeed force sensor [33], melt pressure sensors [31, 34], and thermocouples [31, 34] to monitor the state of the thermoplastic and increase observability. These systems provide an understanding of the rheological state within the hot end, which is critical for stable and repeatable processing [35–37].

System identification is a technique that enables analysis of systems [38, 39], allowing them to be modeled using transient input and output data. Transient modeling by system identification is of significant importance because there are step changes in material extrusion rates (related to stepper motor actuation) with a dynamic pressure response that is affected by the thermoplastic’s properties as well as the hot end and nozzle design [2]. System identification has been performed on roads deposited by material extrusion for control purposes, but not on signals collected during processing [40]. Transient behaviors have implications on the methods used to characterize materials and could affect the measurement of contact pressures [3, 41], die-swell [42, 43] and compressibility [2], which influence the welding between the layers of deposited material [22] and their geometry. Understanding the transient behaviors present during extrusion will also provide support to efforts for sustainability [15, 44], since designers can use these techniques to evaluate their systems relative to others, and to advanced processes that leverage standard material extrusion systems [45].

This paper focuses on applying system identification to data from an instrumented material extrusion printhead, with the intent of modeling the transient response to the start of extrusion. The experimental apparatus is designed to measure motor torque, infeed force, melt pressure, and melt temperature. The response (output) of interest for this work is the melt pressure and the input of interest is the volumetric flow rate. This work looks at two types of system identification models. The simpler model of the two is a first-order model [46, 47] with one pole as follows:

$${M}_{1}\left(s\right)=\frac{{K}_{p}}{1+{T}_{p}s}$$

where ${K}_{p}$ is the static gain, ${T}_{p}$ is the time constant for the pole, and $s$ is the frequency. The second-order model [48] is as follows:

$${M}_{2}\left(s\right)=\frac{{K}_{p}\left(1+{T}_{z}s\right)}{1+2\zeta {T}_{w}s+{\left({T}_{w}s\right)}^{2}}exp\left(-{T}_{d}s\right)$$

where ${T}_{z}$ is the time constant for a process zero, ${T}_{d}$ is the time delay constant, $\zeta$ is the system damping ratio, and ${T}_{w}$ is the time constant. To demonstrate the effect of the model coefficients, Fig. 1 shows the second-order model with varying values of $\zeta$and ${T}_{z}$ with other model parameters (${K}_{p}$ and ${T}_{w}$) held constant as they follow first order behavior.

The effect of the varying experimental factors on the coefficients of both models is investigated in terms of their percent change from their average values across a one-factor-at-a-time (OFAT) experimental design. Principal component analysis (PCA) is used to determine which factors drive the model behavior. The coefficients for the first-order model are also modeled using regression. Lastly, the work estimates the apparent viscosity of the thermoplastic in the hot end by way of the static gain from the first-order model, the known flow rate, and Poiseuille’s equation [49].

2.1 Material

For these experiments, an acrylonitrile butadiene styrene filament (Hatchbox ABS, black) with a diameter of 1.75 mm is used. The material has a recommended processing temperature range of 210°C to 240°C. ABS is used because it is a popular material in both practice and research [8, 50].

2.2 Experimental hardware

The instrumented printhead (see Fig. 2) is mounted on a Creality (Shenzhen, China) Ender 5 Pro. A Micro Swiss (Minneapolis, Minnesota) direct drive extruder replaces the stock extruder. The 0.5 N∙m torque sensor (Dongguan NanLi Co. Ltd., Guangdong, China) and the infeed and melt sensors (SparkFun Electronics, Niwot, CO, Part Number: SEN-14729, combined error of 0.05% of the full-scale value) are included by further modification. To mount the infeed and melt sensors as well as the IR array (HTPA32x32dR2-L5.0/0.85F7.7HiS, Heimann Sensor GMBH, Eltville, Germany), an aluminide plate from i.Materialise (Leuven, Belgium) is used in place of the original Micro Swiss plate. The infeed sensor is attached to the heat sink by an aluminide bracket. Rotation is transmitted from the stepper motor to the torque sensor shaft by way of a timing belt and pulley system.

The system acquires the motor shaft torque, the infeed load, and melt pressure data by the sensors sending voltages, in millivolts, through signal conditioners (OMEGA, model: DRF-LC, Norwalk, CT), which convert the signal to volts, effectively amplifying the signal. After amplification, the signal for each sensor is input to a channel on the data acquisition system, a National Instruments (Austin, TX) NI USB-6210, which is connected to the computer. MATLAB (MathWorks, Natick, MA) is used to control the data acquisition. The sampling rate is 50 kHz with a buffering frequency of 250 Hz. Upon conclusion of data acquisition, the MATLAB script outputs a comma-separated value file for storage and post-processing. The infrared array’s data enters the computer through a separate ethernet cable and uses its own proprietary software. The data is output from the infrared camera software as a text file and post-processed with a MATLAB script. For the purposes of this work, the melt pressure will be the main focus.

The hot end is manufactured in bronze by i.Materialise using a casting process and is a custom design. To enhance the melting of the thermoplastic, the circular inlet lofts from a diameter of 2 mm to a thin slit having a width of 3 mm and a thickness of 1 mm. The thinness of the slit increases the rate of heat transfer. Its wideness allows for the use of a larger diameter load pin, increasing the signal-to-noise ratio for melt pressure measurement relative to Coogan and Kazmer [31]. This geometry also allows for the addition of a window opposing the load pin for viewing the melt. In this slit section, a zinc selenide (ZnSe) lens is held in place by a set screw with a through-hole, allowing the infrared (IR) radiation from the melt to pass through to the IR array. Brass nozzles with varying orifice sizes from McMaster-Carr are used in this work (part numbers 3695N301, 3695N304, 3695N306, Elmhurst, IL, USA). Double the heating power of a typical hot end is generated by two 40 W heater cartridges in this hot end design. A 100 kΩ thermistor is used for hot end temperature control. A pin with a diameter of 3.02 mm transmits the applied stress from the melt inside the hot end to the beam-style load cell. The contact surface on the load cell, to transmit the resulting force to the load cell, is a nylon hex bolt, to minimize heat transfer.

2.3 Experimental factors

The four factors that vary in these experiments are flow rate, hot end temperature, nozzle diameter, and acceleration. Flow rate is the commanded volumetric flow rate at which material is transported into the hot end by the stepper motor. The hot end temperature is the temperature that the controller of the printer is set to target based on the readings of the thermistor in the hot end. The nozzle diameter is the diameter of the nozzle orifice that the thermoplastic passes through immediately before deposition. Acceleration is the rate at which the feed rate used to transport the filament into the hot end changes. The full experimental matrix can be seen in Table 1. This experiment design varies one-factor-at-a-time analysis to provide assessment of the second-order response for each factor independent of the other factors.

Table 1

The design of experiments
Run (description)	Volumetric flow rate [mm³/s]	Temperature [°C]	Nozzle diameter [mm]	Acceleration [mm/s²]
1 (center point)	5	225	0.40	500
2 (low flow rate)	1	225	0.40	500
3 (high flow rate)	10	225	0.40	500
4 (low temperature)	5	200	0.40	500
5 (high temperature)	5	250	0.40	500
6 (small orifice)	5	225	0.25	500
7 (large orifice)	5	225	0.60	500
8 (low acceleration)	5	225	0.40	100
9 (high acceleration)	5	225	0.40	2500

Each experiment starts by allowing the hot end temperature to equilibrate for 2 minutes. Then, for two minutes, a creeping flow of 0.6 mm³/s is performed to ensure a steady process without material degradation. The process then dwells for a minute, without filament retraction. During the experimental flow rate, the material is extruded 75 mm above the print bed for three minutes. Throughout this whole experimental time, the data acquisition system collects data from the torque, infeed, and melt sensors, and the IR array collects images.

2.4 System identification

To perform system identification, the data from the experiments is first loaded into MATLAB. After loading the experimental data, the signal for the volumetric flow rate is generated as a step function that starts at a value of zero and then rises to a value equal to the prescribed volumetric flow rate for that run. Next, a continuous time process model structure is defined with the volumetric flow rate as the input and the observed melt pressure in the hot end as the output using the MATLAB command ‘idproc()’. The models are then fit to the data by a prediction-error minimization using the MATLAB function ‘pem()’. Once the solution converges, the models can be compared, and the final coefficient values are exported for each of the experimental runs of Table 1.

To calculate the pressure from the load cell, the measured voltage is first converted into a force value by multiplying it by a gain value, determined through a calibration experiment where known weights are placed on the load cell and the voltage is recorded. The force is then divided by the cross-sectional area of the load pin, which transmits the force from the melt to the load cell. The pressure drop across the nozzle orifice is estimated by using the following relationship:

$$\varDelta \stackrel{\sim}{P}={\stackrel{\sim}{K}}_{p}Q$$

where $\varDelta \stackrel{\sim}{P}$ is the estimated pressure drop, ${\stackrel{\sim}{K}}_{p}$ is the estimated static gain value, and $Q$ is the flow rate. This relationship is then substituted into the Poiseuille equation [49] as follows:

$$\mu =\frac{{\pi D}^{4}\varDelta \stackrel{\sim}{P}}{128LQ}=\frac{{\pi D}^{4}{\stackrel{\sim}{K}}_{p}}{128L}$$

where $\mu$ is the apparent viscosity at a temperature T, $D$ is the nozzle orifice diameter, and $L$ is the length of the nozzle orifice.

The percent change in the coefficient values for the melt pressure is calculated based on their value relative to the average value of that coefficient for the factor of interest. Regression is performed on the coefficient values for the first-order model using the ‘fitlm()’ function in MATLAB. PCA of the coefficients is performed by implementing the ‘pca()’ function in MATLAB, which outputs the loadings for each coefficient, which are then plotted in a loading plot. The loading plots are designed such that the main factor is directed towards the positive side of the y-axis or ‘north’ so that the relative correlation between the coefficients and the melt pressure can be analyzed.

Figure 3 plots the observed transient shaft, infeed, and nozzle pressure as well as the transient infrared temperature for the experimental run conditions described in Table 1. It is observed that there is typically a small drop in drive force from the measured torque to the infeed load cell, then a larger pressure drop from the infeed load cell to the melt pressure in the hot end just before the nozzle. The torque and nozzle pressure exhibit a smooth response, but the infeed load varies substantially. These behaviors are described in more detail in a recent paper [1], but both the variance in the infeed pressure as well as the pressure drop between the infeed pressure and the nozzle pressure are likely due to dynamic stick-slip behavior and associated wall shear stresses of the driven filament against the inner bore of the hot end. The magnitudes of the pressure responses are dependent on the factor settings indicated in Table 1; their transient behaviors are modeled in detail for the nozzle pressure response.

Before that analysis, however, it is worth noting that the data of Fig. 3 exhibits a slow response relative to the typical process accelerations used in material extrusion additive manufacturing. For example, every start of printing will occur in 0.1 s for an acceleration of 1000 mm/s² when going from 0 to 100 mm/s. Yet, the data of Fig. 3 indicates characteristic time responses on the order of 2 s, indicating that the process may not have reached a steady state at any time during the deposition of a short path. The transient temperature behavior is also worrisome in that the melt temperature of the material being processed varies not only with the cycling of the heater but also with the run conditions (e.g. run 3 with a flow rate of 10 mm³/s). Robust process designs must account for such transient variations.

The primary result of interest is the melt pressure in the nozzle since the melt pressure is highly determinative of the deposition flow rate given compressibility and die swell effects. Figure 4 plots the melt pressure as a function of a step response in flow rate; these results correspond to runs 2, 1, and 3 in Table 1 with flow rates of 1, 5, and 10 mm³/s, respectively, while holding other factors constant. In the top plots, the flow rate and pressure are shown as a function of time. The plots show that, as the flow rate increases, the pressure increases as well. This behavior is expected since the shear stress required to extrude the material increases with shear rate, which is a function of the flow rate. The fitted first- and second-order system identification models are also shown with the standard error measured as a percentage of the mean pressure. The standard error values are on the order of 1%, which indicates that both model forms fit the transient melt pressure response. The second-order model tends to have a lower standard error since the additional model coefficients allow the fitting of more complex behaviors including under-damping and delays beyond the first-order model. As these and later results demonstrate, however, the material extrusion process is typically overdamped or near critically damped with damping ratios, $\zeta$, usually greater than one (see Table 2). For this reason, the first-order model is likely adequate and the focus of further analysis though the second-order model coefficients are reported for completeness.

The bottom left plot of Fig. 4 shows how each of the model coefficients change, relative to its mean, as a function of the flow rate. The static gains, M1.Kp and M2.Kp, are practically superimposed and so only the M2.Kp values are observable. Even though the melt pressure increases as a function of flow rate, the results show that the static gain decreases as a function of flow rate. This decrease is expected given the shear-thinning behavior of polymers and will be investigated in more detail in the discussion section. The two model coefficients that show the most change for the second-order model are the time constant, Tz, and the system damping ratio, $\zeta ,$. This result is also reflected in the loading plot from PCA in the bottom right of Fig. 4, where M2.Tz and M2.Zeta are positively correlated with the flow rate. In this case, the time constant for the process zero is more positively correlated with the flow rate than the system damping ratio. The rest of the coefficients show a lower magnitude and inverse correlation with the flow rate.

Figure 5 provides the transient melt pressure for hot end temperatures of 200, 225, and 250°C respectively while holding other factors constant, corresponding to runs 4, 1, and 5, respectively. These results show that the transient melt pressure decreases as a function of the temperature. This behavior occurs because the temperature increase decreases the material’s viscosity, which leads to a decrease in shear stress, and thus pressure. For this factor, the model standard error values are low, on the order of 1 to 3%, indicating a good fit. For the coefficient percent change as a function of temperature, the coefficient that shows the most change is the system damping ratio, M2.Zeta. This result represents the shift in the damping ratio from the most overdamped behavior with a value of 208.8 (run 4 of Table 2) to a slightly underdamped behavior with a value of 0.878 (run 5 of Table 2). The results conform to expectations given that the viscosity increases at lower melt temperatures and so greatly dampens the transient response. Even at higher melt temperatures, however, the melt does not exhibit significant under-damping or stress overshoot [51], though we have seen similar behaviors in other experiments that we believe are due to lower melt temperatures in the bore of the nozzle due to long dwell times and air cooling of the external nozzle surfaces. Such lower melt temperatures in the bore provide an increased initial flow resistance and nozzle pressure with what appears to be a pressure/stress overshoot. Subsequent melt flow and viscous dissipation then maintains a higher melt temperature and lower steady state pressure. Accordingly, we believe that the material extrusion process typically does not exhibit an underdamped transient response.

In the PCA loadings plot, the relative significance of the melt temperature is clearly exhibited by the magnitude and direction of the model coefficients relative to the size of the melt temperature loading. Specifically, the results show that the static gain, M1.Kp and M2.Kp, are inversely correlated with the melt temperature as consistent with the rheological expectation given the reduction in viscosity with increased melt temperature. The process zero exhibits a slight positive correlation by which the increasing temperature increases the zero, M2.Tz, making the system more stable. The system damping ratio has a negative correlation with the hot end temperature, likely due to the material having both a lower viscosity and more compressibility (lower bulk modulus) at higher temperatures. This reduced damping ratio makes it easier to move more mass of the material, which can lead to a very slight pressure overshoot relative to the steady value within each temperature setpoint.

Figure 6 shows the nozzle pressure behavior as a function of the nozzle diameter corresponding to runs 6, 1, and 7 in Table 1 with respective nozzle orifice diameters of 0.25, 0.4, and 0.6 mm. The steady state pressure is observed to decrease as a function of the nozzle diameter. This behavior is due to the shear stress required for the same volumetric flow rate and temperature decreasing as a function of nozzle diameter, which in turn leads to the slight decrease in the pressure drop across the length of the nozzle orifice. The experimental data for the 0.60 mm nozzle diameter (run 7) demonstrates that the transient nozzle pressure is overdamped in comparison with the smallest nozzle diameters; this behavior is consistent with the system identification values in Table 2 with $\zeta$ values of 2.689 for the 0.60 mm diameter orifice (run 7) and 1.316 for the 0.25 mm diameter orifice (run 6). The modeled static gains, Kp, also vary significantly as consistent with expectations of the Hagen–Poiseuille flow whereby the flow resistance is proportional to D^− 1/4. Specifically, there is an inverse correlation, whereby Kp increases with reduced nozzle diameter orifices. There is also an increase in the process time constant, M1.Tp and M2.Tw, corresponding to a slower system response at larger orifice diameters, which is observable by the shape of the curve in the pressure versus time plot for the 0.60 mm nozzle diameter.

Figure 7 shows the effect of acceleration, the rate at which the feed rate changes from zero to its setpoint value, according to runs 8, 1, and 9 in Table 1 with respective linear acceleration values of 100, 500, and 3000 mm²/s that correspond to volumetric accelerations of 240, 1200, and 7210 mm³/s. The main conclusion is that, for the acceleration values investigated, the transient melt pressure response is not a significant function of the acceleration. The reason is that even the lowest acceleration of 100 mm²/s is high relative to the characteristic process response times around 2 s. Accordingly, lower acceleration settings would be needed to better characterize the effect of acceleration on the transient pressure behavior.

Inspecting the transient plots and system identification model parameters, the value of the steady state pressure is not affected by the acceleration value. This result is expected, as the acceleration should mainly affect the transient response of the pressure. As with the models for the other investigated factors, both the first- and second-order models fit well though there is some cyclic process variation for the 3000 mm²/s acceleration that increases the standard error to nearly 4% of the mean pressure. The coefficient main effects plot shows that the time constant for a process zero and the time constant for the complex pole vary the most with the acceleration. PCA confirms this as well and shows that the system damping ratio does not correlate with the acceleration. The time constant for a delay is shown to be inversely related to the acceleration. The time constants for the process zero and the complex pole increasing with the acceleration makes sense, since the process will reach its steady and stable value faster when the poles and zero are further apart, which would happen when the values for their time constants both increase. The time delay decreasing is due to the material achieving its final flow rate faster as the acceleration increases, thus decreasing the delay between the stepper input and the pressure response, making physical sense.

Table 2 provides the fitted coefficients for the first-order and second-order models and their standard error values expressed as a percentage of the mean pressure. These model coefficient values were used to plot the trends in the main effects plot and the loadings plot from the principal component analysis. Generally, an inspection of the fitted coefficient values shows a relatively narrow range of values with physically meaningful values. The exceptions are the negative value for M2.Tz in run 1 as well as very low values for the characteristic time constant, M2.Tw, in runs 4 and 8. In these instances, the second-order model likely has too much model flexibility for the observed behavior such that there is confounding between the time constants M2.Tz, M2.Tw, and M2.Td. Accordingly, the first-order model is likely more robust and so is the subject for subsequent discussion.

Table 2

Summary of system identification model coefficients
Run	M1, Kp [MPa/ mm³/s]	M1, Tp [s]	M1, standard error [% mean pressure]	M2, Kp [Mpa/ mm³/s]	M2, Tw [s]	M2, Zeta [-]	M2, Td [s]	M2, Tz [s]	M2, standard error [% mean pressure]
1	0.443	1.532	2.28	0.441	0.184	3.463	0.183	-0.038	1.36
2	0.925	2.620	0.83	0.927	0.881	1.716	0.263	0.329	0.63
3	0.319	1.379	2.74	0.322	1.648	1.282	0.281	2.936	1.61
4	0.668	2.632	3.56	0.667	0.006	208.808	0.320	0.291	3.47
5	0.289	0.993	2.37	0.288	0.780	0.878	0.149	0.586	1.11
6	0.515	1.663	2.65	0.527	1.752	1.316	0.157	2.713	1.09
7	0.348	1.404	2.25	0.391	5.415	2.689	0.164	23.465	1.06
8	0.393	1.648	2.95	0.392	0.011	69.341	0.209	0.050	2.73
9	0.389	1.203	3.99	0.386	1.395	0.890	0.000	1.396	3.78

There is a significant need for real-time quality assurance whereby the process performance and product’s critical-to-quality (CTQ) attributes are constantly assessed based on robust process models. Towards this end, we seek to relate the system model identification coefficients to inherent material properties. Only the first-order model is assessed further because it provides good model fidelity and robustness with minimal complexity relative to the second-order model.

The open circle symbols in Fig. 8 represent the modeled coefficients from Table 2 as a function of the factors (flow rate, hot end temperature, and nozzle diameter). The top row of subplots in Fig. 8 provides the static gain of the system, M1.Kp, with the standard error of the fitted model predictions provided by the error bars; the standard errors are very small indicating a very high degree of fit. For all three factors (flow rate, melt temperature, and nozzle orifice diameter), the static gain decreases as a function of the factor in accordance with rheological modeling. For the flow rate and temperature, the shear stress decreases because the viscosity decreases as a function of shear rate (which is a function of flow rate) and the viscosity also decreases as a function of temperature, and shear stress is a function of the viscosity. For the nozzle orifice diameter, higher diameters reduce the flow resistance and so result in a lower static gain at the same melt flow rate and viscosity.

The bottom row of Fig. 8 shows the characteristic time constant, M1.Tp, as a function of the factors. The time constant also decreases as a function of the flow rate, temperature, and nozzle diameter. The time constant decreases with increasing flow rate because the material is transported faster at high flow rates, allowing for the pressure to stabilize sooner. This value decreases as a function of temperature because the rise in temperature leads to a decrease in the viscosity and shear stress at the wall, so there is less resistance to flow and it takes less time to reach equilibrium. The time constant as a function of nozzle diameter tends to decrease due to there being less resistance to flow in the larger nozzle orifice, i.e., lower shear stresses. The time constant is also likely related to the material compressibility, since higher flow rates lead to higher pressures and compresses the material faster than low flow rates, leading to faster pressure stabilization. The decreased resistance to flow that comes with a larger nozzle orifice leads to lower pressures, making the specific volume more stable, and the pressure stabilizes more quickly as a result.

Multiple linear regression was performed as a function of the three factors with a coefficient of determination, R², of 0.87 but all factors having a p-value above 0.05. Viewing the inverse behavior of the model coefficients in Fig. 8 as a function of the flow rate, the multiple linear regression was repeated substituting the inverse of the flow rate for the flow rate. Table 3 shows the results of this regression analysis for the first-order model coefficients. Statistics are included to indicate the significance of each factor with regard to the model parameters as well as the measures of the model significance. With a 95% confidence interval, the inverse of flow rate and the set temperature have a statistically significant effect on the static gain, Kp, and the characteristic response time, Tp. The nozzle orifice diameter has lower statistical significance but is surely physically significant given the Hagen–Poiseuille flow.

Table 3

Regression results summary table
Response & statistics	Factor	Estimate	Standard error	t-statistic	p-value
Kp R²:0.979 F-statistic: 47.2 p-value: 0.00504	Intercept	2.214	0.306	7.226	0.0055
	1/Flow rate [s/mm³]	0.6119	0.061	10.03	0.0021
	Temperature [°C]	-0.00756	0.0013	-5.757	0.0104
	Diameter [mm]	-0.478	0.186	-2.563	0.0829
Tp R²:0.957 F-statistic: 22.3 p-value: 0.0149	Intercept	9.066	1.227	7.391	0.0051
	1/Flow rate [s/mm³]	1.261	0.2443	5.163	0.0141
	Temperature [°C]	-0.0328	0.0053	-6.229	0.0083
	Diameter [mm]	-0.799	0.746	-1.071	0.363

Figure 9 shows the derived apparent viscosity as a function of the shear rate. The viscosity was derived by implementing Eq. (4). The calculated values show the expected trend of decreasing viscosity with increasing shear rate and decreasing viscosity with increasing temperature [52]. This plot demonstrates the potential for these fitted coefficients to be used for rheological property estimation and control. The viscosity predicted by the Cross-Williams Landel-Ferry (Cross-WLF) model is also plotted using equations (5) and (6) with the characterized coefficients of Table 4. The Cross-WLF model is as follows:

$\eta \left(\dot{\gamma },T\right)= {\eta }_{0}/(1+{({\eta }_{0}\dot{\gamma }/{\tau }^{*})}^{1-n})$ ()

where $\eta$ is the shear viscosity, $\dot{\gamma }$ is the shear rate, $T$ is the temperature, ${\tau }^{*}$ is the critical shear stress, $n$ is the power-law index, and ${\eta }_{0}$ is the zero-shear viscosity defined as:

$${\eta }_{0}= {\eta }_{ref}\text{e}\text{x}\text{p}(-({A}_{1}(T- {T}_{ref}))/({A}_{2}+ T -{T}_{ref}\left)\right)$$

${\eta }_{ref}$ is the reference temperature viscosity, ${T}_{ref}$ is the glass transition temperature, and ${A}_{1}$and ${A}_{2}$ are coefficients that are fit to capture the temperature dependency of the viscosity. Please see Table 4 for the model coefficient values as characterized by capillary rheology [53].

Table 4

Cross-WLF model coefficients for ABS.
Coefficient	Value
$n$	0.341
${\tau }^{*}$ [Pa]	41070
${\eta }_{ref}$ [Pa-s]	2.54 x 10¹²
${T}_{ref}$ [K]	365
${A}_{1}$	28.3
${A}_{2}$ [K]	51.6

Table 5 lists the pressures estimated with the system identification model and the resulting apparent viscosity values. The apparent pressures predicted using the fitted static gain overpredict the apparent pressures predicted by the Cross-WLF model. This may be because the static gain values capture variation due the non-isothermal condition within the nozzle, whereas the Cross-WLF model assumes an isothermal condition. Alternatively, the observed melt pressures may be significantly higher than those predicted by the Cross-WLF model characterized with the process material due to internal juncture losses that are not accounted for by Eq. (4).

Table 5

Apparent viscosity calculation values
Flow rate [mm³/s]	Temperature [°C]	Diameter [mm]	Kp [MPa/mm³/s]	Apparent pressure [Pa]	Apparent viscosity, gain [Pa-s]	Apparent viscosity, Cross-WLF [Pa-s]
5	225	0.40	0.443	2214628	696	158
1	225	0.40	0.925	924977	1453	456
10	225	0.40	0.319	3193526	502	100
5	200	0.40	0.668	3338062	1049	273
5	250	0.40	0.289	1446703	454	91
5	225	0.25	0.515	2577252	198	62
5	225	0.60	0.348	1740392	1845	352

The transient behaviors observed in material extrusion are modeled using system identification technique. First-order and second-order process models were implemented. The experiments show that the steady state pressure within the hot end varies as a function of the flow rate, hot end temperature, and the nozzle diameter in a manner consistent with processing and rheology expectations. The steady state pressure did not vary significantly with the implemented acceleration. Generally, the first-order model coefficients were very robust, resulting in model fits with standard errors less than 4% of the mean melt pressure.

The system identification model coefficients were modeled by regression as a function of the flow rate, melt temperature, and nozzle orifice diameter. The apparent viscosity as a function of shear rate is modeled using the Hagen–Poiseuille relationship and the static gain. The observed behavior followed the expected trends though there was a significant shift in the observed behavior with off-line characterization of the material viscosity. While the first-order model is simpler and likely sufficient for many control applications, the second-order model also provides insight into the underlying process physics and damping behavior. It is hypothesized that the characteristic time constant and damping ratio can be used to model the enclosed material volume in the hot end and material compressibility. Such models provide a pathway to understanding the transient behaviors observed in material extrusion and a potential technique for implementing control over these phenomena.

Acknowledgements

This work is supported by the National Science Foundation (NSF) under

GOALI Grant No. #1914651. Any opinions, findings, and conclusions or recommendations expressed in this material are those of the author(s) and do not necessarily reflect the views of the NSF. The authors also acknowledge Stratasys, Ltd. as the University of Massachusetts (UMass) Lowell team’s GOALI partner. Additionally, the authors acknowledge the National Institute of Standards of Technology for providing the nozzles used in the experiments. Lastly, the first author thanks David Rondeau, the technician at UMass Lowell, for his assistance and guidance in machining the hot end and other components

Funding

This work was supported by National Science Foundation (NSF) under GOALI Grant No. #1914651.

Competing Interests

The authors have no relevant financial or non-financial interests to disclose.

Author Contributions

All authors contributed to the study conception and design. Material preparation and data collection was performed by Austin Colon. Data analysis was performed by all authors. Software was handled primarily by David Kazmer. The first draft of the manuscript was written by Austin Colon and all authors commented on previous versions of the manuscript. All authors read and approved the final manuscript.

Colon AR, Kazmer DO, Peterson AM (2023) The dependency chain in material extrusion additive manufacturing: shaft torque, infeed load, melt pressure, and melt temperature. Additive Manuf 103780. http://doi.org/https://doi.org/10.1016/j.addma.2023.103780
Kazmer DO, Colon AR, Peterson AM, Kim SK (2021) Concurrent characterization of compressibility and viscosity in extrusion-based additive manufacturing of acrylonitrile butadiene styrene with fault diagnoses. Additive Manuf 46 Article 102106. http://doi.org/10.1016/j.addma.2021.102106
Kim SK, Kazmer DO, Colon AR, Coogan TJ, Peterson AM (2021) Non-Newtonian modeling of contact pressure in fused filament fabrication. J Rheol 65(1):27–42. http://doi.org/10.1122/8.0000052
Turner BN, Strong R, Gold SA (2014) A review of melt extrusion additive manufacturing processes: I. Process design and modeling. Rapid Prototyp J 20(3):192–204. http://doi.org/10.1108/rpj-01-2013-0012
Gilmer EL, Miller D, Chatham CA, Zawaski C, Fallon JJ, Pekkanen A, Long TE, Williams CB, Bortner MJ (2018) Model analysis of feedstock behavior in fused filament fabrication: enabling rapid materials screening. Polymer 152:51–61. http://doi.org/10.1016/j.polymer.2017.11.068
Vidakis N, Petousis M, Grammatikos S, Papadakis V, Korlos A, Mountakis N (2022) High performance polycarbonate nanocomposites mechanically boosted with titanium carbide in material extrusion additive manufacturing. Nanomaterials 12 Article 1068. http://doi.org/10.3390/nano12071068
Hong Y, Mrinal M, Phan HS, Tran VD, Liu X, Luo C (2022) In-situ observation of the extrusion processes of acrylonitrile butadiene styrene and polylactic acid for material extrusion additive manufacturing. Additive Manuf 49. Article 102507, http://doi.org/10.1016/j.addma.2021.102507
Peterson AM (2019) Review of acrylonitrile butadiene styrene in fused filament fabrication: A plastics engineering-focused perspective. Additive Manuf 27:363–371. http://doi.org/10.1016/j.addma.2019.03.030
Ovlaque P, Bayart M, Soulestin J, Trolez Y, Olivier D, Bujeau B, Charlon S (2022) On the temperature evolution and related crystallinity of polypropylene parts processed via material extrusion. Additive Manuf 58 Article 103065. http://doi.org/10.1016/j.addma.2022.103065
Hsiang Loh G, Pei E, Gonzalez-Gutierrez J, Monzón M (2020) An overview of material extrusion troubleshooting. Appl Sci 10(14) Article 4776. http://doi.org/10.3390/app10144776
Javaid M, Haleem A (2018) Additive manufacturing applications in orthopaedics: A review. J Clin Orthop Trauma 9(3):202–206. http://doi.org/10.1016/j.jcot.2018.04.008
Mohseni M, Hutmacher DW, Castro NJ (2018) Independent Evaluation of Medical-Grade Bioresorbable Filaments for Fused Deposition Modelling/Fused Filament Fabrication of Tissue Engineered Constructs. Polymers 10(1). http://doi.org/10.3390/polym10010040
Koteras R, Wieczorowski M, Znaniecki P, Swojak N (2018) The reliability of reproducing automotive components made by FDM printing in reverse engineering. Archives of Mechanical Technology and Materials 38(1):67–70. http://doi.org/10.2478/amtm-2018-0011
Steenhuis H-J, Pretorius L (2016) Consumer additive manufacturing or 3D printing adoption: an exploratory study. J Manuf Technol Manage 27(7):990–1012. http://doi.org/10.1108/jmtm-01-2016-0002
Chen D, Heyer S, Ibbotson S, Salonitis K, Steingrímsson JG, Thiede S (2015) Direct digital manufacturing: definition, evolution, and sustainability implications. J Clean Prod 107:615–625. http://doi.org/10.1016/j.jclepro.2015.05.009
Ufodike CO, Nzebuka GC (2022) Investigation of thermal evolution and fluid flow in the hot-end of a material extrusion 3D Printer using melting model. Additive Manuf 49 Article 102502. http://doi.org/10.1016/j.addma.2021.102502
Luo C, Wang X, Migler KB, Seppala JE (2020) Upper bound of feed rates in thermoplastic material extrusion additive manufacturing. Additive Manuf 32 Article 101019. http://doi.org/10.1016/j.addma.2019.101019
Go J, Schiffres SN, Stevens AG, Hart AJ (2017) Rate limits of additive manufacturing by fused filament fabrication and guidelines for high-throughput system design. Additive Manuf 16:1–11. http://doi.org/10.1016/j.addma.2017.03.007
Tamburrino F, Graziosi S, Bordegoni M (2019) The influence of slicing parameters on the multi-material adhesion mechanisms of FDM printed parts: an exploratory study. Virtual and Physical Prototyping 14(4):316–332. http://doi.org/10.1080/17452759.2019.1607758
Oleff A, Küster B, Stonis M, Overmeyer L (2021) Process monitoring for material extrusion additive manufacturing: a state-of-the-art review. Progress in Additive Manufacturing 6(4):705–730. http://doi.org/10.1007/s40964-021-00192-4
Serdeczny MP, Comminal R, Mollah MT, Pedersen DB, Spangenberg J (2020) Numerical modeling of the polymer flow through the hot-end in filament-based material extrusion additive manufacturing. Additive Manuf 36 Article 101454. http://doi.org/10.1016/j.addma.2020.101454
Seppala JE, Hoon Han S, Hillgartner KE, Davis CS, Migler KB (2017) Weld formation during material extrusion additive manufacturing. Soft Matter 13(38):6761–6769. http://doi.org/10.1039/c7sm00950j
Peterson AM, Kazmer DO (2022) Predicting mechanical properties of material extrusion additive manufacturing-fabricated structures with limited information. Sci Rep 12 Article 14736. http://doi.org/10.1038/s41598-022-19053-3
Vidakis N, David C, Petousis M, Sagris D, Mountakis N, Moutsopoulou A (2022) The effect of six key process control parameters on the surface roughness, dimensional accuracy, and porosity in material extrusion 3D printing of polylactic acid: Prediction models and optimization supported by robust design analysis. Adv Industrial Manuf Eng 5 Article 100104. http://doi.org/10.1016/j.aime.2022.100104
Serdeczny MP, Comminal R, Mollah MT, Pedersen DB, Spangenberg J (2022) Viscoelastic simulation and optimisation of the polymer flow through the hot-end during filament-based material extrusion additive manufacturing. Virtual and Physical Prototyping 17(2):205–219. http://doi.org/10.1080/17452759.2022.2028522
Osswald TA, Puentes J, Kattinger J (2018) Fused filament fabrication melting model. Additive Manuf 22:51–59. http://doi.org/10.1016/j.addma.2018.04.030
Phan DD, Horner JS, Swain ZR, Beris AN, Mackay ME (2020) Computational fluid dynamics simulation of the melting process in the fused filament fabrication additive manufacturing technique. Additive Manuf 33 Article 101161. http://doi.org/10.1016/j.addma.2020.101161
Kattinger J, Ebinger T, Kurz R, Bonten C (2022) Numerical simulation of the complex flow during material extrusion in fused filament fabrication. Additive Manuf 49 Article 102476. http://doi.org/10.1016/j.addma.2021.102476
Tlegenov Y, Wong YS, Hong GS (2017) A dynamic model for nozzle clog monitoring in fused deposition modelling. Rapid Prototyp J 23(2):391–400. http://doi.org/10.1108/rpj-04-2016-0054
Marion S, Sardo L, Joffre T, Pigeonneau F (2023) First steps of the melting of an amorphous polymer through a hot-end of a material extrusion additive manufacturing. Additive Manuf 65 Article 103435. http://doi.org/10.1016/j.addma.2023.103435
Coogan TJ, Kazmer DO (2019) In-line rheological monitoring of fused deposition modeling. J Rheol 63(1):141–155. http://doi.org/10.1122/1.5054648
Fischer J, Echsel M, Springer P, Refle O (2023) In-line measurement of extrusion force and use for nozzle comparison in filament based additive manufacturing. Progress in Additive Manufacturing. http://doi.org/10.1007/s40964-022-00385-5
Serdeczny MP, Comminal R, Pedersen DB, Spangenberg J (2020) Experimental and analytical study of the polymer melt flow through the hot-end in material extrusion additive manufacturing. Additive Manuf 32 Article 100997. http://doi.org/10.1016/j.addma.2019.100997
Anderegg DA, Bryant HA, Ruffin DC, Skrip SM, Fallon JJ, Gilmer EL, Bortner MJ (2019) In-situ monitoring of polymer flow temperature and pressure in extrusion based additive manufacturing. Additive Manuf 26:76–83. http://doi.org/10.1016/j.addma.2019.01.002
Phan DD, Swain ZR, Mackay ME (2018) Rheological and heat transfer effects in fused filament fabrication. J Rheol 62(5). http://doi.org/https://doi.org/10.1122/1.5022982
Mackay ME (2018) The importance of rheological behavior in the additive manufacturing technique material extrusion. J Rheol 62(6):1549–1561. http://doi.org/10.1122/1.5037687
Das A, Gilmer EL, Biria S, Bortner MJ (2021) Importance of polymer rheology on material extrusion additive manufacturing: correlating process physics to print properties. ACS Appl Polym Mater 3(3):1218–1249. http://doi.org/10.1021/acsapm.0c01228
Anitha M, Prusty SB (2023) System Identification and Analysis of Temperature Control System. IEEE, pp 622–626
Melsa JL (1971) System Identification (Mathematics in Science and Engineering). Academic Press, New York. (in English)
Wu P, Ramani KS, Okwudire CE (2021) Accurate linear and nonlinear model-based feedforward deposition control for material extrusion additive manufacturing. Additive Manuf 48 Article 102389. http://doi.org/10.1016/j.addma.2021.102389
Coogan TJ, Kazmer DO (2019) Modeling of interlayer contact and contact pressure during fused filament fabrication. J Rheol 63(4):655–672. http://doi.org/10.1122/1.5093033
Colon AR, Kazmer DO, Peterson AM, Seppala JE (2023) Characterization of die-swell in thermoplastic material extrusion. Additive Manuf 73 Article 103700. http://doi.org/10.1016/j.addma.2023.103700
Heller BP, Smith DE, Jack DA (2016) Effects of extrudate swell and nozzle geometry on fiber orientation in Fused Filament Fabrication nozzle flow. Additive Manuf 12:252–264. http://doi.org/10.1016/j.addma.2016.06.005
Kazmer D, Peterson AM, Masato D, Colon AR, Krantz J (2023) Strategic cost and sustainability analyses of injection molding and material extrusion additive manufacturing. Polym Eng Sci 63(3):943–958. http://doi.org/10.1002/pen.26256
Kazmer DO, Colon A (2020) Injection printing: additive molding via shell material extrusion and filling. Additive Manuf 36 Article 101469. http://doi.org/10.1016/j.addma.2020.101469
Júnior GA, dos Santos JBM, Barros PR (2009) On Simple Identification Techniques for First-Order plus Time-Delay Systems, IFAC Proceedings Volumes 42 (10) 605–610, http://doi.org/10.3182/20090706-3-FR-2004.00100
Niemann H, Miklos R (2014) A Simple Method for Estimation of Parameters in First order Systems, Journal of Physics: Conference Series 570 (1) 1–1, http://doi.org/10.1088/1742-6596/570/1/012001
Chen L, Li J, Ding R (2011) Identification for the second-order systems based on the step response. Math Comput Model 53(5–6):1074–1083. http://doi.org/10.1016/j.mcm.2010.11.070
O'Carroll C, Sorbie KS (1993) Generalization of the Poiseuille law for one- and two-phase flow in a random capillary network. Phys Rev E Stat Phys Plasmas Fluids Relat Interdiscip Topics 47(5):3467–3476. http://doi.org/10.1103/physreve.47.3467
Seppala JE, Migler KD (2016) Infrared thermography of welding zones produced by polymer extrusion additive manufacturing. Additive Manuf 12:71–76. http://doi.org/10.1016/j.addma.2016.06.007
Divoux T, Barentin C, Manneville S (2011) Stress overshoot in a simple yield stress fluid: An extensive study combining rheology and velocimetry. Soft Matter 7(19). http://doi.org/10.1039/c1sm05740e
Osswald T, Rudolph N (2014) Generalized Newtonian Fluid Models. Polymer Rheology - Fundamentals and Applications. Hanser Publications, Cincinnati, Ohio, pp 59–97
Osswald T, Rudolph N, Rheometry (2014) Polymer Rheology - Fundamentals and Applications. Hanser, pp 187–220

Download PDF

Version 1

posted

You are reading this latest preprint version

Transient modeling of material extrusion by system identification

Status:

Version 1

Abstract

Figures

Highlights

1. Introduction

2. Methodology

2.1 Material

2.2 Experimental hardware

2.3 Experimental factors

2.4 System identification

3. Analysis

4. Results

5. Discussion

6. Conclusions

Declarations

References

Status:

Version 1

Coefficient	Value
\(n\)	0.341
\({\tau }^{*}\) [Pa]	41070
\({\eta }_{ref}\) [Pa-s]	2.54 x 10¹²
\({T}_{ref}\) [K]	365
\({A}_{1}\)	28.3
\({A}_{2}\) [K]	51.6