Exploiting nozzle geometry to predict resolution in extrusion-based bioprinting: mathematical modelling of a power-law fluid

TL;DR

This paper presents a mathematical model to predict the resolution of extrusion-based bioprinting by considering nozzle geometry and material properties.

Contribution

The novelty lies in the development of a user-friendly model that predicts printed filament diameter using nozzle geometry and fluid properties.

Findings

The model accurately predicts filament diameter with an R2 value greater than 0.97 when shear-thinning is considered.

The model was validated using two materials, three nozzle sizes, and two temperature conditions.

The model helps reduce reliance on trial-and-error methods by narrowing down optimal process parameters.

Abstract

Extrusion-based additive manufacturing (AM) is a popular technique used in the fabrication of three-dimensional constructs. Owing to the nonlinear manner in which process parameters affect resolution and printability, the optimal combination remains platform and material specific. This study proposes a user-friendly, adaptable model to predict the diameter of a printed line of material through extrusion-based bioprinting. Exploiting the geometry of an arbitrary, axisymmetric nozzle and assuming a power-law fluid, the model generated determines a relationship between the printed filament diameter and the pressure drop, nozzle travel speed, nozzle geometry and material flow properties. Employing the model prior to printing enables engineers to restrict process parameter space and minimize the dependence on the current print-and-test methodology before an optimal combination of process parameters is determined. Two materials (a poly(vinyl alcohol)-based (PVA-based) hydrogel and Nivea Crème), two temperature conditions and three nozzle sizes were used for model validation, presenting good agreement with model predictions. When the shear-thinning property is included, the coefficient of determination, R2, is greater than 0.97. This model provides context and direction for future optimization-driven design research for this advancing manufacturing technology.

Figures

The extrusion-based bioprinting process, with P the pressure applied to the hydrogel, d the layer thickness and DN, LN, UN and H—the nozzle diameter, length, travel speed and offset, respectively.

A flowchart presenting the current design for the AM workflow. Arrows (red) returning to earlier steps illustrate their current iterative nature.

A schematic of (a) a conical nozzle with an inner diameter Di (inner radius ri), outer diameter Do (outer radius ro), length Lc and an angle θ between the nozzle wall and the upward vertical, and (b) a straight, cylindrical nozzle with a diameter DN (radius RN) and length LN.r denotes the radius of an arbitrary small segment of the conical nozzle of height δz with a pressure drop δP; ΔP denotes the pressure drop between the pressure applied at the top and the pressure at the nozzle tip, and τc and τs the shear stress along the conical and straight nozzle walls, respectively.

A print of the nine lines using Nivea Crème with a 22G conical nozzle, a 5 mm s−1 nozzle travel speed, 0.3 mm nozzle offset and an air pressure of 62×103 Pa. The red lines marked on the first line indicate the three measurements taken before an average for the filament width was taken for each line.

A regression analysis to determine the relationship between filament width (y-axis) and pressure drop (x-axis), with (a) a poly(vinyl alcohol)-based hydrogel and (b)–(d) Nivea Crème, with a 25G, 22G and 20G conical nozzle, respectively. In each case, the top graph exhibits a curve fitting of the form d=ax1/2n, with a a constant coefficient and n the flow behaviour index, listed in Table 1. The residual plots (bottom) present the deviation between the data and fitted curve, with error bars representing 95% (uncorrected) confidence intervals and significance based on FDR-corrected p-values.

Three prints using Nivea Crème with a 22G conical nozzle, a 5 mm s−1 nozzle travel speed, and 0.3 mm nozzle offset, ranging from under- (left) to over-extrusion (right). An air pressure of 58×103 Pa was used on the left, 62×103 Pa in the middle and 80×103 Pa on the right.

A regression analysis to determine the relationship between filament width (y-axis) and pressure drop (x-axis), with (a) a poly(vinyl alcohol)-based hydrogel and (b)-(d) Nivea Crème, with a 25G, 22G and 20G conical nozzle, respectively. In each case, the top graph exhibits a curve fitting of the form d=ax1/2n, with n=1 and a a constant coefficient listed in Table 2. The residual plots (bottom) present the deviation between the data and the fitted curve, with error bars representing 95% (uncorrected) confidence intervals and significance based on FDR-corrected p-values.

A flowchart presenting how the user may implement the model to establish a window of printability and reduce dependence on the current print-and-test methodology.

Tables

Table 1.. Values of the constant coefficient, a, for the flow behaviour and consistency indices, n and K (Pa sn), respectively, corresponding to the optimal non-Newtonian curve fitting of the form d = a(ΔP)1/2n, for a poly(vinyl alcohol)-based hydrogel and Nivea Crème (using a variety of conical nozzles) obtained via regression analysis conducted in MATLAB® and captured in Figures 5a–d, with ΔP the pressure drop and d the printed filament width. The goodness of fit of each fitting is captured by the coefficient of determination, R2.

material and nozzle used	a	n	K (Pa sn)	R 2
poly(vinyl alcohol)-based hydrogel, 25G	8.930 × 10–9	0.4570	2.9596 × 103	0.9940
Nivea Crème, 25G	1.138 × 10–18	0.1694	11.09052	0.9913
Nivea Crème, 22G	1.220 × 10–18	0.1654	10.4465	0.9725
Nivea Crème, 20G	1.203 × 10–18	0.1635	9.2926	0.9932

Table 2.. Values of the constant coefficient, a, for the consistency index, K (Pa sn), corresponding to the optimal Newtonian curve fitting of the form d = a(ΔP)1/2n, where n = 1, for a poly(vinyl alcohol)-based hydrogel and Nivea Crème (using a variety of conical nozzles) obtained via regression analysis conducted in MATLAB® and captured in Figures 5a–d, with ΔP the pressure drop and d the printed filament width. The goodness of fit of each fitting is captured by the coefficient of determination, R2.

material and nozzle used	a	K (Pa sn)	R 2
poly(vinyl alcohol)-based hydrogel, 25G	4.448 × 10–6	2.2388 × 107	0.7532
Nivea Crème, 25G	1.296 × 10–6	1.9014 × 106	0.3146
Nivea Crème, 22G	2.020 × 10–6	1.0917 × 106	0.2920
Nivea Crème, 20G	2.806 × 10–6	7.7960 × 105	0.2883

Table 3.. Values of the flow behaviour and consistency indices, n and K (Pa sn), respectively, for poly(vinyl alcohol)-based hydrogels and Nivea Crème obtained via standard rheological tests.

study	materials	details	n	K (Pa sn)
Paxton et al. [12]	Nivea Crème	temperature: 24°C	0.552	26.1
Tu et al. [40]	Nivea Crème	temperature: unspecified	0.045	867
Jawad et al. [41]	PVA	5% w/w solution, with molecular weight, hydrolysis, and temperature unspecified	0.18	—
Lewandowska et al. [42]	PVA	10% w/w concentration, molecular weight 30–70 kDa, hydrolysis 87–90%, at a temperature of 30°C	0.98	0.037

Equations

Full text

Introduction

Additive manufacturing (AM) is a rapidly advancing technology that presents far greater design freedom compared to traditional manufacturing methods, allowing optimization to drive design and generate more efficient parts [1,2]. Its technologies have been widely employed in tissue engineering and have adopted the name ‘bioprinting’ [3]. Bioprinting involves the assembly of three-dimensional, often patient-specific, tissue-like structures via accurate deposition of a bioink (a formulation compatible with cell culture, sometimes containing encapsulated cells, and suitable for biofabrication) in a layer-by-layer fashion. Bioprinting provides a promising trajectory for tissue replacement and may in time satisfy organ shortages, emphasizing the importance of research into the optimization of its design strategy [4].

Inspired by the established extrusion technique used in fused deposition modelling (a method of polymer AM), extrusion-based bioprinting involves continuous filaments of bioink that are pneumatically or mechanically extruded from a cartridge on to the build platform (as seen in Figure 1). Owing to its compatibility with a wide variety of viscosities and the high cell density that may be achieved within the three-dimensional structures produced, extrusion-based bioprinting is the most popular choice of bioprinting method described within the literature [5].

The extrusion-based bioprinting process, with P the pressure applied to the hydrogel, d the layer thickness and DN, LN, UN and H—the nozzle diameter, length, travel speed and offset, respectively.

Design for AM comprises the entire process chain for fabrication, from assessing the appropriate manufacturing method to preparing the geometric data and validating for a specific AM technology. To date, the design workflow remains fairly fragmented and time-consuming, presenting a sequence of manual steps repeated several times until the requirements for a particular component are met. Figure 2 outlines the part and process-specific design for AM workflow currently undertaken by design engineers (according to Wiberg et al. [6]), where the interdependent nature of part and process design is identified. A more comprehensive framework to guide fabrication can be found in Pradel et al. [7]. Automation of the workflow, with adaptable parametric computer-aided design models that synchronize with feedback loops, is necessary to enable designers to control the geometry and machine settings simultaneously and optimize the system as a whole.

A flowchart presenting the current design for the AM workflow. Arrows (red) returning to earlier steps illustrate their current iterative nature.

Sustainable design is of utmost importance in the manufacturing industry. It is desirable to determine a methodology that generates minimal waste. Simulating the current iterative print-and-test methodology (seen in Figure 2) is a logical step towards achieving this aim. In this study, we focus on the initial preparation of the bioprinter, aiming to restrict process parameter space through prediction of the printed thread diameter and, thus, minimizing the dependence on the current print-and-test methodology to determine the optimal combination.

Although a consistent definition is yet to be established, we shall define the printability of a bioink to be the deviation of a printed construct from the initial three-dimensional model data for a particular application, characterized in terms of extrudability, shape fidelity and structural integrity [8]. Here, extrudability refers to the capability of a bioink to be continuously extruded through a nozzle in a controlled manner; shape fidelity to the difference in cross-sectional area of a filament from that set in the three-dimensional model data; and structural integrity to the capability of a construct to hold its intended structure post-print. Furthermore, when considering the viability of cells encapsulated in a bioink, the success of a construct is no longer just dependent on the printability of the material, but on the proportion of healthy cells remaining post-print [9]. In fact, owing to the complex rheological behaviour observed at the microscale, obtaining consistent results between applications has proven extremely challenging [10]. While the ultimate goal remains to optimize the printability of cell-laden materials, it is necessary to first understand fluids with similar rheological properties. As such, this paper primarily focuses on capturing fluid flow properties independent of cell concentration.

Future developments in the bioprinting field involve the generation of quantitative tools and methods to objectively analyse the performance of different bioinks during the extrusion-based bioprinting process. By simulating the process, any experimental analysis and between-laboratory biases can be greatly reduced, aiding the reproducibility of results across the literature, and increasing time and cost efficiency [11]. Approaches include mathematical models that use data from the rheological properties of a particular bioink as an input to predict the optimal printing environment for desired printability. These in silico models include power-law [12] (or Ostwald de Waele model [13]) and Herschel–Bulkley models [14]. The power-law model takes the form

$$[eqn]$$

where $[eqn]$ denotes shear stress, $[eqn]$ denotes the consistency index (dependent on temperature, air pressure, shear rate and the flow behaviour index [15,16]), $[eqn]$ denotes shear rate and $[eqn]$ denotes the flow behaviour index (with shear-thinning materials identified when $[eqn]$ and Newtonian materials when $[eqn]$ ). We note that a Newtonian fluid presents a linear relationship between shear stress and shear rate, while fluids presenting a deviation from this linearity are termed non-Newtonian. An ideal bioink should exhibit a non-Newtonian behaviour: a shear-thinning property upon extrusion and rapid gelation upon deposition, to aid shape retention [17].

He et al. [3] identified the key process parameters necessary to gain precise control of the printability of a gelatin-alginate composite hydrogel via a series of experiments. They determined that the resolution of a printed line of hydrogel is affected by the nozzle travel speed, $[eqn]$ , nozzle offset, $[eqn]$ and air pressure, $[eqn]$ . Schwab et al. [14] extended this shortlist by identifying a further dependence on the geometry of the nozzle. Their effect on the resolution of a print was examined theoretically (and validated experimentally using a Pluronic F-127 hydrogel) by Suntornnond et al. [18] via a mathematical model derived from the definition of dynamic viscosity, $[eqn]$ . Although their work was based in a platform- and material-specific setting, the agreement between the theoretical and the experimental results emphasizes that mathematical modelling techniques pose a promising alternative to the existing methodology. The model presented in this study provides a generalization of that presented by Suntornnond et al. [18], where a straight nozzle is assumed.

Exploiting the geometry of an arbitrary, axisymmetric nozzle and assuming a power-law fluid, the aim of this study is to present a user-friendly, adaptable mathematical model to determine the printed filament diameter. In particular, a relationship between the printed filament diameter and the pressure drop, nozzle travel speed, nozzle geometry and material flow properties is determined. Examples of a conical and a straight nozzle are presented.

The model results are validated by using two distinctly different materials, a poly(vinyl alcohol)-based (PVA-based) hydrogel and the highly printable Nivea Crème, each using two different temperature conditions and three different conical nozzle sizes.

Employing the model prior to printing enables a ‘window of printability’ to be established, reducing the number of iterations of the current print-and-test methodology before an optimal combination of process parameters is obtained. The implications in the context of cell-laden materials are discussed.

Material and methods

Mathematical model

2.1.

Conservation of mass is an important physical mechanism driving fluid flow; what leaves the nozzle must reach the build platform. As such, the foundation of the model is the concept of flux. We consider a viscous thread of material continuously extruded from a nozzle, with no-slip at the boundary between the fluid and the nozzle wall. We assume the thread maintains a circular cross section from the point of extrusion to the point it reaches the linearly moving build platform and assume the fluid flow to be steady and axisymmetric.

Flux leaving the nozzle

2.1.1.

Considering a small segment of an axisymmetric nozzle with an arbitrary cross-sectional area of height $[eqn]$ and radius $[eqn]$ , the forces parallel to its centreline may be resolved, i.e.

$$[eqn]$$

where $[eqn]$ is the pressure applied perpendicular to the cross section, $[eqn]$ the frictional shear stresses along the nozzle wall and $[eqn]$ the downward pointing unit vector. As such, we generate

$$[eqn]$$

with $[eqn]$ the pressure drop (the difference between the pressure applied to the ink at the top and the pressure at the bottom of the segment), $[eqn]$ the shear stress (along the nozzle wall), $[eqn]$ the angle between the nozzle wall and the upward vertical and $[eqn]$ the radius of the segment [19]. Here, the first term corresponds to the pressure drop across the cross-sectional area of the segment ( $[eqn]$ ), and the second, the component of shear stress parallel to the centreline ( $[eqn]$ ) across its surface area ( $[eqn]$ ). Assuming a power-law fluid in a tube, where the shear rate is given by

$$[eqn]$$

with $[eqn]$ the volumetric flow rate in the nozzle and $[eqn]$ the flow behaviour index [20], equations (1.1) and (2.3) can be substituted into equation (2.2) for $[eqn]$ to obtain

$$[eqn]$$

Integrating over the length of the nozzle, $[eqn]$ , from $[eqn]$ (at the top of the nozzle) to $[eqn]$ (at the nozzle tip), generates

$$[eqn]$$

where $[eqn]$ , with $[eqn]$ the pressure applied to the ink at the top of the nozzle and $[eqn]$ at the nozzle tip (assumed to be atmospheric pressure, with 1 atm = $[eqn]$  Pa). For a particular nozzle geometry, equation (2.5) can be employed to obtain the volumetric flow rate, $[eqn]$ .

Assuming a conical nozzle (seen in Figure 3a), equation (2.5) becomes

A schematic of (a) a conical nozzle with an inner diameter Di (inner radius ri), outer diameter Do (outer radius ro), length Lc and an angle θ between the nozzle wall and the upward vertical, and (b) a straight, cylindrical nozzle with a diameter DN (radius RN) and length LN.r denotes the radius of an arbitrary small segment of the conical nozzle of height δz with a pressure drop δP; ΔP denotes the pressure drop between the pressure applied at the top and the pressure at the nozzle tip, and τc and τs the shear stress along the conical and straight nozzle walls, respectively.

$$[eqn]$$

where $[eqn]$ is the inner radius, $[eqn]$ the outer radius and $[eqn]$ the length of the conical nozzle, in the limits of integration. Solving and rearranging for $[eqn]$ , generates

$$[eqn]$$

using $[eqn]$ . Exploiting the definition of the volumetric flow rate, we obtain an average extrusion speed from the conical nozzle, $[eqn]$ , by dividing through by the cross-sectional area of the nozzle tip

$$[eqn]$$

Following a similar methodology for a straight, cylindrical nozzle (seen in Figure 3b), the volumetric flow rate is given by

$$[eqn]$$

with $[eqn]$ and $[eqn]$ the radius and length of the straight nozzle, respectively, and an extrusion speed, $[eqn]$ , given by

$$[eqn]$$

Flux on the build platform

2.1.2.

The volumetric flow rate on the build platform takes the form

$$[eqn]$$

with $[eqn]$ the printed filament diameter and $[eqn]$ the build platform travel speed. Since flux at the nozzle tip and build platform must balance, equation (2.11) must equal equations (2.7) and (2.9). Equating equations (2.7) and (2.9) to equation (2.11) and rearranging for $[eqn]$ , expressions for the printed filament diameters for conical and straight nozzles, $[eqn]$ and $[eqn]$ , respectively, are generated:

$$[eqn]$$

with $[eqn]$ and $[eqn]$ , the inner diameter, $[eqn]$ , and outer diameter, $[eqn]$ , of the conical nozzle, respectively, and $[eqn]$ the diameter of the straight nozzle. Since an ideal bioink exhibits a shear-thinning property, we expect $[eqn]$ . In the Newtonian limit, where $[eqn]$ and $[eqn]$ , equation (2.12) reduces to

$$[eqn]$$

Note that equation (2.13a) becomes equation (2.13b) when $[eqn]$ . In fact, equation (2.13b) is an expression exploiting Poiseuille’s law [21].

Assuming the nozzle geometry, build platform travel speed, $[eqn]$ , and material properties, $[eqn]$ and $[eqn]$ , are known, equation (2.12) reduces to a proportionality, $[eqn]$ . With the syringe diameter being much greater than the nozzle tip diameter, we assume the pressure drop in the syringe to be negligible compared to that across the nozzle and, therefore, $[eqn]$ becomes, $[eqn]$ , with $[eqn]$ the air pressure applied to the hydrogel at the top of the syringe and $[eqn]$ at the nozzle tip. That is, $[eqn]$ can be assumed as $[eqn]$ in the earlier definition obtained in equation (2.7). We may also replace $[eqn]$ by the nozzle travel speed, $[eqn]$ , without altering its meaning, enabling the use of the model across three-dimensional bioprinters. Whether the build platform travel speed guides the dragging of the filament along the build platform or the nozzle travel speed, the definition of the rate of fluid flow through an arbitrary cross section of the thread on the build platform remains the same.

Experimental methods

2.2.

Two inks were used. The first, a PVA with molecular weight 146–186 kDa, hydrolysis of 99+%, and a $[eqn]$ w/w concentration (Sigma Aldrich, St Louis, MO, USA) and the second, Nivea Crème (Beiersdorf Global AG, Germany), used as purchased. The PVA was dissolved in deionized water via mechanical stirring for 1 hour at $[eqn]$ C. The hotplate was then removed and stirring continued for 1 hour, until the solution reached room temperature [22].

Nivea Crème is an example of a soft colloidal ink, maintaining a constant quality and composition [12]. Owing to its formula, it is temperature stable and able to withstand low and high temperatures without significant changes in its consistency [23]. As such, Nivea Crème is considered a highly printable material for use in bioprinting, presenting an ideal extrudability, shape fidelity and structural integrity. It has even been established as a demonstration ink by bioprinting pioneer REGENHU (Switzerland) [12]. On the other hand, the properties of PVA depend on its molecular weight and concentration as well as on temperature. PVA is a thermoresponsive material, behaving as a Newtonian fluid at room temperature; sub-zero temperatures initiate the non-Newtonian flow necessary for it to hold its shape. Choosing these two materials, a wide range of air pressures and two different build platform temperatures could be considered. This enabled the printing conditions required to achieve printability and emphasise the nonlinear manner to which process parameters affect it.

Four experiments were undertaken using an INKREDIBLE+3D bioprinter (CellInk, Gothenburg, Sweden) with a 5 mm s^−1^ nozzle travel speed. As noted in §1, the experiments focused on process parameters, nozzle travel speed, nozzle offset, air pressure and nozzle geometry, with the aim of determining their effect on print resolution. The temperature of the build platform was also controlled in the PVA-based hydrogel case to enable filament shape retention (not being necessary for the Nivea Crème, the build platform varied with ambient temperature from approximately 21–28°C). For each set-up and material considered, nine lines of hydrogel were continuously printed on to the build platform, with each line measured at three points (start, middle, end) and an average taken, as seen in Figure 4. Note that with the accuracy displayed on the pressure gauge being to the nearest integer, close attention was paid to set it as close to the intended pressure as possible, that is, aiming for as close to $[eqn]$ Pa as opposed to $[eqn]$ Pa.

A print of the nine lines using Nivea Crème with a 22G conical nozzle, a 5 mm s−1 nozzle travel speed, 0.3 mm nozzle offset and an air pressure of 62×103 Pa. The red lines marked on the first line indicate the three measurements taken before an average for the filament width was taken for each line.

In experiment 1, the PVA-based hydrogel was printed using a frozen build platform ( $[eqn]$ C), a 25G conical nozzle (outer diameter: 0.25 mm) and a 0.3 mm nozzle offset, with nine different air pressures considered (10, 15, 20, 25, 30, 35, 40, 45, 50 × 10^3^ Pa). This sub-zero three-dimensional printing methodology is described further in Crolla et al. [22].

The next three experiments, 2(a), 2(b) and 2(c), were conducted with the Nivea Crème, exploiting its highly printable nature to determine the effect of nozzle geometry on shape fidelity, using a build platform at room temperature [12]. For 2(a), a 25G conical nozzle and a 0.1 mm nozzle offset were used, with 30 different air pressures considered ( $[eqn]$ Pa); for 2(b), a 22G conical nozzle (outer diameter: 0.41 mm) and a 0.3 mm nozzle offset were used, with 18 different air pressures considered ( $[eqn]$ Pa); and for 2(c), a 20G conical nozzle (outer diameter: 0.58 mm) and a 0.4 mm nozzle offset were used, with 17 different air pressures considered ( $[eqn]$ Pa). For each air pressure, the print was repeated three times and an average line width was taken. With the nozzle travel speed fixed at 5 mm s^−1^, the nozzle offset was varied to enable adhesion at the range of nozzle tip diameters considered. Therefore, to limit the number of parameters changed across the Nivea Crème experiments, the ratio between nozzle tip radius and nozzle offset (the aspect ratio) was kept approximately the same throughout. Based on measurements taken by Gómez-Blanco et al. [24], the inner diameter remains the same across conical nozzle sizes and is taken to be 4.02 mm.

Each experiment was captured using a Dino-Lite Digital Microscope (AnMo Electronics Corporation, Taiwan) and analysed in ImageJ (NIH, MD, USA). All subsequent statistical analyses were performed in MATLAB $[eqn]$ (MathWorks, MA, USA) v. 23.2.0.2380103 (R2023b). Regression analysis was conducted using the Optimization Toolbox^TM^, statistical hypothesis testing using Student’s t‐test in the Statistics and Machine Learning Toolbox^TM^, and a false discovery rate (FDR) correction made using the Bioinformatics Toolbox^TM^.

Results and discussion

Fundamentally, the air pressure set on a bioprinter is the additional pressure necessary to overcome atmospheric pressure at the nozzle tip and initiate flow (relative to the viscosity of the medium). As such, air pressure and pressure drop are used interchangeably throughout this section. Figures 5a–d presents measurements for the average filament width, $[eqn]$ , for a PVA-based hydrogel and Nivea Crème, respectively, using a variation of conical nozzle sizes and air pressures (also found in electronic supplementary material, Tables S1 and S2). Results from Student’s t‐test (with FDR correction) conducted on the residual between the experimental data and the model prediction are also presented (also found in electronic supplementary material, Tables S3–S6).

A regression analysis to determine the relationship between filament width (y-axis) and pressure drop (x-axis), with (a) a poly(vinyl alcohol)-based hydrogel and (b)–(d) Nivea Crème, with a 25G, 22G and 20G conical nozzle, respectively. In each case, the top graph exhibits a curve fitting of the form d=ax1/2n, with a a constant coefficient and n the flow behaviour index, listed in Table 1. The residual plots (bottom) present the deviation between the data and fitted curve, with error bars representing 95% (uncorrected) confidence intervals and significance based on FDR-corrected p-values.

Parameter fitting

3.1.

An empirical model (equation (2.12a)) was derived, relating the printed filament diameter, $[eqn]$ , to the pressure drop, $[eqn]$ , via a constant of proportionality, $[eqn]$ , incorporating nozzle geometry (conical nozzle, $[eqn]$ , $[eqn]$ and $[eqn]$ , for experimental validation purposes), the hydrogel consistency index, $[eqn]$ , the build platform travel speed, $[eqn]$ (or nozzle travel speed, $[eqn]$ ) and the flow behaviour index, $[eqn]$ . Figures 5a–d captures fittings to the experimental data for the PVA-based hydrogel and Nivea Crème, for the 25G, 22G and 20G conical nozzles, respectively, with functional form

$$[eqn]$$

where

$$[eqn]$$

as presented in equation (2.12a), with $[eqn]$ replaced by $[eqn]$ , since the INKREDIBLE+3D bioprinter relies on a moving nozzle. Table 1 lists the values of the constant coefficient, $[eqn]$ , for the flow behaviour and consistency indices, $[eqn]$ and $[eqn]$ (Pa s $[eqn]$ ), respectively, for each fitting.

Exploiting the flow behaviour index

3.1.1.

Across experiments, the data suggest that an increase in air pressure results in an increase in the average printed filament diameter, moving from an under-extruded state to an over-extruded state via the optimal air pressure, where the filament width, $[eqn]$ , matches the nozzle tip diameter (as seen in Figure 6). In each case, the flow behaviour and consistency indices, $[eqn]$ and $[eqn]$ , respectively, were chosen to optimize the fitting to the data (the $[eqn]$ and $[eqn]$ corresponding to the minimum sum of squared residuals) and, from the residual plots, the data points are seen to lie randomly around the zero line (the fitted line), with an even distribution above and below, indicating that the model captures the data well (with goodness of fit, $[eqn]$ , $[eqn]$ , $[eqn]$ and $[eqn]$ , respectively, as seen in Figure 5). As seen in Table 1, the optimal $[eqn]$ located in each case lies in the range $[eqn]$ , capturing the non-Newtonian, shear-thinning behaviour typical of materials used for extrusion-based bioprinting, with the Nivea Crème presenting a greater shear-thinning property compared to the PVA-based hydrogel, highlighting its highly printable nature [12]. Furthermore, $[eqn]$ and $[eqn]$ remain relatively similar (at an order of magnitude approx. $[eqn]$ ) across the three Nivea Crème experiments, as expected, with differences potentially a result of experimental variation. In fact, since $[eqn]$ , for the average $[eqn]$ across these experiments, taken to be 0.1661, the value of $[eqn]$ corresponding to the best fit in each case remains relatively similar (at an order of magnitude approx. $[eqn]$ ), with $[eqn]$ and $[eqn]$ Pa s $[eqn]$ for the 25G, 22G and 20G nozzles, respectively [16].

Three prints using Nivea Crème with a 22G conical nozzle, a 5 mm s−1 nozzle travel speed, and 0.3 mm nozzle offset, ranging from under- (left) to over-extrusion (right). An air pressure of 58×103 Pa was used on the left, 62×103 Pa in the middle and 80×103 Pa on the right.

Figures 7a–d presents the proportionality fitting for the Newtonian limit when $[eqn]$ (with optimal curve fittings listed in Table 2). Although it captures a similar trend to the data in the case of the PVA-based hydrogel ( $[eqn]$ ), it does not fit the shape closely, lying outside of the data towards the two end points (as seen in Figure 7a). In fact, across all cases, the residual plots reveal the model to match closely in the middle and then deviate from the data towards the two end points, presenting a statistically significant difference between the model prediction and experimental average across all points and indicating a poor fit overall. This reinforces the significance of $[eqn]$ in the proportionality, capturing the non-Newtonian fluid flow behaviour revealed in the experimental data and, as such, obtaining a closer fit. This is reflected in Figure 5, where any statistically significant difference between the model and the experimental data is reduced as a result of implementing the non-Newtonian model.

A regression analysis to determine the relationship between filament width (y-axis) and pressure drop (x-axis), with (a) a poly(vinyl alcohol)-based hydrogel and (b)-(d) Nivea Crème, with a 25G, 22G and 20G conical nozzle, respectively. In each case, the top graph exhibits a curve fitting of the form d=ax1/2n, with n=1 and a a constant coefficient listed in Table 2. The residual plots (bottom) present the deviation between the data and the fitted curve, with error bars representing 95% (uncorrected) confidence intervals and significance based on FDR-corrected p-values.

In the PVA-based hydrogel and Nivea Crème experiments, a total of 27 and 81 filament width measurements were recorded per air pressure, respectively. With fewer measurements in the PVA-based hydrogel case, average values may be skewed by measurement error resulting in wider 95% confidence intervals compared to the Nivea Crème experiments, as seen in Figures 5a and 7a. As such, comparison of experimental averages with model prediction may be unreliable, with any difference lacking statistical significance. Therefore, repeating the PVA-based hydrogel experiments three times, as conducted with the Nivea Crème, may greatly improve its reliability and enable a greater statistical power.

Practical use

3.1.2.

Presenting a good agreement with experiment, equations (2.12a) and (2.12b) can be implemented by the user to establish a window of printability and reduce the number of iterations of the current print-and-test methodology before an optimal combination of process parameters is obtained. Figure 8 presents a flowchart of how the user may apply the model in this way.

A flowchart presenting how the user may implement the model to establish a window of printability and reduce dependence on the current print-and-test methodology.

Since we require the printed diameter, $[eqn]$ , to be as close to the nozzle tip diameter as possible, we may set $[eqn]$ in equation (2.12a) and $[eqn]$ in equation (2.12b), to some small tolerance. Then, fixing the material, $[eqn]$ and $[eqn]$ are specified (provided the temperature and solute concentration are controlled) while selection of a nozzle geometry provides $[eqn]$ , $[eqn]$ and $[eqn]$ , for conical, or $[eqn]$ and $[eqn]$ , for straight nozzles. As such, equations (2.12a) and (2.12b) reduce to

$$[eqn]$$

where $[eqn]$ and

$$[eqn]$$

for conical and straight nozzles, respectively. Therefore, the pressure drop, $[eqn]$ , and build platform travel speed, $[eqn]$ (or nozzle travel speed, $[eqn]$ ), may be varied such that equation (3.3) is satisfied to some small tolerance. In doing so, a window of printability is established, and the current print-and-test methodology is conducted on a reduced parameter space. By generalizing to an arbitrary, axisymmetric nozzle geometry, this methodology generalizes that implemented by Paxton et al. [12], where a straight nozzle is assumed, and enhances the applicability of the tool in the broader context of printing with cell-laden materials.

Implications for cell-laden bioinks

3.1.3.

At present, the principal objective for the model is to optimize printability, specifically, shape fidelity. However, as revealed in §1, in the case of cell-laden bioinks, the complexity of the optimization problem increases twofold. When cell-laden, a suitable trade-off between printability and cell viability has to be determined to achieve optimal mechanical and biological performance of the intended construct [25–28].

The residence time (how long a material particle is exposed to certain shear conditions [12]) and shear forces imposed on the cells need to be controlled. If they are too high, they can cause a loss in immediate and long-term cell viability [29], particularly for those cells travelling along the nozzle wall [12,30]. Temperature control can also be crucial to cell viability. Ouyang et al. [31] explored the biofabrication of embryonic stem cells into cell-laden constructs and discovered that control of the nozzle and chamber temperatures could lead to an increase in cell viability from 55.5 to 90%. In fact, the time period required for assembly and maturation of the cell-laden construct could even exceed the cell survival time [32]. On the other hand, there is evidence to suggest that the encapsulated cells can affect the rheological properties of the bioink and the chemical processes involved and, thus, correspond to a loss in printability, but the extent of these effects still remains elusive [8,14,33,34].

The model is agnostic about whether the material is cell-laden or not. Nevertheless, it provides the engineer with a shortlist of parameters to experiment with, from which they can choose the combination corresponding to the optimal trade-off between printability and cell viability. The model therefore reduces dependency on print-and-test methodologies. This reduction is particularly valuable when printing with cell-laden materials, owing to the costliness of cells [11]. Furthermore, by considering an arbitrary, axisymmetric nozzle, a nozzle geometry more desirable for printing with cell-laden materials may be specified.

Exploring nozzle geometry

3.1.4.

The importance of nozzle selection is emphasized throughout the literature [35]. Mohammed et al. [36] optimized the material and process parameters for a semisolid extrusion-based three-dimensional printing technique involved in the three-dimensional printing of drug products. They concluded that a straight nozzle requires a higher pressure with a greater chance of clogging, while conical nozzles require a much lower pressure, with chances of clogging greatly reduced, providing a much smoother, continuous flow from the nozzle tip. However, throughout the literature, a straight, cylindrical nozzle is typically assumed where empirical modelling techniques are employed [18,37–39]. For example, Suntornnond et al. [18] present a non-Newtonian extension of equation (2.13b), derived from the definition of dynamic viscosity, $[eqn]$ and assume a straight nozzle. They similarly find a direct proportionality between filament width and pressure drop but neglect the flow behaviour index, $[eqn]$ , obtaining $[eqn]$ . For the Pluronic F-127 hydrogel used, a good fit was obtained ( $[eqn]$ ) suggesting the material and its behaviour on the build platform could further affect the success of the fitting.

A straight nozzle was considered in the Nivea Crème experiments, but the minimum pressure necessary to extrude the material was extremely high, reducing the range of air pressures possible to explore on the bioprinter chosen. Furthermore, extrusion from the straight nozzle was relatively sporadic compared to the continuous flow achieved from the conical nozzles. As such, this nozzle geometry was excluded from the experiments.

To adjust to the experiments conducted, we specified a conical nozzle geometry in equation (2.12a), which presented a good agreement with the experiment, with $[eqn]$ . But, ultimately, having considered an arbitrary, axisymmetric nozzle, a more general expression was established which can be adapted to any given axisymmetric nozzle geometry.

Rheometric testing

3.1.5.

Since the model is fitted to the data through the optimization of the flow behaviour index, $[eqn]$ , and consistency index, $[eqn]$ , it provides an efficient alternative to gain an understanding of the flow properties of the material. Typically, a rotational rheometer is used to measure the viscosity of a material as a function of shear rate, enabling the fitting of the power-law model to the linear region of the double logarithmic curve to determine $[eqn]$ and $[eqn]$ as follows:

$$[eqn]$$

Despite standard rheological tests being applied, the flow behaviour and consistency indices for Nivea Crème are seen to vary significantly across the literature (over an order of magnitude, as listed in Table 3), while for the PVA-based hydrogel used in this study, a value for these indices is unreported. These discrepancies suggest a high sensitivity to the rheometric testing conditions. The flow behaviour index determined from the curve fitting of equation (2.12) to experimental data sits within the range established in the literature. This suggests that the model could be capturing an acceptable estimation of the true flow behaviour (with additional information provided in the form of the shape fidelity expected). However, we emphasize that the experimental conditions and fitting methods employed will influence the data generated and, thus, affect the value of $[eqn]$ and $[eqn]$ obtained. Nevertheless, values determined through rheometric testing can be introduced into the optimization problem as the initial guess, reducing the number of iterations before a local minimum is located and increasing the chances of finding the values corresponding to the best fit. To validate these hypotheses, rheometric testing needs to be conducted in the future.

Model limitations and extensions

3.1.6.

As a widely accepted model, the power-law model was used here, but it comes with its limitations. It is typically used to identify shear-thinning materials and, as such, is limited to a particular range of shear rates (generally 10^1^–10^4^ s $[eqn]$ [12]). Tu et al. [43] present a rheological characterization of Nivea Crème, highlighting the presence of a yield stress, $[eqn]$ (the minimum shear stress required to initiate flow). In fact, bioinks typically possess a yield stress, behaving elastically below and as a power-law fluid above this critical stress (a viscoplastic behaviour). With a focus on a power-law framework (when $[eqn]$ ), the model fails to capture this minimum shear stress, predicting flow even at very low shear stresses [44]. As such, for viscoplastic materials, the model may produce misleading results. While the power-law model may be sufficient for the prediction of suitable combinations of printing parameters (where only higher shear rates are relevant), for the purpose of optimizing shape fidelity, the Herschel–Bulkley model may be better suited [12]. To capture this yield stress behaviour, other models, such as the Herschel–Bulkley model, can be readily implemented into the proposed methodology. Georgiou [45] highlights that viscoplastic materials can slide even at shear rates below this critical rate. This suggests that the current no-slip condition between the nozzle wall and the fluid may need relaxing.

Owing to the nonlinear manner in which process parameters affect resolution, the parameters chosen in each case are unlikely to be the optimal combination. For example, on varying the air pressure, it is clear from the experimental results that there is an optimal choice (the point where the filament width, $[eqn]$ , matches the nozzle tip diameter) of which above or below the material is over- or under-extruded, respectively (as seen in Figure 6), but, as the nozzle diameter increases, this optimal choice of air pressure is seen to decrease. Paxton et al. [12] suggest that the nozzle travel speed, $[eqn]$ , should match the extrusion speed from the nozzle tip, but with its dependence on process parameters, such as air pressure, nozzle geometry, viscosity of the material (as captured theoretically in equations (2.8) and (2.10)) and temperature, the extrusion speed is rarely known in an experimental set-up. This emphasizes the necessity to return to the underlying physics of the system to derive a generalized theoretical mathematical model that enables a deeper understanding of how process parameters influence the complex rheological profile presented by a material [46].

Temperature control is crucial when bioprinting with thermoresponsive materials (such as the PVA-based hydrogel used in this study) since their sol–gel transition takes place at a certain temperature. Increased temperatures generally lead to a decrease in the viscosity of a bioink and, thus, correspond to a decrease in $[eqn]$ and an increase in $[eqn]$ ; the fluid characterization approaches the Newtonian limit ( $[eqn]$ ) [47]. A common methodology to capture this behaviour is via the expression of the indices as a function of temperature and is typically in the form of an exponential ratio of temperature (an Arrhenius-type relationship) [15,48]. Although this relationship is deemed suitable from a mathematical point of view, an inverse relationship of $[eqn]$ with temperature, as presented by Abu-Jdayil et al. [49], may better reflect the Newtonian limit at high temperatures. Ultimately, however, the choice of function for $[eqn]$ and $[eqn]$ is heavily influenced by the material and set-up used [50]. Furthermore, the model neglects to capture any thread deformation at the nozzle tip, where the velocity profile adjusts from pipe flow in the nozzle to extensional free-surface flow (e.g. die swell). Although die swell is often observed when extruding with hydrogels, the problem remains complex and minimal research has been conducted to quantify this behaviour [51]. Recent advancements in fused deposition modelling may propose strategies for advancement in extrusion-based bioprinting in the future [52,53].

Conclusion

Although the optimal combination of process parameters still remains platform- and material-specific, quantitative tools are being developed to minimize any between-laboratory biases. A user-friendly, adaptable model has been derived to predict the resolution of a line of material via extrusion-based bioprinting with an arbitrary, axisymmetric nozzle. In particular, the relationship between the filament diameter and the pressure drop, nozzle travel speed, nozzle geometry and flow properties of the material is determined. A PVA-based hydrogel and Nivea Crème were extruded pneumatically using a variety of air pressures and conical nozzle sizes, and the printed filament diameter was measured. Model predictions presented a good agreement with the experimental data, particularly when the shear-thinning property was captured, with quantified errors arising owing to the assumptions made and process parameters considered. Exploitation of the model enables a window of printability to be established, reducing the current dependency on print-and-test methodologies before an optimal combination of process parameters is obtained. By devising mathematical models that exploit the non-Newtonian fluid flow driving the material through the system, we hypothesize that a deeper understanding will be achieved of how process parameters influence this complex rheological profile. The insight gained from undertaking this work has the potential to guide future research in the automation and optimization of design for extrusion-based AM.