Flow and hydrodynamic shear stress inside a printing needle during biofabrication

We present a simple but accurate algorithm to calculate the flow and shear rate profile of shear thinning fluids, as typically used in biofabrication applications, with an arbitrary viscosity-shear rate relationship in a cylindrical nozzle. By interpolating the viscosity with a set of power-law functions, we obtain a mathematically exact piecewise solution to the incompressible Navier-Stokes equation. The algorithm is validated with known solutions for a simplified Carreau-Yasuda fluid, full numerical simulations for a realistic chitosan hydrogel as well as experimental velocity profiles of alginate and chitosan solutions in a microfluidic channel. We implement the algorithm in an easy-to-use Python tool, included as Supplementary Material, to calculate the velocity and shear rate profile during the printing process, depending on the shear thinning behavior of the bioink and printing parameters such as pressure and nozzle size. We confirm that the shear stress varies in an exactly linear fashion, starting from zero at the nozzle center to the maximum shear stress at the wall, independent of the shear thinning properties of the bioink. Finally, we demonstrate how our method can be inverted to obtain rheological bioink parameters in-situ directly before or even during printing from experimentally measured flow rate versus pressure data.


Introduction
Biofabrication, or bioprinting, is a novel technology aimed at applying common 3D printing techniques to fabricate living tissues.In extrusion-based biofabrication, the survival and functionality of printed cells strongly depend on the hydrodynamic stresses that the cells experience during printing [1][2][3][4][5].These stresses arise mainly from viscous shear forces in the printer nozzle and are thus directly related to the flow profile and the viscosity of the bioink [6][7][8][9][10][11] in which the cells are suspended.In an effort to reduce hydrodynamics stresses, shear thinning bioinks have been designed that exhibit a nearly flat velocity profile and correspondingly low shear rates in the nozzle center, in contrast to purely Newtonian liquids that develop a parabolic flow profile with higher shear rates throughout most of the nozzle [12][13][14][15][16][17][18][19].Consequently, with a consistency parameter K i and a dimensionless exponent n i .The i th interval is bounded by the shear rates _ G iÀ 1 and _ G i , and we demand Zð _ gÞ to be continuous at these bounds.This continuity condition together with a set of _ G i uniquely determines the power-law parameters K i and n i in every interpolation interval, as detailed in section S-1.2.
We note that this approach can be applied to any material described in terms of generalized Newtonian fluids, including yield stress fluids.Furthermore, our method includes cell-laden bioinks, as, on the one hand, the presence of cells has been shown to only slightly alter the materials' rheological behavior [5,30].On the other hand, the macroscopic rheology of a cell suspension, determined e. g. via shear rheometry or our capillary rheometry method presented in section 3, can be used as input for our method.

Governing equations
Analogously to the well-known Poiseuille flow of a Newtonian fluid [31, pp. 180 ff.], we assume a stationary, laminar, and pressure driven flow, with the velocity having only an axial component u depending on the radial position r.We consider a cylindrical channel and neglect entrance and exit effects.Applying these flow conditions, the incompressible Navier-Stokes equations reduce to the ordinary differential equation as shown in section S-1.3: Here, the constant pressure gradient G≔ @p @z ¼ Dp L is defined by the pressure drop Δp = p 0 − p L < 0 over a channel segment of length L. For a Newtonian fluid, i. e. Zð _ gÞ ¼ Z, integration of Eq (1) directly yields the well-known linear radial dependency of the shear stress: Inserting these piecewise profiles into the Navier-Stokes equation Eq (1) yields the following system of equations where i denotes the intervals as above: In order to solve this system of equations we assume the axial velocity to be continuously differentiable and the shear rate to be continuous.The flow shall further fulfill a no-slip boundary condition at the channel wall and have its maximum at the channel center.The mathematical solution to the system of equations Eqs (3) and ( 4) is detailed in section S-1.4 and section S-1.5 and can be summarized as follows: the shear rate profile is obtained by integrating Eq (3) over the radial position once.Inserting this solution into Eq (4) yields the velocity profile after another integration over r.Both integrations come along with integration constants that are determined employing the boundary conditions of the flow and the continuity conditions as stated above.

Results
The first equation Eq (3) can be rearranged and integrated once to obtain the shear rate profile in the i th interval: From this, the velocity profile is obtained by integrating over r, which ultimately yields (cf.(S-58)): Here, the newly introduced index k denotes the radial interval that contains the physical boundary of the channel, i. e. R k−1 � A � R k with the channel radius A.
The radial shear stress profile can, similarly to the Newtonian case, be derived from Eq (1), yielding the same linear behavior: This shows that the shear stress profile in a cylindrical channel is independent of the shear thinning properties of the material.Using the solutions for the shear rate (5) and the velocity (6), we derive mathematical expressions for the flow rate as well as the average velocity, shear rate, viscosity, and shear stress.Details of the derivation and the corresponding solutions can be found in (S-63), (S-66), (S-71) and (S-74), respectively (cf.section S-1.6).The flow rate or, equivalently, the average flow velocity determines the printing speed in 3D bioprinting processes.The average shear rate and shear stress can be used to estimate cell damage during printing [2,4] as detailed in section S-1.8 of the S1 File.We discuss the inclusion of possible wall-slip effects [35] in section S-1.9 of the S1 File.

Validation
To validate our method, we implement the presented algorithm in a Python [32] tool, included as S1 File together with an explanatory tutorial in section S-3 and available at https://github.com/sjmuellerbt/CYprofiles.Our tool performs the viscosity interpolation according to section 1.1 for a five-parameter Carreau-Yasuda fluid, given in Eq (9).The radial profiles for velocity, shear rate, viscosity, and shear stress and their respective averaged quantities are calculated after providing the printing parameters, i. e. the nozzle radius and the pressure gradient or an imposed flow rate.
We first validate our algorithm using an exact global mathematical solution for a simplified CY model.Next, we compare our algorithm with Lattice Boltzmann simulations for a general CY model using the open source software package ESPResSo [33,34], for which we extended both the CPU and GPU implementation with several inelastic viscosity models, including the CY model.We finally perform experimental velocity profile measurements in a microfluidic channel and confirm the theoretical prediction of our Lattice Boltzmann simulations.

Validation with global solution
We consider a simplified Carreau-Yassuda (CY) model of the following form where Z0 is the viscosity in the limit of zero shear rate and K is a time constant.For this model, an exact global solution to the NSE Eq (1) can be found as described in(S-79) and (S-81) (cf.section S-1.7).As shown in Fig 2, we find excellent agreement between this exact solution and the calculated profiles using our Python tool.

Validation with Lattice Boltzmann simulations
The general CY model [25] is given by where Z1 is the viscosity in the limit of infinite shear rates and the exponents a 1 and a 2 determine the shape of the transition between the zero-shear Newtonian plateau and the power-law region as well as the power-law behavior.For this general CY fluid a global mathematical solution to the NSE does not exist.We therefore compare our algorithm to Lattice-Boltzmann simulations using realistic bioink and printing parameters for a chitosan hydrogel taken from [26] with the following rheological parameters: Z0 ¼ 5807 Pas, K = 5.33 s, a 1 = 1.35 and a 2 = 0.87.The simulation setup consists of a 5 × 400 × 400 (x × y × z) box with a cylindrical boundary along the x-axis corresponding to a physical radius of A = 100 μm.The flow is periodic in xdirection thus leading to an effectively infinitely long channel.Further details of the Lattice-Boltzmann simulations are given in the S1 File.
The calculated and simulated flow profiles are in excellent agreement (Fig 3) thus validating our algorithm for a general CY fluid.

Validation with experimental flow profile measurements
As experimental proof, we measure the flow profile of an alginate solution along the centerline of a microchannel and compare our findings to Lattice Boltzmann simulations of the same geometry.
We prepare a 2.0% alginate solution by mixing 800 mg of alginate (Grindsted PH 176, Dupont, USA) in 50 ml Dulbecco's phosphate buffered saline under constant stirring overnight at room temperature together with yellow-green fluorescent beads (FluoroSphere carboxylated beads, Invitrogen, diameter: 0.5 μm).The alginate solution is injected under defined pressure into a polymethylmetacrylate microfluidic channel equipped with male mini Luer lock connectors (Darwin Microfluidics, France, internal volume: 8 μl) via a 15 cm long silicon tube (inner diameter: 1 mm).The channel has a length of 58 mm and a quadratic cross section of 190 μm × 190 μm, similar in size to the cross section of a typical printing needle.A square cross section of the channel was chosen to avoid optical distortions that would arise from the curvature of a cylindrical glass capillary in combination with the refractive index differences between glass and alginate.We visualize the flow of alginate using an epifluorescence microscope (DM4, Leica Microsystems, Germany) equipped with a CCD camera (frame rate: 100 Hz, Prosilica GE680, Allied Vision, Germany) and a 100 mW laser diode (473 nm).The microscope is focussed at the mid-section of the channel (height: 95 μm).
We perform measurements at a pressure of 300 kPa, close to actual printing conditions.The maximum flow speed in the center of the channel is around 2 cm s −1 , which is too fast to track the beads between successive frames.Instead, the velocity is estimated from the length of the linear streaks of the beads during exposure, as shown in We perform Lattice Boltzmann simulations of the pressure driven flow of the alginate solution in a square microchannel.The simulation setup consists of a 5 × 400 × 400 (x × y × z) box with plane boundaries in y-and z-direction forming a square channel that corresponds to the 190 μm × 190 μm microfluidic channel used in the experiment.The viscosity parameters were obtained using our capillary rheometry method described in section 3 as Z0 ¼ 3:65 Pas, _ g c ¼ 21:71 s À 1 and α = 0.67, according to equation Eq (10) below.Since we do not know the pressure drop across the mini Luer-lock connectors and tubings, we estimate the pressure gradient from the maximum flow speed measured at the center of the channel.Accordingly, we find that 79% of the total pressure drop of 300 kPa occurs across the 58 mm long channel, while the mini Luer-lock connectors and tubings account for the remaining 21%.
Fig 4b depicts the measured flow profile in comparison to our Lattice Boltzmann simulations.We see excellent agreement of the measured velocity profile and our numerical prediction.Further measurements of alginate hydrogels with different concentrations, as well as a chitosan hydrogel, at various printing pressures and in different channel geometries are included as S1 File.

Inverse application for a capillary rheometer
Not all laboratories working in bioprinting may have access to sophisticated rheometers for measuring the non-linear viscosity of their bioinks.Moreover, bioinks are often highly sensitive fluids with a large batch-to-batch variation, and the sample used for rheometer measurements may not behave in the same way as the sample used for the actual printing process.In this section, we show how our method can be inverted to perform in-situ capillary rheometry measurements using only a bioprinter and a standard laboratory scale.For this, we measure the pressure versus flow rate relationship [35] for a range of discrete pressure values.Using our Python tool, we then extract from this data the non-linear viscosity parameters of the bioink.

Experimental setup
We prepare a 2.5% alginate solution using the protocol described in section 2.3 without the addition of fluorescent beads.We measure the viscosity of the alginate solution at a temperature of 25˚C at shear rates between 0.01 s −1 and 100 s −1 using a cone-plate rheometer (DHR-3, TA-Instruments, USA).Alternatively, we measure the viscosity with a custom-made bioprinter that we use here as a capillary rheometer.A schematic of the experimental setup is shown in Fig 5 .The alginate solution is driven with a defined pressure (K8P electronic pressure regulator, Camozzi Automation, Italy) through a steel needle (21G blunt cannula #9180109-02, B-Braun, Germany, 28 mm length, 551 μm inner diameter).The pressure is increased stepwise from 20 kPa to 200 kPa in steps of 20 kPa.The driving pressure is measured with a pressure transducer (DRMOD-I2C-R10B, B+B Thermo-Technik GmbH, Germany), and the flow rate of the extruded alginate is measured with a precision scale (DI-100, Denver Instrument, USA).
We then fit the zero-shear viscosity, Z0 , the corner shear rate, _ g c , and the power-law shear thinning exponent, α, of a 3-parameter Carreau-Yasuda fluid to match the measured flow rate versus pressure relationship.The viscosity-shear rate relationship is given by which is derived from Eq (9) by omitting the infinite-shear viscosity (Z 1 ¼ 0), and introducing the corner-shear rate _ g c ¼ K À 1 as well as the single exponent α = a 1 = a 2 .Fitting is performed using a Marquard-Levenberg least-squares method implemented in the Python library SciPy, where the squared difference between the measured and the computed flow rate is minimized for each pressure level.The flow rate is computed according to (S-64) with the printing parameters mentioned above and N = 150 interpolation intervals between shear rates of 10 −6 s −1 to 10 8 s −1 .Since the inner diameter of the printer cartridge is large compared to that of the nozzle, we neglect a possible pressure drop along the cartridge.

Results
When measured with a cone-plate rheometer, the viscosity of a 2.5% alginate solution displays a pronounced shear rate dependency (Fig 6a), which is well described by a 3-parameter CY model according to Eq (10).Specifically, at shear rates below the corner shear rate _ g c � 17:8 s À 1 , the viscosity is approximately constant, with Z0 � 7:9 Pas.At shear rates above _ g c , the viscosity decreases according to a power-law with exponent α � 0.74.
When the same 2.5% alginate solution is extruded through a 28 mm long 551 μm diameter capillary, we find an over-proportional increase in flow rate with increasing pressure (Fig 6b).Specifically, a doubling in pressure causes an approximately 10-fold increase in flow rate.This experimentally measured flow rate versus pressure relationship is exactly predicted by our numerical solution (blue line in Fig 6b), adding further support to the validity of our algorithm.
If a rheometer is not available, the above procedure can be inverted to obtain the rheological properties of the bioink as follows: starting from a first guess of the CY parameters, the pressure versus flow rate is computed using our Python tool.Subsequently, the viscosity parameters are refined until the prediction matches with the experimental data as shown in This is likely due to shear rheometers not being able to achieve the large shear rates that occur under realistic printing conditions, while the capillary rheometer method intrinsically accounts for that.
A specific advantage of this capillary rheometry approach is that the experiment can be performed with the very same bioink that is currently in the printing cartridge prior to the actual printing process.Since most bioprinters are pressure controlled, i.e. the bioink is extruded through a printing needle with a constant pressure, the highly non-linear increase of the flow rate with increasing pressure makes it difficult to find the optimal printing parameters and to predict the material shear stresses during the printing process.Our algorithm solves both problems.

Conclusion
We presented a simple yet highly accurate algorithm to calculate the velocity and shear rate profiles for generalized Newtonian fluids, such as shear thinning bioinks, in cylindrical nozzles.For this, an arbitrary experimentally known viscosity-shear rate relation is split into a set  (10).This shear thinning behavior can be predicted (blue line) from an independent capillary rheometry experiment using our Python tool.(b) Flow rate versus pressure relationship of the alginate solution when extruded through a 28 mm long 551 μm diameter capillary (red squares).This relationship follows our numerical solution using 3 fit parameters (blue line).The flow rate versus pressure relationship can similarly be predicted (black line) from the viscosity values obtained from an independent cone-plate rheometer experiment shown in the upper panel, showing significant deviations with increasing pressure.https://doi.org/10.1371/journal.pone.0236371.g006 of continuous intervals described by power-laws.This includes the possibility to predict velocity and shear stress profiles in pure as well as cell-laden bioinks.In each interval, an exact solution for the shear rate and velocity is computed and connected to neighboring intervals to obtain a continuous smooth profile over the entire nozzle diameter.For the shear stress, the linear radial dependency independent of the fluid rheology was confirmed.In addition, the total flow rate as well as the average viscosity, shear rate and shear stress are also found mathematically.
We implemented our method as an easy-to-use Python tool for calculating the velocity and shear rate profiles for a Carreau-Yasuda fluid.To validate this tool, we compared our predictions to a mathematically exact global solution and to Lattice Boltzmann simulations for realistic chitosan hydrogels under typical bioprinting conditions.In both cases, we found excellent agreement.We further measured the velocity profile of an alginate solution in a microfluidic channel and found good agreement with Lattice Boltzmann simulations.
An important experimental application of our theoretical method is capillary rheometry.Here, the flow rate versus pressure relationship for a given hydrogel is obtained using a standard bioprinter.This data can then be fit to our theoretical predictions yielding the corresponding rheological parameters of the bioink.We illustrated this application for alginate and found good agreement with classical rheometer data.
Our method and the accompanying Python implementation provide a fast and simple tool to predict flow rates and shear stresses during bioprinting for a given bioink and thus will help to optimize printing parameters, especially for shear stress-sensitive living cells.

S-1. Theory
Our algorithm starts from an experimentally known viscosity-shear rate relation η( γ) and interpolates it by a series of power-law functions.This interpolation is subsequently used for writing down a similar series of Navier-Stokes equations which are solved first for the shear rate and then for the velocity profile.

S-1.1. Viscosity model
The viscosity-shear rate relationship of the bioink, or any other generalized Newtonian fluid, can be approximated by a continuous, piecewise function as depicted in figure 1a of the main text.In every interval the viscosity-shear rate relation is described by a power-law model with a consistency parameter K i having the physical unit Pa s n i , and a dimensionless exponent n i , according to the literature [1][2][3][4][5][6][7].We note that the shear rate can also be understood as dimensionless quantity, normalized to a constant shear rate of 1 s −1 without changing its numerical value.Doing so, the consistency parameter can be interpreted as a reference viscosity with the more meaningful physical unit Pa s.The i th interval is bounded by the shear rates Γi−1 and Γi .The condition ensures the continuity of (S-1) across the interval boundary Γi (i = 0, . . ., N − 1).Since real fluids usually exhibit Newtonian behavior for zero and infinite shear rates, we take for the power-law exponents in the first and last interval.We note that instead of the power-law interpolation, a linear interpolation would also be possible.However, since most bioinks show power-law shear thinning over a wide range of shear rates, a power-law interpolation is computationally more efficient when logarithmically-spaced shear rate intervals are used, as shown in figure 1a.

S-1.2. Determination of K i and n i
Starting from an experimentally known viscosity-shear rate relation η ( γ) which can be given either as raw rheological data or as a viscosity model with known parameters (e.g.[8,9]) such as the Carreau-Yasuda model, the consistency indices K i and exponents n i in each interpolation interval are determined as follows.The lowest and highest consistency indices are fixed by eq.(S-3) as Since rheological data often spans multiple decades, we choose an equidistant partitioning of the interval [ Γ0 , ΓN−1 ] on a logarithmic scale, as shown in figure 1a.Given the bounds of this interval and the number of interpolated points, the intermediate shear rates are given by (S-5) The parameters of the interpolating power-law functions η i ( γ) are found by inserting the known viscosity values at the interval bounds.Thus, the following system of equations needs to be solved: By division of the two equations, the power-law exponent is found to be Multiplication of (S-6) by (S-7) gives an expression for the consistency index: By inserting a functional form or raw data for η ( γ) into (S-8) and (S-9) the interpolation can be performed in the entire range of shear rates.

S-1.3. Governing equations
The Navier-Stokes equations to determine the flow field u read with the fluid mass density , the pressure gradient ∇p, the viscous stress tensor τ , and an external force term f .The viscous stress tensor is related to the viscosity and the strain rate tensor ε via where the strain rate tensor is defined as Here, ∇ u denotes the dyadic product of the gradient operator and the velocity vector and (∇ u) its transpose.The shear rate can be obtained as invariant of the strain rate tensor, i. e.In the following, we prove that the same holds for an arbitrary generalized Newtonian fluid.The strain rate tensor in cylindrical coordinates reads: where the notation ∂ x = ∂ ∂x denotes a partial spatial derivative with respect to the coordinate x.The "−" signs indicate the symmetric components of the tensor.The underbraced terms vanish due to the flow conditions.Thus, the viscous stress tensor reduces to a single component, (S-20) The components of the NSE yield: This shows that the pressure gradient has only a z-component.By applying the derivative ∂ z again on the remaining z-component of the NSE, we obtain: which shows that the pressure gradient is indeed constant and allows us to define where ∆p = p L − p 0 < 0 is the pressure difference along a channel segment of length L. Applying the flow conditions, the Navier-Stokes equations reduce to the ordinary differential equation ( 1): This equation is however still highly non-linear due to the dependency of η ( γ) on ∂ r u via the shear rate γ (cf.(S-31), ( 4)).
S-1.3.3.Ansatz and boundary conditions.Similar to the piecewise viscosity model in (S-1), we decompose the axial velocity u (r) into intervals: as illustrated in figure 1b.The essential difference between (S-1) and (S-29) is that the interval boundaries are determined by shear rates Γi for the former and by radial positions R i for the latter.While the interval boundaries for the viscosity Γi are an input quantity (see section S-1.2), the radial boundaries are determined a posteriori from the shear rate profile by the condition γ (R i ) = Γi (S-30) as will be shown in (S-42) below.The shear rate as a function of the radial position γ (r) is given by the first derivative of the velocity with respect to the radial position, i. e. γ (r) = − ∂u (r) ∂r , (S-31) and can also be written in the same piecewise manner: As for classical Poiseuille flow, we assume the common case of the velocity monotonically decreasing with the radial position.According to (S-31), the shear rate is therefore always positive.For u (r) to be continuously differentiable and finite at the channel center, the piecewise definitions of the velocity and the shear rate must be equal at the intermediate points, R i , i. e.
The flow shall further fulfill a no-slip boundary condition at the cylindrical channel wall at r = A, i. e.
To ensure the continuous differentiability of the axially symmetric flow field, the flow must have a maximum at the channel center, r = 0. Therefore, the shear rate has to vanish at this point: γ (0) = 0 (S-36)

S-1.4. Solution to the shear rate profile
Inserting the ansatz (S-29) and (S-32) into the Navier-Stokes equation (S-28) yields the following system of equations, ( 3) and ( 4), where i ∈ {0, . . ., N } denotes the intervals as above: (S-38) The first equation (S-37), ( 3), can be rearranged and integrated once to obtain γi where the c i are a set of integration constants that are determined next using the continuity conditions (S-34) and the boundary condition (S-36).
S-1.4.1.Determination of the integration constants of the shear rate profile.The integration constants c i can be shown to be zero using the complete induction proof described in the following.The base clause, (S-39) for i = 0 with the boundary condition (S-36) gives which is only fulfilled if the integration constant vanishes, thus, c 0 = 0. Assuming that c i = 0, c i+1 can be determined using the continuity condition (S-34), where the Γi are given.The equality Γi = γi (R i ) yields an expression for the radial position R i of the interfacial point between γi (r) and γi+1 (r), that can be inserted into the second part of (S-41) using (S-39): Employing the continuity condition of the viscosity model (S-2) gives an expression for the ratio of the consistency parameters, i. e.
Inserting (S-46) into (S-44) finally yields which is equal only if c i+1 = 0 thus completing the proof.With that, the final form for the shear rate profile in the i th interval is obtained as (cf.( 5)): Note that this solution reduces to the simple power-law model solution if the index i is dropped.

S-1.5. Solution to the velocity profile
The velocity profile is obtained by inserting (S-48), ( 5), into the second part of the system of differential equations (S-38), ( 4), and integrating over r: The integration constants ci are determined next using the no-slip boundary condition (S-35) and the continuity conditions for the velocity field (S-33).
S-1.5.1.Determination of the integration constants of the velocity profile.Since the number of intervals of the viscosity model N is independent of the choice of the flow parameters, G and A, and this choice uniquely determines the R i via (S-42), the outer channel boundary R is not necessarily located in the last interval of the velocity ansatz function u N (r).Instead, the radius of the channel lies in the k th interval, i. e.
where 0 < k ≤ N .Intervals with i > k whose boundaries R i lie beyond the channel radius A have no physical significance and are disregarded in the following.Consequently, the no-slip boundary condition applies to the k th interval: The integration constant can therefore easily be found as For i < k, the continuity condition for the velocity field (S-33) can be written as Assuming that ci+1 is known and rearranging this equation for ci yields where the underbraced terms can be identified as the shear rates at the interfacial position which are equal by the continuity conditions (S-34).Hence, Finally, inserting the expression for the known outermost integration constant, ck , the interior integration constants can be determined as Combining (S-50), (S-53) and (S-57), the velocity profile in the i th interval is given by (6):

Calculation of averages
In the following, we derive mathematical expressions for the flow rate as well as the average shear rate, viscosity, and shear stress.The flow rate or, equivalently, the average flow velocity determines the printing speed in 3D bioprinting processes.The average shear rate and shear stress can be used to estimate cell damage during printing [7,11].
S-1.6.1.Average velocity and flow rate.The cross-sectional average of the velocity field is given by (S-59)

S-1.7. Global analytical solution for a simplified CY model
In order to validate the algorithm presented above according to the procedure described in section 2 of the manuscript, we calculate a global mathematical solution to the Navier-Stokes equation (S-28), ( 1), for a simplified Carreau-Yasuda model.In the following, we derive an analytical solution for the flow profiles of a CY model (cf.( 9)) with the following simplification: The viscosity as a function of the shear rate is therefore given as (cf.( 8)) Using the same assumptions as in section S-1.3.1, the NSE yields: After a first integration and rearrangement of the equations one obtains and application of boundary condition (S-36) determines the integration constant as c 1 = 0. Therefore, the shear rate profile is given as: Inserting this result into (S-38) and integrating the resulting equation gives the velocity profile: The second integration is easily found by applying the no-slip boundary condition (S-35).Thus, the velocity profile is given by:

Estimate of force and deformation experienced by flowing cells
In this section, we provide a simple approach to estimate the force sensed by a cell and its resulting deformation during printing.Both quantities depend on the radial position r at which the cell transitions through the channel.Considering, in a first approximation, the cell as a sphere with radius R c , the shear force acting on it is given by the surface integral of the shear stress over the sphere.Due to the linearity of the shear stress with the radial position, and as long the cell is small compared to the channel, the force is obtained as the product of the surface area and the shear stress at the radial position of the sphere center, as detailed next.
With the shear stress given in (S-73), the shear force acting on the cell is obtained as: The result is depicted in figure S-1a.Another critical quantity is the cell deformation, which can be approximated using the shear stress and the mechanical properties of the cell.As a rough estimate, we assume the cell to behave linearly elastic with the stress-strain relationship given as where the strain ε quantifies the relative stretching of the cell.The Young's modulus E is chosen as 1 kPa to 10 kPa to cover the typical range of stiffness for cells [13].This leads to significant deformations as shown in figure S-1b which reiterates the importance of hydrodynamic shear forces in bioprinting.The corresponding flow profile is shown in fig. 3. We note that the assumption of linear elasticity predicts relatively large deformations, which would not be the case for a more realistic, strain-hardening behavior.Wall-slip effects are sometimes reported, especially for fluids exhibiting non-Newtonian behavior [12] or highly hydrophobic channel coatings.The general approach to include a velocity slip at a wall is to allow for a finite tangential velocity at this point or, equivalently, to shift (in the calculations) the channel wall further outwards by a distance known as the slip length.It is thus straightforward to incorporate slip effects into our calculations if the slip length is known.Alternatively, if the slip velocity next to the wall is known instead of the slip length, simply shifting upwards the computed no-slip velocity profile by a constant value represents a very good approximation.The shear rate would stay unchanged, and likewise the viscosity and the shear stress.Therefore, the inclusion of slip effects in our algorithm would be unproblematic.

S-2.1. Additional experiments
In this section, we provide more validation to our algorithm with experimental measurements.Using the same approach as detailed in 2.3 of the manuscript, we performed velocity profile measurements of 2 % alginate at 100, 200, and 300 kPa, of 3 % alginate at 300, 400, and 500 kPa, and 3 % chitosan at 300, and 400 kPa.These measurements as well as the velocity profiles calculated using our Lattice Boltzmann method are depicted in figure S-2, showing good agreement.Since the pressure drop in the connectors and tubings is unknown but depends on the rheology of the hydrogel, we assume a constant pressure drop before the microchannel of 21 % for 2 % alginate, 10 % for 3 % alginate, and 23 % for 3 % chitosan, respectively.Additionally, we measured the flow profile of 3 % alginate in a rectangular microchannel with 1000 µm × 200 µm cross section.Due to limitations of the field of view of the microscope, less than one half of the channel could be focused during the measurements.In figure S-3, the mirrored experimentally measured velocity profile is shown in comparison to our Lattice Boltzmann calculations.Due to connectors and tubing with a diameter of approximately the size of the microchannel, the pressure drop of 48 % is reasonable.

S-2.2. Error quantification
In the following, we present an error calculation for the different flow experiments.First, we calculate an averaged velocity profile for the measurements as well as our calculations by averaging the data in a given d-interval.   in the range of d = 75 µm to 500 µm).The average velocity in each bin is computed by Experimentally measured velocity profiles of 3 % alginate at 400 and 500 kPa in a 1 mm × 200 µm microchannel in comparison to numerical results using the Lattice Boltzmann method.The experimental data is mirrored with respect to the channel center.where N i is the number of data points in the i-th bin.Using this approach, we additionally calculate the standard deviation of the data from the average value as . (S-85) Using the averaged profiles, we calculate the relative error between measurement and Lattice Boltzmann computation as: The averaged profiles with a range of ±σ u,i are shown in figure S-4 for the square microchannel and in figure S-5 for the rectangular microchannel, where we find relative errors in the range of = 3.3 % to 13.5 %, and = 2 %, respectively.

S-3. Computational procedure and user guide
This section describes the structure and usage of the Python classes, found in the Supplementary Material, that implement the presented algorithm for the Carreau-Yasuda model.The first part gives an overview of the four implemented Python classes.
The second part provides a short user guide explaining how to use the classes for flow profile calculations.Note that their usage requires a Python (version 2 or 3) installation [14,15].The use of a Python IDE, e. g.Spyder [16], Thonny [17] or PyCharm [18], is optional but can be advantageous.The coloring in the code examples is the following: blue denotes classes, green denotes variables, and gray means a comment.For the Carreau-Yasuda model [19] the viscosity is given by η where η0 and η∞ are the viscosities in the limit of zero and infinite shear rates and K is a time constant with the unit s.Its inverse, K −1 = γc , is sometimes referred to as corner shear rate and determines the transition to the zero-shear Newtonian plateau.The exponents a 1 and a 2 determine the shape of the transition between the zero-shear Newtonian plateau and the power-law region as well as the power-law behavior.

S-3.1. Overview of Python classes
The tool uses four classes that hold the input parameters, perform the calculations and save or plot output data.They can be found in the file CYprofiles.py.The four classes are: (i) Analytical Viscosity(): instances of this class hold the parameters of the Carreau-Yasuda model and can calculate the viscosity for a given shear rate according to (S-87).
(ii) Interpolation(): instances of this class perform the interpolation of a given Analytical Viscosity() in a provided range of shear rates using the partitioning described in section S-1.2.
(iii) Printing Parameters(): instances of this class hold the printing parameters, i.e. the nozzle radius and the pressure gradient or the flow rate.
(iv) Profiles(): instances of this class perform the calculation of the velocity, shear rate, and viscosity profile for given Interpolation() and Printing Parameters() according to the presented algorithm.If a flow rate is provided, the corresponding pressure gradient is calculated iteratively to match the given flow rate.

S-3.2. User guide
This section is meant to serve as an explanatory tutorial for our Python tool.It will cover the two main steps necessary for calculating a flow profile: the viscosity interpolation according to section S-1.2 and the profile calculation according to the presented algorithm.The following code examples can be found in the file tutorial.py.
Step To perform the interpolation, the range of shear rates to interpolate and the number of (power-law) intervals is required: This section briefly summarizes the Lattice Boltzmann method and the extension we added to the open-source package ESPResSo [20].For an introduction into the Lattice Boltzmann method we refer the interested reader to the book by Krüger et al. [21].The Lattice Boltzmann equation for the multiple relaxation time scheme used in ESPResSo reads: M −1 ω M ij (f j ( x, t) − f eq i ( x, t)) (S-88) It describes the collision and streaming of the population distribution f i (i = 0, . . ., 18) during one time step ∆t.Here, c i are the discretized lattice velocities, M denotes transformation matrix that maps the populations onto moment space, ω is the diagonal relaxation frequency matrix, and f eq i denote the equilibrium population distributions.The relaxation frequency for the shear moments ω S is related to the dynamic viscosity of the fluid via [22] η = c 2 with the fluid mass density and the lattice speed of sound c s .The calculation of the viscosity according to the rheological model requires the local shear rate at each lattice node.Chai et al [22] showed that the local strain rate tensor can be obtained from the populations by ( c i ) α ( c i ) β M −1 ωM ij (f j ( x, t) − f eq i ( x, t)) . (S-90) The shear rate is then obtained as invariant of the strain rate tensor according to (S-13).
From the local shear rate, the viscosity according to the rheological model and the At the boundary of the cylindrical channel a bounce-back algorithm is applied to realize a no-slip boundary condition.The flow is driven by a pressure gradient along the zdirection, which is realized as external force density in the algorithm.

Fig 1 .
Fig 1. Viscosity and flow profile interpolation.(a) The viscosity-shear rate relationship of an arbitrary shear thinning fluid obtained e. g. from a rheometer measurement is interpolated by power-law intervals.The bounds of the intervals (vertical lines) are given by the intermediate shear rates, _ G i .By using a large number of intervals, any arbitrary viscosity-shear rate relationship can be approximated as closely as desired.(b) A long cylinder with uniaxial, stationary flow is used as a model for the flow of a bioink through a printer nozzle.The flow profile is split into radial intervals R i determined implicitly via the intermediate shear rates _ G i .https://doi.org/10.1371/journal.pone.0236371.g001

Fig 2 .
Fig 2. Validation with a mathematical solution.Flow profiles for the simplified CY model: the global mathematical solution and the prediction by our algorithm agree very well.The parameters are N = 1000, η 0 = 100 Pa s, K = 1.0 s and G = −1.95× 10 6 Pa m −1 .https://doi.org/10.1371/journal.pone.0236371.g002 Fig 4a, divided by the exposure time of 7 ms.

Fig 4 .
Fig 4. Validation with experimental flow measurements.Experimental measurement of the flow profile of a 2% alginate solution in a 190 μm × 190 μm microchannel.(a) Example micrograph of the bead tracking procedure.The velocity with respect to the lateral position is obtained as the length (yellow circles and labels) of the streaks divided by the exposure time.(b) The measured flow profile is in excellent agreement with our Lattice Boltzmann simulations.https://doi.org/10.1371/journal.pone.0236371.g004

Fig 5 .
Fig 5. Experimental capillary rheometer setup.Schematic of the experimental setup using a custom-made bioprinter as capillary rheometer: the bioink is driven through a syringe under defined pressure, and the flow rate of the extruded alginate is measured with a precision scale.https://doi.org/10.1371/journal.pone.0236371.g005 Fig 6b.The parameters obtained from the flow-rate versus pressure data (red squares in Fig 6b) are Z0 � 6:8 Pas, _ g c � 27:9 s À 1 , and α � 0.78 and differ only slightly from the parameters extracted from the cone-plate rheometer measurements.Accordingly, also only a slight difference between both parameter sets is seen in the velocity, shear rate and viscosity profiles shown in Fig 7. Also visible in Fig 6b is an increasing deviation of the flow rate versus pressure prediction for the cone-plate rheometer from the measured data with increasing pressure.

Fig 6 .
Fig 6.Comparison between capillary rheometer and cone-plate rheometer results.(a) Viscosity versus shear rate for a 2.5% alginate solution as measured with a cone-plate rheometer (data from 4 independent measurements, red squares) shows the pronounced shear thinning of a CY fluid that is well characterized by 3 fit parameters (black line) according to Eq(10).This shear thinning behavior can be predicted (blue line) from an independent capillary rheometry experiment using our Python tool.(b) Flow rate versus pressure relationship of the alginate solution when extruded through a 28 mm long 551 μm diameter capillary (red squares).This relationship follows our numerical solution using 3 fit parameters (blue line).The flow rate versus pressure relationship can similarly be predicted (black line) from the viscosity values obtained from an independent cone-plate rheometer experiment shown in the upper panel, showing significant deviations with increasing pressure.

Fig 7 .
Fig 7. Alginate flow profile from capillary rheometer and cone-plate rheometer data.Flow profiles of 2.5% alginate hydrogel with a pressure difference of Δp = −10 5 Pa and N = 150.There is only a slight difference between the flow profiles calculated from the viscosity parameters obtained with a cone-plate rheometer (black line) and our capillary rheometer (blue line).https://doi.org/10.1371/journal.pone.0236371.g007

εu
αβ ε αβ .(S-13)For a purely Newtonian fluid, the viscosity in (S-11) would be a constant.S-1.3.1.Flow conditions.Analogously to the well-known Poiseuille flow of a Newtonian fluid[10, pp.180  ff.], we assume a stationary, laminar, and pressure driven flow, with the velocity having only an axial component depending on the radial position.We consider a cylindrical channel and neglect entrance and exit effects.For the following derivation, a cylindrical coordinate system with a radial component r, an azimuthal component φ, and an axial component z, is employed.In these coordinates, the flow conditions read: (r, φ, z) = u e z (S-18) S-1.3.2.Constant pressure gradient.For a purely Newtionan fluid, the flow conditions (S-14)-(S-18) imply a spatially constant pressure gradient throughout the entire channel.

Figure S- 1 .
Figure S-1.(a) The force F acting on cells as radial function according to (S-82).(b)The strain according to (S-83) as a deformation measure for a linear elastic cell.The corresponding flow profile is shown in fig.3.We note that the assumption of linear elasticity predicts relatively large deformations, which would not be the case for a more realistic, strain-hardening behavior.

S- 1 . 9 .
Inclusion of wall-slip effects For the square channel, we choose a bin width of d bin = 7.31 µm (corresponding to N bins = 26 bins in the range of d = −95 µm to 95 µm), for the rectangular channel d bin = 16.3 µm (corresponding to N bins = 26 bins

Figure S- 4 .
Figure S-4.Averaged profiles from figure S-2.The gray area indicates the range of one standard deviation from the mean curve.gives the relative error calculated according to (S-86).

Figure S- 5 .u mm s − 1
Figure S-5.Averaged profiles from figure S-3.The gray area indicates the range of one standard deviation from the mean curve.gives the relative error calculated according to (S-86).

1 #
initialize Interpolation instance interpol= Interpolation(gamma0=gamma0, gammaN=gammaN, Ninterpol=Ninterpol, analytical=analytical) The calculation is then simply performed by executing interpol.calculateinterpolation() and the interpolation can be checked by plotting the calculated data via interpol.plotinterpolation() .To save the viscosity interpolation, one can use interpol.saveinterpolation(file) to save the data in viscosity-shear rate format and interpol.saveinterpolation parameters(file) to save the power-law parameters for all intervals in a file.# initialize Printing Parameters instance printparams = Printing Parameters( flowrate=flowrate, channelRadius=rchannel) S-4.Lattice Boltzmann algorithm for generalized Newtonian fluids local relaxation time according to (S-89) are computed at each lattice node and updated in every time step.In order to ensure simulation stability, we choose the time step globally according to Krüger et al [21, p. 273] as ∆t = c 2 s τ − global relaxation parameter τ = 1, and a reference kinematic viscosity ν .The latter is provided, for instance, by the upper Newtonian viscosity plateau of the corresponding CY model.