Simulations of 3D bioprinting: predicting bioprintability of nanofibrillar inks

The bioinks presented and compared in this paper are based on cellulose nanofibril (CNF) dispersions in water and exhibit distinct differences in flow behaviour. The first bioink, ink 6040, is based on CNF dispersions and alginate, and is a less viscous ink with characteristics governed by surface tension forces. CNF was kindly provided by RISE Bioeconomy (Sweden) and alginate (SLG100, Mw = 150–250kDA, above 60% α − 1− glucuronic acid ) was purchased from FMC Biopolymers (Norway). 3 wt% CNF was mixed with 3 wt% alginate at a dry weight ratio of 60:40. Ink 6040 has been studied within the tissue engineering field and shown promising results as an ink with regard to printing and cell survival [21, 22, 39]. During processing of CNF from pulp, the pulp is subjected to either enzymatic or chemical pretreatment to increase the efficiency of the mechanical fibrillation and high pressure homogenisation. Ink 6040 was prepared from enzymatically pretreated CNF [40]. The second bioink, 4% CNF ink, was used as a comparison to ink 6040 due to its difference in viscoelastic properties. It consisted of a 4 wt% CNF dispersion in water, and is a more viscous bioink with characteristics driven by viscoelastic stresses. Instead of being enzymatically pretreated, the 4 wt% CNF had been pretreated by carboxymethylation [41] which leaves the CNF more charged, and thereby less aggregated and more dispersed in water.

2.1.1. Rheology measurements

All inks were printed with the same pattern (figure 1), nozzle geometry, height of nozzle and printing speed. Tissue engineering scaffolds with a porous architecture are preferred to allow for nutrients and oxygen to reach the cells embedded in the bioink [42]. Therefore, 3D bioprinting of grid-like structures are 3D bioprinted when assessing bioinks and their compatibility with cells [43–45]. To capture the most important features of 3D bioprinted grids, the pattern was designed to include straight paths, corners and cross overs. A conical nozzle with an outlet diameter of $0.42\,\mathrm{mm}$ was used. The distance between the outlet of the nozzle and the printing plate was $0.4\,\mathrm{mm}$ for the first layer (continuous line in figure 1) and $0.8\,\mathrm{mm}$ for the second layer (dotted line in figure 1). The printing speed, which is the movement of the printer head in the X- and Y-direction, was set to $10\,\mathrm{mm}\,{{\rm{s}}}^{-1}$. The flow rate is controlled by adjusting the applied pressure. The pressure was chosen so that a continuous flow was obtained and so that all inks were printed with a flow rate of 90–100 $\,\mathrm{mg}\,{\min }^{-1}$. This applied pressure varied for each ink due to the difference in viscosity and yield stress of the inks. Generally, an ink with a higher yield stress requires a higher pressure force to overcome the fluid stresses and accelerate the fluid ink. The specific flow rate for each ink was determined by weighing the mass of ink dispensed during $30\,{\rm{s}}$ at the pressure required to reach the aimed flow rate. The applied pressure and corresponding flow rate (averaged from four measurements) for each ink are presented in table 1. Directly after printing, the printed grid was photographed for comparison with the results from the corresponding simulation. Images of the printed grid from the top and from the side were used to measure the width (n = 50) and height (n = 20) of the printed lines by ImageJ software [46].

Zoom In Zoom Out Reset image size

Table 1.  Measured flow rates at pressure used while printing.

Ink	Pressure ($\,\mathrm{kPa}$)	Flow rate ($\,{\mathrm{cm}}^{3}\,{\min }^{-1}$)
Ink 6040 (CNF/Alg)	14	0.10
4% CNF ink	40	0.09

The simulation software used is IPS IBOFlow, which is an incompressible, segregated Navier–Stokes solver based on immersed boundary methods [47, 48]. It uses a Cartesian octree grid that can be dynamically refined to get a higher resolution of the fluid-air interface and the flow close to immersed objects. The software has been previously used to e.g. efficiently simulate robot application of sealing material and adhesives on automotive industry [49–52].

In the simulations a computational domain consisting of $10\times 10\times 1$ cubic cells is used, where each side of the cell is $1.0\,\mathrm{mm}$. The dynamic refinements of the fluid-air interface cells is applied four times, making the smallest cell size $0.062\,5\,\mathrm{mm}$. A constant time step of $1.0\times {10}^{-5}\,{\rm{s}}$ is used during all simulations. To mimic the outlet of the nozzle used in experiments, a circular moving immersed inlet condition of diameter $0.42\,\mathrm{mm}$ is used. To simulate printing of bioink, a fixed volumetric flow rate is applied to the inlet condition. The volumetric flow rate used for each bioink is presented in table 1. The circular inlet condition moves according to the printed grid design seen in figure 1 at the prescribed speed of $10\,\mathrm{mm}\,{{\rm{s}}}^{-1}$ in X- and Y-direction.

2.3.1. Rheology model

The rheology of the bioinks in IPS IBOFlow is modelled by a linear PTT-model [53]. This is a full stress tensor based model which captures the viscoelastic behaviour of the solder paste. The derived constitutive equation for the viscoelastic stress can be written as

$$\begin{eqnarray}&&F({\tau }_{{ij}}){\tau }_{{ij}}+\lambda {\mathop{\tau }\limits^{ \triangledown }}_{{ij}}=2\eta {S}_{{ij}},\end{eqnarray} \tag{ 1 }$$

where λ is the relaxation time of the material, η is the viscosity contribution to the elastic stress, ${\tau }_{{ij}}$ is the viscoelastic stress tensor, S_ij is the strain rate tensor and F is the PTT-model stress function. The stress function takes the nonlinear elongational behaviour of the viscoelastic stresses into account and is written as

$$\begin{eqnarray}&&F({\tau }_{{ij}})=1+\displaystyle \frac{\lambda \varepsilon }{\eta }{\tau }_{{kk}},\end{eqnarray} \tag{ 2 }$$

where ε is a parameter related to the elongational behaviour of the flow. In order for the model to work, the material related parameters (λ and η) need to be determined. Oscillation frequency measurements shown in section 2.1.1 give data of both storage modulus, G', and loss modulus, G'', over a wide range of frequencies, ω. Both G' and G'' can also be calculated as a function of the material parameters λ and η. The storage modulus is calculated as

$$\begin{eqnarray}&&G^{\prime} =\displaystyle \frac{\eta \lambda {\omega }^{2}}{1.0+{\omega }^{2}{\lambda }^{2}},\end{eqnarray} \tag{ 3 }$$

and the loss modulus is calculated as

$$\begin{eqnarray}&&G^{\prime\prime} =\displaystyle \frac{\eta \omega }{1.0+{\omega }^{2}{\lambda }^{2}}.\end{eqnarray} \tag{ 4 }$$

The complex modulus, G^*, is calculated as

$$\begin{eqnarray}&&{G}^{* }=\sqrt{G{{\prime} }^{2}+{G}^{{\prime\prime} 2}},\end{eqnarray} \tag{ 5 }$$

are used for fitting of the parameters against the provided experimental data.

2.3.2. Surface tension force model

The surface tension force is approximated by a body force across the fluid-air interface and is modelled using the continuum surface force method by Brackbill et al [54]. It is formulated as

$$\begin{eqnarray}&&{\vec{f}}_{\sigma }=\sigma \kappa {\rm{\nabla }}\alpha ,\end{eqnarray} \tag{ 6 }$$

where κ is the curvature of the liquid–gas interface, σ is the value of the surface tension between the two fluids and ∇α is the gradient of the volume fraction. The curvature is approximated as the divergence of the interface normal,

$$\begin{eqnarray}&&\kappa ={\rm{\nabla }}\cdot \hat{n},\end{eqnarray} \tag{ 7 }$$

where $\hat{n}$ is the interface normal, which is in turn approximated by the normalised volume fraction gradient

$$\begin{eqnarray}&&\hat{n}=\displaystyle \frac{{\rm{\nabla }}\alpha }{| {\rm{\nabla }}\alpha | }.\end{eqnarray} \tag{ 8 }$$

2.3.3. Quasi-dynamic contact angle model

The quasi-dynamic contact angle model takes into account if the three phase contact line, the line where fluid–fluid–solid meet, is advancing or receding. Since there only is dependence on the contact line direction, two basic equations govern this model

$$\begin{eqnarray}&&\theta ={\theta }_{A},\quad {\rm{if}}\quad {\vec{V}}_{{cl}}{\rm{}}\,{\rm{is}}\,{\rm{advancing}},\\ &&\theta ={\theta }_{R},\quad {\rm{if}}\quad {\vec{V}}_{{cl}}{\rm{}}\,{\rm{is}}\,{\rm{receding}},\end{eqnarray} \tag{ 9 }$$

where ${\vec{V}}_{{cl}}$ is the contact line velocity, θ is the resulting dynamic contact angle from the model, θ_A is the advancing contact angle and θ_R is the receding contact angle. The contact angle is implemented as a boundary condition for the continuum surface force method by altering the interface curvature at the three phase contact line [55]. The interface normal is otherwise calculated according to equation (8), but to alter the interface curvature at the three phase contact line, it is instead calculated as

$$\begin{eqnarray}&&{\hat{n}}_{f}={\hat{n}}_{s}\cos (\theta )+{\hat{n}}_{t}\cos (\theta ),\end{eqnarray} \tag{ 10 }$$

where ${\hat{n}}_{s}$ is the normal pointing outward from the solid surface and ${\hat{n}}_{t}$ is the tangential surface normal pointing along the surface, towards the centre of the heavier fluid.

Both of the tested bioinks, ink 6040 and 4% CNF ink showed a solid like behaviour in the oscillation frequency sweeps since the storage modulus was higher than the loss modulus (figure 2). Ink 6040 had a lower moduli and was therefore expected to flow more easily than 4% CNF ink.

Zoom In Zoom Out Reset image size

In figure 3 the experimental data of G^* and the resulting models of the parameter fitting are presented. The resulting values of λ and η for the respective bioinks are summarised in table 2.

Zoom In Zoom Out Reset image size

Table 2.  Material parameters of viscoelastic PTT-models.

Bioink	λ ($\,{\rm{s}}$)	η ($\,\mathrm{Pas}$)
Ink 6040 (CNF/Alg)	3.014	1.048 × 104
4% CNF ink	4.716	4.380 × 104

The rheology models of the complex modulus resembled the experimental data well enough to be used in the simulations. The tuned viscoelastic rheology model is then used in the framework for a 3D bioprinting simulation of the grid structure shown in figure 1. The same flow rate, printing speed, nozzle diameter and layer height is used in the simulation as for the experimental setup described in section 2.2. The results from the simulation are then compared to photos of experimentally printed grid structures. By using 4% CNF ink and ink 6040, which exhibit clear differences in viscoelastic properties, see figure 2, the ability of the framework to simulate different types of bioink was tested. The comparison for the 4% CNF ink, presented in figure 4, show close similarities between simulation and experiment. The uneven curvature of the printed line, as well as the thickening of the printed line before certain corners is accurately captured. The overhanging lines are quite distinct for the more viscous 4% CNF ink, which is seen in the simulations as well.

Zoom In Zoom Out Reset image size

For the less viscous bioink, ink 6040, the flow is more driven by surface tension forces. Therefore the surface tension model and the static contact angle model plays an important role to accurately describe the interface curvature and surface flow. And as seen in figure 5, the smoother surface of the less viscous ink is reproduced in the simulations by the surface tension model. The static contact angle model correctly models the surface flow by governing the motion of the three phase contact line on the substrate. Overall the agreement between experiment and simulation is good, and this can especially be seen on the rounding of corners. Both the experiment and the simulation yields less distinct overhanging lines compared to the more viscous 4% CNF ink.

Zoom In Zoom Out Reset image size

Figure 6 shows the average height and width of the printed line compared to the simulated line. The simulation continuously measures the height and width of the first segment of the grid structure, while the experimental data is from multiple measurement points along the printed line. The standard deviation together with minimum and maximum values are also presented distinguish differences between simulation and experiment more clearly. The results from the two experiments verify the accuracy of the numerical simulation and indicate a well functioning rheology model. The simulation of 4% CNF ink has very accurate mean values for both line height and line width. Minor differences can be seen when comparing the standard deviation and the min/max-values, where the fluctuations in line height and line width differ slightly. The fluctuations in line height are slightly larger, while the fluctuations in line width are smaller. The overall comparison for 4% CNF ink, both visual and measured, indicates an excellent agreement between simulation and experiment. For ink 6040, the measurements of line height and line width indicate a very good agreement as well, with nearly identical mean values between simulation and experiment. The simulation has a smaller standard deviation for both height and width, indicating less size fluctuation along the printed line. This is also seen in the min/max-values comparison, where the simulation has a higher minimum value and a lower maximum value when it comes to printed line width. However, the overall comparison, including both the visual and the measurements, show that the simulation framework is able to capture the characteristics of ink 6040 and produce accurate results.

Zoom In Zoom Out Reset image size

Knowledge of the height and the width of a line printed with a certain bioink is important for determining suitable printing settings. The height will govern the suitable thickness of each 3D printed layer. If the thickness is set higher than the actual printed line height, the nozzle will after a few layers loose contact with the previously printed line causing the bioink to be dispensed inconsistently. In contrast, a too low thickness will cause layers being dispensed into each other, or alternatively hinder the flow of bioink out from the nozzle. The width is used when setting the distance between each printed line. If the intention is to print interconnecting lines to form a printed sheet, setting a distance larger than the line width will print lines with spacing in between while a too low distance will print lines that overlap each other. Thus, the simulation tells us that for printing 4% CNF ink with the flow parameters provided in section 2.2 (10 mm s⁻¹, 100 mg min⁻¹), the layer thickness and line width should be set to $400\,\mu {\rm{m}}$ and $420\,\mu {\rm{m}}$, respectively. On the other hand, for ink 6040 the layer thickness should be set to $230\,\mu {\rm{m}}$ and the line width should be set to $780\,\mu {\rm{m}}$.

The simulation showed good correlation with the experimental data and was therefore further used to study how the printed lines of the bioink are affected by changing the printing speed and nozzle height.

3.3.1. Different printing speeds

The resulting layer thickness and line width presented in the previous section, are valid for a printing speed of $10\,\mathrm{mm}\,{{\rm{s}}}^{-1}$ and a flow rate of 90–100 mg min⁻¹. Assuming a constant flow is kept, with increasing printing speed the line width should theoretically decrease, due to the lower volume dispensed per time unit. This correlation has also been shown experimentally in studies on 3D bioprinting where the line width controls the printing resolution [56]. By simulating, the change in line height, as a result of varying printing speed, can be predicted in addition to line width. Simulation can thereby be used to study how the cross section changes depending on the ink and the printing speed as shown in figures 7 and 8 for ink 6040 and 4% CNF ink.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The ratio of width to height, does not vary proportionally when changing the printing speed. For 4% CNF ink, the ratio calculated from the average line width and line height decreases from 1.5 to 0.9 when increasing the printing speed, and the data for both height and width is shown in figure 9. This result show the difficulty in experimentally determining the optimal layer thickness and line distance to set while printing. Also the standard deviation and min/max-values of the line width and the line height gives information on the printability of the inks at different printing speeds. For ink 6040, the standard deviation in figure 9 is larger at 5 and $20\,\mathrm{mm}\,{{\rm{s}}}^{-1}$ indicating that the printed lines will have a more consistent shape throughout the printing at $10\,\mathrm{mm}\,{{\rm{s}}}^{-1}$ and thus result in better printability. The prediction for the best suited printing speed of 4% CNF ink is more inconclusive, but it probably lies somewhere in between 10 and $20\,\mathrm{mm}\,{{\rm{s}}}^{-1}$. This is because the standard deviation of printed line height is lowest at $20\,\mathrm{mm}\,{{\rm{s}}}^{-1}$, but the printed line width has the lowest standard deviation at $10\,\mathrm{mm}\,{{\rm{s}}}^{-1}$.

Zoom In Zoom Out Reset image size

The inks are intended for printing with cells for construction of synthetic tissues. There are several factors in the bioprinting process which may influence the cell survival; mixing cells with the bioink, osmolarity of the bioink, shear forces within the nozzle, diameter of the nozzle outlet, temperature, and, access to nutrients and oxygen. While dispensing the bioink, which is the simulated bioprinting step, there is a transition point from flowing vertically out from the nozzle to being dragged along the fixed printing plate. During this transition there is a build up of viscoelastic stresses on the outside of the nozzle tip, which may influence the cell survival. Therefore it is valuable that simulation can be used to predict these stresses, as shown in figure 10 where the distribution of viscoelastic stresses is visualised within a printed line. For both inks, the stresses are highest closest to the nozzle outlet and decrease further from the outlet along the printed path. Comparing the viscoelastic stresses between the two bioinks, the stresses are higher in the 4% CNF ink. This is an important factor for the characteristics of a 4% CNF ink printed line, whereas the surface tension forces are of greater importance for the characteristics of a ink 6040 printed line.

Zoom In Zoom Out Reset image size

3.3.2. Different nozzle heights

The nozzle height is the distance between the outlet of the printer head nozzle and the printing plate upon which ink is being dispensed. The measurements of height and width of simulated printed lines is presented in figure 11. It is seen that decreasing or increasing the nozzle height with $0.1\,\mathrm{mm}$, to $0.3\,\mathrm{mm}$ and $0.5\,\mathrm{mm}$, respectively, does influence the line resolution, but not to a very great extent. For 4% CNF ink there is a trend towards better line resolution with fewer fluctuations as the nozzle height is increased. Both the standard deviation and the difference between the minimum and maximum values are very low for a nozzle height of $0.5\,\mathrm{mm}$. The same cannot be said for ink 6040, where the results suggest that using a lower nozzle height of either 0.3 or $0.4\,\mathrm{mm}$ would be favourable. At the highest nozzle height, the difference between the minimum and the maximum value becomes considerably larger, and an increase in standard deviation is observed as well.

Zoom In Zoom Out Reset image size

Another factor which the nozzle height may influence is the printability of bioinks containing cells, since the distribution of viscoelastic stresses will change as shown by figure 12. For both inks, the viscoelastic stresses are clearly decreasing with increasing nozzle height, although a stronger trend is seen for 4% CNF ink. In the way of viewing changes to both nozzle height and printing speeds as tools to increase the printability of bioinks containing cells, by using them to reduce the viscoelastic stresses, the printing speed showed to have greater effect on the printability than the nozzle height.

Zoom In Zoom Out Reset image size

The 3D bioprinting simulation in IPS IBOFlow is controlled through a Lua script. The first part of the script contains the setup of the simulation, which involves reading options from a separate options file, importing nozzle definitions and defining the boundary conditions used in the simulation. The boundary condition used for all boundaries, except the negative Z boundary, is an outlet boundary condition. In the negative Z direction a no-slip wall boundary condition is utilised, and both advancing and receding contact angles are prescribed at this boundary for the dynamic contact angle model. The path of the nozzle and the speed at which it follows it, is setup after the nozzle definitions has been imported, and it consists of multiple position points, which describes the printing path, and the nozzle speed. After these steps, the simulation is setup. The last part of the script contains the initialisation of the nozzle, activation of dynamic grid refinements and execution of the simulation.

The nozzle definition is kept in a separate file and read by the main script. It contains definitions of both the nozzle and the bioink. The nozzle definitions are shape, which is usually a cylinder, and diameter, which varies from nozzle to nozzle. The flow out of the nozzle is also defined, as either mass flow rate or volumetric flow rate. The bioink is defined by entering density, rheology model and parameters of the rheology model. A viscoelastic rheology model is used for 3D bioprinting and the parameters of the model are explained in section 2.3.1.

The options related to the simulation is stored in an options file. It contains options related to time settings, such as time step length, simulations end time and how often to write VTK output, which contains the data of the simulation at that time step. It also contains options for the surrounding fluid in the computational domain, which is the density and viscosity of air. VOF related options such as surface tension and type of dynamic contact angle model is also set here. Options related to the computational domain is also found here, by defining the number of base grid cells in each direction, and the length of a base grid cell.