Lattice Boltzmann Multicomponent Model for Direct-Writing Printing
Michele Monteferrante, Andrea Montessori, Sauro Succi, Dario Pisignano, Marco Lauricella
LLattice Boltzmann Multicomponent Model for Direct-Writing Printing
Michele Monteferrante, a) Andrea Montessori, Sauro Succi,
1, 2, 3
Dario Pisignano,
4, 5 and Marco Lauricella b) Istituto per le Applicazioni del Calcolo CNR, Via dei Taurini 19, 00185 Rome,Italy Center for Life Nano Science @ La Sapienza, Istituto Italiano di Tecnologia, 00161 Roma,Italy Harvard Institute for Applied Computational Science, Cambridge, Massachusetts,United States Dipartimento di Fisica, Università di Pisa, Largo B. Pontecorvo 16 3, 56127 Pisa,Italy NEST, Istituto Nanoscienze-CNR, Piazza S. Silvestro 12, 56127 Pisa, Italy
We introduce a mesoscale approach for the simulation of multicomponent flows to model the direct-writingprinting process, along with the early stage of ink deposition. As an application scenario, alginate solutionsat different concentrations are numerically investigated alongside processing parameters, such as apparentviscosity, extrusion rate, and print head velocity. The present approach offers useful insights on the inkrheological effects upon printed products, susceptible to geometric accuracy and shear stress, by manufacturingprocesses such as the direct-writing printing for complex photonic circuitry, bio-scaffold fabrication, and tissueengineering.
I. INTRODUCTION
In the last decade, 3D printing processes have gained enormous attention as tool for additive manufacturing in manyfields of science and engineering. The major success of 3D printing is mostly due to the digital process control, whichoffers remarkable flexibility in terms of patterning through material deposition. By exploiting a computer-controlledlayer-by-layer fabrication technique, soft materials are utilised in fused deposition modelling and in extrusion direct-writing bio-plotters that apply a pressure gradient to fluids, possibly generating architected matter with qualitativelynew properties . As a consequence, this set of technologies is nowadays exploited in a wide variety of applications,such as tissue engineering (e.g., bio-compatible scaffolds), microsystems (lab-on-chip), microelectronics (sensors), andaerospace structures (aircraft engine bracket), to name a few . The vast potential of additive manufacturing requires,however, an unprecedented control over several aspects of the soft materials involved in the 3D printing process. Theirdynamics, composition, structure, function and rheology are all key elements, which severely affect the features of thefinally produced parts and structures.In this framework, efficient numerical simulations might offer a crucial help to understand the relevant, and in-terplaying, characteristics of the fluids and experimental setups, similarly to other manufacturing processes wherecomputational tools have been successfully applied in the last two decades (e.g. electrospinning , electrospray ). Inthe direct-writing printing context, numerical simulations could be used to maximise the printability of a given ink,avoid process failures and anticipate the microstructural properties of the products, by providing a set of observables(e.g., flow rate, stress tensor), which are often difficult to access experimentally. Printability is usually studied em-pirically, and a priori criteria are not available that, given an ink and a prototype model, allow the success of themanufacturing process to be reliably predicted. Moreover, the challenging characteristic length scale of the process, i.e. the diameter of the print head nozzle used to deposit materials, significantly constrain the possible choices withinavailable computational methods. For instance, microscopic techniques, such as molecular dynamics, are generallyunable to access length and time scales of experimental relevance, for want of computing resources. On the otherhand, finite volume or finite element methods may also become computationally expensive in the presence of movingboundary conditions, such as the ones describing the moving print head.Mesoscale techniques, and particularly the Lattice Boltzmann (LB) method offer an appealing alternative toboth methods above, eventually striking an optimal balance between the two. Indeed LB moves noise-free discreteprobability distributions along force-free (straight) trajectories and represents the effect of molecular collisions throughthe relaxation towards a suitable lattice local-equilibrium. Once the lattice symmetry and the local equilibria aresuitably designed, the scheme can be shown to reproduce quantitatively the Navier-Stokes equations of fluid flows.The result is a very elegant and efficient computational scheme, featuring outstanding amenability to parallel imple-mentations also in the presence of strong nonlinearities and complex boundary conditions . a) First author b) Corresponding author: [email protected] a r X i v : . [ phy s i c s . f l u - dyn ] F e b In this work, we open a new route for predicting 3D printability, developing the regularised version of the ColourGradient (CG) Lattice Boltzmann (LB) model to account for the Non-Newtonian rheological behaviour, typicalof 3D printed pseudo-plastic inks. These systems endure a largely varying apparent viscosity, depending on the shearrate . The regularised version of the LB method mitigates issues related to both low and high viscosities , the firstthreatening numerical stability, while the latter undermining the very hydrodynamic limit of the LB scheme. Further,the shift in the nozzle position during actual 3D printing processes is also included. As a practical application, wefocus on sodium alginate solutions, which are widely used in direct-writing printing to manufacture scaffolds for cellcultures and tissue regeneration . The printing accuracy is discussed in terms of a Parameter Optimization Index(POI) , which is predicted in terms of the numerical inputs. II. MODEL DETAILSA. Regularized colour-gradient lattice Boltzmann model
The regularized CG method for multicomponent-multiphase systems provides computationally efficient access tocapillary number regimes relevant to microfluidics .The general form of CG LB equations writes as follows: f ki ( ~x + ~c i ∆ t, t + ∆ t ) = f ki ( ~x, t ) + Ω ki [ f ki ( ~x, t )] , (1)where k is the colour component, i is the index spans the lattice discrete directions, and Ω ki denotes the collisionoperator of the Colour Gradient model. In Eq. 1, f ki , represents the probability of finding a particle of the k − th component at position ~x and time t with discrete lattice velocity ~c i . For the sake of simplicity, we adopt the standardD2Q9 lattice, and only two colours are assumed in the system, i.e. yellow and blue, standing for the dense ink and air,respectively. In the following, i = 0 denotes the resting population with zero velocity, while i ∈ [1 , . . . ,
8] represent thedirections at angle θ = ( i − π/ | ~c i | = 1 for i ∈ [1 , , , | ~c i | = √ i ∈ [2 , , ,
8] inlattice units, assumed ∆ x LB = 1 and ∆ t LB = 1. A comparison among the regularized version of the Colour Gradient(CG) lattice Boltzmann (LB) model and other LB diffuse interface approaches (e.g., pseudo-potential model, freeenergy model) was discussed by S. Leclaire and coworkers . The density, ρ k , of the k − th fluid component is assessedas the zeroth moment of the distribution functions ρ k ( ~x, t ) = X i f ki ( ~x, t ) , (2)while the total momentum, ρ~u , is defined by the first order moment: ρ~u = X k X i f ki ( ~x, t ) ~c i , (3)being ρ the sum of the two component densities. The collision operator, Ω ki , results from the combination of threesub-operators: Ω ki = (Ω ki ) (3) [(Ω ki ) (1) + (Ω ki ) (2) ] . (4)The first, (Ω ki ) (1) , is the standard BGK operator:(Ω ki ) (1) = − τ [ f ki ( ~x, t ) − f k,eqi ( ~x, t )] , (5)where τ is a relaxation time setting the numerical viscosity of the mixture (see below) and f k,eqi ( ~x ) is a modifiedequilibrium distribution function: f k,eqi ( ~x, ρ k , ~u ) = ρ k (cid:16) φ ki + w i h ~c i · ~uc s + ( ~c i · ~u ) c s − ( ~u ) c s i(cid:17) , (6)with c s the lattice sound speed and w i the weights of the standard D2Q9 lattice : w = 4 / w , , , = 1 / w , , , = 1 /
36. Here, the coefficients, φ ki , read : φ ki = α k , i = 0 , (1 − α k ) / , i = 1 , , , , (1 − α k ) / , i = 2 , , , , (7)and are tuned to simulate systems with different density ratio γ : γ = ρ Y ρ B = 1 − α B − α Y , (8)with the apexes Y and B standing for yellow and blue fluid component, respectively. The partial pressure of k − thcomponent reads: p k = 35 ρ k (1 − φ k ) . (9)The second operator, (Ω ki ) (2) , called perturbation operator, generates the interfacial tension and has the form:(Ω ki ) (2) = A |∇ ρ N | h w i ( ~F cg · ~c i ) | ~F cg | − B i i , (10)where ~F cg denotes the colour gradient force, reading: ~F cg = ρ B ρ ∇ (cid:16) ρ Y ρ (cid:17) − ρ Y ρ ∇ (cid:16) ρ B ρ (cid:17) . (11)As observed in Refs , the gradient for an arbitrary observable χ can be obtained by the second-order isotropiccentral scheme : ∇ χ ( ~x, t ) = 1 c s X i w i χ ( ~x + ~c i , t ) ~c i (12)In Eq. 10, the B i coefficients depend on the lattice (taken: B = − / B , , , = 2 / B , , , = 5 /
108 fromRef. ), whereas A is a free parameter modeling the surface tension, σ , that is: : σ = 2 τ A, (13)where τ is the effective relaxation time. The recoloring operator (Ω ki ) (3) is necessary since the perturbation operatoralone does not guarantee the phase separation:(Ω Yi ) (3) = ρ Y ρ f ∗ i + β ρ Y ρ B ρ cos( ϕ i ) X k f eqi ( ~x, ρ k , ~u = 0) (14)(Ω Bi ) (3) = ρ B ρ f ∗ i − β ρ Y ρ B ρ cos( ϕ i ) X k f eqi ( ~x, ρ k , ~u ) (15)Here, β is a parameter tuning the thickness of the diffuse interface, f ∗ i is the post collision total density in the latticedirection i , f eqi = P k f k,eqi and finally: cos( ϕ i ) = ~c i · ∇ ρ N | ~c i ||∇ ρ N | . (16)The kinematic viscosity, ν , is assessed as : 1 ν = ρ Y ρ ν Y + ρ B ρ ν B , (17)being ν Y and ν B the density of the yellow and blue component, respectively. In order to model the wettability on thedifferent walls in the system, see Fig. 1, and compute the gradients of ρ k by Eq. 12 also close to the boundaries, weset a fictitious fluid density, ρ ks , for the two components on all the wall nodes . The fictitious densities are estimatedby the extrapolation of the color function at neighboring fluid lattice nodes by the formula: ρ ks ( ~x, t ) = ζ k ( ~x, t ) X i w i ρ k ( ~x + ~c i , t ) w i s ( ~x + ~c i , t ) , (18)where ζ k ( ~x, t ) is a positive parameter tuning the affinity of the different walls, see Fig. 1, at the position ~x for a givenfluid component, and s is a switch function taking value one if the site at ~x + ~c i is a fluid and is zero otherwise. Notethe present strategy can be interpreted as a simplified version of the approaches reported in Refs . Although lessaccurate in reproducing the contact angle, our approach is a simple procedure to model hydrophobicity ( ζ k <
1) orhydrophilicity ( ζ k >
1) of the walls as given in Fig. 1.Implying the Einstein convention for summation over Greek indices (see Appendix A of Ref. ), the regularizationstep reads : f k,regi ( ~x, t ) = f k,eqi ( ~x, ρ k , u ) + w i c s Q iαβ Π neq,kαβ , (19)where Q iαβ = ( c iα c iβ − c s δ αβ ) and Π neq,kαβ = ( P i f ki c iα c iβ ) − ( P i f k,eqi c iα c iβ ) with α, β denoting Cartesian directionsand δ the Kronecker delta. Note that Eq. 19 consists of a projection of a distribution functions, f ki , onto the set ofHermite basis. In doing so, we obtain a set of filtered distribution functions, f k,regi , which depends only on the firstand second macroscopic hydrodynamic moments without higher-order non-equilibrium information often referred asghost moments . It was shown that the procedure provides general benefits in terms stability in the BGKLB scheme, which can be decisive in the case of low-viscosity simulations. Hence, the regularized distributions, f k,regi are used in Eq. 1.The hydrodynamic limit of Eq. 1 is found to converge to a set of equations for the conservation of mass andlinear momentum: ∂ρ∂t + ∇ · ρ~u = 0 (20) ∂ρ~u∂t + ∇ · ρ~u~u = −∇ p + ∇ · [ ρν ( ∇ ~u + ∇ ~u T )] ++ ∇ · Σ (21)where p = P k p k is the pressure, ν = c s ( τ − /
2) is the kinematic viscosity of the mixture, Σ = − τ P k P i (cid:0) Ω ki (cid:1) (2) ~c i ~c i is the stress tensor of the curved interface, and τ is a time tuning the relaxation of the fluid flow towards its localequilibrium used in the collision operator, (Ω ki ) (1) , of Eq. 5. At each time step before the collision in Eq. 1, all thepopulations, f ki , are filtered by applying the regularization step .It is worth to highlight that the LB approach avoids two potential and serious issues of computational physics indiscretizing the Navier-Stokes equations of continuum fluid mechanics: strong non-linearity and complex geometrywithin a time-dependent formulation. In particular, the discretization of the non-linear term, ∇ · ρ~u~u , in the Navier-Stokes equations requires the non-locally derivative approximations over the adjacent space in numerical techniquessuch as finite-difference methods and finite element methods. In contrast, the LB method disentangles the non-localityand non-linearity of the problem, since the non-linearity is treated locally (collision step of Eq. 1), and the non-localityis treated linearly (streaming step of Eq. 1) as a shift in memory over the adjacent nodes. Thus, it turns out thatthe LB approach is a very attractive computational bargain to high-performance computing on parallel architectures,including GPUs . B. Extension to non-Newtonian flow and moving print head
To model the non Newtonian fluids, the model is extended in similarity with the approach reported in Refs .The extension consists essentially of determining the local value of the relaxation time, in such a way that the desiredlocal value of the viscosity is recovered to match a constitutive equation for the stress tensor . We assume thatthe shear-thinning model introduced originally by M. Cross (in the following referred to as Cross model) adequatelydescribes the ink viscosity. Note that the Cross model was already employed to describe the Non-Newtonian behaviorof alginate solutions . However, other possible models can be freely adopted without losing the applicability of thepresent approach. The Cross model states: µ ( ˙ γ ) = µ ∞ + µ − µ ∞ λ ˙ γ ) n , (22)where µ is the apparent dynamic viscosity, n the flow index ( n < µ the zero shear viscosity, µ ∞ the asymptotic value, and λ the retardation time at which the shear-thinning starts. In the following, the yellowdense component (the ink) is assumed to show a Non-Newtonian behaviour while the blue component (air) is aNewtonian fluid. Far from the interface, the stress tensor and the strain tensor are mainly represented by the k-thcomponent so that Π αβ ∼ Π kαβ and Γ αβ ∼ Γ kαβ , respectively. Following the literature , the stress tensor Π αβ relateswith the strain tensor Γ αβ by the relation Γ αβ = − (1 / ρτ c s )Π αβ , where the stress tensor Π αβ = P i ( f i − f eqi ) ~c iα ~c iβ .Thus, in the yellow fluid bulk the last relation can be rewritten as:˙ γ Y ( | Π Yαβ | ) = | Π Yαβ | ρ Y τ Y ( ˙ γ Y ) c s , (23)where the stress tensor of the yellow component reads Π Yαβ = P i (cid:16) f Yi − f Y,eqi (cid:17) ~c iα ~c iβ , and the tensor norms arecomputed as ˙ γ Y = 2 | Γ Yαβ | = 2 qP αβ Γ Yαβ Γ Yαβ and | Π Yαβ | = qP Π Yαβ Π Yαβ with the relaxation parameter τ Y settingthe kinematic viscosity of the yellow fluid, ν Y = c s ( τ Y − / µ ( ˙ γ ) = ρν ( ˙ γ ) and µ ( ˙ γ ) = ρc s ( τ ( ˙ γ ) − / τ Y ( ˙ γ y ), requested in Eq. 23, can be obtained by Eq.22 as: τ Y ( ˙ γ Y ) = 12 + 1 c s [ ν ∞ ,Y + ν ,Y − ν ∞ ,Y λ ˙ γ Y ) n ] . (24)Inserting Eq. 24 in Eq. 23 provides an implicit equation in the variable ˙ γ Y , which is solved iteratively, performingseveral iterations as long as the value of ˙ γ Y is not converged below a given threshold. If close to the interface, τ iscomputed from the interpolated value of the viscosity by Eq. 17. A similar approach was exploited by Pontrelli etal. to model a pseudo-plastic single-phase fluid, and it was validated by comparison with the analytical solution ofPoiseuille flow with the power-law model.Since the print head moves during the process, we needed a particular treatment of the boundary conditions ofthe nozzle walls and fluid nodes around the nozzle. Inspired to the trailblazing work by Antony Ladd , we define atemplate of solid nodes with an internal reference system, which is translated along the time evolution by integratingan equation of motion. In accordance with the formulation proposed by F. Jansen and J. Harting , only the exteriorregions are filled with fluid, whereas the interior parts of the nozzle is considered solid nodes.Denoting with f ∗ i ( ~x b , t ) the post collisional distribution at the boundary position ~x b hitting the nozzle wall in thedirection ~c i and located in the middle position along the link connecting the solid node ~x s from the boundary fluidnode ~x b , we exploit a simple generalization of the halfway bounce-back rule . Hence, the streaming step proceedsas: f k ¯ i ( ~x b , t + 1) = f ki ( ~x b , t ) − w i ρ kw ~υ nozzle · ~c i c s (25)where ¯ i is the lattice direction − ~c i . The symbol ρ kw in Eq. 25 denotes the density at the wall position, which isobtained by a third order interpolation in the direction − ~c i .Because the print head moves over the lattice nodes, it happens that a subset of fluid boundary nodes in front ofthe moving nozzle crosses its surface, becoming solid nodes. Similarly, a subset of interior nodes on the surface isreleased at the back of the nozzle. The two distinct events require the destruction and the creation of fluid nodes,respectively. Following the previous strategy reported in the literature, whenever a fluid node changes to solid, thefluid is deleted without transferring its linear momentum to the nozzle beneath the nozzle infinite mass hypothesis.In the creation fluid node step, as the nozzle leaves a lattice site, new fluid populations are initialized from theequilibrium distributions, f eq,ki (¯ ρ k , ~υ nozzle ), for the two k -th components with the velocity of the nozzle wall, ~υ nozzle and the k − th fluid density taken as its average value, ¯ ρ k , computed over the neighbouring fluid nodes .As a first approximation, the solid-fluid interactions are accounted only for the part concerning the effect of themoving nozzle on the surrounding fluid and not viceversa, which is equivalent to assume that the motion of the printhead is fully controlled by the digital system of the 3D printer (nozzle infinite mass hypothesis). Hence, a constantvelocity ~υ nozzle of the nozzle (print head) is taken as an input parameter to describe the linear motion of the nozzle,and a drag force is added close to the deposition zone in order to model the friction between the ink and the collector.Inside the nozzle, the yellow component is inserted with constant velocity ~u ink with respect to the internal referenceframe on the print head. Considering that the nozzle reference frame is moving with ~υ nozzle , the total fluid velocityinserted at the inlet reads: ~u Inlet = ~u ink + ~υ nozzle . (26)Hence, ~u Inlet replaces ~υ nozzle in Eq. 25, which is used to model the fluid inlet inside the print head.We also used a Dirichlet boundary conditions in our work to maintain the pressure (density) of the blue fluid (air)constant. For the Dirichlet condition, the anti-bounce back scheme is used instead for constant pressure (densities)boundaries: f k ¯ i ( ~x b , t + 1) = − f k ∗ i ( ~x b , t ) ++2 ρ kw n φ ki + w i h ( ~c i · ~u w ) c s − ~u w c s io , (27)where ρ kw is the imposed density at the open boundary and ~u w is the velocity at half-way point estimated by a secondorder interpolation along the direction − ~c i .The drag force modelling the friction between the ink and the collector reads: F d ( ~x, t ) = − γρ Y ( ~x, t ) ~u ( ~x, t ) , (28)where γ is the friction coefficient tuning the drag force. This force is turned on at grid points which are closer than 4lattice units from the deposition wall. The drag force is added in Eq. 5 by the exact difference method proposed byKupershtokh and coworkers . III. SYSTEM SETUP
It is worth to remark the main simplifying assumptions adopted in the simulations. In the present paper, the Crossmodel is adopted to describe the non-Newtonian fluid, albeit any other rheological relation can be adopted, obviouslyin the context of pseudo-plastic models (e.g., Carreau Model, Herschel-Bulkley model, etc.). Further, the adhesionproperty of the fluid onto the deposition surface can be relevant. The contact angle is set equal to 90 ◦ , correspondingto the neutral affinity of the ink to the surface (neither hydrophilic nor hydrophobic). Nonetheless, different contactangle values can be investigated by the present model. Finally, we exploit a two-dimensional description of the systemto reduce both the computational cost and the wall-clock time necessary for each simulation. As a consequence,the lateral shear effect of the slender fluid filament along the third dimension is neglected in the two-dimensionaldescription. Nonetheless, the comparison with experimental data in Section IV will show as the two-dimensionalapproximation does not invalidate the agreement between the numerical results and the experimental counterpart. Inother words, the validation of the present model is not prevented from the dimension reduction.The alginate concentration in water is taken in the range 1% −
3% w/v, allowing the investigation of the fluidcharacteristics for optimal printing. Alginate solutions are shear-thinning non-Newtonian fluids . Thus, theviscosity decreases to smaller value as the shear rate increases. This dependence of the viscosity on the shear ratemakes the whole printing process complex, since on the one side, low viscosity reduces shear forces, thus speeding upprinting, but it can also reduce both resolution and accuracy . The apparent viscosity (0.2-3.0 Pa · s) observed inalginate solutions represents a good compromise between the above criteria . The parameters of the non Newtonianfluid are taken by . Sarker and Chen have also investigated the rheology of the alginate solutions, although someparameters as λ and n were not explicitly given. The rheological parameters of the yellow fluid modelling the inkare reported in Tab. I. The kinematic viscosity of the blue fluid modelling the air is set ν B = 1 . · − m /s,corresponding to the kinematic viscosity at 25 ◦ C.To model the fluid-air system we simulate two fluids with a density ratio γ = 842 . ≈ the water/air density ratioat 25 ◦ C), while the surface tension is set at the typical value, σ = 50 . · − N/m . The simulation box consist of240 ×
880 lattice nodes, the nozzle diameter of the channels was fixed at d = 60 lattice nodes, and the same value wasassigned to the distance of the nozzle by the deposition surface corresponding to 60 lattice nodes. This configurationis characteristic of 3D bio-plotters . The system is initialised with the nozzle filled up with ink and located on theright side of Fig. 1, that also shows graphically the various boundary conditions used in the simulations. These wereset periodic along the x − axis, while along the y − axis, the bottom side is a no-slip wall and the top boundary outsidethe nozzle is set to a constant blue (air) density ρ B .The ink velocity inside the nozzle in the internal reference frame is obtained from the mass flow rate given inliterature , by noting that: v ink = ψρπ ( d/ (29)As reported in , for a nozzle diameter of 0 . ψ are between 7.7 mg/s and 27 mg/s.Setting ψ =14 mg/s and the ink fluid velocity v ink by Eq. 29, the nozzle velocity, v nozzle , was varied in range from0 . v ink to 1 . v ink . Hence, υ inlet is computed as: ~u Inlet = ~v ink + ~υ nozzle . ν (m s) ν ∞ (m s) λ (s) n ν air ( m s)1%-25 ◦ C 15 . − .
11 10 − .
16 10 − .
552 10 − ◦ C 66 . − .
69 10 − .
96 10 − .
552 10 − ◦ C 463 . − .
64 10 − . − .
552 10 − ◦ C 186 . − .
14 10 − . − .
552 10 − ν (LB) ν ∞ (LB) λ (LB) n ν air (LB)1%-25 ◦ C 0.134 0.0181 3 .
310 10 − . ◦ C 0.430 0.0045 8 .
325 10 − . ◦ C 1.991 0.0156 131 . − . ◦ C 0.802 0.0221 21 .
37 10 − . Figure 1. Representation of the system alongside with the different treatment of boundary conditions.
Denoting by the subscripts LB and MKS the physical observable in lattice and MKS system of units respectively, weadopted the following rules to convert the lattice units ∆ x LB , ∆ tx LB , ∆ m LB in the corresponding physical quantities∆ x MKS , ∆ tx MKS , ∆ m MKS . Assuming ∆ x LB , ∆ tx LB , ∆ m LB equal to one, the lattice conversion rules are:∆ x MKS = d MKS d LB = 3 . ¯3 · − m (30)∆ m MKS = (∆ x MKS ) ρ air MKS ρ B LB = 4 . · − kg (31)∆ t MKS = (∆ x MKS ) ν B LB ν B MKS = [4 . − . · − s (32)In particular, if we assume the nozzle diameter d MKS = 0 . · − m from Ref. corresponding to d LB = 60 latticenodes, the ∆ x MKS is determined by Eq. 30, while ∆ m MKS is obtained by fixing the air density in lattice units ρ B LB equal to one. Since the kinematic viscosity of the air is always lower than the corresponding value in the ink, ∆ t MKS is obtained by fixing τ B and, thus, ν B LB is also determined by the relation ν BLB = c s ( τ B − . τ B LB and ν B LB were fixed depending on the cases under investigations (see Tab. I ), so that ∆ t MKS spans over the range reported in Eq. 32. However, τ B ∈ [0 . , . . . , .
54] is taken sufficiently far from the limitingvalue 0 . . IV. RESULTS AND DISCUSSION
In 3D printing, the ultimate printability of a given prototype depends both on the printer device and on the physicalproperties of the ink fluid. In order to assess the quality of the print process with respect to tunable parameters, weintroduce a POI following : P OI = accuracy/theoretical shear stress . In fact, it was found that the shear stresscan be minimised by manipulating printing parameters , since it is proportional to the inlet pressure p and inverselyproportional to the nozzle diameter d . Hence, assuming the accuracy to scale inversely with the thickness (height h )of a single printed thread, the POI is written as : P OI ∝ dh p (33)We include the coverage ratio (that is, the ratio of the sectional area of the ink thread, A i , to the area of the rectanglecircumscribing the thread, A r , as in the middle panel of Fig. 2a) in the definition of the accuracy and obtain: P OI ∝ A i A r dh ∆ p m (34)where: we consider the shear stress as proportional to the variation of the pressure in the ink, that reaches a maximumpressure variation ∆ p m (variation with respect to the equilibrium pressure p ) outside the nozzle at completion of thedeposition process; h is taken as the largest height value of the thread behind the moving nozzle.Hence, all the POI values are normalised to a reference value, in order compare the results with the differentparameters set. We take the reference POI as the largest value, max i P OI i , corresponding to a perfect coverage ratio,( A i /A r ) = 1, across all the simulations.As a result, the normalized POI reads: P OI normi = P OI i max i { P OI i } . Nonetheless, it is worth highlighting that the POIis here aimed to determine the process quality in the context of alginate-type inks used for manufacturing applicationsin the context, among others, of bio-scaffolds. Hence, the POI index involves the maximum pressure observed in thesimulation to monitor the shear forces in the fluid. For other applications, such as manufacturing processes withpolymeric inks, the shear stress can be relatively less important. Thus, other indexes of printing quality could bemainly focused on the geometrical precision rather than the shear forces in the fluid, for instance, in the contexts ofnano-printing or electrode fabrication .For all the simulations, we stopped the run as the nozzle covers six time the nozzle diameter d from the latticeposition where the deposited ink first touches the collector, thus allowing the geometry of the printed thread to becompletely developed in high-resolution printing .In top panel of Fig. 2a, a set of snapshots are reported, showing the ink mass densities map ( ρ Y ) at the end of thesimulation for three different velocity conditions of the nozzle υ nozzle = { . , . , . } υ ink and alginate solution withthe lowest concentration 1%. The shape of the deposited ink is found to be strongly dependent on υ nozzle , overflowingbeyond the travelled length of six diameters for υ nozzle = 0 . υ ink such to provide a poor printing quality. That a lowdispensing speed compared to υ ink , providing a surplus of ink compared to the space spanned by the nozzle, woulddecrease the printing accuracy is in agreement with previous results, both experimental and numerical .The height of the printed thread in Fig. 3a is also clearly dependent on the parameter υ nozzle . In order to studythe dynamic evolution, we probe the maximum height, h , of the fluid thread behind the moving nozzle, investigatingwhether a stationary condition is reached along the deposition process.The thread height, h , as a function of the distance travelled by the moving nozzle is reported in Fig. 3. In boththe case of 1% and 3% concentrations, the height of the thread for lowest nozzle velocities, υ nozzle = 0 . υ ink doesnot show any asymptotic trend, confirming that the printing precision is deteriorated by a low dispensing velocity.We also observe that the asymptotic values in h decrease as the nozzle velocity increases for the 1% concentration,highlighting a clear, monotonic trend, which is in agreement with the experimental observations reported by Webband Doyle . It is worth to highlight as the numerical results trace qualitatively the experimental trend, although themodel has been implemented in the two-dimensional framework, endorsing the validity of the dimension reduction.The height evolution at high velocity and concentration 3% w/v shows a second jump, namely the thickness obtainedat higher values of υ nozzle overcomes the value measured with lower velocities υ nozzle = [0 . , . , . υ ink , providinga non monotonic trend. This suggests that, at least at high ink concentration, an optimum operating value exists forthe dispensing velocity, compared to the ink delivery rate, which minimises the thread height. Further, the sequence ofasymptotes is found to be monotonic also in the case of 2% solutions at 25 ◦ C, while the sequence with concentration3% w/v and 40 ◦ C shows the same non-monotonic trend already observed at 25 ◦ C. The non monotonic hight trendsobserved for 3% alginate concentrations (both at 25 ◦ C and 40 ◦ C) is produced by the presence of irregularities inthe thread shapes as the one represented in Fig. 4. These irregularities manifest for large nozzle velocites and 3%alginate concetrations are explained by viscous effect (see the discussion below) and determine also the behaviourof the coverage ratio. The coverage ratio, reported for the four cases in Figure 2(b), allows a similar classification,showing a maximum around υ nozzle ≈ . υ ink for solutions with 3% concentration. Then, the coverage ratio decreasesat higher nozzle velocity values due to the irregular shape of the deposited ink as reported in the ink density map ofFig. 2 (panel a). The distinct irregular signature in the thread decreases the coverage ration at υ nozzle = 1 . υ ink . Figure 2. (a) Yellow fluid densities for 1% alginate solution. From top to bottom the nozzle velocities υ nozzle are equal to1 . υ ink , 1 . υ ink , and 0 . υ ink . The snapshots are taken as the nozzle has travelled six times, 6 d , the nozzle diameter fromthe lattice position where the ink touches the deposition substrate. In the middle panel, we also sketch thread height ( h ), thearea of the ink thread ( A i ), and the area of the rectangle circumscribing the thread ( A r ) defining the coverage ratio, A i /A r .(b)Coverage ratio obtained by different printing parameters. In Fig. 4, the fluid viscosity of the mixture is reported for the case 3% w/v at 25 ◦ C after the ink is deposited. As afirst, we observe in Fig. 4 overall a low viscosity in the extruded fluid part which is the result the shear rate enforcedamong the moving nozzle and the substrate. In all the simulations, we observe a detachment point of the ink from thesubstrate. In particular, the shearing force produces the detachment point just after the ink reaches the substrate.Later, the detachment point remains visible as an irregular blob in the thread (see Fig. 4). Then, the tread reabsorbsthe blob under the action of the capillary pressure. Hence, the rheological behaviour of non-Newtonian inks play acentral role in this process. In particular, the relation between shear-rate and the viscosity tunes the magnitude ofthe transmitted nozzle movements to the deposited ink, biasing both the thread shape and the quality of the final0
Figure 3. Thread height as a function of the nozzle displacement, for a 1% (a) and a 3% w/v alginate solution. Temperature= 25 ◦ C products. Further, it is observed the presence of a low viscosity close to the wall of the nozzle(see Fig. 4), whichis consistent with the Poiseuille flow as a consequence of the larger velocity gradient, ∂υ/∂x , close to the no-slipboundaries. On the other hand, the viscosity profile shows a high peak in the middle of the nozzle (corresponding tothe lowest velocity gradient point), which can be relevant for the cell viability in bio-inks. Indeed, in the context ofcell culture applications, shear stress is essential to control cell viability, which may be compromised by the impactforce generated by high gradient in viscosity within the nozzle channel .The P OI normi values, for thread heights less than two times the nozzle diameter, are assessed and shown in Fig 5.Increasing the alginate concentration results in higher
P OI normi values, which is mainly due to the higher coverageratios A i /A r alongside with smaller variation in the ink pressure. In particular, for 3% w/v concentrations, the ratio∆ p m /p is found in the range from 0.1 to 0.2 for slow nozzle velocity υ nozzle < . υ ink . On the contrary, ∆ p m / ∆ p isalways larger than 0.4 for 1% and 2% w/v and for all the nozzle velocities, providing lower P OI normi values. Thus,the
P OI normi shows a peak at moderate nozzle speed in the range 0.6-0.8 υ ink at high alginate concentration 3 %w/v, due to the simultaneous concourse of high coverage ratios and small ink pressure variations.Finally, in some cases, small porosity was noted in the ink fluid (see the upper panel of Fig. 2). In order to addressthe question, the POI values were reconsidered, taking into account the porosity. Indeed, since A i represents the area1 Figure 4. Kinematic viscosity map of the mixture, ν , for the case 3% w/v at 25 ◦ C, in LB units, after the ink is deposited.The dashed white line highlights the fluid interface. The ink on the right side features higher viscosity than in the contact linewith the substrate.Figure 5. Parameter optimization index, reported for thread heights less than two times the nozzle diameter, for differentalginate concentrations and temperatures. covered by the ink, the area decreases as the porosity increases in the fluid, whenever the porous are excluded in the A i assessment. The POI results are practically unaffected by this new definition, with A i values always changing lessthan 1%. Consequently, the porosity in the trend-line does not bias the features of printed material in the presentsimulations, and its effect can be reasonably neglected. V. CONCLUSIONS
Summarising, we have introduced a multi-component model of non-Newtonian inks through a regularised versionof the colour gradient LB model and used it to simulate the 2D printing process as a function of a number of design2parameters. The model allows to calculate the shear stresses during the printing process of non Newtonian inks,directly accessible by simulations, that is very important to control the cell viability in bio-inks. The print accuracywas quantitatively analysed using the same indexes used in experimental studies . The impact of the pseudo-plasticrheology on the printing accuracy was investigated for a set of solutions at different alginate concentration. Systematicinvestigations of processes are enabled on a broad viscosity range, providing a useful tool to probe the dynamics ofthe forces acting inside and on the ink during additive manufacturing. In real systems, shear thinning fluids areusually employed in order to favour the throughput of the device (small viscosity at high shear rates) and to obtain astable and regular thread at the end of deposition (high viscosity at small shear rates). However the accuracy of thedeposited threads is deteriorated for high viscosity alginate concentration at high print speed since these compositesfavour the transmission of the inertia of the fluid impacting the substrate which produce irregularities in the threads. DATA AVAILABILITY
Data available on request from the authors.
ACKNOWLEDGMENTS
The research leading to these results has received funding from MIUR under the project “3D-Phys” (PRIN2017PHRM8X), and from the European Research Council under the European Union’s Horizon 2020 FrameworkProgramme (No. FP/2014- 2020)/ERC Grant Agreement No. 739964 ("COPMAT").
REFERENCES R. L. Truby and J. A. Lewis, Nature , 371 (2016). T. D. Ngo, A. Kashani, G. Imbalzano, K. T. Nguyen, and D. Hui, Composites Part B: Engineering , 172 (2018). R. D. Farahani, M. Dubé, and D. Therriault, Advanced Materials , 5794 (2016). E. Axpe and M. L. Oyen, International journal of molecular sciences , 1976 (2016). B. Utela, D. Storti, R. Anderson, and M. Ganter, Journal of Manufacturing Processes , 96 (2008). M. Lauricella, S. Succi, E. Zussman, D. Pisignano, and A. L. Yarin, Reviews of Modern Physics , 035004 (2020). A. M. Gañán-Calvo, J. M. López-Herrera, M. A. Herrada, A. Ramos, and J. M. Montanero, Journal of Aerosol Science , 32 (2018). S. Succi,
The lattice Boltzmann equation: for complex states of flowing matter (Oxford University Press, 2018). T. Krüger, H. Kusumaatmaja, A. Kuzmin, O. Shardt, G. Silva, and E. M. Viggen, Springer International Publishing , 978 (2017). H. Huang, M. Sukop, and X. Lu,
Multiphase lattice Boltzmann methods: Theory and application (John Wiley & Sons, 2015). R. Benzi, S. Succi, and M. Vergassola, Physics Reports , 145 (1992). M. M. Cross, Journal of colloid science , 417 (1965). J. Latt and B. Chopard, Mathematics and Computers in Simulation , 165 (2006). B. Webb and B. J. Doyle, Bioprinting , 8 (2017). A. Montessori, M. Lauricella, M. La Rocca, S. Succi, E. Stolovicki, R. Ziblat, and D. Weitz, Computers & Fluids , 33 (2018). S. Leclaire, A. Parmigiani, O. Malaspinas, B. Chopard, and J. Latt, Physical Review E , 033306 (2017). S. Leclaire, M. Reggio, and J.-Y. Trépanier, Journal of Computational Physics , 318 (2013). S. Leclaire, A. Parmigiani, B. Chopard, and J. Latt, International Journal of Modern Physics C , 1750085 (2017). Z. Wen, Q. Li, Y. Yu, and K. H. Luo, Physical Review E , 023301 (2019). S. Saito, Y. Abe, and K. Koyama, Physical Review E , 013317 (2017). H. Liu, A. J. Valocchi, and Q. Kang, Physical Review E , 046309 (2012). T. Reis and T. Phillips, Journal of Physics A: Mathematical and Theoretical , 4033 (2007). M. Latva-Kokko and D. H. Rothman, Physical Review E , 046701 (2005). T. Akai, B. Bijeljic, and M. J. Blunt, Advances in water resources , 56 (2018). S. Leclaire, K. Abahri, R. Belarbi, and R. Bennacer, International Journal for Numerical Methods in Fluids , 451 (2016). R. Zhang, X. Shan, and H. Chen, Physical Review E , 046703 (2006). C. G. Coreixas,
High-order extension of the recursive regularized lattice Boltzmann method , Ph.D. thesis (2018). A. Montessori, G. Falcucci, P. Prestininzi, M. La Rocca, and S. Succi, Physical Review E , 053317 (2014). J. Latt,
Hydrodynamic limit of lattice Boltzmann equations , Ph.D. thesis, University of Geneva (2007). M. Lauricella, S. Melchionna, A. Montessori, D. Pisignano, G. Pontrelli, and S. Succi, Physical Review E , 033308 (2018). G. Pontrelli, S. Ubertini, and S. Succi, Journal of Statistical Mechanics: Theory and Experiment , P06005 (2009). S. Gabbanelli, G. Drazer, and J. Koplik, Physical review E , 046312 (2005). E. Aharonov and D. H. Rothman, Geophysical Research Letters , 679 (1993). O. Malaspinas, G. Courbebaisse, and M. Deville, International Journal of Modern Physics C , 1939 (2007). R. Ouared and B. Chopard, Journal of statistical physics , 209 (2005). B. S. Roopa and S. Bhattacharya, International journal of food science & technology , 2583 (2009). A. J. Ladd, Journal of fluid mechanics , 285 (1994). F. Jansen and J. Harting, Physical Review E , 046707 (2011). A. Ladd and R. Verberg, Journal of statistical physics , 1191 (2001). C. K. Aidun, Y. Lu, and E.-J. Ding, Journal of Fluid Mechanics , 287 (1998). I. Ginzburg, Advances in Water Resources , 1196 (2005). A. Kupershtokh, D. Medvedev, and D. Karpov, Computers & Mathematics with Applications , 965 (2009). M. Sarker and X. Chen, Journal of Manufacturing Science and Engineering (2017). P. Del Gaudio, P. Colombo, G. Colombo, P. Russo, and F. Sonvico, International journal of pharmaceutics , 1 (2005). Y. He, F. Yang, H. M. Zhao, Q. Gao, B. Xia, and J. Z. Fub, Scientific Reports , 29977 (2016). J. Ventrici de Souza, Y. Liu, S. Wang, P. Dörig, T. L. Kuhl, J. Frommer, and G.-y. Liu, The Journal of Physical Chemistry B , 956(2018). D. Ye, Y. Ding, Y. Duan, J. Su, Z. Yin, and Y. A. Huang, Small , 1703521 (2018). M. Wei, F. Zhang, W. Wang, P. Alexandridis, C. Zhou, and G. Wu, Journal of Power Sources , 134 (2017). Y. Jin, W. Chai, and Y. Huang, Materials Science and Engineering: C , 313 (2017). Z. Zhang, Y. Jin, J. Yin, C. Xu, R. Xiong, K. Christensen, B. R. Ringeisen, D. B. Chrisey, and Y. Huang, Applied Physics Reviews ,041304 (2018). J.-F. Agassant, F. Pigeonneau, L. Sardo, and M. Vincent, Additive Manufacturing , 100794 (2019). J. Shi, B. Wu, S. Li, J. Song, B. Song, and W. F. Lu, Biomedical Physics & Engineering Express , 045028 (2018). F. Lee, C. Iliescu, F. Yu, and H. Yu, in