Bioink Temperature Inﬂuence on Shear Stress, Pressure and Velocity Using Computational Simulation

: Bioinks are usually cell-laden hydrogels widely studied in bioprinting performing experimental tests to tune their rheological properties, thus increasing research time and development costs. Computational Fluids Dynamics (CFD) is a powerful tool that can minimize iterations and costs simulating the material behavior using parametric changes in rheological properties under testing. Additionally, most bioinks have speciﬁc functionalities and their properties might widely change with temperature. Therefore, commercial bioinks are an excellent way to standardize bioprinting process, but they are not analyzed in detail. Therefore, the objective of this work is to study how three temperatures of the Cellink Bioink inﬂuence shear stress pressure and velocity through computational simulation. A comparison of three conical nozzles (20, 22, and 25G) for each temperature has been performed. The results show that shear stress, pressure, and velocity vary in negligible ranges for all combinations. Although these ranges are small and deﬁne a good thermo-responsive bioink, they do not generate a ﬁlament on the air and make drops during extrusion. In conclusion, this bioink provides a very stable behavior with low shear stress, but other bioprinting parameters must be set up to get a stable ﬁlament width. A.D.-P. and J.B.P.; Formal analysis, J.C.G.-B.; Investigation, J.C.G.-B., E.M.-S. and J.B.P.; Methodology, J.C.G.-B. and E.M.-S.; Project administration, J.C.G.-B. and J.B.P.; Software, J.C.G.-B., E.M.-S., A.C.M. and M.M.; Supervision, J.B.P.; Writing—original draft, J.C.G.-B.; Writing—review and editing,


Introduction
Additive manufacturing technology is currently contributing with many possibilities to tissue engineering. In this sense, 2D structures created by standard procedures of tissue engineering can evolve into complex 3D structures using bioprinting [1]. Specifically, bioprinting could produce these complex structures by superposing biomaterial layers with several biological compounds that finally can generate artificial tissues and organs [2]. Bioprinting can also minimize rejection risk when patient's cells are used in the creation of autologous tissues and/or organs [1]. Bioprinting is usually divided into four main technologies: micro-extrusion, inkjet, laser-assisted, and stereolithography [3]. However, properties such as versatility, printing speed, and the possibility of using high viscous materials with a high cell density make micro-extrusion the most used bioprinting technique [1,3,4].
Because of their high importance in the bioprinting process, several studies analyzed how different materials affect cellular survival [5,6], printability [7-9], curing or cross-linking [9][10][11][12], and shape fidelity [13,14]. The bioinks used in bioprinting are usually cell-laden hydrogels with these parameters are essential to determine whether cells survive. Despite using the same mathematical model, none of them reported the numerical methods used, so the reproducibility of their research is hindered. Furthermore, three main differences can be noticed between both works. Firstly, they provide different levels of detail for the mesh. Liravi et al. [40] used a very specific mesh generation with a remeshing procedure developed by Wilkes et al. [42] to properly generate droplets. On the contrary, Samanipour et al. [41] did not even mention the mesh they generated for their simulations. Secondly, the characterization of the material is done under different fluid assumptions. Liravi et al. [40] used a Carreau-Yasuda potential model to fit the rheological data of a non-Newtonian material, while Samanipour et al. [41] used a constant value for viscosity for simplicity. The former method is used to fit all data, at the expense of using a complex equation, while a simpler Potential Law cannot fit very low or very high shear rate values. Lastly, regarding materials, Liravi et al. [40] used several concentrations of a polysiloxane-based hydrogel extruded into the air. Samanipour et al. [41] used a GelMA hydrogel extruded into oil. Therefore, all previous simulations analyze the behavior of non-commercial bioinks and, as far as the authors know, no previous computational analysis has been performed for Cellink Bioink. This bioink is made of alginate with nanocellulose fibers and some authors have experimentally demonstrated the adequate bioprinting properties of this type of bioink (good rheological behavior, cellular viability or mechanical response) [43][44][45]. Among those authors, Müller et al. have also performed additional CFD computational simulations to study pressure and shear stress inside a needle [46].
Hence, the objective of this work is to analyze the Cellink Bioink behavior while bioprinting using computational simulations. Specifically, the pressure, velocity, and shear stress of this bioink are studied using three temperatures (15,25, and 37 • C) and three conical tips' geometries (20,22, and 25G) as inlet parameters. Additionally, an analysis of the bioink filament volume is performed to measure its stability undero nine different combinations of temperatures and geometries in simple conditions.

Computational Model
Three geometrical models were created and simulated in COMSOL Multiphysics 5.4a (COMSOL Inc., Burlington, MA, USA, 2018) through a 2D axisymmetric model and a Two-Phase Flow level set interface. Commercial 20G, 22G, and 25G conical tip geometries ( Figure 1) were selected for simulations. A 22G conical tip is the recommended tip for Cellink Bioink by Cellink in its bioprinting manual [46], while 20 and 25G were selected because they are the bigger and smaller most used sizes after 22G. All conical tips were modelled after experimental measurements using a caliper. The geometrical model and its measurements can be seen in Figure 2. left, where X is 0.30, 0.20, and 0.13 mm for 20, 22, and 25G, respectively. The conical tip was modelled using a trapezium and a rectangle, while the air was composed of a rectangle and a trapezium. Two different domains were considered in the geometrical model. The first domain was related to the nozzle where the hydrogel was placed, and the second domain was the outside of the nozzle and corresponds to the air where the bioprinting material was ejected (Figure 2 right). Models were meshed using COMSOL-optimized mesh generation for fluid dynamics with 2D triangular elements obtaining a total of 6353 elements with sizes ranging from 0.015 to 0.335 mm. The average skewness mesh quality is 0.9121 with a minimum value of 0.5715.     Level Set method is a Eulerian transport method capable of capturing the interface of two fluids and the changes of the interface because of motion. This method couples the track of the interface with a fluids dynamic of each fluid, expressed by incompressible Navier-Stokes equations where ρ is the density, u is the speed of the fluid, p denotes the pressure, I is the identity matrix, µ is the fluid viscosity, g is the gravity, and F st is the surface tension force calculated as where σ is the surface tension, φ is the contour line of the gas-liquid two-phase flow interface, and δ is the Dirac delta function formulized as follows Initially, the flow is assumed to be laminar, but a posterior verification is needed after simulation by calculating the Reynolds number.
The material used for the simulations was the commercially available Cellink Bioink, composed of alginate and nanocellulose fibers, with a density equal to 1000 kg/m 3 . Its surface tension was measured and calculated using a KRUSS G20/DSA10 drop shape analyzer, obtaining a value of 55.8 mN/m. Additionally, its viscosity/shear rate at 15, 25, and 37 • C was provided by Cellink and fitted to a simple viscosity Potential Law where µ is the dynamic viscosity (Pa·s), m is the fluid consistency index, . γ is the shear rate (s −1 ), and n is the flow behavior index. The selection of a simpler viscosity law pursues minimizing the computational cost, as long as the shear rate is within the proper values. The fitting parameters and goodness of the fit are listed in Table 1 and in Figure 3. All bioink data were introduced in COMSOL as a user defined material and air data was obtained from COMSOL material library.
In the Level Set method, Cellink Bioink is expressed by φ = 0, the air is expressed by φ = 1, and the level set interface (transition area between both materials) is expressed by φ = 0.5. The level set equation can be seen as the volume percentage of liquid in the gas-liquid two-phase flow [47].
Density and viscosity of the involved material can change in the interface following the level-set functions ρ = ρ air + (ρ bioink − ρ air )φ (6) µ = µ air + (µ bioink − µ air )φ (7) The track of the level set interface is described in the following equation when it moves under the velocity field u ∂φ ∂t + u·∇φ = 0 (8) Only the normal component of the velocity is needed because the level-set method considers the interface movement to be normal to itself. Therefore, Equation (8)   All bioink data were introduced in COMSOL as a user defined material and air data was obtained from COMSOL material library.
In the Level Set method, Cellink Bioink is expressed by = 0, the air is expressed by = 1, and the level set interface (transition area between both materials) is expressed by = 0.5. The level set equation can be seen as the volume percentage of liquid in the gas-liquid two-phase flow [47].
Density and viscosity of the involved material can change in the interface following the level-set functions The track of the level set interface is described in the following equation when it moves under the velocity field Only the normal component of the velocity is needed because the level-set method considers the interface movement to be normal to itself. Therefore, Equation (8) can be reformulated as follows The Level Set method needs its function to be a distance function within all simulation to depict the interface. To assure a correct depict of the interface, numerical reinitialization and stabilization terms are added to the level set equation. Therefore, the non-conservative reinitialized level set function can be formulated as ∂φ ∂t where φ is the contour line of the gas-liquid two-phase flow interface, γ is the reinitialization parameter (approximately the maximum value of the velocity field), and is the interface thickness parameter (usually half of the mesh size in the region). Here, the bioprinting material corresponds to the domain where φ < 0.5, and air corresponds to the domain where φ > 0.5. Boundary conditions were set as shown in Figure 2, they were stablished taking into account that all simulations are a 2D axisymmetric simulation and to make possible a realistic shear-thinning non-Newtonian flow. In this sense, the laminar flow inlet was set as a full developed flow, described with the following equations to assure tangential flow component on the boundary is zero (Equation (11)) and to set the real inlet value (Equation (12)) u − (u·n)n = 0 (11) [−pI where p is the user average pressure (N/m 2 ), in our simulations 15 kPa, as recommended by Cellink in their bioprinting manuals for Cellink Bioink [46], K is the viscous stress (N/m 2 ), and p in is the real inlet pressure (N/m 2 ). The outlet boundary condition was set as an open boundary condition with no normal stress to allow bioprinting material to fill this domain and air to leave. Additionally, wall boundaries were set with a non-slip condition, which means that the velocity on the walls is zero, and level set interface in the wall is described by Because the simulation is an 2D axisymmetric simulation, a symmetry boundary condition was set, as shown in Figure 2. It was modelled as a combination of Neumann and Dirichlet boundary conditions, allowing flow to not penetrate this boundary, and vanishing shear stress for an incompressible flow

Simulation
Nine simulations of 10 s with a 1 ms step were carried out for each geometry and temperature. Each simulation was composed of two study steps: Phase Initialization and Time-Dependent. The Phase Initialization step is in charge of obtaining all initial values of the Level Set method on every mesh element. In this sense, Phase Initialization is solved in COMSOL using the distance to the initial interface, D wi (m), and initializes the level set variable φ to ensure a smoothly variation between 0 and 1 (maximum and minimum values) and to minimize numerical instabilities. Then, this initial level set value is translated to the Time-Dependent step using the following expressions for the two different domains in domains initially filled with bioprinting material, and in domains initially filled with air, where φ is the domain reference (volumetric fraction), and is the interface thickness (m). For the Phase Initialization step, a stationary solver was used to calculate the level set initial values that later were used as t = 0 values by the time-dependent solver of the Time-Dependent step. Both the Phase initialization and the Time-Dependent steps used a Newton non-linear method in a fully coupled solver with a Parallel Direct sparse Solver (PARDISO). The Newton non-linear method is in charge of the successive iterative calculation of all coupled Fluids Dynamics and Level Set formulas described before. This method evaluates all non-linear expressions and the Jacobean on each iteration and assures that the calculation error between successive iterations is below the user tolerance (set at 10 −6 ). If the error is higher than tolerance, the damping factor is automatically changed, and the time step is reduced. The time step is calculated using a Backward Euler BDF method which is known for its stability and is the most recommended one for fluid transportation in COMSOL.
The PARDISO is a solver based on LU decomposition that tries to improve the sequential and parallel sparse numerical factorization performance. This method can use all computer cores to perform parallel calculations, reducing the simulation time to the detriment of computer expenditure.

Results and Discussion
Reynolds numbers were calculated for all simulations obtaining a maximum value of the order of 10 −4 , which means that the initial assumption of laminar flow is correct. In the same way, the shear rate of all simulations has been checked, ranging from 0.1 to 155 1/s. Therefore, the Cellink Bioink behavior curve can be modelled using the Potential Law.
All simulation errors were measured in the velocity field. The error between successive time steps was below 10 −6 in the 18th step with a later average error of 10 −13 . To achieve these low errors, a very small-time step (approximately 10 −6 s) is calculated by the BDF method, ranging the simulation time between 3 and 10 h.

Outlet Pressure
Simulated outlet pressure was measured using a line probe at the very end of the nozzle, specifically in the initial gas-liquid interface. Figure 4 shows all pressures for 20, 22, and 25G conical tips at 15, 25, and 37 • C.  It might be expected Cellink Bioink at 37 °C to have the lowest pressure for each one of the different conical tips. This supposition is based in the viscosity analysis, as a low viscosity fluid is supposed to provoke lower inner pressures when compared to a higher-viscosity fluid for the same conical tip. Specifically, in bioprinting materials, other authors as Bartnikowski et al. [48] have checked that hydrogels' viscosity usually decreases when temperature increases. However, there are some cases, such as Pluronic, that an increment in temperature provokes an increase in viscosity. This specific material contains (poly (ethylene oxide))100-(poly (propylene oxide))65-(poly (ethylene oxide))100. At a low temperature, it is an individual block of copolymers but when temperature increases, its internal structure changes and forms micelles that increase its viscosity [49]. In our case, Cellink Bioink is composed by alginate and nanocellulose fibers. Its behavior was analyzed in several works [26,43,45,50] and they all determine that the normal behavior of this kind of bioink is that the Three out of nine simulations present a low-pressure peak in their temporal evolution of maximum outlet pressure (25 and 37 • C for 20G and 37 • C for 22G). Additionally, the outlet pressure of 22G conical tip at 25 • C decreases at the end of the simulation, which implies that a low-pressure peak is being formed beyond the 10 s, similarly to the peaks in the other simulations. In this sense, to properly analyze and compare the pressure behavior of all simulations, a set of relevant points are identified. These points are: (1) the maximum value before the low-pressure peak, (2) the minimum pressure value, (3) the stable pressure value just after the peak, and (4) the pressure value at the end of the simulation. All plots were analyzed to identify these key points when available (for plots without peak, only the value at the end of the simulation is provided) and their pressure values are shown in Table 2. It might be expected Cellink Bioink at 37 • C to have the lowest pressure for each one of the different conical tips. This supposition is based in the viscosity analysis, as a low viscosity fluid is supposed to provoke lower inner pressures when compared to a higher-viscosity fluid for the same conical tip. Specifically, in bioprinting materials, other authors as Bartnikowski et al. [48] have checked that hydrogels' viscosity usually decreases when temperature increases. However, there are some cases, such as Pluronic, that an increment in temperature provokes an increase in viscosity. This specific material contains (poly (ethylene oxide)) 100 -(poly (propylene oxide)) 65 -(poly (ethylene oxide)) 100 . At a low temperature, it is an individual block of copolymers but when temperature increases, its internal structure changes and forms micelles that increase its viscosity [49]. In our case, Cellink Bioink is composed by alginate and nanocellulose fibers. Its behavior was analyzed in several works [26,43,45,50] and they all determine that the normal behavior of this kind of bioink is that the higher the temperature is, the lower the viscosity. Therefore, considering Table 2 and Figure 4, it is clear that this expected behavior is not obtained by all but one of the analyzed conical tips. More specifically, 22G conical tip is the only geometry where outlet pressure is higher at 15 • C than at 25 and 37 • C. This unexpected behavior is caused by the different inlet volumetric flow in all simulations. In this sense, the total extruded volume of Cellink Bioink along the time is shown in Figure 5. Additionally, the total volume values referring to the key points defined before can be seen in Table 3. Therefore, when pressure is selected as simulation inlet parameter, the volumetric flow cannot be further controlled. It depends on the fluid viscosity and the inner geometry where the fluid flows.
Therefore, two main considerations can be extracted from the pressure results. On the one hand, the geometry can change the outlet pressure value. In general, the bigger the outlet diameter is (the smaller the conical tip gauge), the lower the pressure is ( Figure 4). On the other hand, the temperature has a slight influence on pressure for the same geometry. Despite different bioink temperatures leading to different pressures, the average differences are around 600 Pa, so they can be considered negligible.
In order to properly compare our pressure results with those from bibliography, we have selected the values not considering excluding the low-pressure peak values. Thus, our maximum pressure values vary between 1143 and 1757 Pa. As far as the authors know, Reid et al. [33] have performed the only study that analyzed this pressure under experimental settings. They found that pressure ranged from 101 to 107 kPa using "a fluid with similar properties to blood", such as the bioink, but further information is not provided. They focused on the geometrical optimization of the inner nozzle geometry using different shapes and lengths and set up the inlet flow at 0.1 mm 3 /s. In this regard, the lack of information about the used bioink makes difficult to compare results. However, it seems that their high-pressure values are caused by the very small outlet diameter (60 µm) used.  Therefore, two main considerations can be extracted from the pressure results. On the one hand, the geometry can change the outlet pressure value. In general, the bigger the outlet diameter is (the smaller the conical tip gauge), the lower the pressure is ( Figure 4). On the other hand, the temperature has a slight influence on pressure for the same geometry. Despite different bioink temperatures  Although Lee et al.
[51] determined that alginate bioinks are not thermo-responsive and the viscosity variation is not noticeable, our results (Table 3) show that temperature has an important influence on bioprinting procedure. For a specific conical tip, an increase in temperature provokes the bioink viscosity to decrease, therefore generating a higher flow rate. Controlling the amount of extruded material is one of the key points in a successful bioprinting procedure. Most commercial bioprinters are pneumatic-driven and their main limitation is that the bioink flow cannot be precisely controlled [52]. To create any tissue, it is necessary to deposit a determined amount of material in a certain point, so, if temperature influences the bioink flow, other printing parameters, such as XY-speed, must be carefully selected to achieve proper results.
As mentioned before, there are three simulations with low-pressure peaks and a fourth one where a peak is foreseen. Comparing Figures 4 and 5, the low-pressure peaks are produced at the exact time an abrupt change in volumetric flow is observed. Those peaks are produced at different times depending on the geometry and the temperature: 2.31, 3.14 and 6.86 s for 20G at 37 • C, 20G at 25 • C, and 22G at 37 • C, respectively. At those instants, an initial droplet is formed (Figure 6), which would eventually fall, taking into account that the cross-area is reduced in such way that it can be considered that the droplet will fall. This droplet creation also depends on the distance to the printing substrate. In the light of the presented results, it seems that none of the simulated temperature/geometry configurations at 15 kPa will form a stable filament but droplets. Simulation time (t = 10 s) is not enough to even generate the initial droplet at the lowest temperature (15 • C) or lowest inner diameter (25G). The four images in Figure 6 show the simulations where a droplet is generated and drops (A, B, C) or is about to drop (D). Table 3 shows that the bioink volume to produce a falling droplet varies depending on the geometry. In this sense, the total volume of Cellink Bioink hanging from the conical tip is around 106 mm 2 on average for 20G, and 90 mm 2 for 22G. The formation of droplets depends on factors such as the volumetric flow, the viscosity of the material, the cross-section area and, in particular for hydrogels, the distance between the nozzle tip and the printing substrate. In our simulations, droplets are formed because two main reasons. First, both 20 and 22G have the same wall thickness (0.2 mm, Figure 2), but the total cross area of the conical tip outlet is 0.50 and 0.38 mm 2 for 20 and 22G, respectively. Since the total droplet volume before falling is dependent on the cross-area in contact as well as the surface tension and this is constant, the cross-area limits the maximum volume of the droplet. Secondly, the distance to the printing substrate, defined as "h" by He et al. [9], usually varies from µm to few mm in a standard bioprinting procedure, but we have set a larger h in our simulations to check whether Cellink Bioink is capable of creating a stable filament on the air. Therefore, in Figure 6A,B, it is shown the simulations at the exact time (3.14 and 2.31 s, respectively) when a visible reduction in the cross-area near the conical tip of 20G at 25 • C and 37 • C that will cause the droplet fall. Additionally, in Figure 6D, the falling droplet of 22G at 37 • C in t = 6.86 s is represented, while Figure 6C shows the 22G at 25 • C simulation in t = 10 s, in which the total extruded volume is not enough to generate a falling droplet. According to our results, none of the analyzed combinations can generate a stable filament flow on the air, due to either the generation of droplets or an insufficient extruded volume. Therefore, other printing parameters related to these printability and shape fidelity features, such as h and XY-plane speed, should be correctly selected to create such stable filament flow.  Finally, Figure 4 shows an increase in pressure since 3.01 and 4.10 s in 20G at 37 °C and 20G at 25 °C, respectively, until the end of the simulations. This final increase in pressure is provoked by the material accumulation on the printing substrate. This accumulation reaches de conical tip and causes the rise of pressure to continue to the extruding material. Therefore, to the existing extrusion pressure, the force needed to push away the already extruded material must be added.

Outlet Velocity
Simulated maximum outlet velocity was measured using the same probe as in pressure, at the very end of the nozzle. Figure 7 shows the maximum velocity for 20, 22, and 25G conical tips for 15, Finally, Figure 4 shows an increase in pressure since 3.01 and 4.10 s in 20G at 37 • C and 20G at 25 • C, respectively, until the end of the simulations. This final increase in pressure is provoked by the material accumulation on the printing substrate. This accumulation reaches de conical tip and causes the rise of pressure to continue to the extruding material. Therefore, to the existing extrusion pressure, the force needed to push away the already extruded material must be added.

Outlet Velocity
Simulated maximum outlet velocity was measured using the same probe as in pressure, at the very end of the nozzle. Figure 7 shows the maximum velocity for 20, 22, and 25G conical tips for 15, 25, and 37 • C along the time. Additionally, Table 4 presents all velocity values in the previously defined key points. It can be observed that there are high-velocity peaks produced at the same time as the low-pressure peaks. The velocity increment is produced by the droplet material pulling the bioink outside the nozzle. When the droplet is completely separated from the nozzle, the velocity peak disappears, and the previous value is restored. It would be expected that the outlet velocity depends on the outlet geometry for a certain inlet flow, as defined by continuum equation. In this sense, and as explained before, an inlet boundary condition of 15 kPa makes the bioink inlet flow to be dependent on its viscosity and the nozzle geometry. This dependency results from a variability in the flow for all simulations, as can be seen in Figure 5. Similar to the findings in the pressure, there are two important considerations for velocity. Firstly, for the same geometry, velocity changes with temperature. As explained before, the lower the temperature is, the lower the viscosity and the higher the extruded volume. An increase in volumetric Additionally, Table 4 presents all velocity values in the previously defined key points. It can be observed that there are high-velocity peaks produced at the same time as the low-pressure peaks. The velocity increment is produced by the droplet material pulling the bioink outside the nozzle. When the droplet is completely separated from the nozzle, the velocity peak disappears, and the previous value is restored. It would be expected that the outlet velocity depends on the outlet geometry for a certain inlet flow, as defined by continuum equation. In this sense, and as explained before, an inlet boundary condition of 15 kPa makes the bioink inlet flow to be dependent on its viscosity and the nozzle geometry. This dependency results from a variability in the flow for all simulations, as can be seen in Figure 5. Similar to the findings in the pressure, there are two important considerations for velocity. Firstly, for the same geometry, velocity changes with temperature. As explained before, the lower the temperature is, the lower the viscosity and the higher the extruded volume. An increase in volumetric flow leads to an increase in velocity according to continuum equation when cross-area remains constant. Therefore, temperature changes the extrusion velocity. Lastly, the geometry has also an important influence on velocity, and this influence is higher when the conical tip gauge increases. As shown in Figure 7 the influence of temperature on velocity is higher at lower gauges, while at higher gauges the effect of temperature is reduced and the cross-section area reduction plays a major role.
Regarding to other authors velocities, they obtained different results: velocities equal to 5.50 and 7.20 cm/s for a 60 µm diameter conical and needle tip, respectively [33], or 36.70 cm/s for a 28G gauge needle [36]. Again, differences between their results and our simulations might be explained by differences in viscosities, which cannot be assured, as authors do not provide this parameter for their bioink. Nevertheless, our results confirm that differences in cross-area have an important influence on velocities.

Shear Stress
Simulated shear stress is measured in the whole inner extruder domain using a surface probe. Maximum shear rate values are presented in Figure 8 and Table 5. important influence on velocity, and this influence is higher when the conical tip gauge increases. As shown in Figure 7 the influence of temperature on velocity is higher at lower gauges, while at higher gauges the effect of temperature is reduced and the cross-section area reduction plays a major role.
Regarding to other authors velocities, they obtained different results: velocities equal to 5.50 and 7.20 cm/s for a 60 μm diameter conical and needle tip, respectively [33], or 36.70 cm/s for a 28G gauge needle [36]. Again, differences between their results and our simulations might be explained by differences in viscosities, which cannot be assured, as authors do not provide this parameter for their bioink. Nevertheless, our results confirm that differences in cross-area have an important influence on velocities.

Shear Stress
Simulated shear stress is measured in the whole inner extruder domain using a surface probe. Maximum shear rate values are presented in Figure 8 and Table 5. Since shear stress depends on shear rate, and shear rate depends on velocity, it might be expected that shear stress behaves similarly to velocity. In any case, the analysis of shear stress must be done carefully for non-Newtonians fluids, as the variations in viscosity with shear rate might change the shear stress behavior. As can be seen in Figure 8, both shear stress and velocity have a similar temporal behavior.
Results show that shear stress varies with temperature and geometry. For a defined geometry, the shear stress decreases when temperature increases, and the temperature influence is reduced when the conical tip gauge increases. Additionally, high-shear stress peaks appear at the same time as pressure or velocity peaks. Nevertheless, an odd behavior of 25G at 15 °C, and 25G at 25 °C can be found, with an increment of shear rate at the end of the simulation, which might be caused by the high viscosity of the material at those temperatures where the bioink has a low velocity (shear rate). Since shear stress depends on shear rate, and shear rate depends on velocity, it might be expected that shear stress behaves similarly to velocity. In any case, the analysis of shear stress must be done carefully for non-Newtonians fluids, as the variations in viscosity with shear rate might change the shear stress behavior. As can be seen in Figure 8, both shear stress and velocity have a similar temporal behavior.
Results show that shear stress varies with temperature and geometry. For a defined geometry, the shear stress decreases when temperature increases, and the temperature influence is reduced when the conical tip gauge increases. Additionally, high-shear stress peaks appear at the same time as pressure or velocity peaks. Nevertheless, an odd behavior of 25G at 15 • C, and 25G at 25 • C can be found, with an increment of shear rate at the end of the simulation, which might be caused by the high viscosity of the material at those temperatures where the bioink has a low velocity (shear rate). The measurement of shear stress is done in the whole domain, so the maximum value may not be found in a fixed position. However, Liu et al. [30] observed that maximum shear stress is placed at the very tip of the conical tip nozzle. In this regard, the shear stress distribution in all simulations is very similar to the 22G at 37 • C (Figure 9).  The measurement of shear stress is done in the whole domain, so the maximum value may not be found in a fixed position. However, Liu et al. [30] observed that maximum shear stress is placed at the very tip of the conical tip nozzle. In this regard, the shear stress distribution in all simulations is very similar to the 22G at 37 °C (Figure 9). Many authors agree that shear stress is one of the main parameters to analyze and control in order to reduce cell death during the bioprinting process. Similarly to the finding of Yan et al. [28], our results specifically confirm that increasing pressures or increasing outlet diameter provoke wall shear stress to increase. On the contrary, our simulations do not coincide with Liu et al. [30], as we obtain lower shear stress with higher viscosity of the bioink. This different behavior might be caused by the different inlet used in our simulations and their experiments, because Liu et al. used a constant mass flow, while we have used a constant pressure. Müller et al. [45] also studied the shear stress using several nozzles with a very similar bioink, but in the same way as with Liu et al., their results are not directly comparable to ours due to the different geometries and inlet pressures. Hence, previous experimental tests are hardly comparable due to different boundary conditions or scarce definition of bioprinting parameters. On the other hand, shear stress can be found ranging from approximately 200 Pa to 20 kPa. In this sense, our results are quite similar to those obtained by Li et al. [37]. Specifically, those related to an inlet flow of 0.015 mL/s using an alginate hydrogel (0.41 Pa·s).
According to Blaeser et al. [29], shear stresses below 5 kPa do not have a harmful effect on cells (cellular viability over 96%). Although our results show high variability in shear stress values, all of them remain below 5 kPa. Therefore, low damage effect on cells is expected, which allows for high cellular viability. Many authors agree that shear stress is one of the main parameters to analyze and control in order to reduce cell death during the bioprinting process. Similarly to the finding of Yan et al. [28], our results specifically confirm that increasing pressures or increasing outlet diameter provoke wall shear stress to increase. On the contrary, our simulations do not coincide with Liu et al. [30], as we obtain lower shear stress with higher viscosity of the bioink. This different behavior might be caused by the different inlet used in our simulations and their experiments, because Liu et al. used a constant mass flow, while we have used a constant pressure. Müller et al. [45] also studied the shear stress using several nozzles with a very similar bioink, but in the same way as with Liu et al., their results are not directly comparable to ours due to the different geometries and inlet pressures. Hence, previous experimental tests are hardly comparable due to different boundary conditions or scarce definition of bioprinting parameters. On the other hand, shear stress can be found ranging from approximately 200 Pa to 20 kPa. In this sense, our results are quite similar to those obtained by Li et al. [37]. Specifically, those related to an inlet flow of 0.015 mL/s using an alginate hydrogel (0.41 Pa·s).
According to Blaeser et al. [29], shear stresses below 5 kPa do not have a harmful effect on cells (cellular viability over 96%). Although our results show high variability in shear stress values, all of them remain below 5 kPa. Therefore, low damage effect on cells is expected, which allows for high cellular viability.

Conclusions
In this work, several simulations have been done to study the impact of temperature and geometry in a commercial bioink: Cellink Bioink, its rheological data and inlet pressures were provided by the company bioprinting protocols. The simulation results demonstrated the suitability of this commercial bioink to be used in micro-extrusion bioprinting techniques, regardless of the temperature or the conical tip used (15,25, and 37 • C and 20, 22, and 25G). However, it is recommended to use this bioink at 37 • C, not with the aim of having a minimum shear stress, but in order to obtain the higher volumetric flow. A higher volumetric flow leads to higher bioprinting speed, so cells are under pressure for a shorter time. Additionally, shear stress obtained from simulations forecasts a proper cellular viability at all temperatures, according to previous studies for values under 5 kPa.
Despite the suitability of this studied bioink, simulations have been performed with several simplifications, such as not-defined wall friction, the use of a bioink without cells inside, and no interaction between the bioink filament and the printing substrate. Additionally, only conical nozzle geometries have been simulated. As future works, a comparative study of how conical and needle tips geometries should deeply analyze the effect of shear stress under different temperatures of the bioink. Furthermore, other bioprinting settings, such as different height (h) or XY-plane speed, should be also included in future simulations to analyze the bioink filament in dynamic conditions during its deposition on the printing substrate. Finally, experimental tests using this commercial bioink with a 3D bioprinter should be performed to validate simulation results with the actual behavior of bioinks.