Investigation of cell-viability in the bioprinting process

IN

DEGREE PROJECT MECHANICAL ENGINEERING, SECOND CYCLE, 30 CREDITS

STOCKHOLM SWEDEN 2016 ,

Investigation of cell-viability in the bioprinting process

VARUNA DHARMADASA

Investigation of cell-viability in the bioprinting process

Varuna Dharmadasa

Supervised by:

Lisa Prahl Wittberg, Assoc. Prof.

Karl H˚ akansson, PhD Examiner:

Fredrik Lundell, Assoc. Prof.

Department of Mechanics Kungliga Tekniska H¨ ogskolan

Sweden

20/09/16

Table of Contents

1 Introduction 3

1.1 Background . . . . 3

2 Bioink Fluid characteristics 5 2.1 Rheological properties . . . . 5

2.2 Cell properties . . . . 9

3 Case setup 10 3.1 Geometry . . . . 10

3.2 Boundary conditions . . . . 12

3.3 Transient and Steady-state solver comparison . . . . 12

3.4 Grid study . . . . 16

3.4.1 Procedure of creating the grid . . . . 16

3.4.2 Simulation setup . . . . 18

3.4.3 Results . . . . 19

4 Post-processing 26 4.1 Modelling the cells in the flow . . . . 26

4.2 Cell deformation model in 2D . . . . 29

5 Results & discussion 32 5.1 Case 1: Near-wall region . . . . 32

5.2 Case 2: Outside the Near-wall region . . . . 37

5.3 Parameter study . . . . 44

6 Conclusion & Future work 47

Abstract

The human society has throughout history been faced with great chal- lenges. These challenges however hold the opportunity for us humans to learn and grow as a species in whole. As the human population is increas- ing, more attention has been focused to the medical field to deal with the challenge of curing and treating people in larger scales at a faster rate.

A particular challenge today is to meet the high demand in organ trans- plants. The number of human donors is scarce relative to the demand, and the transplantation is never guaranteed to be successful. Therefore allot of research is being conducted regarding the potential of 3D-bioprinting.

3D-bioprinting is an interesting field with a lot of potential where the

ultimate goal is to produce human organs for transplantation with the

use of a 3D-printer. However, there are still many cases in which the cell

viability in the bioprinting process is significantly low. If the reason is

biological or mechanical due to the strains in the flow through the bio-

printer is sometimes unclear. Here presented is an investigation on the

fluid stresses present in the nozzle of the bioprinter. This is done by simu-

lating the flow through the nozzle tip using CFD software and calculating

the principal stresses on the cells in the post processing step. By using

simple elastic deformation models the total area strain is calculated along

the particle track of a cell to predict how the cells may deform throughout

its particle track in the nozzle. It is found that the fluid stresses present

in the converging nozzle considered in this case are significant, and cannot

be excluded as a prime reason for the death of the cells in the bioprinting

process. Due to the non-newtonian character of the bioink considered in

this case, the cells close to the wall experience principal stresses signifi-

cantly higher than in the mainflow. Generally the character of the stresses

experienced by the cells along their particle tracks is observed to be highly

exponential, thus it proposed for future work to investigate how much the

maximum magnitude of the stresses at the outlet can be decreased by

considering shorter nozzle tips.

1 Introduction

1.1 Background

The application of technology within the realms of biology has so far made a remarkable impact on society and continues to do so at a faster pace. Today many people around the world are in need of a transplant of some sort, wether it be kidneys, liver, lungs or heart. Unfortunately failure of important organs such as those mentioned above have fatal consequences for many. Organ transplants are therefore in high demand, but the number of human donors are relatively scarce. Therefore searching for options outside relatives and human donors in general becomes more and more relevant. Using animal transplants for humans is not a new idea, however it is a difficult task since the outcome is many times a rejection from the human immune system. Therefore it is of great need to find new alternative ways to produce organ translpants. 3D-bioprinting is a very interesting and relatively new field that may contribute in this area.

3D-bioprinting is the process of producing a structure on which stem cells can grown and eventually become an organ. This is done using a 3D-printer, and similar to how ink is used for normal paper printers, a bioink is needed to do 3D-bioprinting. The Bioink consists of the ink itself, and the stem cells.

The bioink has to be of good printability whilst simultaneously being a good environment for the cells to survive and reproduce in, this is known as the biofabrication window and is discussed later.

Figure 1: The stages in 3D-bioprinting, figure from [9]

3D-bioprinting is a promising field. However it has been noted that the cell

viability during the bioprinting process is an issue that appears in some cases

without an evident explanation. It is suspected that the mechanical strains

on the stem cells may be what is causing the fatal outcome for the cells, but

biological reasons cannot be excluded. In this thesis an investigation will be

performed to see how the fluid mechanics of the bioprinting process may affect

stem cells in terms of stresses. First simulations of the flow of the bioink (without

the cells) through the nozzle tip,which is assumed to be the critical part of the

bioprinter, is carried out. Then a lagrangian particle tracking method called

DPM-model (Discrete phase method) is employed in FLUENT to inject and

track the cells in the flow. The cells are modeled as inert solid spheres. Four

cells are then injected into the flow from different radial positions at the inlet,

the first being very close to the wall (without interacting with the wall) and the

fourth being in the mainflow. Along the particle tracks of the cells the shear-

and strain-rates of the flow field are extracted to calculate the shear- and normal

stresses along the same tracks. Once the stress tensor is known at every point

along the way of the cells paricle track, the principal stresses are calculated. A

simple elastic cell deformation model is then employed to visualize how the cells

may deform in the flow. Finally a parameter study is done where important

parameters such as the elasticity of the cell and Poisson’s ratio (used in the cell

deformation model) is varied to see how much the results are affected. The stem

cell considered in this report is the induced Pluripotent Stem Cell (iPSC).

2 Bioink Fluid characteristics

In order to predict the behavior of the IPS cells in the bioprinting process, the carrier fluid has to be modeled correctly. In this case the carrier fluid is a bioink consisting of nanofibrillated cellulose (NFC) and alginate. Within the field of bioprionting a successful bioink must have qualities that meet the values within what is known as the biofabrication window. This basically means that both the biological requirement for the ink to be a environment in which the cell can survive and reproduce, and the requirement for the ink to be of good printability to be fulfilled. Combining the shear thinning viscous properties of the NFC with the fast cross-linking properties of the alginate allowing the ink to hold its shape after printing as it has a high viscosity at ”zero” shear-rates, a successful bioink can be produced [8]. This ink exhibits a shear thinning non- newtonian behavior meaning that the the viscosity of the fluid decreases as the shear rate increases. The bioink mainly considered in this study is one with proportions of 60% nanofibrillated cellulose and 40% alginate, comparison of the results will then be made with a bionk with NFC/alginate proportions of 80 : 20. The water content of the resulting bioink solution is 97.5%. The bioink here described will from here on be referred to as Ink6040.

2.1 Rheological properties

Data from rheological experiments conducted at Chalmers institute of technol-

ogy illustrate the shear thinning properties of Ink6040 in the figure below (Note

that the scale of the figure is logarithmic). For very low values of shear-rates,

Ink6040 has a viscosity of 62 200 Pa s, and for shear rates at 10

s

experimental

data shows that the ink has a viscosity of 0.3 Pa s. To describe the shear-thinning

properties of Ink6040 in a continuous manner a curve was approximated to the

datapoints.

Figure 2: Ink6040 viscosity variation with shear rate (data from private com- munication, plots of data can be found in [8]).

The relationship between the viscosity of Ink6040 and the shear rate it is exposed to can be described by

η = 109.73 · ˙γ

. (1)

Note however that experimental data only exists for shear rates up to 10

s

, thus the above relationship may be inaccurate for shear rates above this value.

Note that water has a constant viscosity of 8.9 × 10

Pa s, and since Ink6040

consists of 97.5% water, it is unreasonable that the viscosity of Ink6040 would

go below the viscosity of water for high shear rate values.

Figure 3: Ink6040 viscosity compared to water viscosity.

The decreasing manner of the viscosity with shear rate depicted in the fig- ure above is a property of shear-thinning fluids. An example of a shear-thinning fluid would be modern day paint. When the paint is applied to a surface using a brush, shear is applied which causes the viscosity of the paint to decrease and make the paint thin out on the surface evenly. As the brushing motion is stopped, the paint regains its viscosity therefore does not drip. In the analysis below, it is assumed that the viscosity of Ink6040 varies according to the rela- tionship given above until shear rate values reach 10

s

(where the curve ap- proximation of Ink6040’s viscosity intersects water’s viscosity, see figure above), and for shear rates greater than this value the viscosity is assumed to be that of water (8.9 × 10

Pa s).

The shear stress variation due to shear rates up to 10

s

can then be approximated using the relation above as

τ = η · ˙γ = 109.73 · ˙γ

· ˙γ = 109.73 · ˙γ

. (2)

Since Ink6040 is of pseudoplastic character, the shear thinning of the ink

increases as it is subjected to stress, which also means a decrease in the viscosity

of the ink. The Shear stress, due to the assumption of Ink6040 having the same

viscosity as water for shear rates above 10

s

10

[

], increases at a higher rate

after this value and reaches 10 kPa at shear rates just above 10

s

. In regions

of the bioink where the shear rates are high (higher than 10

s

, the shear

stress will be approximated by the relation

τ = η

· ˙γ = 8.94 × 10

· ˙γ (3) In the sections to come it is investigated how the shear rate varies in the bioink when flowing through the converging nozzle.

Figure 4: Shear stress dependancy on shear rate for Ink6040. Note that the abrupt change in the derivative of the shear stress is due to the viscous assump- tion regarding the viscosity of water for shear rates above 10

s

.

Since Ink6040 consists of 97.5% water, it is modeled as fluid with the density of 998.2 kg m

with the shear thinning properties explained above in FLUENT.

Among the viscosity models in FLUENT, the model that best describes the viscous behavior of Ink6040 is the non-newtonian power law given by

η = k · ˙γ

e

, (4)

where the consistency index k and the paramenter n can be found from (1) to

be 109.73 and 0.154 respectively. Note that the parameter n is a measure of

how much the fluid deviates from being Newtonian (i.e. if n = 1 the fluid is

Newtonian). In this study the viscosity of Ink6040 is assumed to be temperature

independant, therfore T

= 0. In FLUENT it is also required to state the

upper and lower limit of the viscosity of Ink6040, from the experimental data

above and the assumptions regarding the water of viscosity discussed before,

η

= 62 200 Pa s and η

= 8.94 × 10

Pa s.

2.2 Cell properties

Depending on the cell type, the size and shape of the cell will vary. However, different cell types may also imply different mechanical characteristics. This could mean that different cells behave or deform differently when exposed to the same stress. Therefore it is important to properyly define the cell type in order to model the impact from the fluid in which it is suspended. As stated before the cell type of interest in this study is the induced pluripotent stem cells (iPSC).

Pluripotent stem cells are stem cells that have the potential to differentiate into any type of cell ranging from lung cells to blood cells to epidermal cells (any cell that is present in the formation of the embryo). The difference between iPSC’s and regular pluripotent stem cells is that the former are pluripotent stem cells that have been artificially derived by forcing an expression of specific genes and transcription factors in adult cells. Below is a table with the characteristics of the iPSC that will be used in this study. The density of the iPSC is computed by approximating the shape of the cell as a sphere with a diameter of 12 µm, and a weight of 10

kg.

Cell parameter Value Diameter d

10 to 15 µm

Density ρ

1105 kg m

Stiffness E 0.9 to 1.3 kPa

Note that both the cell diameter and the cell stiffness vary from cell to cell

within the range given in the table above, however a diameter of 12 µm and a

stiffness of 1.3 kPa will be used in the main study. A parameter study is the

conducted in which the stiffness of the cell is varied to see how the results may

vary.

Figure 5: Illustration of the different ways bioink can be extruded through a microextrusion bioprinter [9].

3 Case setup

The most frequently used bioprinting technology today include inkjet, microex- trusion and laser-assisted technology. In this case a microextrusion bioprinter is considered. The microextrusion bioprinter extrudes continuous beads of the bioink with the aid of either pneumatic or mechanical dispensing system (see Figure 5 above). The part of the bioprinter which is relevant in terms of fluid mechanics, is the syringe through which the bioink is extruded. In this section it is described how the geometry of the syringe is modeled. Next the boundary- conditions applied on the geometry is discussed. The flow through the syringe is assumed to be of steady state, to validate this theory a comparison between the transient and the steady-state solver was made which is also presented in this section. lastly, the methodology behind creating and choosing the appropriate grid for the simulations is presented.

3.1 Geometry

The syringe through which the bioink is extruded consists of a tube with a converging nozzle at the end (see Figure 6). The converging nozzle is suspected to be the critical part through which the mechanical stress on the cells steadily increase and cause fatal injury to the cells. Therefore the converging nozzle is of most concern. However, the bioink is extruded from upstream the nozzle in the tube either pneumatically or mechanically, where it then flows into the nozzle and accelerates. Therefore the bioink already has a certain velocity profile when entering the nozzle. Thus the tube cannot be neglected when modelling the geometry of the syringe.

The dimensions of the converging nozzle is available through measurement

of the actual nozzle which is interchangeable of the printerhead. It is assumed

that the bioink is extruded one nozzle length upstream the end of the tube

(where nozzle is attached). Therefore only this part of the tube will be mod-

eled, although it may extend further up. For the converging nozzle, the length

l

= l

= 31mm

l

= 31mm D

= 5mm

D

= 0.41mm

Figure 6: Schematic showing the dimensions of the modeled syringe (not to scale). l

is the length of the converging noz- zle, and l

is assumed to be of same length as the nozzle.

is measured to be 31 mm with a inlet diameter of 5 mm and an outlet diameter

of 0.41 mm.

∆P = 20 kPa

P

V

V

= V

P

Figure 7: Schematic of how the bound- ary conditions are set for the tube and the converging nozzle respectively.

The inlet condition for the tube is set to 20 kPa above operating condi- tions (atmospheric pressure), and at- mospheric pressure at the outlet. From this simulation the velocity profile at the outlet is then extracted and used as inlet condition for the converging noz- zle. The outlet condition for the nozzle is atmospheric pressure.

3.2 Boundary conditions

As stated in the beginning of the section, the converging nozzle is of most concern. If however the pressure boundary condition is applied directly at the inlet of the converging nozzle, the volume flow rate will be unreasonably small (since the ink is accelerated from more or less zero velocity). Therefore it would be more reasonable to apply the pressure boundary condition at the inlet of the assumed tube so that the flow already has a velocity when it approaches the converging nozzle. In order to reduce the computational cost, the tube is separated from the converging nozzle. This way it is possible to generate a coarser mesh for the tube, since only the tube outlet velocity is of interest, and a fine mesh for the converging nozzle since this is assumed to be the critical part on which this thesis is focused on.

By simulating the flow through the tube with a pressure inlet condition, one can get a reasonable estimate of the velocity profile at the inlet of the nozzle.

The velocity profile at the outlet of the tube is then extracted and applied as the inlet boundary condition for the nozzle in FLUENT. The pressure inlet condition at the inlet of the tube is set to 20 kPa above atmospheric pressure as this setting produced best bioprinting results at Chalmers. The pressure outlet condition at the outlet of the nozzle is set to atmospheric conditions of 101 325 Pa (default in FLUENT).

3.3 Transient and Steady-state solver comparison

For the comparison of the results between the transient and the steady-state

solver in FLUENT, only the converging nozzle will be studied since this part

is of main concern. When the bioink is flowing through the nozzle, the flow

is for the majority of the time likely to be in a steady state compared to the

initial transient phase when the bioink is first extruded through the nozzle due

to the backpressure. Therefore it is reasonable to assume that most of the dam-

age inflicted on the cells during the flow probably occurs in the steady-state

phase of the flow. In order to compute the stress on these cells numerically,

one can therefore simulate the flow in steady-state directly in FLUENT. How- ever the steady-state solver in the software neglects higher order time terms whereas the transient solver does not, hence the steady state approached with the time-marching solver will be more accurate. It is also interesting to inves- tigate whether or not the inflow velocity profile has a significant effect on the transient solution compared to the steady state solution. Hence a study is done to determine how consistent the results using the transient solver are with the steady-state solver in FLUENT.

As the flow is assumed to be incompressible, it is simulated in FLUENT by solving the mass conservation and momentum equations. The mass conservation law is given by

∂ρ

∂t + ∇ · (ρ~ v) = 0,

note that as the flow is modeled to be incompressible, the first term dissap- pears and we end up ∇ · ~ v = 0. The momentum equations being solved for the flow is given by

∂(ρ~ v)

∂t + ∇ · (ρ~ v~ v) = −∇P + ∇ · (¯ τ ) + ρ~ ¯ g + ~ F ,

where P is the static pressure, ρ~ g is the gravitational force and ~ F is any external body forces (for example the effect of the cells on the continuous fluid).

The stress ¯ τ is defined as ¯

¯ ¯ τ = µ h

(∇~ v + ∇~ v

) − 2 3 ∇ · ~vI i

,

where µ is the viscosity and the second term represents the effects of volume dilatation, but since the flow in this case is incompressible the above equation simplifies to

¯ ¯ τ = µ h

(∇~ v + ∇~ v

) i .

Note also here that for Ink6040 is non-newtonian, therefore the viscosity µ will be given according to the relation discussed in previous sections where it is a function of ∇~ v + ∇~ v

.

The flow of Ink6040 through the nozzle is a low speed incompressible flow, therefore the pressure-based solver is used both for the steady state and transient case. When using the pressure-based solver a numerical algorithm is required to derive a pressure equation (or pressure correction) from a combination of the continuity and momentum equations. For the steady state solver, the default SIMPLE (Semi-Implicit Method for Pressure-Linked Equations) algorithm is picked as the pressure-velocity coupling algorithm. However for the transient case, the PISO (Pressure-Implicit with splitting of Operators) algorithm is used since the flow is initially unsteady.

As can be seen from the results below (Figure 8 and 9), the results from the

transient solver converge to the results of the steady-state solver. The time step

was calculated to be ∆t = 2.5 × 10

s by dividing the smallest dimension of the smallest element in the mesh by the max velocity at the inlet. The velocity profiles converge much smoother at the outlet with the transient solver compared to halfway down the nozzle. However at t = 1.25 × 10

s the transient results both at the outlet and halfway down the nozzle have converged to the steady- state solver results.

-1.5 -1 -0.5 0 0.5 1 1.5

[m] _×10^-3

-8 -7 -6 -5 -4 -3 -2 -1 0

[m/s]

Velocity profiles halfway down the nozzle

Steady state solution Transient solution at t=2.58e-03 Transient solution at t=4.35e-03 Transient solution at t=1.25e-02

Figure 8: Comparison between velocity profiles halfways through the nozzle

between the steady-state and transient solver at different times.

-3 -2 -1 0 1 2 3

[m] _×10^-4

-250 -200 -150 -100 -50 0

[m/s]

Velocity profiles at outlet

Steady state solution Transient solution at t=2.58e-03 Transient solution at t=4.35e-03 Transient solution at t=1.25e-02

Figure 9: Comparison between velocity profiles at outlet of nozzle between the

steady-state and transient solver at different times.

3.4 Grid study

In order to make reliable simulations, it is important that the results are not de- pendant on the grid domain. Therefore a grid study was conducted investigating how the velocity profiles vary with increasing number of nodes and different grid strategies. The velocity profiles were compared for the different grids halfway down and at the outlet of the converging nozzle. Due to symmetry, only half of the nozzle is modeled. Note that the grid study below only concerns the nozzle, and not the tube part.

3.4.1 Procedure of creating the grid

ICEM CFD was used to create the grids. The region of the flow near the wall is of particular interest since this is a critical region with high shear rates and therefore also the most harmful region for the IPS cells. The dimensions of the nozzle are such that it has a height of 31 mm, a inlet radius of 2.5 mm and an outlet radius of 0.205 mm. In order to create a grid that is more resolved near the walls compared to the center of the nozzle, an ”ogrid” blocking strategy was implemented. This way it is possible to tailor the number of nodes and mesh laws in the region near the wall such that the grid becomes denser closer to the wall. Initially, three grids consisting of 226 000, 725 000 and 1 646 725 nodes respectively, were created. All three grids were created using the strategy stated above. A mesh law in the flow direction was implemented to increase the number of nodes towards the outlet since higher accuracy is needed as the flow accelerates within smaller regions.

Next will be discussed how these simulations were carried out and the results

that were obtained. The software used for simulations is FLUENT. Note that

Fluent is an unstructured solver, thus the grid is saved as unstructured and

exported to FLUENT as .msh extension file. Since the geomtry of the nozzle is

rather simple, the elements however align quite structured and the normal faces

of the elements are aligned with the flow direction.

Figure 10: The coarse mesh of the outlet of the nozzle with a total of 226 000 nodes.

Figure 11: The medium mesh in the study consisting of 725 000 nodes.

Figure 12: The finest mesh in the study. This mesh has 1 646 750 nodes.

3.4.2 Simulation setup

The simulation of the flow of Ink6040 through the nozzle is modeled to be lami- nar in FLUENT. The operating pressure conditions in FLUENT is kept default, 101 325 Pa. According to data from the bioprinting process with the Ink6040, it was assesed that a pressure of 20 kPa for the microextrusion bioprinter worked the best. Thus the pressure boundary conditions were set to 20 kPa and 0 kPa at the inlet and oulet respectively, relative to the operating conditions. No- slip conditions were set on the walls of the nozzle and symmetry condition on the symmetry plane. The symmetry condition sets the normal velocity at the symmetry plane to zero as well as the normal gradients of any variable at the symmetry plane.

The fluid to be simulated to flow through the nozzle should posses similar properties to that of the Ink6040. Since the ink is 97.5% water, the density of the fluid was set to be constant at 998.2 [kg/m

. The bioink is non-Newtonian, and data regarding viscosity variation with shear rate is available, thus fluid is modeled in Fluent using the non-Newtonian power law for viscosity given by

η = k ˙γ

e

.

Our case is temperature independent, thus T

= 0. For Ink6040 the consis-

tency index k and power-law index n is found to be 109.73 and 0.154 respectively

from curve approximation to experimental data.

The flow is simulated in steady-state conditions with the pressure-based solver since the flow is a low-speed incompressible flow. When using the pressure- based solver a equation for the pressure or pressure correction needs to be de- rived from a combination of the momentum and continuity equations, the nu- merical algorithm for how this pressure or pressure correction is to be derived is refered to in FLUENT as Pressure-velocity coupling scheme. The default setting SIMPLE was used since the flow is steady. In hindsight SIMPLEC (SIMPLE- Consistent) could have been used for faster convergence since the model was simply laminar.

3.4.3 Results

Results are presented for halfway down the nozzle and at the outlet of the nozzle.

Although the continuity and momentum residuals were set to 10

, the differ- ences between the results seemed to increase with increasing grid size (number of nodes), which is opposite to the trend sought for when doing a grid-study.

Due to the inlet pressure boundary condition, the fluid at the inlet is accelerated from more or less a zero velocity (very small). The nozzle geometry accelerates the flow, however outlet velocity of the fluid will still be small (order of 10

).

Therefore the convergence criteria of 10

was not enough and had to be in- creased to 10

. A monitor at the center point of the outlet was implemented to ensure convergence. At both halfway down and at the outlet of the nozzle the velocity profiles are plotted and compared pairwise, i.e. the relative deviation between the coarse mesh and the medium mesh is computed, and the relative deviation between the medium mesh and the finer mesh is computed separately.

Note that the plots for the relative deviation are cut-off at 30% for a more clear view of how the deviation varies.

As can be seen from the results below (Figures 13-14), the relative deviation between the velocity profiles decreases as the number of nodes are increased. The difference between the medium mesh and the finer mesh is almost nonexistent around the majority of the nozzle at the outlet, but spikes up to about 4 − 5%

in the edges where the boundary layer is present. Thus, one can ask if it is worth the computational time to achieve such small difference in results between the medium and coarse mesh (almost double the time). However it may be reasonable to produce a mesh coarser around the center and finer at the edges than the medium mesh.

The computational effort for the 1 646 750 nodes grid is more than double that of the medium grid, and thus the 725 000 nodes grid is more favorable since the majority of the relative error as we stated earlier is around 1%. To fine tune the relative deviation of the velocity profile of the 725 000 nodes grid at the edges, a new mesh with 849 000 nodes was created.

The new grid consisting of 849 000 nodes is identical to the grid with 725 000 nodes in the center region, but has around 35% more nodes near the edges.

Also a geometric law was introduced near the edges to have the node density

-1 -0.5 0 0.5 1

Outlet position [m] _×10^-3

-2.5 -2 -1.5 -1 -0.5 0

Velocity [m/s]

×10^-3Velocity profiles midway for 226 000 and 725 000 nodes 226 000 nodes 725 000 nodes

-1 -0.5 0 0.5 1

Outlet position [m] _×10^-3

0 5 10 15 20 25 30

Relative percentage error

Relative error between 226 000 and 725 000 nodes

Relative percentage error 1% line

-1 -0.5 0 0.5 1

Outlet position [m] _×10^-3

-2.5 -2 -1.5 -1 -0.5 0

Velocity [m/s]

×10^-3Velocity profiles midway for 725 000 and 1 646 750 nodes

725 000 nodes 1 646 750

-1 -0.5 0 0.5 1

Outlet position [m] _×10^-3

0 5 10 15 20 25 30

Relative percentage error

Relative error between 725 000 and 1 646 750 nodes

Relative percentage error 1% line

Figure 13: Comparison between the velocity profiles halfway down the nozzle

for the coarse and medium fine mesh (above) and the medimum fine and the

finest mesh (below). Note that the relative error plot is cut-off.

-2 -1 0 1 2

Outlet position [m] _×10^-4

-0.1 -0.08 -0.06 -0.04 -0.02 0

Velocity [m/s]

Velocity profiles at outlet for 226 000 and 725 000 nodes

226 000 nodes 725 000 nodes

-2 -1 0 1 2

Outlet position [m] _×10^-4

0 5 10 15 20 25 30

Relative percentage error

Relative error between 226 000 and 725 000 nodes

Relative percentage error 1% line

-2 -1 0 1 2

Outlet position [m] _×10^-4

-0.1 -0.08 -0.06 -0.04 -0.02 0

Velocity [m/s]

Velocity profiles at outlet for 725 000 and 1 646 750 nodes

725 000 nodes 1 646 750

-2 -1 0 1 2

Outlet position [m] _×10-4

0 5 10 15 20 25 30

Relative percentage error

Relative error between 725 000 and 1 646 750 nodes

Relative percentage error 1% line

Figure 14: Comparison of velocity profiles at the outlet of the nozzle for the

coarse and the finest mesh (above) and the medium fine and the finest mesh

(below). Note that the relative error plot is cut-off

increase toward the edge. The aim with this grid is to hopefully damp out the spontaneous spikes of relative deviation near the edges when compared to the 1 646 750, without increasing the computational effort too much. A comparison of the 725 000 nodes and the new 847 000 nodes mesh can be seen below.

Figure 15: A view of the original mesh at the outlet of the nozzle with a total

of 725 000 nodes.

Figure 16: A view of the mesh with 847 000 nodes. Note the increase node density at the edges.

Similar to before, results are presented comparing the velocity profiles midway

along the nozzle and at the outlet of the nozzle.

-1 -0.5 0 0.5 1

Outlet position [m] _×10^-3

-2 -1.5 -1 -0.5 0

Velocity [m/s]

×10^-3 Velocity profiles midway

725 000 nodes 847 000 nodes 1 646 750 nodes

-1 -0.5 0 0.5 1

Outlet position [m] _×10^-3

0 5 10 15 20

Relative percentage error

Relative error between velocity profiles

Rel. error between 725 000 and 1 646 750 nodes Rel. error between 847 000 and 1 646 750 nodes 1% line

Figure 17: Velocity profile comparison in the crossection halfway down the

nozzle for the different meshes.

-2 -1 0 1 2

Outlet position [m] _×10^-4

-0.1 -0.08 -0.06 -0.04 -0.02 0

Velocity [m/s]

Velocity profiles at outlet

725 000 nodes 847 000 nodes 1 646 750 nodes

-2 -1 0 1 2

Outlet position [m] _×10^-4

0 5 10 15 20

Relative percentage error

Relative error between velocity profiles

Rel. error between 725 000 and 1 646 750 nodes Rel. error between 847 000 and 1 646 750 nodes 1% line

Figure 18: Comparison between the velocity profiles at the outlet of the nozzle for the different meshes.

Studying Figure 17, we can see that the spikes of 8% near the edges is succesfully damped. However the spikes nearest to the wall have increased to around 5%.

This can be because the new mesh might have even better accuracy near the edges than the 1 646 750 nodes mesh due to the geometric law, meaning that this deviation could be more a measure of how much the finest mesh differs from the 847 000 nodes grid and not the other way around. Anyways a higher deviation near very close to the wall would be to prefer since future analysis will probably involve fluid particles not right at the wall, but further in towards the center.

In figure 18, the relative deviation spikes at the edges are damped to almost

half its value with the original 725 000 nodes mesh when using the new mesh

with 847 000 nodes. The spikes of the newer mesh stretch out to around 2% at

the edge.

4 Post-processing

4.1 Modelling the cells in the flow

In previous sections it has been described how the simulation is setup to model the flow of Ink6040 through the nozzle. However, IPS cells have not yet been incorporated in the flow. The flow of Ink6040 together with IPS cells through the nozzle of the bio-printer can be described as a multiphase flow (a flow in which there is more than one material present with different properties). In this case the flow consists of two phases, the Ink6040 which is the carrier or continuous phase and the stem cells which is the discrete phase dispersed in the continuous Ink6040 phase. A important property when modeling a multiphase flow is the volume fraction of the discrete phase (α

), given by

α

= N V

V , (5)

where N is the number of particles (or cells in this case), V

is the volume of the particle and V is the total volume of the phases combined. It is known that the nanocellulose mixture contains 20 million cells/ml, i.e.

=

= 2 × 10

. According to data received from the bio technology department at Chalmers Institute of Technology, the diameter of the IPS cells ranges from 10 to 15 µm . Assuming the average diameter of a IPS cell to be 12 µm and that the shape of the cells are spherical, the volume of one stem cell is calculated to be V

=

3 p

6

= 9.05 × 10

m

. Inserting the values into (5) the volume fraction is calculated to be

α

= N V

V = 2 × 10

9.05 × 10

= 0.0181. (6) According to FLUENT guidelines, multiphase flows with volume fractions less than 10 − 12% are best modeled with the DPM (discrete phase model) in which the particles are tracked in a lagrangian frame of reference. Note also that since the volume fraction of the discrete phase is around 1 − 2%, impact of the discrete phase on the continuous phase cannot be disregarded without further investigation. Since the cell diameter is of small order compared to the inlet of the nozzle, it’s inertia might be so small that flow simply carries it along the streamlines, in which case it is enough to model the flow as one-way coupled between the continuous and discrete phase. A theoretical value for whether the particle track of a certain particle in a flow is dominated by the inertia of the particle or the streamlines of the flow, is the Stokes number. The Stokes number for flow through acceleration nozzle [12] is given by

St = ρ

C

d

U 9µD

, (7)

where ρ

is the density of the particle, C

is the slip correction factor (assuming

it to be equal to unity in this case [4], d

the diameter of the particle, U the

average velocity of the flow at the inlet, µ the dynamic viscosity of the fluid and finally D

is the diameter of the inlet nozzle. With the assumption of one IPS cell weighing 10

kg (Approximate weight of embryonic stem cells which are very similar to stem cells) and computing its volume for an average IPS cell of 12 × 10

m diameter, its density is calculated to be approximately 1105 kg m

. From previous simulations the volume flow (Q rate through the inlet is found to be approximately 1.5 × 10

m

s

, thus the average velocity at the inlet is 1.53 m s

(U = Q/A

). The majority of the flow at the inlet has a dynamic viscosity of 0.3 Pa s and the inlet diameter D

is as stated in previous sections 5 × 10

m. Substituting these values into (7) Stokes number is calculated to be

St = 1105 · 1 · (12 × 10

)

· 1.53

9 · 0.3 · 0.005 = 0.000018  1. (8) Since the Stokes number in this case is much smaller than unity, it can be assumed that the particle tracks of the cells are more or less governed by the streamlines of the flow.

The IPS cells is best suited to be modeled as inert particles in the DPM model, meaning that these do not evaporate if the particle temperature reaches vaporization temperatures. When the density of the continuous phase is near to the density of the discrete phase it is recommended in the FLUENT guide- lines that both virtual mass force and the pressure gradient force are enabled.

Consider the fluid surrounding the particle, the force required to accelerate this fluid is the virtual mass force and the pressure gradient force is simply the force that arises due to the pressure gradient in the fluid.

The DPM model introduces the discrete phase to the continuous phase by injecting a single particle from a specific coordinate or multiple particles through a surface or in group. The most convenient option in this case is injection through the inlet surface, however in FLUENT this means that one particle will be released from each face on the inlet. Since the mesh is such that it is finer near the edges, the distribution of particles will be non-uniformly increasing radially, not to mention that the cell density will exceed that of the real Ink6040 mixture.

In order to track the particles and see for example if they come in contact with the wall or where they end up at the outlet, a stochastic distribution of the cells were made using MATLAB. A code was written to randomly distribute points within a rectangle. This rectangle had the dimension of half the inlet of the nozzle. The points that lie inside the radius of the inlet minus the radius of a particle (a particle can only be so close to the wall), were then selected.

Since MATLAB rounds down microscale value to zero, all data was scaled up to three orders of magnitude and then scaled down again when exporting the data. Depending on how close the area fraction of the particles at the inlet were of 1.8%, the number of points were either increased or deacreased.

A injection file was then created in MATLAB (see figure where the coordi- nates of the particles were given, their diameter and massflow per particle track.

For convenient tracking and post processing all particle diameters were set to

the average cell diameter of 12 µm. Mass flow per particle stream is a way of

simulating the flow of multiple particles, however all particles with the same

initial position is in the DPM-model assumed to have the same particle track.

Changing the mass flow rate does not change the particle track.

Figure 19: Stochastic distribution of cells in MATLAB.

Since the volume fraction theory and the Stokes number in this case opposes

each other in the sense that according to the former, two-way coupling should

be considered as α

≈ 10

, and the latter suggests that the particles will follow

the fluid streamlines due to the very small Stokes number. To be sure a study

was conducted where the particle distribution mentioned above was injected

every 0.001 s (starting from a steady-state flow) until particles started escaping the nozzle. Then, the particles in the flow field were deleted. If the particles would have had an impact on the flow field, injecting a single particle at a certain point at this stage would show different results compared to a injection of a single particle from the same location before the study. The strain-rate along the particle track of a single particle injected at the same location was investigated before and after the continuous injections. The strain-rate did vary along the particle track, showing that the continuous injections did have a slight impact on the flow field, however the changes were around below 1%.

Therefore it was decided to hereafter consider one-way coupling in steady state flow between the phases.

4.2 Cell deformation model in 2D

Presented below is a model for predicting the cell deformation in shear and extensional flow combined. The model is based on the assumption that shear and strain rates in the fluid produces stresses that directly act on the cell. The 2D cell deformation is then predicted by a viscoelastic creep model based on experimental findings [11].

Due to the converging geometry of the nozzle, both shear rates and strain rates will be present in the flow. For low mass flow rates, strain rate values are not significant, however in the case here studied the mass-flow rate is approxi- mately 1.5 × 10

m

s

with strain rates of significant magnitude. Using the DPM model in FLUENT, the particle track of a modeled cell can be predicted, along which both shear rates and strain rates can be extracted. Thus, it is possible to compute the shear and normal stress variation in the fluid (Ink6040) along the particle track. However, modelling the cell deformation due to both shear- and normal-stress can be quite complicated, therefore using the theory of principal stresses simplifies the problem.

Consider the 2D stress tensor

T =  σ

τ

τ

σ

 ,

where τ

= η( ˙γ

) · ˙γ

according to section 2.1 and ˙γ

=

j

+

i

, which means the stress tensor is symmetric. Due to the symmetric property of the stress tensor it is always diagonizable, in other words, one can find eigenvalues λ

and λ

such that det|τ

− λδ

| = 0. The (orthogonal) eigenvectors e

and e

corresponding to the eigenvalues span up the principal axes. If the coordinate system of the stress tensor is rotated so that it coincides with the principal axes, the stress tensor takes form of

T = ¯ σ

0 0 σ

 .

Evidently, when the coordinate system coincides with the principal axes the

shear stresses disappear. The principal stresses are in this case the eigenvalues

of the tensor, and make up the diagonal of the new tensor.

Assuming that the shear and normal stresses act directly on a fluid element, using the knowledge above, it is possible to rotate the fluid element and the original coordinate system such that only normal stresses are present and acting on the element. Assume also that the shape of the fluid element is spherical to emulate the cell geometry. Cells are known to have viscoelastic characteristics.

This means that the deformation of the cells exponentially reaches a max value when subject to a constant stress after a short period of time. From the article mentioned above [11], the peak value is obtained after 5 seconds, however the simulations done using the DPM model show that a particle released at the inlet flows through the nozzle in around 0.02 seconds. Therefore it can be assumed that the deformation of the cell is linear elastic for this very short period of time.

To model the deformation of the cell in 2D it is then assumed that deformation of the axes due to the stress in the respective other directions are related by Poisson’s ratio ν.

y

x

¯

e

e ¯

σ

σ

y

x

σ

σ

E

E

Figure 20: Cell deformation model. σ

and σ

are the principal stresses acting along eigenvectors ¯ e

and ¯ e

respectively.

According to the article cited in the beginning

11

, the deformation of a cell subject to a constant stress is given by

u(t) = 2σh

3E



1 +  τ

τ

− 1

 e



H(t), (9)

E

= τ

τ

E

. (10)

Here, σ is the constant applied stress, h

initial cell height (cell diameter in this

case), E

and E

instantaneous and relaxed modulus respectiveley, lastly τ

and τ

are the stress and creep relaxation time constants respectively. H(t) is

the step function for the applied constant stress.

y

x

σ

σ

y

x

Figure 21: Cell deformation model. Modelling the axes of the cell as linear- elastic and related to each other with Poisson’s ratio.

As the stresses on the cell act for a very short time, the cell deforms linearly elastic. This means that applying a step function is unnecessary. Incorporating Poisson’s ratio into the above equation for small times, one arrives at:

u

= 2(σ

− νσ

)h

3E

, (11)

u

= 2(σ

− νσ

)h

3E

, (12)

for the deformations in the two different directions. It has been reported

that human red blood cells can survive area strains up to 60% (cell-viability is

reported to be ≈ 0 at this strain) if the strain is impulsive and acts under a

very small time [6]. Therefore 60% is here set as the upper limit for the cell

deformation, and if the area strain is calculated to be beyond this, the cell is

assumed to have died from fatal injuries.

x y

∂v

∂y

< 0

∂v

∂y

> 0

Figure 22: Schematic of the near wall regions inside the noz- zle (gray region). In the main flow (white region), the veloc- ity gradients are positive whilst they are negative in the near wall region.

5 Results & discussion

The no-slip condition at the wall together with the converging geometry of the nozzle gives rise to a region near the wall where the velocity gradient is negative.

This region coincides with what in literature is known as the boundary layer, however since no DNS is performed in this study the boundary layer may not be fully resolved, therefore this region will be referred to as the near-wall region.

Outside this region, the velocity gradients are positive. The investigation of the deformation of the cells will be conducted for two cases, inside the near-wall region and outside the near-wall region (in the main flow).

As stated above two different cases are studied to understand the effect of the flow on the cells. The first case to be considered is the near-wall region. In this case, two sub-cases are investigated: the flow effects on a cell as close to the wall as possible at the outlet (without coming in contact with the wall), and a cell that arrives right at the edge of the near-wall region at the outlet. For the second case the region outside the near-wall region is considered. Here the effects on a cell released such that it arrives at the edge of the near-wall region (from outside the near-wall region) and a cell that at the outlet appears well into the main flow at the outlet.

5.1 Case 1: Near-wall region

The no-slip condition at the wall causes the magnitude of shear-rate values to

increase with decreasing distance to the wall. In order to get an estimate of the

greatest impact this particular flow can have on a cell, it is best to follow a cell

released as close to the wall as possible. The behavior of a cell when in contact

with the wall is complicated and does not produce any relevant information,

therefore the setup is such that if the cell comes in contact with the wall it

becomes trapped. The closest a cell can be injected from the wall at the inlet

without being trapped is found to be 2.5 × d

through trial and error (d

is the

particle diameter, 12 µm). To observe the minimum flow impact on a cell inside

the near-wall region, a cell is inject at the inlet so that it’s particle track ends

right at the edge of the near-wall region. This injection position was found to

be approximately 12.5 × d

from the wall at the inlet. From here on the cell

released at the wall is referred to as cell 1, and the cell ending up at the edge

of the near-wall region as cell 2. Below is a figure with the shear-rate variation along the particle track for a particle released 2.5 × d

and 12.5 × d

from the wall at the inlet respectively.

0 0.01 0.02 0.03 0.04

Particle path length [m]

10⁰ 10² 10⁴ 10⁶ 10⁸

Strainrate[1/s]

∂u/∂x

0 0.01 0.02 0.03 0.04

Particle path length [m]

10^-2 10⁰ 10² 10⁴ 10⁶

Shearrate[1/s]

∂u/∂y

0 0.01 0.02 0.03 0.04

Particle path length [m]

-10⁸ -10⁶ -10⁴ -10²

Shearrate[1/s]

∂v/∂x

0 0.01 0.02 0.03 0.04

Particle path length [m]

-10⁸ -10⁶ -10⁴ -10² -10⁰

Strainrate[1/s]

∂v/∂y

Cell 2 Cell 1 (Closest to wall)

Figure 23: Strain- and shear-rates along the particle track for a particle released 2.5 diameters from the wall at the inlet (near the wall) and for a particle released at approximately 12.5 diameters from the wall at the inlet which ends up at the border inside the near wall-region.

Note that the length of the nozzle is 0.031 m, and as the converging half angle of the nozzle is small (≈ 4

), the total path length travelled by the cells from the inlet to the outlet will be approximately the same (0.031 m), making it easier to compare the results for the different cells. From a first glance at the plots above, one can notice that the shear- and strain-rate variation experienced by the cells is much steadier for cell 1 compared to cell 2. This is reasonable as cell 1 is as close to the wall as possible without touching the wall and away from the border of the near-wall region. However it can be seen that for cell 2, which is at all times closer to- and ends up right at the border of the near-wall region at the outlet, that the variation of the shear- and strain-rates are allot more unsteady.

Since

< 0 along the whole path length for both cells, the particle track

for both cells is at all times inside the near-wall region for both cases. In the

same plot there is a sudden fluctuation of

for cell 2 which is opposite to

the general decreasing trend, this is probably due to the fact that as the cell

approaches the outlet it comes closes to the near-wall region border and thus for

a very brief moment the strain rate experienced by cell 2 is increased. Initially there is a rapid decrease in

for cell 2, which is not as apparent for cell 1. As cell 1 is injected very close to the wall, there is no great relative difference in what the cell experiences in the

field. Cell 2 on the other hand is injected further away from the wall, the streamline (and thus the cell) therefore carries the cell vertically straight down initially just before the converging effect of the wall pushes the streamline away from the wall (this happens very quickly), thus cell 2 is carried closer to the wall which is why there is an initial decrease in

for cell 2. The ”push” from the wall on cell 2 can be observed in the plot for

as there is an initial rapid increase in the strain-rate experienced by the cell.

For the shear-rate plots of

and

, note that the sign convention is due to how the coordinate system is defined for this case. For example, when cell 2 is injected into the flow at the inlet it travels straight down until it experiences the converging effect of the wall, therefore the velocity in the x-direction increases in the streamwise direction initially, however the streamwise direction is in the negative y-direction. Hence, the initial values for the cells in the shear-rate plot

∂u

∂y

are negative for both cells (more so for cell 2) which is why these values do not show in the logarithmic scale. The near encounter with the near-wall region border can also be seen in the plot for

where values are missing around (0.025 m) of the particle path length. Due to the fact that the velocity gradient decreases along the nozzle wall inside the near-wall region the cell will be compressed in the y-direction. Similarly as

> 0 inside the near-wall region, the cell will elongate in the x-direction. Generally it can be observed that the shear- and strain rates are greater in magnitude for cell 1 (closest to the wall) than cell 2 (at the border of the near-wall region).

Below are the corresponding fluid stresses along the particle tracks of the

two cells.

0 0.02 0.04 Particle path length [m]

10² 10³ 10⁴

Normalstressσxx[Pa]

σ_xx

0 0.02 0.04

Particle path length [m]

-10⁵ -10⁴ -10³ -10²

Shearstressτxy[Pa]

τ_xy

0 0.02 0.04

Particle path length [m]

-10⁴ -10³ -10²

Normalstressσyy[Pa]

σ_yy

Cell 2 Cell 1

Figure 24: Normal and shear-stress variation along the particle track for two particles that end up as near the wall as possible and at the border of the near wall region respectively.

To eliminate the shear-stress and simplify the cell deformation prediction,

the principal stresses are calculated according to section 2.5. It is understood

from the velocity gradient variation inside the near-wall region that the maxi-

mum principal stress acts along the principal axis closest to the x-axis and the

minimum principal stress along the principal axis 90

anticlockwise from the

maximum stress direction.

0 0.02 0.04 Particle path length [m]

10² 10³ 10⁴ 10⁵

Stress[Pa]

Principal stress in x^′-direction

0 0.02 0.04

Particle path length [m]

-10⁵ -10⁴ -10³ -10²

Stress[Pa]

Principal stress in y^′-direction

Cell 2 Cell 1

Figure 25: Principal stress variation along the particle track for the two particles (ending up near the wall and at the border of the near wall-region at the outlet respectively).

In this simplified cell-deformation model the principal stresses are assumed to act on a sphere in a undeformed state. Below is a table showing the predicted strains on the axes at different time steps and the resulting area change of cell 1 closest to the wall. Note that if the axis strains is less than -1 or if the total area change is greater than 60% the result is reported as unphysical. The cell therefore is assumed to have ruptured at this stage.

Time [s] Distance covered by cell 1 [m] Total area change [%]

0.0 0.0 -15.5%

0.024 0.013 -26.5%

0.0265 0.018 -36.5%

0.0275 0.021 -48.5%

0.0278 0.0230 -58.6%

0.0278-0.0284 0.023-0.031 <-60%

In the table below are the strain and total area change data along the particle

track for cell 2 at the edge of the near-wall region.

Time [s] Distance covered by cell 2 [m] Total area change [%]

0.0 0.0 -10.8%

0.0092 0.018 -22.3%

0.0099 0.027 -34.1%

0.00996 0.0285 -45.8%

0.00998 0.0301 -59.9%

0.009980-0.009988 0.0301-0.031 <-60%

Although the total area change is below 60% for a greater portion of the particle track time compared to the previous case, the results of the predicted cell deformation goes beyond the assumed limit of 60% during the final instants.

Note that cell 2 comes very close to the outlet, the difference in time between where the cell has deformed −45.8%, to when it reaches the outlet is not resolved in MATLAB.

5.2 Case 2: Outside the Near-wall region

As discussed in the beginning of the section, gradients change outside the near- wall region. More specifically,

> 0 causing the cells to deform in a different manner. Two subcases are also considered here, the first one being just outside the near-wall region in the main flow (cell 3), and the second cell well into the main flow region (cell 4). Outside the near-wall region the cell is assumed to elongate in the y-direction and become compressed in the x-direction due to the velocity gradients in this region.

To find the position at the inlet where the cell assumes ends up in the main

flow just outside the near-wall region, single cell injections were made near

injection position of the last case (12.5 × d

from the wall). If the injected cell

assumes high values of

> 0 at the outlet, it is assumed to be outside the near-

wall region. From trial and error a cell injected approximately 25 × d

from the

wall at the inlet, turned out to fulfill the criteria just mentioned and is assumed

to be right outside the near-wall region at the outlet of the nozzle. Since the

nozzle is converging. Along the particle track values of gradients with opposite

signs of the ones just mentioned can be found, suggesting that the particle is

inside the near-wall region, however once the particle approaches the outlet

the particle is certainly outside the near-wall layer. Due to the transitioning

phenomena of the particle tracks and the change in signs of the gradients, for

this case the results are not plotted in logarithmic scale as it would have excluded

a significant amount of the results.