Numerical Approaches to Complex Fluids

We are surrounded by a variety of fluids in our everyday life. Besides water and air, it is common to deal with fluids with a peculiar behaviour such as gel, mayonnaise, ketchup and toothpaste. While simple fluids made by identical molecules show a linear response to the applied forces, complex fluids with a microstructure, such as suspensions, may show a very complex response. In this chapter, we introduce numerical approaches for complex fluids focusing on the way the additional stress, due to the presence of a microstructure, is modelled and how rigid and deformable intrusions can be simulated.


Newtonian and Non-Newtonian Rheology
The macroscopic rheological behaviour of a viscous fluid is well characterised in a Couette flow, i.e. the flow between two parallel walls of area A and at distance b: the upper wall moving at constant (low) velocity U 0 and the lower at rest. To keep the upper wall moving at constant velocity we need to apply a force F which is proportional to the wall area: F ∝ A. Therefore it is more general to consider the stress τ = F /A instead of the force F itself. In a Newtonian fluid the shear stress is proportional to the velocity of the upper wall and to the inverse of the wall distance b, i.e. τ ∝ U 0 /b. This linear response defines Newtonian fluids, such as air, water, oil and many others. Note that, in a simple Couette flow the ratio U 0 /b equals the wall-normal derivative of the velocity profile and the shear (deformation) rate: du/dy =γ = U 0 /b. Thus, for a Newtonian fluid we can express the law relating the applied force with the response, i.e. the shear stress τ with the shear rateγ , as τ = μγ , (1.1) where the proportionality coefficient μ is called dynamic viscosity with dimension P a s in the SI. Many Newtonian fluids exist, each with a different value of viscosity, and therefore flowing at different velocity when subject to the same stress. The viscosity coefficient of a Newtonian fluid does not depend on the shear rate, but may vary with the temperature. Indeed, the viscosity usually increases with temperature in gases, while it decreases in liquids. This behaviour is related to the effect of the temperature on the molecular structure of the fluid, but this is outside the scope of present chapter and the reader is refereed to specialised textbooks.
Fluids that exhibit a non-linear behaviour between the shear stress τ and the shear rateγ are called non-Newtonian and fluids whose response does not depend explicitly on time but only on the present shear rate are denoted generalised Newtonian fluids. In particular, when the shear stress increases more than linearly with the shear rate, the fluid is called dilatant or shear-thickening, whereas in the case of opposite behaviour, i.e. when the shear stress increases less than linearly with the shear rate, the fluid is called pseudoplastic or shear-thinning. Examples of typical profiles of the shear stress τ as a function of the shear rateγ for Newtonian, shear-thickening and shear-thinning fluids are shown in the left panel of Fig. 1.1. The ratio of the applied stress and the resulting deformation rate is the so-called apparent effective viscosity μ e = τ/γ : it increases withγ for shear-thickening fluids, while it reduces for shear-thinning ones, which means that the fluidity of shear-thickening fluids reduces increasing the shear rate, while the opposite is true for shear-thinning fluids. Examples of shear-thinning fluids are ketchup, mayonnaise and toothpaste, while corn-starch water mixtures and dense non-colloidal suspensions usually exhibit a shear-thickening behaviour. Note that, sometimes, the same fluids can have plastic or elastic responses depending on the flow configuration.
Complex fluids may behave as solids, with a finite deformation, when the applied stress is below a certain threshold τ 0 , while for stresses above it, they start flowing as liquids. These fluids are called yield stress or Bingham fluids: when the applied stress exceeds the so-called yield stress, τ 0 , these fluids can exhibit a linear relation between stress and deformation similar to Newtonian fluids or a pseudoplastic response. These macroscopic behaviours are related to changes of the microscopic structure of the fluid, and indeed these fluids are constituted by a Newtonian fluid with one or more suspended phases, such as fibres, polymers, trapped fluids (emulsions). From a qualitative point of view, the material hardly flows and deforms when the connections and interactions between the phases constituting the microstructure are intense. Changing the level of the stress τ applied on these complex fluids may either strengthen, weaken or break these interactions, thus altering their microstructure, and eventually reflecting in their nonlinear rheological behaviour.
In order to describe complex fluids, we need a relation as in Eq. (1.1) between the applied stress, τ , and the deformation rate,γ . A relation that can be used to summarise the behaviours previously described for complex fluids is the Herschel-Bulkley formula τ = τ 0 + Kγ n , (1.2) where τ 0 is the yield stress, n the flow index and K the fluid consistency index. A Newtonian behaviour is recovered when τ 0 = 0, n = 1 and K = μ, while values of the flow index above and below unity, n > 1 and n < 1, denote shear-thickening and shear-thinning fluids, respectively. Finally, yield-stress fluids are characterised by a finite non-zero value of the yield stress τ 0 . The consistency index K measures how strong the fluid responds to the imposed deformation rate. However, the consistency index has the same dimension of a dynamic viscosity only when n = 1, and in general its dimension is a function of n, so that it is not possible to compare different values of K for fluids with different flow indexes n. The fluid discussed so far are inelastic, since the stress is just a function of the present value of the deformation rate, i.e. τ = τ (γ ), and not on the previous history of the deformation rate (no memory effects). Another important class of non-Newtonian fluids, which cannot be described by the Herschel-Bulkley formula, is that of viscoelastic fluids. These materials have property similar to both a viscous liquid and an elastic solid. Indeed, the deformation is not anymore permanent, as in usual fluids, and depends on both viscous and elastic contributions. When a constant stress τ is applied the deformation of a viscoelastic fluid increases with time, but when the applied stress is removed, the fluid tends to recover its original configuration (similarly to elastic solids). Polymer solutions usually experience a viscoelastic behaviour, another culinary example being pizza dough: when softly pressed it deforms, but when the pressure is removed the original shape is recovered. However, if the dough is strongly deformed, we can rearrange it in a new stable configuration similarly to what happens in fluids. Memory and elastic effect are difficult to model, and typically require information about the microstructure deformation.
In some applications, complex fluids can be successfully modelled just by considering that their response is related to the memory of the deformation rate history; in other words, they have a time-dependent viscosity if exposed to a constant value of the shear rate. Two main kind of such fluids can be identified: thixotropic fluids whose effective viscosity decreases with the accumulated strain and rheopectic fluids, whose effective viscosity increases with the accumulated strain. A classic example of a fluid characterised by a thixotropic behaviour is painting whose apparent viscosity increases when the deformation rate reduces in order to better adhere to a surface. Rheopectic fluids are less common, and an example is the synovial fluid in our knees, whose property facilitates the absorption of shocks. Thixotropic and rheopectic fluids are usually modelled by a timedependent viscosity, function of a scalar parameter that represents the evolution of their microstructure.

Inelastic Shear-Thinning/Thickening Fluids
Shear-thinning and shear-thickening are possibly the simplest non-Newtonian behaviours of fluids, when the viscosity μ decreases and increases under shear, i.e. μ = μ(γ ); these behaviours are only rarely observed in pure materials, but can often occur in suspensions. Despite its simplicity, this behaviour is able to capture the main effects induced by a microstructure in many applications. Several models have been developed to describe these fluids, an example of shear-thinning model being the Carreau law, usually used to describe generalised fluid where viscosity depends upon shear rate. The model is able to properly describe pseudoplastic fluid viscosity for many engineering application [3], and assumes an isotropic viscosity proportional to some power of the shear rate [4]: In the previous relation, μ is the viscosity, μ 0 and μ ∞ the zero and infinite shear rate viscosities, λ the relaxation time and n < 1 the power index; the second invariant of the strain-rate tensorγ can be determined asγ = 2 S ij : S ij , where S ij = which reproduces a monotonic increase of the viscosity with the local shear rate for n > 1. The constant M is called the consistency index and indicates the slope of the viscosity profile. More details on the Carreau and power-law models can be found in Ref. [4]. From a numerical point of view, implementation of a shear dependent viscosity is often straightforward; however, high variations of viscosity may result in significantly time step constraint when explicit schemes are used, and disrupt the solution technique usually used to solve implicitly the viscous terms in the momentum equation. Indeed, the diffusive term cannot be reduced to a constant coefficient Laplace operator since the viscosity is now a function of space. Dodd and Ferrante [5] have introduced a splitting operator technique able to overcome this drawback, initially derived for the pressure Poisson equation; however, this splitting approach can easily be extended to the Helmholtz equation resulting from an implicit (or semi-implicit) integration of the diffusive terms as well. In particular, the viscosity is split in a constant part and in a space-varying component, i.e. μ (x) = μ 0 +μ (x), and the resulting diffusive term split consequently in a constant coefficients operator that can be treated implicitly, and in a variable coefficients operator which can be treated explicitly.

Viscoelastic Fluids
Viscoelasticity is the property of materials that exhibit both viscous and elastic characteristics when undergoing deformation. Unlike purely elastic substances, a viscoelastic substance has an elastic component and a viscous component, and the latter gives the substance a strain rate dependence on time. and a purely elastic spring connected in series, the Kelvin-Voigt model ( Fig. 1.2a), made by a Newtonian damper and a Hookean elastic spring connected in parallel, and the standard linear solid model, which combines the Maxwell model and a Hookean spring in parallel. In 1950 Oldroyd proposed a famous viscoelastic model [6], often called Oldroyd-B model ( Fig. 1.2b), where the fluid is assumed to consist of dumbbells, beads connected elastic springs. In a frame-independent form, it can be expressed in terms of the upper-convected derivative of the stress tensor where τ ij is the stress tensor, λ the relaxation time, η m the material viscosity and S ij the rate of strain tensor. Although the model provides good approximations of viscoelastic fluids in shear flow, it has an unphysical singularity in extensional flow, where the dumbbells are infinitely stretched [7]. In order to overcome this problem, the finite elastic non-linear elastic (FENE) model has been proposed; it consists of a sequence of beads with non-linear springs, with forces governed by the inverse Langevin function. Subsequently, the finite elastic non-linear extensibility-Peterlin (FENE-P) model has been developed, by extending the dumbbell version of the FENE model and assuming the Peterlin statistical closure for the restoring force. The model is suited for numerical simulations, since it removes the need of statistical averaging at each grid point at any instant in time, and because the polymer suspension is treated as a continuum and its dynamics represented by an evolution equation of the phase-averaged configuration tensor C ij , a symmetric second-order tensor defined as C ij =< q i q j >, where q i are the components of the end-to-end vector for a polymer molecule. The evolution of the polymer conformation is governed by the balance of stretching and restoring forces in an Eulerian framework, such that the transport equation for the conformation tensor can be expressed as where τ ij is the polymeric stress tensor, defined as with L the maximum polymer extensibility, δ ij the Kronecker delta and λ the polymer relaxation time. A non-dimensional number can be defined based on the polymer relaxation time λ, which is usually called Weissenberg number W e, and is defined as The previous transport equation is a balance between the advection of the configuration tensor on the left-hand side, and the stretching and relaxation of the polymer, represented by the first two terms and the last one on the right-hand side, respectively. Polymer stresses result from the action of polymer molecules to keep their configuration close to the highest entropic state, i.e. the coiled configuration (see Refs. [3,8]). The polymer stress is then added to the momentum equation, the Navier-Stokes equation for an incompressible flow in the case of polymer solutions. The numerical solution of Eq. (1.6) is cumbersome, and many researchers showed that the numerical solution of a viscoelastic fluid is unstable, especially in the case of high Weissenberg numbers, since any disturbance amplifies over time [9][10][11]. Indeed, the numerical solution of this equation can easily diverge and lead to the numerical breakdown since it is an advection equation without any diffusion term [12]. One of the earliest solution to this problem has been to introduce a global artificial diffusivity (AD) to the transport equation of the conformation tensor [11,13,14] by adding to the right-hand side of Eq. (1.6) the term k where k is a coefficient. Subsequently, global AD was replaced by local AD, where the diffusion is applied only to locations where the tensor C ij experiences a loss of positiveness. Recently, researchers started to use high-order weighted essentially non-oscillatory (WENO) schemes [15] for the advection terms in the equation. WENO scheme are non-linear finite-volume or finite-difference methods which can numerically approximate solutions of hyperbolic conservation laws and other convection dominated problems with high-order accuracy in smooth regions and essentially non-oscillatory transition for solution discontinuities. Apart from that, the governing differential equations can be solved on a staggered grid using a second-order central finite-difference scheme. This methodology has been proved to work properly by Sugiyama et al. [16] and also successfully used in Refs. [17][18][19]. A comprehensive review on the properties of different numerical schemes for the advection terms is reported by Min et al. [9].
An alternative methodology to overcome such problems is the so-called logrepresentation of the conformation tensor that ensures the positive-definiteness of the tensor C ij , even at high Weissenberg number [20][21][22][23]; this consists in solving equivalent transport equations for A ij = log C ij , instead of the ones for the conformation tensor C ij . Following the notation used in Ref. [23], we write A = log C = R log DR T , where D is a diagonal matrix containing the eigenvalues of C and R an orthogonal matrix containing the eigenvectors of C. First, we define a decompose of the velocity gradient such that (∇u) T = + B + NC −1 ; note that, and N are antisymmetric and that B is traceless, symmetric and commutes with C. Next, we introduce four new matrices, M = R T (∇u) T R, = R T R, B = R T BR and N = R T NR, and rewrite the decomposition of the velocity gradient as M = + B + ND −1 . Note that, in order to ensure a unique decomposition B is diagonal, while and N are antisymmetric matrices. N and can then be found by satisfying the equations and Finally, the original transport equation for C is rewritten into an equivalent one for A where e A = RDR T and e −A = RD −1 R T .

Plastic Effects
Viscoplasticity is a theory in continuum mechanics that describes the rate-dependent inelastic behaviour of solids. Rate dependence in this context means that the deformation of the material depends on the rate at which loads are applied.
The first viscoplastic rheological model based on yield stress (stress at which a material begins to deform plastically) was proposed by Schwedoff [24] as a plastic viscoelastic version of the Maxwell model: whereε is the rate of deformation, η m the solid viscosity and τ 0 the yield stress. The previous model states that when the stress τ is less than the yield stress τ 0 , the material is completely solid, and the rate of deformation is zero, while when the stress is greater than the yield value, it behaves as a fluid. Note that, at steady state, we obtain τ = τ 0 + η mε . Bingham [25] proposed a similar model: which can be rewritten as (1.14) Bingham model is exactly equivalent to the steady case of the one proposed by Schwedoff for positive rates of deformation. In 1947, Oldroyd modified the Bingham model and proposed the following constitutive equation [26]: which combines the yielding criterion with a linear Hookean elastic behaviour before yielding and a viscous behaviour after yielding. Differently from the previously described models, here when the stress is less than the yield value, the material is not completely rigid. The numerical simulation of a Bingham fluid is not a straightforward task, because of the mathematical non-smoothness of the model and the indeterminacy of the stress tensor below the yield stress threshold [27]. Two kind of solution methods has been proposed in the literature, the regularisation approach [28][29][30][31][32] and the augmented Lagrangian [33][34][35][36][37][38][39][40] algorithm. The former solution method consists in modifying the constitutive equation in order to avoid the numerical and mathematical complexities, while the second consists in solving the whole problem as a minimisation of a functional with a step descent Uzawa algorithm [41]. In other words, the former method consists in solving modified equations which are computationally more permissive, while the second solves the actual yield stress model, but is computationally much more expensive. Among the first category of regularised approaches, in 1987, Papanastasiou [29] developed a modified constitutive relation for Bingham plastics whose main feature is that the tracking of the yield surfaces is completely eliminated. The model assumes where M is a constant that, when chosen sufficiently big, provides a quick stress growth even at relatively low strain rates. This behaviour is consistent with materials in their practically unyielded state, i.e. plastic material that exhibits little or no deformation up to a certain level of stress determined by the yield stress. Due to the fast growing stress, this model has been sometimes used to represent fluids that exhibit extreme shear-thickening behaviour.
Motivated by experimental observations, where yield-stress fluid have an elastic response, Saramito [42,43] combined the Bingham and Oldroyd models, and proposed a model for elastoviscoplastic fluids ( Fig. 1.2c) where the total stress is again σ = ηε + τ . While Schwedoff proposed a rigid behaviour when |τ | ≤ τ 0 and Oldroyd a change of model when reaching the yield value, Saramito assures a continuous change from a solid to a fluid behaviour of the material. The mechanical model is composed by a friction element inserted in the Oldroyd viscoelastic model: at stresses below the yield stress, the friction element remains rigid, and the whole system predicts only recoverable Kelvin-Voigt viscoelastic deformation due to a spring and a viscous element η in parallel. Note that, the elastic behaviour τ = με is expressed in differential form and that μ = η m /λ is the elasticity of the spring. As soon as the strain energy exceeds the level required by the von Mises criterion [44], the friction element breaks allowing deformation of another viscous element (η m ), and the material is described by the Oldroyd viscoelastic model. After expanding the time derivative in the previous equation, the general Saramito model can be written as where τ d ij = τ ij − 1 N τ kk δ ij is the deviatoric part of τ ij , with N = 2 or 3 the dimension of the problem at hand, and δ ij the Kronecker delta. Note that for yield stress τ 0 = 0, the Oldroyd-B model is recovered. A non-dimensional number can be defined based on the field stress τ 0 , which is usually called Bingham number Bn, and is defined as The yield stress value of certain materials, for example, liquid metals, is a function of the temperature [45,46]. Indeed, while in crystal solids the yielding involves bond switch in an orderly manner, in metallic glasses this should be determined by bond breakage [47,48]. By computing separately the mechanical and thermal energies that are required for bond breakage, a simple relation between the yield stress and the temperature can be obtained: where T is the ambient temperature, ρ the density, M the molar mass and T g the glass transition temperature. Guan et al. [49] used molecular dynamic simulations and found that the yield strength and the temperature are well correlated through a simple expression where T 0 and τ 0 are viscosity-dependent, normalised constants. The numerical solution of Eq. (1.18), similarly to Eq. (1.6), may be cumbersome. The use of high-order WENO schemes for the advection terms in the equation is suggested to have high-order accuracy in smooth regions and essentially nonoscillatory transition for solution discontinuities [50,51]. The previously discussed log-representation of the equation can be used as well.

Fluid-Structure Interaction
A fully Eulerian formulation of a fluid structure problem can be obtained with a technique similar to the one discussed in the previous sections. Indeed, we can consider fluid and solid motion governed by the conservation of momentum and the incompressibility constraint: where the suffixes f and s are used to distinguish the fluid and solid phase. In the previous set of equations, σ ij is the Cauchy stress tensor. The kinematic and dynamic interactions between the fluid and solid phases are determined by enforcing the continuity of the velocity and traction force at the interface between the two phases where n i denotes the normal vector. The problem at hand can be solved numerically by using the so-called one-continuum formulation [52], where only one set of equations is solved over the whole domain. This is achieved by introducing a monolithic velocity vector field u i valid everywhere obtained by a volume averaging procedure [53,54], i.e.
where φ s is an indicator function expressing the local solid volume fraction. Thus, we can write the stress in a mixture form as A fully Eulerian formulation is obtained after properly defining the fluid and solid Cauchy stress, with examples given in [16][17][18]55].

Microscopic Approaches
In this section we will discuss approaches used to perform interface-resolved simulations of the intrusions defining the microstructures and thus at the origin of the non-Newtonian behaviours described above. We will consider rigid and deformable particles, as well as two-fluid systems. Indeed, recent developments in computational power and efficient numerical algorithms have allowed the scientific community to numerically resolve the microstructure of suspensions in fluids.

Eulerian/Lagrangian Methods
Eulerian/Lagrangian methods are often used to simulate suspension in fluids, and are also called immersed boundary methods (IBM). The main feature of this method is that the numerical grid does not need to conform to the geometry of the object, which is instead replaced by a body force distribution f that mimics the effect of the body on the fluid and restores the desired velocity boundary values on the immersed surfaces. To do that, two separate grids coexist, the Eulerian fixed grid where the flow is solved, and the Lagrangian grid representing the moving immersed boundary (see Fig. 1.3a); a singular force distribution at the Lagrangian positions is first determined and then applied to the flow equations in the Eulerian frame via a regularised Dirac delta function. The primary advantage of the IB method is associated with the simplification of the grid generation task: indeed, grid complexity and quality are not significantly affected by the complexity of the geometry. The advantage of the IB method becomes eminently clear for flows with moving boundaries, where the process of generating a new grid at each time step is avoided, because the grid remains stationary and non-deforming. A drawback of this approach is that the grid lines are not aligned with the body surface, so in order to obtain the required resolution, higher number of grid points may be required. Many IBMs have been created so far, which differ in the way the immersed boundary force is computed [56][57][58][59][60]. The different methods are often grouped in two categories, continuous and direct forcing: in the first approach the forcing is incorporated into the continuous equations before discretisation, whereas in the second approach the forcing is introduced after the equations are discretised. The continuous forcing approach is very attractive for flows with immersed deforming boundaries, whereas the direct one is more commonly used to simulate rigid boundaries. The original IB method was developed by Peskin [61] for the coupled simulation of blood flow and muscle contraction in a beating heart and is generally suitable for flows with immersed elastic boundaries. The IB is represented by a set of elastic fibres and the location of these fibres is tracked in a Lagrangian fashion by a collection of massless points moving with the local fluid velocity, i.e. the coordinate X k of the k-th Lagrangian point is governed by the equation where u is the local fluid velocity. The stress (denoted by F ) is related to deformation of these elastic fibres by a constitutive law, such as the Hooke's law, and the effect of the IB on the surrounding fluid is captured by transmitting the fibre stress to the fluid through a localised forcing term in the momentum equations where δ is the Dirac delta function. Because the location of the fibres does not generally coincide with the nodal points of the Cartesian grid, the forcing is distributed over a band of cells around each Lagrangian point and added on the momentum equations of the surrounding nodes. Thus, the sharp delta function is replaced by a smoother distribution function, denoted here by d, so that the forcing at any grid point x i,j is given by The fibre velocity in Eq. (1.25) is also obtained through the use of the same smooth function. The choice of the distribution function d is a key ingredient in this method, and several different distribution functions have been derived and employed in the past [61][62][63][64].
In the same spirit, Goldstein et al. [65] developed another model, called feedback forcing, to simulate the flow around rigid and moving bodies, where the effect of the body on the surrounding flow is modelled through a forcing term of the form where the coefficients α and β are selected to best enforce the boundary condition at the immersed solid boundary, whose velocity is V . The above relation is a feedback to the velocity difference u − V and behaves in such a way to enforce u = V on the immersed boundary. Indeed, the first term on the right-hand side of the equation tends to annihilate the difference between u and V , whereas the second term can be interpreted as the resistance opposed by the surface element to assume a velocity u different from V . In an unsteady flow the magnitude of α must be large enough so that the restoring force can react with a frequency which is bigger than any frequency in the flow; however, big values of α and β render the forcing equation stiff and its time integration requires very small time steps. The method has been used to simulate flexible filaments as well [66,67]. Even if the original intent behind Eq. (1.28) is to provide feedback control of the velocity near the surface, from a physical point of view it can also be interpreted as a damped oscillator [68] with frequency 1/ (2π ) √ α and damping coefficient −β/ 2 √ α .

Immersed Boundary Methods for Suspensions of Rigid Particles
Uhlmann [56] proposed a computationally efficient numerical method based on the IBM to simulate suspension of rigid particles. Firstly, the surface of the immersed surface delimiting the body is discretised using N markers, called Lagrangian points X; note that, in general they do not correspond to the grid nodes x. The solution of the incompressible Navier-Stokes is based on the fractional-step method [69]. Indeed, a simple prediction step is first performed, without taking into account the immersed object. The obtained velocity field u * is then interpolated (with an interpolator operator I) onto the embedded geometry , The values of U * are used to determine a distribution of singular forces along the boundary of that restore the prescribed boundary values U as (1.30) The force field defined over is then transformed into a body force distribution applied to the Eulerian grid using a convolution operator C The momentum conservation equation is then solved again with the computed volume force field added as a source term and the time advancement step is completed with the usual solution of the pressure Poisson equation and the projection step where velocity and pressure are corrected to ensure mass conservation. Note that, this procedure is common to the most modern IB methods, and the step that defines the method is the way in which the operators I and C are built: in particular, the interpolation and spreading operations are based on the regularised Dirac delta function by Roma et al. [70], which extends over three grid cells in all coordinate directions. The desired velocity U at a location X on the interface between the fluid and the immersed boundary is given by the rigid-body motion of the solid object: (1.32) where r = X − x c is the position vector relative to the particle centroid, u c is the translational velocity of the particle centroid and ω c is the angular velocity of the particle. The translational and angular velocities of a particle are described by the Newton-Euler equations, which for a sphere reduce to and In the previous two equations, Eqs. (1.33) and (1.34), ρ p is the density of the particle, V p its volume (4/3πR 3 for a sphere with radius R), τ the fluid stress tensor, n the outward-pointing unit normal at the surface δV p of the particle, g the gravitational acceleration and I p the moment of inertia of the particle (2/5ρ p V p R 2 for a solid sphere). F c and T c represent the force and torque acting on the particle as a result of collisions and contact with other particles or solid walls. Breugem [59] proposed two major improvements to the method discussed above. The first is the so-called multidirect forcing scheme. The use of a regularised Dirac delta function for the interpolation and spreading operations results in a diffuse distribution of the IBM force around the interface and because of that, the influence region of neighbouring Lagrangian points overlaps. Eulerian grid points in the overlap region are used to enforce the boundary value multiple times; thus, the resulting forcing is perturbed and the final distribution of the IBM force may not properly enforce the desired boundary condition. The remedy for this problem is to iteratively determine the IBM forces on the relevant Eulerian grid points such that they collectively enforce the desired boundary condition at the different Lagrangian points [71,72]. The second improvement suggested by Breugem is the inward retraction of Lagrangian grid. The delta function of Roma et al. [70] has a width of three grid cells, and because of that the (outer) radius of the particle actually increases from R to R + 3 x/2; this effect results in an increase of the particle drag force, which is partially balanced by an overall permeability of the particle due to a non-perfect boundary condition. As shown by Breugem, the first inaccuracy (increase in the effective radius) is stronger than the other one, and the suggested solution is to slightly retract the Lagrangian points from the surface towards the interior of the particle [73,74].
In Eqs. (1.33) and (1.34), F c and T c are the force and torque acting on the particle as a result of collisions and contact with other particles or solid walls. A recent model for particle-particle and particle-wall interactions in interfaceresolved simulations of particle-laden flows is described by Costa et al. [75]. The model consists of three different interactions: long-and short-range hydrodynamic interactions, and solid-solid contact. The long-range interactions are directly obtained by the immersed boundary method, while the short-range ones are based on a lubrication model employed when the gap between particles is below the grid size, and is based on asymptotic expansions of the analytical solution for canonical lubrication interactions between spheres in the Stokes regime. Roughness effects can be accounted for as well. This correction is applied until the particles reach contact when a linear soft-sphere collision model is used. Note that, the approach described above can be extended to particles of different shapes [76], and that alternative collision models can be found in the literature, for example, Refs. [77,78].

Front-Tracking Methods for Suspensions of Deformable Droplets
The so-called front-tracking method is an evolution of the immersed boundary method used to simulate viscous, incompressible, immiscible two-fluid systems, first developed by Unverdi [79] and Tryggvason [80]. In such multiphase problems, the density and viscosity fields of each fluid remain constant, but they are discontinuous across the interface. In order to avoid numerical diffusion or oscillations problems close to the jump, these fluid properties are not advected directly, but instead a Lagrangian grid is created to describe the boundary between the different fluids, which is then moved with the fluid velocity. Therefore, at every time step it is necessary to reset the fluid properties and to do so, an indicator function F (x) is also introduced, equal to 1 inside one fluid and 0 in the other one. This function is constructed from the known position of the interface and is used to evaluate the proper values of density and viscosity at each grid point: where the suffixes 1 and 2 indicate the two fluids. The jump in the indicator function carried by the interface is spread to the grid points nearest to the interface, in order to ensure that the fluid properties change smoothly across the interface. This generates a grid-gradient field which is zero everywhere except near the interface and has a finite thickness of the order of the mesh size. The spreading of the jump onto the grid is done in such a way that the volume integral of the gradient is conserved, i.e. if G (x) is the gradient of the indicator function evaluated at a stationary grid point x, and D is a distribution function, such as the one introduced by Peskin [81], then 36) where N k is the unit normal vector to the interface element of area S k whose centroid is at X k . The indicator function is found everywhere by solving the following Poisson equation: where the right-hand side is computed by simple numerical differentiation.
Since the fluid velocities are computed on the fixed grid and the front moves with the fluid velocities, the velocity of the interface points must be found by interpolation; thus, similarly to Eq. (1.36), to interpolate the velocity on the interface Lagrangian points we use where the sum is now over the points on the stationary grid in the vicinity of the considered k-th Lagrangian point. Finally, the new position of the interface is found by solving a simple advection equation As the front moves, it deforms and stretches, and the resolution along some parts of the interface can become inadequate or overly crowded. To maintain accuracy, either additional elements must be added when the separation among points becomes too large or points must be redistributed to maintain adequate resolution.
In those calculations where the surface tension σ is needed, the magnitude of the surface tension force is obtained from the local curvature K of the interface: F k = σ K k N k S k . This force is then distributed onto the grid as (1.40) Note that, alternative ways have been proposed to calculate the surface tension without having to find the curvature, but only the local tangent defined by the Lagrangian points on the interface [80].

Eulerian/Eulerian Methods
When dealing with moving and deformable boundaries, an alternative approach to the IB front-tracking methods previously discussed are the so-called front-capturing methods, which are fully Eulerian and handle topology changes automatically. A strong advantage of such methods is that they are easier to parallelised and typically achieve higher efficiency than their Lagrangian counterpart. However, interactions between approaching particles and droplets are difficult to control and may depend on the resolution adopted. Eulerian interface representations include essentially the volume of fluid (VOF) [82] and level-set (LS) [83][84][85] methods, their variants and combination. The VOF method defines different fluids with a discontinuous colour function, and its main advantage is an intrinsic mass conservation; however, it suffers from an inaccurate computation of the interface properties, such as normals and curvatures [86,87]. Contrary to the VOF, the LS method prescribes the interface through a (Lipschitz-)continuous function which usually takes the form of the signed distance to the interface. Thus, normals and curvatures can be readily and accurately computed, while mass loss/gain may occur since the LS function has no volume information. Furthermore, this approach requires a procedure to reshape the LS into a distance function, i.e. the reinitialisation step. More recently, researchers have started to develop coupled VOF-LS methods [88] in order to overcome the disadvantages of both the techniques.

Volume of Fluids
We introduce an indicator (or colour) function H to identify a given fluid so that H = 1 in the region occupied by fluid 1 and H = 0 in fluid 2. Considering that the fluid is transported with the flow velocity, we update H in the Eulerian framework by the following advection equation: where u is the velocity vector. The cell-averaged value of the indicator function is defined as the volume fraction or volume of fluid (VOF) function within a cell Thus, the VOF function assumes values 0 ≤ φ ≤ 1 (see Fig. 1.3b). Combining the two previous equations, we obtain the advection equation of the VOF function in the divergence form: In a conventional VOF method, the interface separating different fluids is piecewise reconstructed for each cell by straight line segments, which are then used to calculate the numerical fluxes necessary to update the VOF function. This geometric reconstruction effectively eliminates the numerical diffusion that smears out the compactness of the transition layer of the interface. Different methodologies have been proposed to accurately recover the exact surface geometry from the discretised VOF function: the simple line interface calculation (SLIC) method [89], the piecewise linear interface calculation (PLIC) [90,91], the latter being further modified by several authors [92][93][94][95][96]. Another technique is the tangent of hyperbola for interface capturing (THINC) method [97]: this avoids the explicit geometric reconstruction by using a continuous sigmoid function rather than the Heaviside function, thus allowing a completely algebraic description of the interface and enabling the computation of the numerical flux. An improvement was proposed by combining the original THINC method with the first-order upwind scheme in the so-called THINC/WLIC (THINC/weighted linear interface capturing) method [98].
Recently, the method has been further developed in the multi-dimensional THINC (MTHINC) method where the fully multi-dimensional hyperbolic tangent function is used to reconstruct the interface [99,100]. The numerical fluxes can be directly evaluated by integrating the hyperbolic tangent function which also prevents the numerical diffusion that smears out the interface transition layer. Another advantage of the method is that the normal vector, curvature and approximate delta function can be directly obtained from the derivatives of the function, thus the standard smoothing or convolution techniques used in conventional VOF methods is not required. The unit normal vector is defined as n = m/|∇φ|, being m the gradient of the VOF function, i.e. m = ∇φ, which can be computed using the usual Young's approach [90,91]. For example, in the 2D case, first the values of the derivative at the four cell corners indexed as i ± 1 2 , j ± 1 2 are calculated by the VOF function on its surroundings, for example, 44) and then averaged to find the cell-centre value (1.45) The curvature k is then found by taking the divergence of the normal vector k = −∇ · n. (1.46) The surface tension force f = σ knδ can be computed using the continuum surface force (CSF) model [101], where the 1D approximate delta function δ is directly approximated by δ ≈ |∇φ|. Thus, we obtain (1.47) Finally, the mixture density and dynamic viscosity are simply averaged in terms of the VOF function (similarly to Eq. (1.35)): Due to the non-uniformity of the density, the Poisson equation used to enforce a divergence-free velocity field becomes 49) which is in an equation with variable coefficients. In order to utilise an efficient FFT-based pressure solver with constant coefficients [5,102], we use the following splitting of the pressure term [103]: where ρ 0 is a constant density equal to the lowest density between the two phase. With this splitting, the Poisson equation can be rewritten as Note that, the correction step of the fractional-step method needs to be modified accordingly and that the term 2p n − p n−1 is a linear extrapolation consistent with using a second-order scheme to integrate in time the momentum equation, for example, Adams-Bashforth. The methodology can be extended to the Runge-Kutta schemes by following Ref. [104].

Level-Set Method
The level-set function φ approximates the signed distance from the interface, thus φ = 0 denotes the interface and φ > 0 or φ < 0 the two different fluids separated by it (see Fig. 1.3c). The motion of the interface is governed by the following transport equation (formally similar to Eq. (1.41)): where u is the flow velocity field. The equation is closed and allows to solve the system of equations in a fully Eulerian fashion. Notwithstanding the formal simplicity of the equation, its numerical solution is challenging. The time integration is often performed by a three-stage total variation diminishing third-order Runge-Kutta scheme [105,106], while the advection term of the equation is usually discretised with one of the following schemes: the high-order upstreamcentral (HOUC) scheme [107], the weighted essentially non-oscillatory (WENO) scheme [108], the semi-Lagrangian scheme [109] or the semi-jet scheme [105]. Quantitative comparisons of these schemes in various test cases can be found in Ref. [110]. As reported by Ge et al. [106], for flows involving moderate deformations, the HOUC scheme is usually sufficient and the most efficient, while for more complex flows, the WENO or the semi-Lagrangian/jet schemes combined with grid refinement should be used. Note that, Eq. (1.52) does not need to be solved over the entire computational domain, but only near the φ = 0 value, where the interface is located and the normal and curvature are needed. Thus, in the so-called narrow band approach [111,112], the level set is computed and stored only around the interface, and fast computation and low memory usage may be achieved. Although the level-set function is initialised to be a signed distance, this property is lost with time, causing numerical issues in the evaluation of the normal and curvature [83]. These issue requires an additional treatment in order to reshape the level-set function φ into a distance function, i.e. |∇φ| = 1. This is usually performed by converting it into a time-dependent Hamilton-Jacobi equation [83] ∂φ ∂T φ 0 being the level-set field before redistancing, T a pseudo-time and S(φ 0 ) the mollified sign function of the original level set. When the steady state solution of the equation is reached, the zero level-set contour is unaltered, while the rest of the field has re-obtained the property of being a signed distance function. In practice, this equation is iterated only for few steps towards its steady state every certain number of time steps. Although an alternative approach exists, i.e. the fast marching method (FMM) [84], the reinitialisation procedure allows the use of high-order schemes and is easy to parallelised; thus, it has been a much more popular choice.
Finally, the unit normal vector n and the local mean curvature κ can be simply computed directly from the LS function as follows: n = ∇φ |∇φ| and κ = −∇ · n, (1.54) and the body force f due to surface tension, included in the momentum equation of the Navier-Stokes equation, is expressed as being δ the delta function and σ the surface tension. The density and viscosity vary across the fluid interface, and can be expressed in a mixture form (similarly to Eq. (1.35)) as where H (φ) is the regularised Heaviside function defined such that it is zero inside one fluid and unity in the other one. In order to utilise an efficient FFT-based pressure solver with constant coefficients, the technique described by Dodd and Ferrante [5,102] and discussed in the previous section can be easily used.
We conclude by noting that, although easy to implement, CSF effectively introduces an artificial spreading of the interface. An alternative approach is the ghost fluid method (GFM) [113][114][115], which provides a finite-difference discretisation of the gradient operator even if the stencil includes shocks.

Phase-Field Methods
From a macroscopic point of view, the interface between two immiscible fluids can be usually assumed to be sharp, since its thickness is of the order of few nanometres. This aspect motivated the numerical methods discussed above, that, however, need to overcome the difficulties associated with evolving a discontinuous front through Eulerian fields. An alternative approach has been therefore proposed and successfully used for different applications, the so-called phase-field (or diffuse interface) method, where the interface is modelled as a thin layer across which the fluid properties continuously change avoiding discontinuous fields. The original idea can be attributed to Van der Waals who suggested that the interface thickness in a binary mixture is determined by the balance of counteracting weakly nonlocal terms in the free energy: these, on the one hand, provide diffusion and, on the other hand, a sharpening of the interface. Cahn and Hilliard [116] used this approach in the context of phase separation problems deriving an evolution equation for the concentration field. The thermodynamic consistency of the coupled Cahn-Hilliard/Navier-Stokes model [117] and its ability to handle topological changes are the main reasons that justify the increasing use of phase-field methods [118][119][120][121][122]. The fundamental variable of phase-field models is a scalar (x i , t) which represents the relative quantity of one of the two phases, and whose extreme values, = ±1, correspond to the two pure fluids. Since the two fluids are immiscible, the flow domain is essentially divided in subdomains containing only one of the two phases, separated by an extremely thin interface where the two substances mix, −1 < < 1. The equation describing the evolution of the phases is where G = δF /δ is the thermodynamic chemical potential defined as the functional derivative of the free energy F with respect to , and M is a proportionality constant called mobility. The left-hand side of Eq. (1.57) represents the convective transport of , whereas the right-hand side the driving force from the chemical potential, ensuring phase separation with the exception of the thin layer constituted by the interface. The mobility M has no intuitive physical meaning and is related to the time scale of the interface dynamics. In the limit of vanishing mobility, we recover pure convection neglecting the inner interface dynamics, whereas for infinite mobility the interface will reach equilibrium immediately and do not vary. Following the seminal work of Cahn and Hilliard [116], it is possible to write the free energy in the following form: with σ the surface tension and the interface thickness. The form of the first term between brackets in Eq. (1.58) guarantees the presence of two minima of F for the two pure fluids at = ±1. The gradient for the free energy drives therefore the system towards phase separation. However, the second term in the brackets is proportional to the square of |∇ |. Hence, when reducing the interface thickness, the free energy increases because of the sharpening of the gradient of . The combination of these two terms tends to create an equilibrium interface of finite thickness, mimicking what originally assumed by Van der Waals (see, for example, [120] for more details). Using Eq. (1.58), the chemical potential G can be explicitly written as which can be directly inserted in Eq. (1.57) to compute the phase field. Eq. (1.57) needs to be coupled with the Navier-Stokes equations (convection) which, in turn, are influenced by the phase field because of surface tension. The surface tension across the relatively thick interface between the two phases can be expressed in terms of the free energy (and chemical potential). From the consideration that, without dissipation, the total energy of a systems needs to be conserved it is possible to find the expression for the momentum forcing arising from the interface (surface tension) [123], (1.60) Considering Eqs.(1.57), (1.59) and (1.60), together with the incompressible Navier-Stokes system, we can write the complete Navier-Stokes/Cahn-Hilliard model for the immiscible fluid dynamics. In the following it is presented directly in dimensionless form: (1.63) Here, the equations have been made dimensionless using the typical flow scales, L for length, U for velocity, μ and ρ for dynamic viscosity and density. We have also introduced the Cahn number Cn = /L as the dimensionless measure of the interface thickness, the Weber number, We = ρU 2 L/σ as the ratio between convection and surface tension effects, the dimensionless mobility M * = 3 Mσ/(2 √ 2UL 2 ) measuring the (dimensionless) mobility intensity and the Reynolds number Re = ρU L/μ. When using the phase-field method it is important to correctly set the Cahn Cn and mobility M * numbers. While the other dimensionless numbers used above are classical and depend only on the macroscopic properties of the system, the values of Cn and M * are set also by numerical considerations.
Since the physical thickness of the interface is usually on the nanometre scale, realistic values of the Cahn number for micro-and macro-scale applications may range between 10 −3 and 10 −10 . As the interface needs to be resolved using 4 ÷ 5 control points for accurate numerical simulations, the use of the realistically tiny values of the Cahn number is unpractical in simulations and an artificial thickening of the interface is necessary. Clearly, a trade-off is necessary when deciding the numerical interface thickness; to this end, it has been proved that the dynamics is not affected by the artificial thickening if the typical interface size is well below the smallest flow length scales. In other words, the Cahn number should be small, yet it can be much larger than the real one, for example, around 10 −2 ; this is the so-called sharp-interface limit, [124]. In this limit, any further decrease of the Cahn number does not produce appreciable differences in the macroscopic dynamics.
More difficult is to precisely define or even measure the mobility M [123]. Moreover, it is crucial to determine how the mobility depends on the interface thickness because of its artificial magnification necessary in simulations. Recently, it has been found using matching asymptotic expansions that the optimal scaling to recover the sharp-interface limit is M * = α Ca 2 (α an order one constant), see Ref. [120] for more details.
The main advantage of phase-field methods is in their strict thermodynamic derivation and the ability to automatically handle interface topology changes, the former opening also possibilities for direct simulations of phase changes. The main drawback is represented by the high computational costs since it is important to use at least 4 ÷ 5 points within the interface thickness which, in turn, needs to be smaller than the flow length scales. The Navier-Stokes/Cahn-Hilliard model has recently been adopted in several studies of immiscible fluids both in laminar and turbulent conditions, for example, Refs. [122,[125][126][127]. Besides the Cahn-Hilliard formulation, different phase-field methods have been proposed, based on alternative forms of the free energy, see, for example, [128] for compressible flows with phase change for cavitation problems. Similar methods exist also in the Lattice-Boltzmann framework, for more details the reader is referred to [129][130][131] and the references therein.

Other Approaches
In many applications the suspended objects are much smaller than the smallest hydrodynamical scale; in these cases the so-called point particle method can be used, in which the particle are treated as Lagrangian points moving with the local flow velocity. When the volume fraction of the particles is small enough hydrodynamic interactions and collisions among particles can be effectively neglected, while for large values of the particle-to-fluid density ratio, i.e. for significant mass loads, the momentum exchange between the two phases is still significant and must be accounted for. In general, the motion of the particles is described by a set of ordinary differential equations for the particle velocity and position, with the velocity depending on the Stokes drag and the buoyancy force [132] but also added mass and history force can become important [133,134]; note that, the Stokes drag coefficient is sometimes corrected with empirical correlations when the particle Reynolds number does not remain small. The interested readers are referred to Refs. [135][136][137][138] for more details. Other approaches have been proposed to model the feedback of the point particle on the flow, which should be independent of the resolutions used; details can be found in Refs. [139][140][141][142].
In the recent years, various alternative methods have been proposed which are modifications and/or combinations of the methods previously discussed. A commonly used approach is the force coupling method: this is a bridge between methods for Stokes flow and for low-Reynolds number conditions and is based on a low-order, finite force multipole representation of the effect the particles have on the surrounding fluid flow. In particular, the full Navier-Stokes equations are solved with an additional force density: the force monopole corresponds to the force the particle transmits to the fluid if it were replaced by a rigid particle of mass with the same density as the fluid, while the force dipole is a combination of a symmetric stresslet and a torque that acts on the fluid: the torque is set in a similar manner as the force monopole in terms of the angular momentum of the displaced fluid, while the stresslet is chosen to ensure that the average rate of strain within the particle is zero. The aim is to create a flow outside the particle that matches the actual flow within a short distance from the surface. Note that, the fluid inside the particle volume is an active part of the simulation and satisfies the same integral moments as a rigid particle. The interested reader is referred to the detailed review in Ref. [143] and the references therein.
Another technique is the so-called volume of fluid tensorial penalty method: the method is based on the one-fluid formulation modified for dealing with particle flows. In particular, the fluid and solid phases are treated as two fluids with different rheological property and distinguished by a phase function, such as the solid volume fraction. The solid particle behaviour is recovered in the Navier-Stokes equations by properly decomposing of the stress tensor: in particular, the stress tensor is rewritten in a way to separate the compression, tearing, shearing and rotation contributions [144]. The decomposition is used to separate the stress components operating in a viscous flow and to facilitate the implementation of a numerical penalty method. In the solid phase it is imposed that the local flow admits no shearing and tearing while preserving a constant rotation; these flow constraints are implicitly transmitted to the particle sub-domain as they are solved with the fluid motion. Note that, the previous viscous penalty method is formally equivalent to choosing a viscosity much bigger than the fluid one, similarly to what previously done when discussing plastic effects. However, on a discrete point of view, the two formulations are not equivalent: when the dynamic viscosity is used to impose the solid behaviour, a first-order convergence in space is usually obtained and a rasterisation effect can be produced at the particle-fluid interface; on the other hand, when the viscous stress tensor splitting is used, a more accurate fluid-solid interface is used, thus reducing the rasterisation effect and a second-order convergence can be achieved. To solve the unsteady Navier-Stokes equations together with the incompressibility and solid constraints, the augmented Lagrangian method can be applied. The interested readers are referred to Refs. [145][146][147].

Conclusions
In this chapter, we have introduced present numerical methods for complex fluids. In the first part, we have discussed continuum approaches to viscoplastic and viscoelastic fluids, whereas interface-resolved methods for the simulations of suspensions of rigid and deformable particles, droplets, bubbles have been presented in the second part.
Continuum approaches require rheological models and constitutive equations for the additional stresses due to the fluid microstructure. These are derived theoretically or from experimental data; however, recently, the fast development of computational resources has enabled us to also resolve the microstructure with numerical simulations of the type described in the second part of this chapter. Although important results are continuously reported on the behaviour of viscoplastic and viscoelastic fluids, both in laminar and turbulent flows, these are somewhat restricted to relatively simple geometries. It is therefore relevant to explore the performance of the current numerical algorithms in more complex geometry and thus come closer to industrial applications and natural flows as, for example, avalanches. In addition, new and more sophisticated models are continuously proposed, also thanks to the development of new experimental techniques able to probe fluids subject to time varying shear or stress. From a numerical point of view, these will pose new challenges which we may not be able to tackle with the tools currently available. One such example is the isotropic kinematic hardening idea, based on the concept that the material yield stress builds up and evolves in time together with the flow field, where the steady state yield stress is determined via the back stress modulus (a new material parameter) and the deformation of microstructure (a hidden internal dimensionless evolution variable). Another important point for industrial applications is represented by the case where the fluid stress induces a breakage in a solid structure, i.e. hydraulic fracturing. Once a break-up occurs, the fragment behaves as a suspended particle that may contribute to additional solid breakage events. Despite its relevance in applications, the (automatic) numerical handling of this complex dynamics is still challenging.
Examples of simulations of multiphase and multicomponent flows are numerous and ever increasing, so that the approaches reported here constitute necessarily a quick and limited overview. Two directions are seen here as emerging and of potential interest. On one side the study of intrusions in a complex fluid, ubiquitous in applications, and, on the other side, the need to improve short-range nearcontact interactions among particles/droplets. Indeed, when particles/bubbles are significantly larger than the colloids or macromolecules providing plastic and elastic effects, it is reasonable to model the complex fluid as intrusions in a non-Newtonian matrix, thus combining the two types of approaches presented here, denoted as macro and microscopic.
Although the approaches discussed here fully resolve the surface of the suspended phase, rigid or deformable, solid or fluid, models are necessary for the shortrange interactions: (i) from a numerical point of view, the grid resolution inevitably becomes coarse when two objects approach each other, and (ii) microscale chemical and physical effects determine different interactions, for example, repulsive or attractive forces, slippage, Marangoni effects and depletion forces associated with the microstructure of the suspending phase. These become therefore truly multiscale problems, ranging from nanoscale surface interactions to millimetre/centimetre flows. Coupling of the numerical approaches mentioned here with nanoscale simulations, for example, molecular dynamics, is, therefore, a very active research area.
We would like to conclude by saying that we believe it is fundamental to always be aware of the limitations of the simulations one wishes to perform. This allows the researcher to design configurations and propose, though sometimes unphysical, numerical experiments that can unveil important physics, something which may not be possible in laboratory experiments.