An IB Method for Non-Newtonian-Fluid Flexible-Structure Interactions in Three-Dimensions

Problems involving fluid flexible-structure interactions (FFSI) are ubiquitous in engineering and sciences. Peskin’s immersed boundary (IB) method is the first framework for modeling and simulation of such problems. This paper addresses a three-dimensional extension of the IB framework for non-Newtonian fluids which include power-law fluid, Oldroyd-B fluid, and FENE-P fluid. The motion of the non-Newtonian fluids are modelled by the lattice Boltzmann equations (D3Q19 model). The differential constitutive equations of Oldroyd-B and FENE-P fluids are solved by the D3Q7 model. Numerical results indicate that the new method is first-order accurate and conditionally stable. To show the capability of the new method, it is tested on three FFSI toy problems: a power-law fluid past a flexible sheet fixed at its midline, a flexible sheet being flapped periodically at its midline in an Oldroyd-B fluid, and a flexible sheet being rotated at one edge in a FENE-P fluid.


Introduction
Phenomena of fluid flexible-structure interactions (FFSI) are ubiquitous in various fields of engineering and sciences. For instances, clothes moving in a washing machine, booms floating in ocean preventing oil pollution, blood flowing in deformable blood vessels, to name just a few. Due to complexity of such problems, mathematical modeling and computer simulation of FFSI problems are challenging. The first methodology for modeling and simulation of FFSI problems is probably the immersed boundary (IB) method introduced by Peskin [Peskin (1972[Peskin ( , 1973]. Since the birth of the IB method, numerous methods for FFSI problems are developped. These include immersed boundary method [Iaccarino and Verzicco (2003); Mittal and Iaccarino (2005)], Arbitrary Lagrangian Eulerian (ALE) [Hughes, Liu and Zimmermann (1981); Donea, Giuliani and Halleux (1982); Yang, Sun, Wang et al. (2016)], the lattice Boltzmann [Lallemand and Luo (2003)], fictitious domain [Glowinski, Pan and Periaux (1994a,b)], front tracking [Glimm, Grove, Li et al. (1998)], immersed interface [Leveque and Li (1994); LeVeque and Li (1997); Li and Lai (2001)], blob-projection [Cortez (2000)], phase field [Sun, Xu and Zhang (2014); Wick (2016); Zheng and Karniadakis (2016); Mokbel, Abels and Aland (2018)], immersed John et al. (2018)] are developping an immersed-boundary lattice-Boltzmann method for viscoelastic fluids including Oldroyd-B and FENE-CR fluids. In their work the motion equations of the solid are solved for by finite difference or finite element methods. In this paper we discuss our recent development of a 3D IB methods for three non-Newtonian fluids: power-law, Oldrod-B [Oldroyd (1950)], and FENE [Peterlin (1961)]. These models can describe many non-Newtonian fluids encountered in real-world problems including blood and polymeric fluids. Some preliminary results was reported in a short letter [Zhu (2018)]. Details and more tests of the new method, including the integration of the three non-Newtonian models are presented in the current paper. In our new method, the lattice Boltzmann equations, the D3Q19 model [Qian (1990); SY and GD (1998)], are used to model the fluid flow. The power-law, the Oldroyd-B model, and the FENE-P model are used to model constitutive equations for various non-Newtonian fluids. The formal one (power-law) is incorporated into the lattice Boltzmann D3Q19 model by an algebraic approach; the latter two are numerically solved by the lattice Boltzmann D3Q7 model [Malaspinas, Fiétier and Deville (2010)]. The deformable structure is modelled by elastic fibers which can be stretched, compressed, and bent. The fluid-structure-interaction is modelled by the Dirac delta function, as in Peskin's original immersed boundary method. Note that in our method, the three non-Newtonian models are integrated seamlessly via a model parameter. Selection of a specific model is done by setting a specific numerical value to the parameter. This is a unique feature of our hybrid method thanks to the advantages of using lattice Boltzmann approach for modeling both fluid motions and the corresponding constitutive equations. To examine the new method and demonstrate its capability, we consider three FFSI toy problems-a power-law fluid flow passes a deformable sheet tethered at the middleline, a flexible rectangular sheet is flapped sinusoidally at its midline in a stationary Oldroyd-B fluid, and a flexible rectangular sheet is rotated at one edge in a stationary FENE-P fluid. The remaining article is as follows. Section 2 gives the complete mathematical formulation of the IB method for non-Newtonian fluids. Section 3 discusses details of the numerical methods for the mathematical formulation. Section 4 addresses the verification and validation of the new method and its implementation. Section 5 reports some main simulation results for each of the three test problems. Section 6 concludes the paper with a summary.

Mathematical formulation
The immersed boundary formulation, a nonlinear system of differential-integral equations, describing the motion of a generic flexible structure in a non-Newtonian fluid (power-law, Oldroyd-B or FENE-P) may be written as follows. The equations are listed first and then followed by explanation. Eqs.
(1) and (2) are the classic Navier-Stokes equations for a viscous incompressible fluid where u denotes the fluid velocity, ρ denotes the density, p the pressure. The symbol f ib denotes the Eulerian force applied by the immersed structure to the fluid, b f denotes other external forces exerted on the fluid, for example, the gravity. The letter D = (∇u + ∇u T )/2 denotes the strain rate tensor, η the fluid dynamic viscosity, and Π the viscoelastic stress tensor. Note that the two equations are valid for both Newtonian and non-Newtonian fluids. For Newtonian fluids, η is constant and Π is zero. For non-Newtonian fluids obeying power-law, η is given by Eq.
(3) whereγ is the shear rate, D ij is the ij th component of strain rate tensor, constants η 0 and n are model parameters describing the properties of the non-Newtonian fluid. Note that when n = 1 the fluid becomes Newtonian. For non-Newtonian fluids such as polymeric fluids, the viscoelastic stress Π may be computed from the conformation tensor C (Eq. (4)), where κ p and µ p are the relaxation time and dynamical viscosity of the polymer, respectively. The polymer conformation tensor C is governed by Eq. (4) , i.e., the FENE-P model [Peterlin (1961)]. The variables a and b are given by Eq.
(5) where r p is a model parameter related to the maximum length of polymer molecules which is permitted and tr(C) is the trace of the conformation tensor C. Note that variable a is a nonlinear function of the conformation tensor C and the FENE-P model is a nonlinear system of hyperbolic partial differential equations. When a = b = 1, FENE-P model reduces to the popular linear Oldroyd-B model [Oldroyd (1950)]. Note that the Navier-Stokes equations are parabolic and elliptic in nature, but the FENE-P model (Eq. (4)) is hyperbolic. They are coupled through fluid velocity u and viscoelastic force Π.
The immersed boundary Eulerian force density f ib in Eqs. (1) and (2) can be computed via Eq. (6), where α is the Lagrangian coordinate of the immersed structure. The function δ is the classic Dirac delta function. The symbol X is Lagrangian position of the structure. The function F is the corresponding Lagrangian force density, which is computed from the elastic potential energies (Eq. (8)) of the structure. In Eq. (8), the first integral represents the contribution from stretching/compression (E s ) and the second one represents the contribution from bending (E b ). The quantities K b and K s are bending and stretching coefficients of the elastic constitutive fibers of the deformable structure. Their numerical values are related to the Youngs' modulus and Poisson ratio of the structure [Strychalski, Copos, Lewis et al. (2015)]. The velocity of the immersed structure U (α, t) is interpolated from the velocity of surrounding fluid by Eq. (10). Note that this equation by definition dictates that the immersed structure must follow the motion of fluid because of fluid viscosity. That is, the no-slip boundary condition is enforced on the fluid-solid interface.
There are three major dimensionless ratios in the FFSI problems involving non-Newtonian fluids: Reynolds number Re = UcLc νp+νf , structure flexure modulusK b = Kb ρcU 2 c L 4 c , and Weissenberg number for polymeric fluid W i = κ pγ (or exponent n for power-law fluid). In above definition, ρ c is the characteristic fluid mass density, U c is the characteristic flow speed, L c is the characteristic length of the immersed structure, κ p is the relaxation time of polymer,γ is the characteristic flow shear rate, ν p is the polymer kinematic viscosity, and ν f is the fluid kinematic viscosity. Note that Re measures the ratio of inertial force and viscous force,K b measures the ratio of the elastic force and inertial force, W i measures the ratio of viscoelastic force and viscous force. The index of power-law n < 1 corresponds to shear-thinning fluid, n > 1 corresponds to shear-thickening fluid, and n = 1 corresponds to Newtonian fluid.

Numerical methods
As we can see, the immersed boundary formulation for the FFSI problems of non-Newtonian fluids is a nonlinear system of integral and partial differential equations (PDE). The PDEs are of mixed type (hyperbolic, parabolic, and elliptic). Numerical solution of this kind of hybrid system is challenging. Here we choose to use the lattice Boltzmann approach [Wolf-Gladrow (2000); Guo and Shu (2013); Huang, Sukop and Lu (2015); Succi (2018)] for this system. The lattice Boltzmann approach treats the incompressible flows as slightly compressible flows (which are governed by hyperbolic PDEs). This kind of artificial compressibility approach for the flow equations is consistent with the hyperbolic nature of the constitutive equations of the non-Newtonian fluids such as FENE-P fluids. The details of the discretizations for the flow equations, constitutive equations, and immersed boundary equations are given as below. The lattice Boltzmann equations, the D3Q19 model [Qian (1990); SY and GD (1998)], are used to solve numerically the viscous incompressible flow equations (Eqs. (1) and (2)) for non-Newtonian fluids. Carefully chosen 19 velocities (with three speeds) are used to discretize the particle velocity space ξ, ξ i (i = 0, 1, ..., 18): Use g i (x, t) to denote the single-particle velocity distribution function along the direction of ξ i (i = 0, 1, ..., 18). The lattice Boltzmann equation (LBE) advancing g i (x, t) one-timestep forward along the direction ξ i can be written as Here The weight w i is given as follows: Here c s = c/ √ 3 is the model sound-speed and c is the model lattice-speed: The function g (eq) i is the discretization of the equilibrium distribution function: The method introduced in Guo et al. [Guo, Zheng and Shi (2002)] is applied to treat the external force term. Note that the lattice Boltzmann equation (Eq. (11)) can be regarded as an explicit second-order accurate discretization in space and time of the flow equations for viscous non-Newtonian fluids by a Lagrangian approach [Wolf-Gladrow (2000)]. It is equivalent to a second-order finite difference scheme for the viscous incompressible Navier-Stokes equations [Junk (2001)].
For non-Newtonian fluids modelled by power law, the constitutive equation is incorporated into the lattice Boltzmann flow model (the D3Q19) algebraically by Eq. (3). The shear rate may be calculated as follows:γ = 2D ij D ij , D ij (i = 1, 2, 3; j = 1, 2, 3) is k . Here g k is the velocity distribution function along k th direction, g (eq) k is the equilibrium distribution function, and ξ ki (k = 0, 1, ..., 18; i = 1, 2, 3) is the i th component of the k th discrete direction. See [Zhu, Yu, Liu et al. (2017)] for details. For non-Newtonian fluids whose constitutive laws obeying the FENE-P model (including the Oldroyd-B model), the lattice Boltzmann D3Q7 model for advection-diffusion equations [Malaspinas, Fiétier and Deville (2010)] is applied for numerical solutions of their constitutive equations. It is coupled with the lattice Boltzmann model for flow equations through the flow velocity u and viscoelastic force Π. Lattice Boltzmann particles in the D3Q7 model are allowed to move along six discrete directions ζ i , i = 1, 2, ..., 6 at a node, where Particles can also stay at the node ζ 0 = (0, 0, 0). For a given component of the configuration tensor C αβ , at a given node x, along a given direction ζ i , i = 0, 1, 2, ..., 6, the single-particle velocity distribution function q αβi is evolved according to Where the relaxation time χ is related to the diffusivity constant κ p by χ = 8κp+1 2 . The ratio κ p /µ p is set to be a very small number e.g. 10 −6 , here µ p is polymer dynamical viscosity. The equilibrium distribution q (eq) αβi = w i C αβ (1 + ζikuk c 2 l ), where ζ ik is the k th (k = 1, 2, 3) component of velocity ζ i (i = 0, 1, 2, ..., 6). Function ψ αβ = − 1 κp (aC αβ − bδ αβ ) + C αi ∂uβ ∂xi + C iβ ∂uα ∂xi . The αβ th component of the conformation tensor is computed by C αβ = i=6 i=0 q αβi + 1 2 ψ αβ . The viscoelastic force in Eq. (1) is computed by ∇ · Π = ∇ · ( µp κp (aC − bI)) and the spatial derivative ∂ xi , i = 1, 2, 3 is discretized by the central difference scheme. The macroscopic fluid mass density ρ(x, t) and momentum (hence velocity) ρu(x, t) can be computed from g i (x, t) at each node by Let integer m denote time step: g m = g(x, ξ, m), X m (α) = X(α, m), u m = u(x, m), p m = p(x, m), ρ m = ρ(x, m). Let the flexible structure be discretized by a set of elastic fibers with Lagrangian coordinate α 2 . Let α 2 = k 2 ∆α 2 , where k 2 is an integer. Let a fiber be discretized by a set of points with Lagrangian coordinate α 1 . Let α 1 = k 1 ∆α 1 , where k 1 is an integer. Then the elastic potential energy is discretized by: where the first term corresponds to the stretching/compression energy and the second term corresponds to the bending energy. The Lagrangian force density at node with index l, (F) l , l = 1, 2, ..., N f , is given by Here δ kl is the Kronecker δ: δ kl = 1 if k = l and δ kl = 0 if k = l. The integral equations in the immersed boundary formulation are computed by the trapezoidal rule: Here Γ denotes the immersed structure, Ω denotes the fluid domain. And ∆x, ∆y, and ∆z denote the meshwidth in x, y, and z directions. The Dirac δ function is discretized by δ d : Here φ(r) = 0.25(1 + cos(0.5πr)) for |r| ≤ 2 and is 0 otherwise. For other choices of φ(r) see Peskin et al. [Peskin and McQueen (1996)]. The motion equation of the structure is discretized by For clarity the algorithm of the IB formulation for non-Newtonian fluids is summarized as follows. Suppose all variables are known at time step t (an integer), the procedure for updating all of the variables for next time step t+1 is as follows. 0) Initialization of all variables; 1) Advance the LBE (Eq. (11)) for flow (D3Q19 model) from t to t+1 using f ib and Π from time step t; compute the new fluid velocity u, velocity gradient ∇u, and mass density ρ; 2) Advance the LBE (Eq. (13)) for the constitutive equations (the D3Q7 model) from t to t+1 using u and C from time step t; compute the viscoelastic force ∇ · Π from the newly updated conformation tensor C, 3) Compute the structure velocity U from the fluid velocity u by Eq. (19); 4) Update the structure position by its velocity via Eq. (21); 5) Compute the immersed boundary force exerted by the structure to the fluid using its new configuration via Eq. (17); 6) Convert the Lagrangian force to Eulerian force by Eq. (18); 7) Compute the new equilibrium distribution functions of the q (eq) αβ and g (eq) using the newly obtained fluid velocity u, conformation tensor C, and mass density ρ; 8) Go to 1).
Note that the algorithm for non-Newtonian fluids can be easily combined with algorithm for Newtonian fluids using the D3Q19 model. Therefore, the new method may be implemented in one computer program with a single model parameter m p (an integer) to switch between Newtonian and non-Newtonian fluids. The option m p = 0 selects Newtonian fluid. The code will bypass the D3Q7 model and set Π = 0, n = 1. The option m p = 1 selects the power-law fluid. The code will bypass the D3Q7 model and set Π = 0. The option m p = 2 selects the Oldroyd-B fluid. The code will execute the D3Q7 model with a = b = 1 and set n = 1. The option m p = 3 selects the FENE-P fluid. The code will execute the D3Q7 model and set n = 1. Thus the Newtonian, power-law, Oldroyd-B, and FENE-P fluids are seamlessly integrated together in the new IB method via the lattice Boltzmann approach and they can be implemented in a single computer code.

Verification and validation
The numerical methods involved and their implementations used in the work have been verified and validated in different settings. The lattice Boltzmann method (D3Q19) and its implementation have been verified and validated in Zhu et al. [Zhu, Tretheway, Petzold et al. (2005)]. The lattice-Boltzmann immersed-boundary method (LB-IB) with its implementation for Newtonian fluid flows have been verified and validated in Zhu et al. [Zhu, He, Wang et al. (2011b)]. The LB-IB for power-law fluid is verified and validated in Zhu et al. [Zhu, Yu, Liu et al. (2017)]. The preliminary verification and validation of the LB-IB method for Oldroyd-B/FENE-P have been reported in a short letter [Zhu (2018)].
In this paper, the newly developped LB-IB method for polymeric flows are further tested on two new FFSI toy problems: a flexible sheet being flapped periodically at the middle and being rotated constantly at one edge in stationary Oldroy-B and FENE-P fluids in three dimensions. Many simulations with various dimensionless parameters indicate that the method is conditionally stable. Mesh refinement studies indicate the method is first-order accurate, which is consistent with the IB framework in general.

Test problems
In this section we consider three FFSI model problems. I) a power-law fluid flows around a flexible rectangular sheet fixed at the midline in a three-dimensional rectangular domain; II) a flexible rectangular sheet is flapped sideways (left and right) at the midline sinusoidally in a 3D rectangular box full of an Oldroyd-B; III) a flexible rectangular sheet is rotated constantly at one edge with a constant speed in a 3D rectangular box full of a FENE-P fluid. In case I, the structure is initially stationary. The flow passes around it and causes it to bend and get aligned with the flow. No-slip boundary condition is applied on the top, bottom, front, and rear rigid walls. Constant velocity is specified at the inlet and outlet boundaries (in x-axis). In cases II and III, the fluid is initially stationary and the structure is forced to move. The active motions of the structures drive the fluid flow. Periodic boundary condition is applied along all directions of the computational domain. In all cases, the x− axis points from left to right; the y− axis points from front to rear; the z− axis points from bottom to top. The sheets are composed of two groups of elastic fibers which can be compressed, stretched, and bent. The two groups of fibers are cross-linked and are orthogonal to each other initially. This type of FFSI problem possesses three significant dimensionless parameters: flow Reynolds number Re, structure bending modulusK b , and fluid Weissenberg number W i (or exponent n for power-law fluid). Numerous simulations on the three model problems using various combinations of these parameters are performed to test the capability of the new method. Some representative simulation results are reported below for each case. Case I) An elastic sheet with aspect ratio 1:2 (width versus length) is placed initially on the y − z plane (i.e., vertical) in the middle of the box (in x, y, and z directions). Its midline is fixed in a power-law fluid flow and the sheet is free to move otherwise. This problem was intensively studied by Zhu et al. [Zhu, Yu, Liu et al. (2017)]. More simulations are performed and some typical results from different combination of values of Re,K b , n are shown here. The left panel in Fig. 1 shows the position and shape of the sheet at several equally distributed time instants from a simulation with Re = 80,Kb = 0.0001, n = 0.5. The left most is the initial position (vertical and flat). The right most is the final position (curved). Starting from the initial configuration, the sheet moves and deforms with the flow, and finally sets down to a quasi-steady state. Note that the position/shape of the sheet partially overlap for some instants. The right panel in Fig. 1 shows streamlines of the flow field from a simulation with Re = 120,Kb = 0.05, n = 0.6. The gray surface is the position and shape of the sheet. The streamlines start from the lower half plane of the inlet (using tens of uniformly spaced seeds). Notice the twisted curves behind the sheet. These streamlines come from the lower half plane at inlet and move to the upper half plane after past the sheet. This reveals the complicated flow patterns right behind the sheet. Case II) A deformable sheet of the same aspect ratio is initially placed the same way as in case I. Its midline is now forced to flap sinusoidally (on the x − y plane along x-direction) in an Oldroyd-B fluid. The x-coordinate of the leading edge is given by . Here x(t) is the x-coordinate of the sheet middle-line, A denotes prescribed amplitude of flapping, P denotes flapping period, and t denotes time. Some typical results from different combination of values of Re,K b , W i are given below. In Fig. 2 the left panel plots the shape of the sheet at several equally separated instants from a simulation with Re = 40,Kb = 0.005, W i = 0.1. The left most is the initial position (vertical and flat). The rest is the shape (not physical position) of the sheet at several time instants within a period. Note that the sheet physical position (x-coordinate) is shifted horizontally by a constant for the purpose of displaying the shapes at multiple instants on the same panel. Starting from the initial position, the sheet moves and deforms with the flapping midline. Spontaneous motion of the sheet along y and z directions are not seen. The right panel in Fig. 2 shows streamlines of the flow field from a simulation with Re = 60,Kb = 0.004, W i = 0.1. The gray surface (partially buried in the curves) at the center is the position and shape of the sheet. All streamlines start from a vertical plane (parallel to the initial position of the sheet) near the left boundary of the computational domain (using 25 uniformly spaced seeds). Notice the colors of the curves denote the velocity magnitude. The twisted curves around the sheet indicates the complexity of flow patterns in the vicinity of the structure. Case III) An elastic sheet with aspect ratio of 1:4 (width versus length) is initially put on a horizontal plane of the x and y axes in the middle of the box in all three directions (x, y, and z). Its right edge is rotated constantly and periodically with a period P (on the y − z plane anticlockwise) in a still FENE-P fluid. Again, the structure is not restricted elsewhere and is allowed to move freely in other directions. Some typical results from different combination of values of Re,K b , W i are displayed here. In Fig. 3 the left panel shows the shape/position of the sheet at a few equally distributed instants from a simulation with Re = 10,Kb = 0.005, W i = 1.0. The top most is the initial position (horizontal and flat). The remaining is the shape (not physical position) of the sheet at several time instants within a period. Note that the sheet vertical position is shifted down by a constant for the same purpose as in case II. We see that as the right edge is being rotated, the rest of the sheet follows the rotation motion. Due to flexibility (instead of being rigid), the leading edge (initially straight) deforms into a curve and the entire sheet deforms into a helical structure. As time goes by, more helical structures appear and they appear to move downstream (from righ to left) in an animation. Interestingly the entire sheet moves forward towards right boundary slowly. The right panel in Fig. 3 shows streamlines of the flow field from a simulation with Re = 10,Kb = 0.005, W i = 1.0. The gray surface (partially blocked by the curves) is the position and shape of the twisted sheet. All streamlines start from a horizontal plane parallel and close to the initial position of the sheet (20 uniformly spaced seeds). The colors of the curves denote the velocity magnitude. The tightly wound curves around the sheet indicate a rotating flow near the sheet. The existing immersed boundary (IB) framework has been extended in three dimensions to FFSI problems involving non-Newtonian fluids. The fluids may be power-law, Oldroyd-B, or FENE-P. The viscous incompressible Navier-Stokes equations for the flow and the constitutive equations for the fluid (Oldroyd-B and FENE-P) are simultaneously solved with the lattice Boltzmann approaches by the D3Q19 and the D3Q7 models, respectively. The power-law is incorporated algebraically into the lattice Boltzmann flow model. The new method is tested on three FFSI toy problems: deformable sheets interacting with power-law, Oldroyd-B, and FENE-P fluids in three dimensions. Our results show that the new IB method is first-order accurate and conditionally stable.