Zhe Li
email: zhe.li@ec-nantes.fr

Julien Favier

Umberto D'ortona

Sébastien Poncet

An immersed boundary-lattice Boltzmann method for single-and multi-component fluid flows

Keywords: Immersed boundary method, Lattice Boltzmann method, Multi-component fluid flow, Wetting conditions

The paper presents a numerical method to simulate single-and multi-component fluid flows around moving/deformable solid boundaries, based on the coupling of Immersed Boundary (IB) and Lattice Boltzmann (LB) methods. The fluid domain is simulated with LB method using the single relaxation time BGK model, in which an interparticle potential model is applied for multi-component fluid flows. The IB-related force is directly calculated with the interpolated definition of the fluid macroscopic velocity on the Lagrangian points that define the immersed solid boundary. The present IB-LB method can better ensure the no-slip solid boundary condition, thanks to an improved spreading operator. The proposed method is validated through several 2D/3D single-and multicomponent fluid test cases with a particular emphasis on wetting conditions on solid wall. Finally, a 3D two-fluid application case is given to show the feasibility of modeling the fluid transport via a cluster of beating cilia.

Introduction

The numerical simulation of phenomena involving moving/deformable boundaries in multi-component flows is attracting more and more the interest of physicists and engineers, especially in recent biofluidic applications. In this context, the scope of this work is to develop a numerical tool to simulate the beating of epithelial cilia to transport the mucus in airways [START_REF] Ikegami | Tubulin polyglutamylation is essential for airway ciliary function through the regulation of beating asymmetry[END_REF]. Epithelial cells have whiplike appendages extending from their surface, designed to move the surrounding fluid. The numerical challenges thus include the simulation of slender flexible structures in large deformations, in a two-component fluid flow environment: the periciliary fluid whose properties are close to those of water and the mucus [START_REF] Button | A Periciliary Brush Promotes the Lung Health by Separating the Mucus Layer from Airway Epithelia[END_REF]. Note that although this work mainly aims at simulating the transport of mucus by ciliary motion, a wide spectrum of applications can be found in biomechanics (inner ear, sperm cells), aerodynamics (control of boundary layers around airfoils [START_REF] Favier | Passive separation control using a self-adaptive hairy coating[END_REF]), or animal propulsion studies [START_REF] Fadlun | Combined immersed-boundary finite-difference methods for three-dimensional complex flow simulations[END_REF].

In the present work, the numerical simulation involving mucus and periciliary fluid is tackled via an interparticle potential Multi-Component Multi-Phase (MCMP) Lattice Boltzmann (LB) model [START_REF] Porter | Multicomponent interparticle-potential lattice Boltzmann model for fluids with large viscosity ratios[END_REF]. Based on the original Shan-Chen's model [START_REF] Shan | Lattice Boltzmann model for simulating flows with multiple phases and components[END_REF][START_REF] Shan | Multicomponent Lattice-Boltzmann Model with Interparticle Interaction[END_REF], this MCMP scheme applies an Explicit-Forcing (EF) scheme [START_REF] He | A novel thermal model for the lattice Boltzmann method in incompressible limit[END_REF] in order to take into account external force effects, instead of modifying the equilibrium velocity as in [START_REF] Shan | Lattice Boltzmann model for simulating flows with multiple phases and components[END_REF][START_REF] Shan | Multicomponent Lattice-Boltzmann Model with Interparticle Interaction[END_REF]. As shown in the work of Porter et al. [START_REF] Porter | Multicomponent interparticle-potential lattice Boltzmann model for fluids with large viscosity ratios[END_REF], the EF interparticle potential LB model adopted in the present paper can allow one to reduce the magnitude of spurious currents near the fluid-fluid interface, as well as to increase the density or viscosity ratio, compared to the original SC model. More details about the thermodynamic consistency of different forcing schemes can be found in the works of Huang et al. [START_REF] Huang | Forcing term in single-phase and Shan-Chen-type multiphase lattice Boltzmann models[END_REF] and Li et al. [START_REF] Li | Forcing scheme in pseudopotential lattice Boltzmann model for multiphase flows[END_REF].

As an alternative to conventional numerical methods, LB method has been widely applied for simulating fluid flow problems with diverse physical phenomena, because of its inherent advantages, such as the microkinetic level for dealing with the fluid properties and the simplicity for parallel implementation, etc. Besides, one important feature of LB method is that it can fully recover the Navier-Stokes (NS) equations at the macroscopic scale [START_REF] Chen | Recovery of the Navier-Stokes equations using a lattice-gas Boltzmann method[END_REF][START_REF] Qian | Lattice BGK Models for Navier-Stokes Equation[END_REF].

To incorporate moving solid boundaries in the fluid domain, we choose to couple the LB method with an Immersed Boundary (IB) method. Originally introduced by Peskin [START_REF] Peskin | Flow patterns around heart valves: a numerical method[END_REF] for simulating blood flow in the heart, IB method is now widely used and several variants have been proposed within the framework of NS solvers [START_REF] Fadlun | Combined immersed-boundary finite-difference methods for three-dimensional complex flow simulations[END_REF][START_REF] Zhu | Simulation of a flapping flexible filament in a flowing soap film by the immersed boundary method[END_REF][START_REF] Uhlmann | An immersed boundary method with direct forcing for the simulation of particulate flows[END_REF][START_REF] Pinelli | Immersed-boundary methods for general finite-difference and finite-volume Navier-Stokes solvers[END_REF]. Using IB method allows one to easily introduce moving solid boundaries in fluid flows simulated on a fixed Cartesian grid, which is the case for the present LB solver.

When considering body force term in the fluid domain, one should redefine or correct the macroscopic velocity in the framework of LB method for both single-and multi-component fluid cases in order to fully recover NS equations [START_REF] Shan | Multicomponent Lattice-Boltzmann Model with Interparticle Interaction[END_REF][START_REF] Guo | Discrete lattice effects on the forcing term in the lattice Boltzmann method[END_REF]. In such case, the macroscopic fluid velocity is calculated as the sum of two terms: one is related to the distribution functions of LB model, and the other is the body force-related term. The main idea of IB method is to calculate an appropriate body force-related term so that the fluid can have a desired macroscopic velocity on the Lagrangian points defining the solid wall. Hence, a straightforward method is to first calculate the IB-related force at these Lagrangian points, based on the definition of the macroscopic fluid velocity, and then distribute it onto the neighboring Eulerian nodes with the spreading operator, following Peskin [START_REF] Peskin | Flow patterns around heart valves: a numerical method[END_REF][START_REF] Peskin | The immersed boundary method[END_REF]. Such kind of IB-LB method, like the one of Chen et al. [START_REF] Chen | Immersed boundary method based lattice Boltzmann method to simulate 2D and 3D complex geometry flows[END_REF], can be categorized as explicit IB-LB algorithm, which is now known to have issues in ensuring the no-slip solid boundary condition [START_REF] Wu | Implicit velocity correction-based immersed boundarylattice Boltzmann method and its applications[END_REF]. To overcome this drawback, Wu and Shu [START_REF] Wu | Implicit velocity correction-based immersed boundarylattice Boltzmann method and its applications[END_REF] proposed an implicit velocitycorrection IB-LB method, which relies on resolving a system of equations at each time-step in order to calculate the IB-related force that ensures the no-slip boundary condition.

In the present paper, an IB-LB method is proposed to simulate single-and multi-component fluid flows in the presence of fixed or moving solid boundaries. The macroscopic fluid velocity is split into two parts, and the IB-related force is directly obtained from the definition of the macroscopic fluid velocity at the Lagrangian points. However, we follow the method proposed by Pinelli et al. [START_REF] Pinelli | Immersed-boundary methods for general finite-difference and finite-volume Navier-Stokes solvers[END_REF] to amend the spreading operator in order to improve the reciprocity of interpolation-spreading operations. As a consequence, with such spreading operator, the no-slip solid boundary condition can be better ensured, comparing to the ordinary explicit IB-LB coupling approach. Similar to the method proposed by Wu and Shu [START_REF] Wu | Implicit velocity correction-based immersed boundarylattice Boltzmann method and its applications[END_REF], the resolution of a linear system of equations is also required, but only when the positions of the Lagrangian points change with time. It means that in the case of static solid boundary cases the resolution of the system will be carried out only once at the beginning of the simulation, which turns out to be an advantage of the present method. In addition, it may be worth noting that the present IB-LB method appears to be formally faster than the one proposed by Favier et al. [START_REF] Favier | A Lattice Boltzmann -Immersed Boundary method to simulate the fluid interaction with moving and slender flexible objects[END_REF], since only one LB calculation is required within one time-step, avoiding the need of an extra prediction step.

Besides, another novelty of the present work lies in the extension of the proposed IB-LB method to multi-component fluid cases. The key point is to calculate the IB-related force by means of the definition of the macroscopic velocity for each fluid component. Once these forces are obtained, the spreading procedure is the same as in the single-component fluid case. Furthermore, with the proposed IB-LB coupling method one can take into account different wetting properties of solid wall, by adding the IB-related force after considering the fluidsolid adhesion force model proposed by Martys and Chen [START_REF] Martys | Simulation of multicomponent fluids in complex three-dimensional geometries by the lattice Boltzmann method[END_REF]. For the sake of clarity, we summarize the advantages of the proposed IB-LB coupling method as follows:

• Moving and curved solid boundary condition can be easily incorporated into the fluid simulation. The extension to multi-component cases is straightforward and different wetting conditions can be imposed on curved walls, which is relatively difficult to handle using the classical bounce-back method

• IB-related force is incorporated by means of the definition of macroscopic velocity, in such a way that no extra prediction sub-step is required at each time-step

• The no-slip solid boundary condition can be enhanced with an improved spreading operator. For static IB condition, the resolution of a linear system of equations is required only once at the beginning of numerical simulation, because of the explicit feature of the proposed method To our knowledge, the only attempt to combine multi-phase flow and immersed boundary was done by Shao et al. [START_REF] Shao | Development of an immersed boundaryphase field-lattice Boltzmann method for Neumann boundary condition to study contact line dynamics[END_REF]. One of the interest of this article is the development of a method to implement Neumann boundary condition in the frame of immersed boundary method. Then they apply this boundary condition to the contact line of a droplet on a wall. We successfully reproduce the contact line dynamics, as Shao et al. [START_REF] Shao | Development of an immersed boundaryphase field-lattice Boltzmann method for Neumann boundary condition to study contact line dynamics[END_REF] did, but using a different interparticle potential MCMP model, a different IB method implementation, and more importantly, Dirichlet boundary conditions appear to be sufficient in our numerical framework.

The article is organized as follows. Section 2 first gives the basic equations of the LB method for single-and multi-component models. Then the main idea of the proposed IB-LB method is shown in the second part of this section, followed by a brief introduction of the approach to define the spreading operator. In Section 3, several 2D/3D numerical test cases are shown to validate the proposed IB-LB method. A 3D two-fluid case on the fluid transport by a cluster of beating cilia is presented in Section 4 in order to show the feasibility of simulating mucus flows in airways. Finally, the conclusions are drawn in Section 5.

Mathematical formulations and numerical methods

2.1. Lattice Boltzmann method for multi-component fluid flow 2.1.1. Basic equations of lattice Boltzmann method (single-component)

In LB method, the fluid status is calculated or updated in time by resolving the discrete LB equation given by [START_REF] Shan | Kinetic theory representation of hydrodynamics: a way beyond the Navier-Stokes equation[END_REF]:

f i (x + e i ∆t, t + ∆t) = f i (x, t) - ∆t τ f i (x, t) -f (eq) i (x, t) (1)
where f i (x, t) denotes the distribution function at the site x and the time t, in the i th direction of the used D2Q9 lattice for 2D cases or D3Q19 lattice for 3D cases, as shown in Figure 1, and e i is the i th discrete velocity vector. For example, the set of discrete velocity vectors of D2Q9 lattice can be obtained in the following way:

        e 0 . . . e i . . . e 8         = ∆x ∆t 0 +1 0 -1 0 +1 -1 -1 +1 0 0 +1 0 -1 +1 +1 -1 -1 ⊤ e x e y (2)
where ∆x is the lattice spacing and ∆t is the time step, e x and e y are two unit vectors in x-and y-directions, respectively. Equation (1) allows one to calculate the distribution function f i (x, t+ ∆t) at the next time step, so that one can update the macroscopic state of the singlecomponent fluid: the density ρ (mass per unit volume) and the velocity vector u:

       ρ = i f i ρu = i e i f i (3)
In addition, the single relaxation time Bhatnagar-Gross-Krook (BGK) [START_REF] Bhatnagar | A model for collision processes in gases. I: small amplitude processes in charged and neutral one-component system[END_REF] collision model is adopted in Equation (1). τ is the relaxation time, which is related to the kinetic viscosity ν by ν = c 2 s (τ -0.5∆t) with c s = 1/ √ 3 being the speed of sound. In the present work, the isothermal LB model is applied, which gives the fluid pressure as p = c 2 s ρ. Here we adopt the classical normalization so that ∆x = 1, ∆t = 1 in the LB calculation. In Equation (1), f (eq) i (x, t) is called the discrete equilibrium distribution function which can be obtained by Hermite series expansion of the Maxwell-Boltzmann equilibrium distribution [START_REF] Shan | Kinetic theory representation of hydrodynamics: a way beyond the Navier-Stokes equation[END_REF]:

f (eq) i = ρω i 1 + e i • u c 2 s + (e i • u) 2 2c 4 s - u 2 2c 2 s ( 4
)
where the weight coefficients ω i are ω 0 = 4/9, ω 1-4 = 1/9 and ω 5-8 = 1/36 for D2Q9 lattice, and ω 0 = 1/3, ω 1-6 = 1/18 and ω 7-18 = 1/36 for D3Q19 lattice [START_REF] Qian | Lattice BGK Models for Navier-Stokes Equation[END_REF]. When considering the body force effects in the fluid system, one has to modify the discrete LB equation (1) by adding a body force term [START_REF] Guo | Discrete lattice effects on the forcing term in the lattice Boltzmann method[END_REF]:

f i (x + e i ∆t, t + ∆t) = f i (x, t) - ∆t τ f i (x, t) -f (eq) i (x, t) + ∆tF i (x, t) (5)
where the discrete force term F i is given as:

F i = 1 - ∆t 2τ ω i e i -u c 2 s + e i • u c 4 s e i • F (6)
in which F denotes the body force per unit volume. The definition of the macroscopic velocity has to be accordingly changed so that the numerical model can recover the Navier-Stokes equations [START_REF] Guo | Discrete lattice effects on the forcing term in the lattice Boltzmann method[END_REF]:

ρu = i e i f i + ∆t 2 F (7)
It is worth noting that Equation [START_REF] Shan | Multicomponent Lattice-Boltzmann Model with Interparticle Interaction[END_REF] offers a relation between the macroscopic fluid velocity and the body force at any time. The proposed IB-LB method is based on calculating the force related to the immersed body with the help of relation [START_REF] Shan | Multicomponent Lattice-Boltzmann Model with Interparticle Interaction[END_REF]. If no other body force is considered, then we have F = F IB . More details about the IB-LB coupling algorithm will be given subsequently.

Explicit forcing interparticle potential multi-component LB model

In the adopted multi-component LB model [START_REF] Porter | Multicomponent interparticle-potential lattice Boltzmann model for fluids with large viscosity ratios[END_REF], each fluid component obeys the discrete LB equation written as:

f σ i (x + e i ∆t, t + ∆t) = f σ i (x, t) - ∆t τ σ f σ i (x, t) -f σ(eq) i (x, t) + ∆tF σ i (x, t) (8)
where f σ i (x, t) denotes the distribution function for the σ th component, and τ σ is the relaxation time of the single relaxation time BGK model. The equilibrium distribution function f σ(eq) i is calculated as:

f σ(eq) i = ρ σ ω i 1 + e i • u eq σ c 2 s + (e i • u eq σ ) 2 2c 4 s - u eq σ • u eq σ 2c 2 s ( 9
)
in which ρ σ denotes the density of the σ th component. The equilibrium velocity u eq σ is the same for each component and equal to the common velocity u ′ designed to ensure the conservation of total momentum [START_REF] Porter | Multicomponent interparticle-potential lattice Boltzmann model for fluids with large viscosity ratios[END_REF]:

u eq σ = u ′ = σ i e i f σ i /τ σ σ i f σ i /τ σ (10)
Furthermore, in Equation ( 8), F σ i is the explicit forcing term related to the total body force F σ per unit volume exerted on the σ th component, which is calculated by [START_REF] Porter | Multicomponent interparticle-potential lattice Boltzmann model for fluids with large viscosity ratios[END_REF]:

F σ i = 1 - ∆t τ σ F σ • (e i -u eq σ ) ρ σ c 2 s f σ(eq) i (11)
Based on the methodology of Martys and Chen [START_REF] Martys | Simulation of multicomponent fluids in complex three-dimensional geometries by the lattice Boltzmann method[END_REF], we propose an IB-LB coupling algorithm which relies on adding an extra IB-related force in the total body force vector F σ in Equation ( 11), in such a way that, at the sites near the solid wall, F σ will include, if no other body forces are considered, the Shan-Chen-type fluid-fluid cohesion force F SC σ , the solid-fluid adhesion force F AD σ and the IB-related force

F IB σ , i.e. F σ = F SC σ + F AD σ + F IB σ .
As given in [START_REF] Martys | Simulation of multicomponent fluids in complex three-dimensional geometries by the lattice Boltzmann method[END_REF] and [START_REF] Huang | Proposed approximation for contact angles in Shan-and-Chen-type multicomponent multiphase lattice Boltzmann models[END_REF], the fluid-fluid cohesion force acting on the σ th fluid can be calculated by:

F SC σ (x, t) = -G coh ρ σ (x, t) i ω i ρ σ ′ (x + e i ∆t, t)e i ( 12
)
where G coh is the parameter that controls the strength of the cohesion force, σ ′ represents another fluid component different from σ, and i denotes the indice of lattice velocities for the local site x. Notice that Equation ( 12) takes into account only the i th direction whose neighbouring node x + e i ∆t is located inside the fluid domain.

Similarly, the solid-fluid adhesion force acting on the σ th fluid is calculated by:

F AD σ (x, t) = -G σ adh ρ σ (x, t) i ω i s(x + e i ∆t, t)e i (13)
where s(x + e i ∆t, t) represents an indicator function that equals 1 or 0 for a solid or a fluid neighbouring node, respectively. By changing the parameter G σ adh , one can adjust the strength of the interaction force between the solid wall and each fluid component. As presented in [START_REF] Martys | Simulation of multicomponent fluids in complex three-dimensional geometries by the lattice Boltzmann method[END_REF] and [START_REF] Huang | Proposed approximation for contact angles in Shan-and-Chen-type multicomponent multiphase lattice Boltzmann models[END_REF], for the solid wall, a positive G σ adh simulates a non-wetting fluid and G σ adh should be negative for a wetting fluid.

At last, the IB-related force F IB σ will be calculated with the proposed IB-LB method so that the macroscopic velocity of each fluid component can have the desired value on the Lagrangian points where the immersed boundary condition is imposed. The definition of the macroscopic velocity u σ of each component is given as follows [START_REF] Porter | Multicomponent interparticle-potential lattice Boltzmann model for fluids with large viscosity ratios[END_REF]:

ρ σ u σ = i e i f σ i + ∆t 2 F σ (14)
where ρ σ = i f σ i denotes the density of the σ th fluid component.

The proposed immersed boundary method

In the present work, the Immersed Boundary (IB) method [START_REF] Peskin | Flow patterns around heart valves: a numerical method[END_REF][START_REF] Peskin | The immersed boundary method[END_REF] is applied to impose solid boundary conditions in single-and multi-component fluid flows simulated by LB method.

Mostly, the idea of IB method is to add an extra body force F IB in the fluid equations so that the macroscopic fluid velocity can be approximately equal to the velocity of the immersed solid boundary. In addition, a fixed Eulerian mesh is often adopted for the fluid domain, whereas the solid boundary is usually presented by a set of Lagrangian points, which can move in space, as shown in Figure 2.

Like the conventional IB method, the present approach proceeds in two main steps:

• Interpolation -Interpolate the fluid information from the Eulerian nodes to the Lagrangian points and then calculate the IB-related body force at the Lagrangian points;

• Spreading -Spread the obtained force from the Lagrangian points to the neighbouring Eulerian nodes and accomplish the resolution of the fluid equations at the Eulerian nodes with the IB-related body force.

Eulerian node

Lagrangian point

Immersed boundary

Kernel support

Interpolation stage

As presented previously, in single-and multi-component LB models, the definition of macroscopic velocity, Equations ( 7) and ( 14), provide a relation between the fluid velocity and the body force. The key point of the proposed method is to calculate the IB-related body force vector by ensuring this definition or relation at the Lagrangian points, as shown in Figure 2.

For instance in single-component case, at the Lagrangian points (see Figure 2), the velocity-force relation can be written as:

I [ρ] s u s = I i e i f i s + ∆t 2 F s ( 15
)
where u s denotes the fluid velocity at the s th Lagrangian point, which equals the velocity of solid boundary in order to satisfy the no-slip condition, e.g. for a stationary solid boundary, we have u s = 0. In addition, F s is the body force per unit volume at this Lagrangian point, and if no other body force is considered, F s = F IB s . Finally, I [•] s represents the interpolation operator defined as:

φ(x s , t) = I [φ(x, t)] s = φ(x, t)δ(x -x s ) dx ≃ j∈Ds φ(x j , t) δ(x j -x s )(∆x) d (16)
which provides the interpolated value of a certain variable φ(x, t) at the Lagrangian point x s = x(X s , t) with X s being its time-independent Lagrangian label. Notice that this interpolation operator makes use of the sampling property of the Dirac delta function δ(x). In Equation ( 16), φ(x j , t) is the variable's value at the j th Eulerian node located inside the support domain D s of the Lagrangian point x s , ∆x is the spacing of the used uniform lattice and d denotes the number of dimensions of the investigated problem. Finally δ(x jx s ) is a mollifier or a smooth approximation to the Dirac delta function. In the present work, we adopt the mollifier proposed by Roma et al. [START_REF] Roma | An adaptive version of the immersed boundary method[END_REF], for instance in 2D case:

δ(x j -x s ) = 1 ∆x δx ( |x j -x s | ∆x ) 1 ∆y δy ( |y j -y s | ∆y ) (17)
with:

δx (r) = δy (r) =            1 3 1 + -3r 2 + 1 0 ≤ r ≤ 0.5 1 6 5 -3r --3(1 -r) 2 + 1 0.5 ≤ r ≤ 1.5 0 otherwise ( 18
)
where r denotes the non-dimensional distance between the Lagrangian point and the Eulerian node in each direction, i.e. r = |x s -x j |/∆x or r = |y s -y j |/∆y. Consequently, we can calculate the IB-related body force F IB s in Equation ( 15), since u s is equal to the velocity of solid boundary when the no-slip condition is applied. Now, we can apply a similar procedure for multi-component fluid cases. Let's consider the macroscopic velocity definition [START_REF] Zhu | Simulation of a flapping flexible filament in a flowing soap film by the immersed boundary method[END_REF] for the σ th component at the s th Lagrangian point:

I [ρ σ ] s u σ,s = I i e i f σ i s + I ∆t 2 (F SC σ + F AD σ ) s + ∆t 2 F IB σ,s (19)
where u σ,s and F IB σ,s denote the macroscopic velocity of the σ th fluid component and the IB-related force at the s th Lagrangian point, respectively.

At this interpolation stage of the proposed algorithm (see Sub-Section 2.3, Algorithm 2, Step 7) for calculating the IB-related force F IB σ,s at the instant t = t n+1 , all the distribution functions have already been updated to t n+1 , hence we have the density, fluid-fluid cohesion force, and fluid-solid adhesion force fields at t n+1 , so that we can calculate the three interpolation terms in Equation [START_REF] Chen | Immersed boundary method based lattice Boltzmann method to simulate 2D and 3D complex geometry flows[END_REF]. Meanwhile, u σ,s has the appropriate value required to impose the no-slip condition. Here we impose that all fluid components have the same velocity at the Lagrangian points, with which we can calculate the individual IB-related force for each component by means of Equation [START_REF] Chen | Immersed boundary method based lattice Boltzmann method to simulate 2D and 3D complex geometry flows[END_REF].

Spreading stage

Once the IB-related body force F IB σ,s is obtained, we can spread this force from the Lagrangian points onto the j th Eulerian node by means of the classical spreading operator S[•] j proposed by Peskin [START_REF] Peskin | The immersed boundary method[END_REF], for instance in the 2D multicomponent case: [START_REF] Wu | Implicit velocity correction-based immersed boundarylattice Boltzmann method and its applications[END_REF] where ∆s denotes the arc length of the s th segment of the immersed boundary Γ IB , of which the center is located at x s , inside the support domain D j of the j th Eulerian node. Notice that we apply the notations F IB σ,s = F IB σ (x s , t) and δjs = δ(x jx s ) for the sake of simplicity.

F IB σ (x j , t) = S F IB σ,s j = Γ IB F IB σ (x s , t)δ(x j -x s ) ds ≃ s∈Dj F IB σ,s δjs ∆s
Pinelli et al. [START_REF] Pinelli | Immersed-boundary methods for general finite-difference and finite-volume Navier-Stokes solvers[END_REF] proposed to modify Equation ( 20) by introducing a characteristic strip width ǫ s = ǫ(x s , t) related to the local dilation coefficients of the mollifier δjs = δ(x jx s ), in such a way that the new spreading operator S * [•] j gets:

F IB σ (x j , t) = S * F IB σ,s j ≃ s∈Dj F IB σ,s δjs ǫ s ∆s (21)
If the immersed boundary Γ IB is discretized into N s segments, N s equations will be required to determine the correct value of the local strip width ǫ s . More details about the procedure to obtain ǫ s can be found in the work of Pinelli et al. [START_REF] Pinelli | Immersed-boundary methods for general finite-difference and finite-volume Navier-Stokes solvers[END_REF] and Favier et al. [START_REF] Favier | A Lattice Boltzmann -Immersed Boundary method to simulate the fluid interaction with moving and slender flexible objects[END_REF]. Briefly, the main idea of the procedure is to calculate ǫ s such that the interpolation-spreading operations are reciprocal, i.e. φ(x s , t) = I S * [φ(x s , t)] j s

, which provides one equation for the s th segment:

φ s = j∈Ds δjs (∆x) d s ′ ∈Dj φ s ′ δjs ′ ǫ s ′ ∆s ′ = j∈Ds δjs (∆x) d N s ′ s ′ φ s ′ δjs ′ ǫ s ′ ∆s ′ (22)
where δjs ′ = δ(x jx s ′ ) and s ′ is another solid segment indice different from s. Regrouping all the N s equations provides a system of equations in matrix form:

Φ = Adiag(ǫ)Φ (23)
where

Φ = [φ 1 , φ 2 , • • • , φ Ns ] ⊤ and ǫ = [ǫ 1 , ǫ 2 , • • • , ǫ Ns ] ⊤
. The element of matrix A at the k th row and the l th column is:

A(k, l) = j∈D k (∆x) d δkj δlj ∆s l (24)
where k ∈ [1, N s ] and l ∈ [1, N s ] with N s being the total number of solid segments. Finally, as presented in [START_REF] Pinelli | Immersed-boundary methods for general finite-difference and finite-volume Navier-Stokes solvers[END_REF] and [START_REF] Favier | A Lattice Boltzmann -Immersed Boundary method to simulate the fluid interaction with moving and slender flexible objects[END_REF], the value of each ǫ s can be obtained by resolving the system of equations:

Aǫ = 1 (25)
where

1 = [1, 1, • • • , 1] ⊤ .
Notice that such spreading procedure is the same for single-and multicomponent cases. Hence, once the IB-related body force F IB s or F IB σ,s is obtained at the Lagrangian points, we will spread the force onto the Eulerian nodes with this spreading procedure proposed by Pinelli et al. [START_REF] Pinelli | Immersed-boundary methods for general finite-difference and finite-volume Navier-Stokes solvers[END_REF].

Overview of the proposed IB-LB coupling algorithm

In this sub-section, an overview of the proposed IB-LB coupling algorithm is given to show step by step how to couple the immersed boundary and the lattice Boltzmann methods for the numerical simulation of fluid flows in the presence of fixed or moving solid boundaries. Although the coupling procedures are very similar for single-and multi-component cases, two separate algorithms are provided for the sake of clarity.

First of all, let's suppose that all the variables are already known at the current moment t, such as the densities ρ(x, t) or ρ σ (x, t), the fluid velocity vector u(x, t), all the distribution functions f i (x, t) or f σ i (x, t) and the force vectors F(x, t) or F σ (x, t). Hence the objective is to update the fluid status from the moment t to the next time step t + ∆t.

Algorithms 1 and 2 show the time advancing procedures for single-and multi-component fluid cases, respectively.

Algorithm 1: IB-LB coupling algorithm for single-component fluid cases

Require: ρ(x, t), u(x, t) and F(x, t) from the previous time-step 1: Calculate f (eq) i (x, t) with Eq. (4) 2: Calculate F i (x, t) with Eq. ( 6) 3: Collision step of LB method with Eq. ( 5) 4: Streaming step of LB method to obtain the distribution functions f i (x, t + ∆t) at the next time-step 5: Calculate the fluid mass density ρ(x, t + ∆t) with Eq. (3) 6: Interpolation stage (15) of IB method to obtain F IB (x s , t + ∆t) at the Lagrangian points with the imposed solid velocity u s (x s , t + ∆t) 7: Spreading stage [START_REF] Favier | A Lattice Boltzmann -Immersed Boundary method to simulate the fluid interaction with moving and slender flexible objects[END_REF] of IB method to obtain F IB (x j , t + ∆t) at the Eulerian nodes 8: Calculate the fluid macroscopic velocity u(x, t + ∆t) with Eq. ( 7) 9: Update the fluid status to the moment t + ∆t and go to Step 1 for the next time-step

Validation of the proposed IB-LB method

3.1. Stationary undeformable object: the wake behind a stationary circular cylinder at Re = 30 The capability of the present IB-LB coupling method is now evaluated by simulating the wake behind a stationary circular cylinder of diameter D. The Reynolds number defined by Re = U 0 D/ν, with U 0 the inlet velocity, is fixed at Re = 30, such that the wake is steady and symmetric with two large recirculation bubbles as shown in Figure 3.

Notice that Figure 3 presents only a rectangular zone around the fixed cylinder. The computational domain Ω f is discretized with a 1000 × 500 uniform lattice. The center of the cylinder is placed at the centerline that is parallel Algorithm 2: IB-LB coupling algorithm for multi-component fluid cases Require: ρ σ (x, t), f σ i (x, t) and F σ (x, t) from the previous time-step 1: Calculate the equilibrium velocity u eq σ (x, t) with Eq. ( 10) 2: Calculate f σ(eq) i (x, t) with Eq. ( 9) 3: Collision step of LB method with Eq. ( 8) 4: Streaming step of LB method to obtain f σ i (x, t + ∆t) 5: Calculate ρ σ (x, t + ∆t) = i f σ i (x, t + ∆t) for each component 6: Calculate F SC σ (x, t + ∆t) with Eq. ( 12) and F AD σ (x, t + ∆t) with Eq. ( 13) using the newly updated density fields 7: Interpolation stage of IB method [START_REF] Chen | Immersed boundary method based lattice Boltzmann method to simulate 2D and 3D complex geometry flows[END_REF] to obtain F IB σ (x s , t + ∆t) at the Lagrangian points 8: Spreading stage of IB method ( 21) to obtain F IB σ (x j , t + ∆t) at the Eulerian nodes 9: Calculate the total body force F σ (x, t + ∆t) for the σ th component 10: Calculate the macroscopic velocity u σ (x, t + ∆t) with Eq. ( 14) 11: Update the fluid status of each component to the moment t + ∆t and go to Step 1 for the next time-step to x-direction, and at a distance of 8D from the inlet with D = 32∆x. The cylinder is discretized with a certain number of IB segments, of which the length is equal to the uniform lattice spacing, which means ∆s ≃ ∆x = ∆y. The inlet velocity and the outlet density or pressure conditions are imposed by means of the method presented by Zou and He [START_REF] Zou | On pressure and velocity boundary conditions for the lattice Boltzmann BGK model[END_REF]. At the top and bottom boundaries, Neumann boundary conditions (zero-gradient condition) [START_REF] Lou | Evaluation of outflow boundary conditions for two-phase lattice Boltzmann equation[END_REF] are imposed in order to have a free-slip condition as in [START_REF] Pinelli | Immersed-boundary methods for general finite-difference and finite-volume Navier-Stokes solvers[END_REF]. Besides, the initial fluid field is given as:

∀(x, y) ∈ Ω f      u x (x, y, 0) = U 0 u y (x, y, 0) = 0 ρ(x, y, 0) = 1 (26)
The present results are compared to previous experimental and numerical data available in the literature, as shown in Table 1. Lengths a and b are the distance between the back of the cylinder and the vortex core and the distance between the two vortex cores respectively, as shown in Figure 3. It is known to be a very sensitive test case for numerical methods. The dimensionless lengths a/D and b/D are particularly well predicted, in good agreement with the experimental results of Coutanceau and Bouard [START_REF] Coutanceau | Experimental determination of the main features of the viscous flow in the wake of a circular cylinder in uniform translation. Part 1. Steady flow[END_REF].

Furthermore, the drag and lift coefficients are calculated with:

       C D = F f →s • e x 0.5ρU 2 0 D C L = F f →s • e y 0.5ρU 2 0 D (27)
where e x and e y are two unit vectors in x-and y-directions, respectively, the fluid density ρ = 1, and F f →s denotes the fluid-induced force vector exerted on the solid cylinder, which can be calculated by:

F f →s = - s∈Γ IB F IB s ǫ s ∆s ( 28
)
where Γ IB denotes the immersed circular boundary, F IB s is the IB body force vector per unit area at the s th segment of Γ IB . ǫ s is the s th solution of Equation [START_REF] Bhatnagar | A model for collision processes in gases. I: small amplitude processes in charged and neutral one-component system[END_REF], which can be considered as the local thickness of the segment [START_REF] Pinelli | Immersed-boundary methods for general finite-difference and finite-volume Navier-Stokes solvers[END_REF][START_REF] Favier | A Lattice Boltzmann -Immersed Boundary method to simulate the fluid interaction with moving and slender flexible objects[END_REF], and finally ∆s is the length of the segment.

In the converged state, the lift coefficient C L is -6.7 × 10 -4 confirming the symmetric nature of the wake, and as presented in Table 1, the drag coefficient C D has a value of 1.79 which agrees with the one obtained by Pinelli et al. [START_REF] Pinelli | Immersed-boundary methods for general finite-difference and finite-volume Navier-Stokes solvers[END_REF] using a Navier-Stokes solver coupled with an IB method.

As presented previously in the introduction, thanks to the spreading operator S * [•] j given by Equation ( 21), the present IB-LB method can better ensure no-slip boundary condition, comparing to the ordinary explicit IB-LB method using the spreading operator S[•] j given by Equation [START_REF] Wu | Implicit velocity correction-based immersed boundarylattice Boltzmann method and its applications[END_REF]. Such improvement is shown in Figure 4 by comparing the streamlines in the numerical simulations carried out with the ordinary explicit and the present IB-LB methods. One can Present method 0.59 0.55 1.79 Pinelli et al. [START_REF] Pinelli | Immersed-boundary methods for general finite-difference and finite-volume Navier-Stokes solvers[END_REF] 0.56 0.52 1.80

Exp.

Coutanceau and Bouard [START_REF] Coutanceau | Experimental determination of the main features of the viscous flow in the wake of a circular cylinder in uniform translation. Part 1. Steady flow[END_REF] 0.54 0.54 -Tritton [START_REF] Tritton | Experiments on the flow past a circular cylinder at low Reynolds numbers[END_REF] --1.74 observe that for the ordinary explicit IB-LB method, the streamlines cross the cylinder, especially around the leading point, which means that a penetration flow occurs at such place. In contrast, for the proposed IB-LB method, since no streamline crosses the IB cylinder, there is no obvious penetration flow. Hence the no-slip solid boundary condition can be better ensured with the present IB-LB method than the ordinary explicit one.

Moving undeformable object

Impulsively started 2D flat plate

We consider now the case of a two-dimensional infinitesimally thin flat and undeformable plate of height h. The plate is suddenly accelerated from rest to a constant velocity U 0 in the direction normal to its surface in a fluid at rest, as shown in Figure 5-a. The Reynolds number is fixed at Re h = U 0 h/ν = 1000 to enable direct comparisons with three previous simulations: Koumoutsakos and Shields [START_REF] Koumoutsakos | Simulations of the viscous flow normal to an impulsively started and uniformly accelerated flat plate[END_REF] using a vortex-particle method, Mittal et al. [START_REF] Mittal | A versatile sharp interface immersed boundary method for incompressible flows with complex boundaries[END_REF] using a directforcing method and recently Pinelli et al. [START_REF] Pinelli | Immersed-boundary methods for general finite-difference and finite-volume Navier-Stokes solvers[END_REF], who developed a Navier-Stokes solver coupled with an IB method.

In the present simulation, the computational domain Ω f of size 15h × 9h is discretized with a 900∆x × 540∆y uniform grid (∆x = ∆y). This 2D plate is discretized with 60 segments (∆s = ∆x) and initially placed in the middle of the domain Ω f of which the four boundaries are all treated with the Bounce-Back rule.

Figure 5-a shows the snapshots of the vorticity magnitude field at three moments: t * = tU 0 /h = 0.5, 1.0 and 2.0. One can observe that, while the plate keeps moving from the left to the right of the domain, two vortices are symmetrically generated near the two ends of the plate. Meanwhile, a recirculation bubble appears and develops behind this 2D moving plate. To validate the proposed method, we plot in Figure 5-b the evolution in time of the computed non-dimensional bubble length s/h that is the ratio between the length of the reverse flow region s measured on the centerline behind the plate and the plate height h. A fairly close agreement is found with the previous results of Koumoutsakos and Shields [START_REF] Koumoutsakos | Simulations of the viscous flow normal to an impulsively started and uniformly accelerated flat plate[END_REF], Mittal et al. [START_REF] Mittal | A versatile sharp interface immersed boundary method for incompressible flows with complex boundaries[END_REF] and Pinelli et al. [START_REF] Pinelli | Immersed-boundary methods for general finite-difference and finite-volume Navier-Stokes solvers[END_REF], which shows the ability of this new forcing algorithm to handle infinitesimally thin moving undeformable bodies.

Oscillating 2D plate

To examine the convergence order of the proposed IB-LB method, we carry out a numerical test case where a 2D infinitely thin plate smoothly oscillates in a square fluid domain (Figure 6). This fluid domain has four periodic boundaries of length L and is initially at rest, i.e. ρ(x, y, 0) = 1 and u x (x, y, 0) = u y (x, y, 0) = 0. The length of the plate l = L/5 and M is its central point. To have a sufficiently smooth fluid field [START_REF] Griffith | On the order of accuracy of the immersed boundary method: Higher order convergence rates for sufficiently smooth problems[END_REF], we impose a translational movement (no rotation) for the plate with x(M ) = A m cos(ωt) + x(O) where x(M ) and x(O) denote the x-direction coordinates of the points M and O, respectively. In addition, the 2D plate oscillates with an amplitude A m = l/2 and a frequency ω = 0.5πU 0 /A m with U 0 being the characteristic velocity determined by the definition of the Reynolds number Re l = U 0 l/ν.

In the present test case, we choose to keep the Reynolds number at Re l = 20 and the kinetic viscosity ν = 0.05 for the LB computation, while changing the lattice resolution. A series of simulations have been carried out with different discretization points for the square's side L: N = 101, 201, 301, 401, 501, 601, 701, 801, 1001, 1601, from which we consider the solution using 1601 points as the reference in order to calculate the L 2 -norm of the numerical error defined as:

L 2 -Error = 1 N 2 101 (u num x -u ref x ) 2 + u num y -u ref y 2 ( 29
)
where N 2 101 = 101 2 represents the total number of the measurement points possessing the same physical positions (at the coarsest lattice nodes) for all the lattice resolutions. Figure 7 shows the decreasing error while increasing the discretization points N 2 . One can observe that the proposed IB-LB coupling scheme tends to be 2 nd order accurate when the physical lattice spacing becomes sufficiently small. This means that the present IB-LB method does not decrease the order of accuracy of LB method. Here one has to note that it is possible for the present method to be 1 st order accurate in some particular cases where the velocity field is not sufficiently smooth, e.g. the first order spatial derivative of the fluid velocity is not continuous across the IB [START_REF] Peskin | The immersed boundary method[END_REF][START_REF] Griffith | On the order of accuracy of the immersed boundary method: Higher order convergence rates for sufficiently smooth problems[END_REF].

Wetting properties for a solid wall

It is possible for a solid wall to get in contact with a wetting or a non-wetting fluid. When such two kinds of fluid come into contact with each other at a solid wall, one will have a typical two fluids-solid wall interaction model in which the contact angle is usually considered as the key parameter to characterize different wetting properties. It is also challenging for numerical method to predict it accurately.

To assess the capability of the present method to handle different wetting conditions at a solid wall, three numerical test cases have been carried out involving two fluids-solid wall interaction model.

2D case: static contact angles under different wetting conditions

Based on the work of Martys and Chen [START_REF] Martys | Simulation of multicomponent fluids in complex three-dimensional geometries by the lattice Boltzmann method[END_REF], Huang et al. [START_REF] Huang | Proposed approximation for contact angles in Shan-and-Chen-type multicomponent multiphase lattice Boltzmann models[END_REF] proposed a theoretical equation for calculating the static contact angle in the two fluidssolid wall interaction model, which is given by:

cos θ 1 = G 2 adh -G 1 adh G coh ρ 1 -ρ 2 2 (30)
where θ 1 denotes the contact angle defined in Figure 8, G 1 adh and G 2 adh are the two adhesion force parameters in Equation ( 13 respectively. In addition, G coh is the fluid-fluid cohesion force parameter of the Shan-Chen-type force model in Equation [START_REF] Qian | Lattice BGK Models for Navier-Stokes Equation[END_REF]. In the droplet, ρ 1 is the equilibrium main density of fluid 1 and ρ 2 is the associated dissolved density of fluid 2 [START_REF] Huang | Proposed approximation for contact angles in Shan-and-Chen-type multicomponent multiphase lattice Boltzmann models[END_REF]. In the following numerical test cases, we choose G coh = 1.8, ρ 1 = 2.05, ρ 2 = 0.05, and G 1 adh = -G 2 adh . The value of G 1 adh or G 2 adh varies according to different wetting conditions.

It is worth mentioning that one should impose the values of G 1 adh and G 2 adh not only for the Eulerian nodes above the IB interface Γ IB but also for the ones under Γ IB . This is because one needs to solve the fluid equations under the immersed boundary Γ IB in the framework of IB method, although Γ IB represents the perimeter of a solid object. In this work, the fluid under Γ IB is set to be the same as the ambient fluid, i.e. we have the fluid 2 under the IB solid wall, as shown in Figure 8. Furthermore, we use G 1 adh = G 2 adh = 0 for the fluid nodes under Γ IB , which corresponds to a neutral wetting condition.

Besides, in this two-dimensional stationary test case, one sets τ 1 = τ 2 = 1, as also used by Huang et al. [START_REF] Huang | Proposed approximation for contact angles in Shan-and-Chen-type multicomponent multiphase lattice Boltzmann models[END_REF]. The whole fluid domain Ω f is a 160∆x × 100∆y rectangle discretized with a uniform (∆x = ∆y) Eulerian lattice. The Bounce-Back rule is applied at the top and the bottom of Ω f in order to ensure no-slip boundary conditions. Meanwhile, the left and right boundaries are treated by periodic conditions. The fluid-solid interface is represented by an immersed boundary line Γ IB as shown in Figure 8, on which 160 Lagrangian points are uniformly distributed. Hence one has ∆s = ∆x = ∆y. The initial form of this droplet of fluid 1 is a semi-circle of diameter D = 50∆x with ρ 1 (x, y, 0) = 2.05 and ρ 2 (x, y, 0) = 0.05 inside the semi-circle surrounded by the ambient fluid 2 where ρ 1 (x, y, 0) = 0.05 and ρ 2 (x, y, 0) = 2.05. Figure 9 shows the results after adh and G 2 adh given in Table 2. To obtain the contact angle from the static density field, we adopt the method used by Huang et al. [START_REF] Huang | Proposed approximation for contact angles in Shan-and-Chen-type multicomponent multiphase lattice Boltzmann models[END_REF], which is based on measuring the height of droplet H and the length of contact area L, as shown in Figure 8. The contact angle θ 1 can be calculated by:

       R = 4H 2 + L 2 8H tan θ 1 = L 2(R -H) (31)
where R represents the estimated radius of the static droplet.

However, when measuring H and L from the density field, one can not exactly determine the contact area, due to the characteristic of the immersed boundary method. Hence, we propose to measure the values of H and L, which are formally written as H = H 0 ± ∆H and L = L 0 ± ∆L where H 0 and L 0 denote the measured values of H and L, respectively, and ∆H and ∆L their estimated errors. In this case, we choose ∆H = ∆L ≃ 2∆x. Figure 10 shows the comparison between the present numerical results and the ones obtained by Equation [START_REF] Coutanceau | Experimental determination of the main features of the viscous flow in the wake of a circular cylinder in uniform translation. Part 1. Steady flow[END_REF] proposed by Huang et al. [START_REF] Huang | Proposed approximation for contact angles in Shan-and-Chen-type multicomponent multiphase lattice Boltzmann models[END_REF]. As one can observe in Figure 10, the present results are in good agreement with the theoretical ones.

2D case: contact line on stationary circular cylinders

This test case is performed in order to show the capability of the proposed method to handle different wetting conditions on curved walls. It was used by Shao et al. [START_REF] Shao | Development of an immersed boundaryphase field-lattice Boltzmann method for Neumann boundary condition to study contact line dynamics[END_REF] to validate their immersed boundary-phase field LB coupling method. In this test case, a circular cylinder of radius 40∆x is fixed at the center of the fluid domain of size 200∆x × 200∆y (∆x = ∆y). Periodic boundary conditions are applied at the left and right boundaries and neutral wetting condition is used at the upper and lower walls. Initially, the fluid-fluid interface is located at the horizontal middle line of fluid domain, for the fluid nodes both inside and outside the circular cylinder. The lower region is set as fluid 1, and the initial densities are the same as in the previous test case. Notice that the solid-fluid adhesion force needs to be zero for the fluid nodes inside the cylinder, i.e. the adhesion coefficient G 1 adh = G 2 adh = 0 for the inner fluid nodes. Whereas for the outer fluid nodes, three combinations of G 1 adh and G 2 adh are used to impose three different wetting conditions, which are 11 shows the numerical results of the proposed method, which appear to be in very good agreement with the desired contact angles θ ref .

G 1 adh = -G 2 adh = -0.45 for θ 1 ref = 60 • , G 1 adh = -G 2 adh = 0 for θ 1 ref = 90 • and G 1 adh = -G 2 adh = 0.45 for θ 1 ref = 120 • . Figure

3D case: breaking of a thin fluid film on a hydrophobic planar wall

The numerical solver has been extended to 3D to simulate the breaking of a thin fluid film, with the same physical parameters as in the previous test cases. At rest, two fluids with the same density and viscosity are stratified in a small bounded box of size L × L × H as shown in Figure 12, which is discretized with a regular (∆x = ∆y = ∆z) lattice with L = 200∆x = 200∆y and H = 24∆z. The top and bottom are modeled as solid walls by means of the bounce-back (BB) rule, and the four lateral surfaces are treated by periodic boundary conditions. The IB-solid wall is a plane surface constituted of 200×200 elemental surfaces of area ∆s x × ∆s y with ∆s x ≃ ∆x and ∆s y ≃ ∆y. Under zero-gravity condition, fluid 1 is confined between fluid 2 and the IB-solid wall with a given height h = 5∆z. The lower region is set as ambient fluid 2 and G 1 adh = G 2 adh = 0 for the fluid nodes under Γ IB . At the time t = 0, a small circular hole of diameter D = 10∆x is set at the center of the fluid 1 film. This academic flow configuration is well detailed in the monograph of De Gennes et al. [START_REF] De Gennes | Gouttes, bulles, perles et ondes[END_REF]. Experimentally, as presented in [START_REF] De Gennes | Gouttes, bulles, perles et ondes[END_REF], when t > 0, it is observed that the diameter D of the hole increases with time t with a constant velocity V D . Besides, this hole of the fluid 1 film is surrounded by a flanged edge that collects the fluid from the center.

In this test case, hydrophobic wetting condition is applied for the IB wall surface with the coefficients of adhesion force G 1 adh = 0.4 and G 2 adh = -0.4. For comparison, a calculation under the same numerical conditions have been performed with a classical bounce-back approach for the bottom wall. The goal is here twofold: to show the capability of the new solver to take into account wetting conditions in the 3D case without using any particular boundary condition for the triple contact line and to link the values of the adhesion coefficients to real wetting conditions for future numerical simulations.

Since the configuration of such test case is symmetric with respect to the central axis in z-direction, the numerical results obtained with the proposed IB-LB coupling method (left) and the classical bounce-back rule (right) are pre- sented in Figure 13 in order to better illustrate the comparison between the two numerical approaches. One can observe that this 3D film breaks progressively as time goes on, because of the existence of the initial hole situated at the center. A good agreement can be found between the results of the IB-LB coupling method and the classical bounce-back rule, in terms of the hole's growing diameter and the shape of its flanged edge.

To estimate the growing diameter for the fluid 1 hole, we provide the evolution in time of the hole's diameter of the present IB-LB method in Figure 14, compared with the one obtained by means of the classical bounce-back rule. As one can observe, the two results are in very good agreement.

Flow induced by a cluster of beating cilia

In the last application case, we apply the proposed IB-LB method to simulate the transport of fluid by a cluster of beating cilia in order to show the feasibility for the present method to handle moving solid boundary condition in a two-fluid environment.

As shown in Figure 15, a cluster of 9 cilia are settled regularly at the center of the bottom of the fluid domain including two immiscible fluids, which are initially stratified in z-direction under a zero-gravity condition. No-slip and free-slip boundary conditions are applied at the bottom and the top surfaces of the fluid domain, respectively. The 4 lateral sides are treated with periodic boundary conditions. The whole system is initially at rest, and the two fluids have the following properties: ρ mucus /ρ water = 1 and ν mucus /ν water = 20. To discretize the fluid domain of size 3l × 3l × 1.6l, we apply a uniform lattice of 60∆x × 60∆y × 32∆z with ∆x = ∆y = ∆z. Each of the 9 cilia is considered as an infinitely thin fiber of length l, which is discretized with 20 segments (∆s ≃ ∆x). All cilia move together in a synchronized way under the beating-motion law governed by a 1D transport equation proposed by Chatelin and Poncet [START_REF] Chatelin | A hybrid grid-particle method for moving bodies in 3D Stokes flow with variable viscosity[END_REF]. As shown in Figure 16, there are two stroke phases for the beating cilium, which are the recovery and effective stroke phases [START_REF] Sanderson | Ciliary activity of cultured rabbit tracheal epithelium: beat pattern and metachrony[END_REF]. Such two kinds of beating phases are imposed with a time-varying coefficient in the 1D transport equation of Chatelin and Poncet [START_REF] Chatelin | A hybrid grid-particle method for moving bodies in 3D Stokes flow with variable viscosity[END_REF]. As the cilia keep beating in the x-z plane, an oscillating fluid flow is induced due to the no-slip solid boundary condition imposed by the proposed IB-LB method. Figure 17 shows the periodically oscillating fluid velocities in x-direction at the three points A, B and C as shown in Figure 15.

Because of the asymmetric beating motion of the cilia, a mean fluid flow is generated in positive x-direction. To verify this, we calculate the mean velocities in x-direction at the A, B and C points, based on the selected period as shown in Figure 17, and finally we obtain the time-averaged velocities: x-direction confirm that the mean flow is in the positive x-direction, which is the direction of the effective stroke phase of the beating cilia.

Conclusions and perspectives

In the present work, an improved direct-forcing Immersed Boundary-Lattice Boltzmann (IB-LB) coupling method has been proposed and validated on several configurations involving moving objects, deforming geometries and multicomponent flows. An important feature of the present numerical framework is that it can be applied for simulating single-and multi-component fluid flows in the presence of fixed/moving solid boundaries. No-slip boundary conditions on a fixed cylinder are better ensured with the present method, comparing to the ordinary explicit one. Meanwhile, thanks to the one-step explicit time integration, an extra prediction step is not required. Moreover, for static solid boundary, the resolution of a linear system of equations is needed only once at the beginning of the simulation. The present method is found to be close to 2 nd order accurate by using a test case where the fluid velocity field is sufficiently smooth. In addition, the proposed method has been adopted in several 2D and 3D multi-component fluid cases, where good agreements have been found with the results published in the literature. Finally, the present numerical solver has been successfully tested on a 3D configuration involving a cluster of beating cilia in a two-fluid flow. This shows the capability of the present solver to simulate the transport of mucus by cilia motion in airways.

Figure 1 :

 1 Figure 1: The two lattices used in the present work: (a) D2Q9 for 2D cases; (b) D3Q19 for 3D cases.

Figure 2 :

 2 Figure 2: A solid boundary immersed in the fluid domain.

Figure 3 :

 3 Figure 3: Stationary wake flow behind a fixed circular cylinder at Re = 30: the velocity magnitude field normalized by the inlet velocity U 0 .

Figure 4 :

 4 Figure 4: Streamlines in the simulation results of the fixed cylinder test case using: (a) the ordinary explicit IB-LB method; (b) the proposed improved IB-LB method.

  y, T ) and u num y = u num y (x, y, T ) are the numerical velocity components in x-and y-directions at the time t = T with T = 2π/ω being the period of the imposed oscillation. Finally, u ref x = u ref x (x, y, T ) and u ref y = u ref y (x, y, T ) represent the reference velocity components obtained by using the most refined lattice N 2 = 1601 2 .

Figure 5 :

 5 Figure 5: Impulsively started flat plate: (a) the vorticity magnitude fields at three different instants (dashed lines are the contour lines for u • ex = U 0 ); (b) the evolution in time of the dimensionless recirculation bubble length s/h.

Figure 6 :

 6 Figure 6: Configuration of the test case: oscillating 2D infinitely thin plate.

) for fluid 1 and fluid 2 ,Figure 7 :

 27 Figure 7: The order of accuracy of the present IB-LB coupling method for the oscillating 2D plate.

Figure 8 :

 8 Figure 8: Configuration of the wetting of a droplet (fluid 1) over a wall Γ IB treated by IB method. θ 1 is the contact angle in this two fluids-solid wall interaction model.

Figure 9 :

 9 Figure 9: Spreading of a droplet on an IB-solid wall under 5 different wetting properties defined in Table2.

Figure 10 :

 10 Figure 10: The contact angle θ 1 under different wetting conditions imposed by varying G 2 adh

Figure 11 :

 11 Figure 11: The contact lines under 3 different wetting conditions on stationary cylinder.

Figure 12 :

 12 Figure 12: Breaking of a 3D fluid film: (a) the general configuration of the numerical simulation; (b) the section of the fluid domain sliced by the plane ABCD.

Figure 13 :

 13 Figure 13: Numerical results of the breaking of a 3D fluid film on a hydrophobic wall (the sliced plane ABCD in Figure 12), which are obtained with the proposed IB-LB coupling method (left) and the classical bounce-back rule (right) at four instants in time: (a) t = 500∆t; (b) t = 1000∆t; (c) t = 1500∆t; (d) t = 2000∆t (the black dashed line represents the IB-solid wall, and the grey rectangular is the bounce-back wall).

Figure 14 :

 14 Figure 14: Evolution in time of the growing diameter of the hole under hydrophobic wetting condition: G 1 adh = 0.4, G 2 adh = -0.4 and G coh = 1.8.

Figure 15 :

 15 Figure 15: Configuration of the application case "Flow induced by a cluster of beating cilia".

Figure 16 :

 16 Figure 16: Two stroke phases for the beating cilium: (a) recovery stroke phase; (b) effective stroke phase.

u A x = 2 . 2 ×Figure 17 :

 2217 Figure 17: Time history of the velocities in x-direction at the points A, B and C in Figure 15.

Table 1 :

 1 Comparisons with previous results of the wake parameters and drag coefficient for the stationary wake flow around a fixed circular cylinder at Re = 30 .

Table 2 :

 2

	5 different wetting conditions according to different values of G 1 adh and G 2 adh
	contact angle (θ 1 )	150 •	120 • 90 •	60 •	30 •
	G 1 adh G 2 adh	0.7794 -0.7794 -0.45 0.45	0 0	-0.45 -0.7794 0.45 0.7794

Acknowledgment

The financial support of CNRS through the PEPS 2014 program "Appels projets exploratoires 2014" is greatly acknowledged.