An Interface-Fitted Fictitious Domain Finite Element Method for the Simulation of Neutrally Buoyant Particles in Plane Shear Flow

: In this paper, an interface-ﬁtted ﬁctitious domain ﬁnite element method is developed for the simulation of ﬂuid–rigid particle interaction problems in cases of rotated particles with small displacement, where an interface-ﬁtted mesh is employed for the discrete scheme to capture the ﬂuid–rigid particle interface accurately, thereby improving the solution accuracy near the interface. Moreover, a linearization and decoupling process is presented to release the constraint between velocities of ﬂuid and rigid particles in the ﬁnite element space, and to make the developed numerical method easy to be implemented. Our numerical experiments are carried out using two different moving interface-ﬁtted meshes; one is obtained by a rotational arbitrary Lagrangian–Eulerian (ALE) mapping, and the other one through a local smoothing process among interface-cut elements. A uniﬁed velocity is deﬁned in the entire domain based on the ﬁctitious domain method, making it easier to develop an interface-ﬁtted mesh generation algorithm in a ﬁxed domain. Both show that the proposed method has a good performance in accuracy for simulating a neutrally buoyant particle in plane shear ﬂow. This approach can be easily extended to ﬂuid–structure interaction problems involving ﬂuids in different states and structures in different shapes with large displacements or deformations.


Introduction
Fluid-structure interaction (FSI) is a crucial research field in current engineering applications, see for example [1][2][3][4][5][6][7]. Particle motion in shear flow is a typical representative and has been a hot topic among scholars, with widespread applications in slurry processing [8], colloid separation [9], and biological cell motion [10]. As early as 1922, Jeffery [11] studied the motion of elliptical particles in creeping flow and provided an analytical expression for the relationship between the angular velocity and the angle. Ding and Aidun [12] performed direct numerical simulations using the discrete Boltzmann method to study the motion of elliptical particles in simple shear flow and found that the behavior of particles at low Reynolds numbers agrees well with periodic solutions of Jeffery. Lundell et al. [13] also conducted numerical studies on the motion of an inertial ellipsoid in a creeping linear shear flow of a Newtonian fluid and found that the particle motion is similar to Jeffery's solution but with the addition of an orbit drift. Furthermore, Pasquino et al. [14] used both experimental and numerical approaches to study the migration and chaining behavior of noncolloidal spheres in a worm-like micellar, viscoelastic solution within the shear flow.
The fictitious domain method is an effective approach for solving FSI problems, which extends the fictitious fluid into solid particles and imposes constraints on their motions in the fluid. Initially proposed by Saul'ev in 1963 [15], the fictitious domain method was used to extend partial differential equations from complex domains to regular ones. Glowinski et al. [16] analyzed the flow around rigid bodies using the fictitious domain method with the Lagrange multiplier in the case of known solid motions. Furthermore, under the assumption that rigid body motions are not known in advance, they directly simulated the flow around moving rigid bodies in two-dimensional and three-dimensional Newtonian and non-Newtonian viscous fluids using the fictitious domain method with the distributed Lagrange multipliers in order to impose constraints on the motion of solid bodies [17]. The distributed Lagrange multipliers-based fictitious domain method for solving fluid-structure interaction problems can also be found in [3,18]. Hwang et al. [19] proposed a novel finite element scheme for direct simulation of the simple shear flow of an inertia-free particle suspension in a Newtonian fluid, applying rigid body conditions on particle boundaries as traction forces through Lagrange multipliers. Yu et al. [4] use the distributed Lagrange multipliers-based fictitious domain method in the simulation of the motion of a spherical particle in a deterministic lateral displacement device. Wang et al. [20] proposed a one-field fictitious domain method for solving general FSI problems, where only fluid variables (velocity and pressure) are defined in the whole domain, while structural variables are no longer defined. This approach not only improves computational efficiency but also maintains the generality and robustness of the distributed Lagrange multiplierbased fictitious domain method.
Nevertheless, the fictitious domain method typically employs a uniformly fixed grid in the entire domain while a moving structural grid on the top of it, which can lead to computational errors near the moving interface as it cannot effectively track the immersed boundary. For example, Auricchio et al. [21] have mentioned that the lack of global regularity in the solution of moving interface problems would cause suboptimal convergence of standard finite element formulation if the computational mesh near the interface is not appropriately modified. Considering the inherent limitations of the fictitious domain method on fixed grids, it seems necessary to introduce grid modification around the interfaces to accurately track moving interfaces and improve computational accuracy. For instance, Wan et al. [22] presented a fictitious boundary method for the simulation of particulate flows, which incorporates the motion of solids and interfaces as additional constraints to the governing Navier-Stokes equations, thereby extending the fluid domain to the entire region. In order to address the issue of low accuracy in boundary approximation, they dynamically relocated the mesh through a special partial differential equation to capture the region near the surface of the moving particles with high accuracy.
In this paper, we develop an interface-fitted fictitious domain method for the simulation of fluid-particle interaction problems. On the one hand, compared with the classical fictitious domain method, our solution accuracy is improved since an interface-fitted mesh is used to track the motion of particles. We note that additional computational efforts to generate the body-fitted mesh are required in this proposed method, and two interface-fitted mesh generation algorithms are given in cases of rotated particles with small displacement. On the other hand, compared to the classical arbitrary Lagrangian-Eulerian (ALE) method, we eliminate the requirement of a fixed number of vertices/elements in the fluid mesh by defining a unified unknown velocity over the entire domain, and instead need to generate an interface-fitted grid in the fixed entire domain. We note that both the classical ALE method and the proposed method can deal with FSI problems in cases of particles with small displacement. But due to the elimination of the requirement of a fixed number of vertices/elements in the fluid, the proposed method is expected to be capable of dealing with FSI problems in cases of particles with large displacement without any remeshing and interpolation.
The rest of this paper is structured as follows. In Section 2, we introduce the model problem of neutral buoyancy particles moving in a planar shear flow. Then, we present the weak formulation based on the fictitious domain method and its ALE description, followed by temporal and spatial discretization, as detailed in Section 3. Numerical results of circular and elliptical particles in planar shear flows are presented in Section 4. The conclusions are reported in Section 5.

Model Problems
Let Ω ∈ R 2 be a rectangular domain filled with a Newtonian viscous incompressible fluid and a rigid particle, as depicted in Figure 1. The bottom, right, top and left boundary of Ω are denoted by Γ 1 , Γ 2 , Γ 3 and Γ 4 , respectively. The distance between Γ 2 and Γ 4 (resp. Γ 1 and Γ 3 ) is denoted by L (resp. D). For any time t ∈ [0, T], the interface between fluid and rigid particle is denoted by Γ t , which divide the domain Ω into two parts: the fluid domain Ω t f and the particle domain Ω t s . For the fluid, the governing equations are given by the following Navier-Stokes equations [17], where du f dt = ∂u f ∂t + (u f · ∇)u f , u f is the fluid velocity, g is the gravity, σ f = −pI + 2µ f D(u f ) is the stress tensor with D(u f ) = 1 2 ∇u f + (∇u f ) T , p is the fluid pressure, ρ f and µ f are the fluid density and viscosity, respectively. The initial condition is defined as while the boundary conditions are given as follows, and where g 0 = (U 0 /2, 0) T , n f is the outward unit normal vector to the boundary of the fluid domain. Equations (4) and (5) indicate that the fluid flow field to be studied exhibits characteristics of a shear flow. Such a characteristic arises when velocities of adjacent fluid layers that are parallel to each other are different, which can be attributed to the reinforced boundary conditions. Equations (6) and (7) imply that the fluid velocity is periodic in the x 1 direction with the period L, which serves to restore the scenario of unbounded shear flow in the horizontal direction as described in Jeffery's solution.
The rigid particle motion is defined by the following Euler-Newton's equations [17] for any t ∈ (0, T], where M = Ω 0 s ρ s dx, ρ s and are the mass, density, and moment of inertia of the rigid particle, respectively. G = (G 1 , G 2 ) T is the center of mass, U is the velocity of the center of mass, and, ω and θ are the angular velocity and inclination angle of the particle, respectively. Moreover, the hydrodynamical force and torque of the immersed rigid particle are defined as T The initial conditions for (8)-(11) are defined as where U 0 , ω 0 , G 0 and θ 0 are given initial values. On the fluid-particle interface Γ t , the following interface condition is imposed, which indicates that there is no relative sliding between the fluid and the rigid body.
In summary, all symbols involved in the above model descriptions and their SI units are also listed in Table 1. Table 1. Symbol descriptions and their SI units.

Symbols
Description Units Viscosity of fluid Pa · s −1 ρ s Density of particle kg · m −2 M Mass of particle kg U Velocity of center of mass m · s −1 G Position of center of mass m Table 1. Cont.

Symbols Description Units
I G Moment of inertia of particle kg · m 2 ω Angular velocity of particle rad · s −1 θ Inclination angle of particle rad F Hydrodynamic force of particle N T Hydrodynamic torque of particle N · m

Weak Formulation Based on the Fictitious Domain Method
In this section, we introduce a unified velocity field and present a weak formulation of the aforementioned model problem based on the fictitious domain method.
Introduce Sobolev spaces , For any v f , V , ξ ∈ V f , multiplying both sides of (1) by v f , and conducting integration by parts on the stiffness term, we have where we can rewrite the above boundary integral term as follows by using (8), (9), (12) and (13), Thus, for any (v f , V , ξ) ∈ V f and q f ∈ Q f , we obtain the following weak formulation of fluid part, Next, we introduce the particle velocity u s , which is defined by It is easy to check that and For any V ∈ R 2 and ξ ∈ R, Then for any V ∈ R 2 and ξ ∈ R, we can attain Furthermore, we define a unified velocity which is defined over the entire domain Ω. Introduce Sobolev spaces Note that du dt = ∂u ∂t + (u · ∇)u. Thus, by combining (19), (20) and (28), we obtain the final weak formulation of the fluid-particle interaction problem: for any (v, V , ξ) ∈ V and q f ∈ Q f , find (u, U, ω) ∈ V and p f ∈ Q f such that subjecting to the following initial conditions where For the cases of neutrally buoyant particles, i.e., ρ s = ρ f , Equation (32) can be simply rewritten as

Interface-Fitted Mesh and ALE Mapping
Although the rectangular domain Ω is independent of time in which the unified velocity u is defined, the fluid-particle interface moves along in time. Thus, we need to construct an interface-fitted mesh in Ω in order to capture the surface of immersed particle, and therefore improve the accuracy of numerical simulation.
Let N be a positive integer, ∆t = T/N be the time step and t n = n∆t, n = 0, . . . , N. Denote by T n h a triangulation of Ω n := Ω t n , where the superscript n indicates the mesh is interface-fitted with Γ t n . That is, for any edge e in T n h ,e ∩ Γ t n h = ∅, wheree denotes the interior of e. For any subdomain D, T n h (D) denotes the restriction of T n h on D. Generally, constructing such interface-fitted meshes and frequent remeshing can lead to prohibitive CPU time in the simulation of FSI problems, especially when the displacement or deformation of the structure is large. In our fictitious domain finite element scheme, the unified unknown velocity is defined in the fixed and connected domain Ω, which releases the requirement that the number of vertices in the fluid subdomain needs to be fixed when an ALE method is used. Although the number of vertices of T n h needs to keep fixed, it is much easier to generate an interface-fitted mesh for the entire fixed domain Ω. In fact, some existing moving mesh strategies, such as those in [23][24][25][26], can be employed for generating such interface-fitted background meshes with some slight modifications.
First, in following the idea presented in [26], one can construct an interface-fitted mesh and the corresponding discrete ALE mapping in cases of rotated particles with small displacement, which is called moving mesh Algorithm 1 in this paper, as illustrated in Figure 2 and described below.
Algorithm 1: Moving mesh by a rotational ALE mapping.
1. Draw an artificial circular buffer zone, Ω t r f , in the fluid domain Ω t f to enclose the immersed particle Ω t s and share the same axis of rotation. Denote Triangulate Ω 0 f with T 0 h that fits Γ 0 and Γ 0 rs , i.e., for any edge e ∈ T 0 h , let e ∩ Γ 0 (Ω n−1 r f ) together with the particle at the same angular velocity and about the same axis of rotation to obtain a rotating mesh T n h (Ω n r f ). 4. Locally shift the nodes on the artificial interface Γ rs associated with T n h (Ω n r f ) in order to find and match with the fixed coordinates of the closest nodes on Γ rs associated with T n h (Ω n s f ). Thus, the rotating and conforming mesh T n h (Ω n r f ), further, the total mesh T n h = T n h (Ω n r f ) ∪ T n h (Ω n s f ) are obtained.
We note that the discrete ALE-based mesh velocity, w h , only needs to be computed in Ω r f , since the mesh T n h (Ω n s f ) is fixed. Algorithm 1 guarantees a rigid particle's vertex of T 0 h (resp. T 0 h (Ω r f )) is always the same rigid particle's vertex of T n h (resp. T n h (Ω r f )) all the time, establishing a one-to-one mapping (denoted by X n h ) between vertices of T n h (Ω r f ) and of T 0 h (Ω r f ). Next, we present another moving mesh algorithm called Algorithm 2, which is similar to the method proposed in [23].
We remark that Algorithm 2 may deliver such a interface-fitted mesh T n h whose rigid particle's (resp. fluid) vertices at one time step may be fluid (resp. rigid particle's) vertices at another time step (e.g., see the bottom-right part of Figure 3), but still form a one-to-one discrete ALE mapping between vertices of T n h and of T 0 h that is still denoted by X n h .

Algorithm 2:
Moving mesh by local smoothing among interface-cut elements.

1.
Triangulate Ω into a rectangular mesh T h with a given mesh size h x (resp. h y ) in x-axis (resp. y-axis) direction that is independent of time t. 2. Inspect all vertices of T h in the column-wise lexicographic order, and divide them into fluid-(black), interface-(red), and rigid particle's (blue) vertex sets associated with the current position of fluid-particle interface Γ t n , where the red interface vertices are those vertices closest to Γ t n , as shown in the top-left part of Figure 3. 3. Move all red interface vertices onto the fluid-particle interface to set them as the closest intersection points between Γ t n and T h , respectively. 4. Cut each quadrilateral of T h with the slash diagonal to obtain a triangular mesh of Ω, T h , where some edges of T h may intersect with the interface, as shown in the top-right part of Figure 3. 5. Locally adjust interface-cut elements of T h using vertex smoothing and edge swapping techniques [24] to enhance the mesh quality, simultaneously, ensure the edges of each element are all aligned with the fluid-particle interface Γ t n . Thus, an interface-fitted mesh T n h is generated, as depicted in the bottom-left part of Figure 3. Once the interface-fitted mesh T n h is obtained, the discrete ALE velocity w n h , which is a piecewise linear function associated with T n h , can be computed by for any vertex P of T n h . Moreover, the discrete ALE material derivative is defined as Associated with the interface-fitted mesh T n h , the subdomains Ω n f and Ω n s are approximated by Ω n f ,h and Ω n s,h , respectively.

A Interface-Fitted Fictitious Domain Finite Element Method
Define the following finite element spaces , where P k denotes the set of polynomials with order less than or equal to k. In Algorithm 1, since only the mesh T n h (Ω n r f ) is updated at each time step, we define u n, * h to be a piecewise P 2 finite element interpolation of u n h in Ω n r f by letting for any vertex P in T n+1 h (Ω n+1 r f ). The discrete ALE velocity w n+1 h is calculated by (39) only in Ω n+1 r f . In order to describe the method in a unified formulation, we extend w n+1 h from Ω r f to the whole domain with zero, and let u n, * h Ω n s f = u n h . In Algorithm 2, we define u n, * h to be a piecewise P 2 finite element interpolation of u n h in terms of the same formula (42) but now for each vertex P ∈ T n+1 h . The discrete ALE velocity w n+1 h is defined in the whole domain and can be calculated by (39). Then, using the backward Euler scheme, an interface-fitted fictitious domain finite element method for solving the fluid-particle interaction problem can be written as follows: and q n+1 f ,h ∈ Q n+1 f ,h , subjecting to initial values where u 0 h is a P 2 finite element approximation of u 0 . Note that the mesh T n+1 h depends on u n+1 h . To release the constraint between v h and (V h , ξ h ) in the definition of V n h , we employ the following Newton's linearization and decoupling process.

Numerical Experiments
In this section, we shall validate the performance of the proposed interface-fitted fictitious domain finite element method by numerically solving three Stokes flow-rigid particle interaction problems with different particle shapes. Two moving mesh algorithms given in Section 3.2 are both tested. We remark that no remeshing is performed in all tests, instead, we only need to move some vertices to their new positions in order to construct an interface-fitted mesh at each time step. On the other hand, since the Stokes flow is governed by the transient Stokes equations that is obtained by simply dropping the convection term, (u f · ∇)u f , from Navier-Stokes Equations (1) and (2) , densities of the fluid and rigid particle are the same as ρ f = ρ s = 1 for cases of neutrally buoyant particles. Moreover, the fluid viscosity µ f = 1, and the shear rate of the fluid G s = 1. Therefore, the velocity on Γ 1 and Γ 3 , defined as boundary conditions in (4) and (5), satisfies U 0 /2 = G s D/2 = 1.

Circular Particle
First of all, we consider the case of a circular particle with radius r = 0.15. The mass center of the circular particle is initially located at (1, 1). The spatial mesh size is h = 1/40 and h = 1/60, and the time step is ∆t = 0.05.
It is known that the relationship between the angular speed and the rotating angle of elliptical particles in an unbounded shear Stokes flow is as follows [11] The theoretical angular velocity given by the above equation is related to the ratio of r a and r b , independent of their specific values. Nevertheless, the values of r a and r b used in the numerical experiment should be as small as possible as Jeffery's solution is derivated under the assumption that the range of the domain is much larger than the particle size so that to minimize the influence of solid particles on the fluid. Thus, the angular speed of a circular particle is about −G s /2 = −0.5 in the case of r a = r b = 0.15. Figure 4 shows snapshots of the velocity field and the pressure field surrounding the circular particle obtained by moving-mesh Algorithms 1 and 2 with h = 1/40 and h = 1/60, respectively. It can be observed that the pressure near the fluid-particle interface is higher on the upstream side and lower on the downstream side, since the particle locates at the center of the plane shear flow.
The plot of the angular speed versus the angular displacement is reported in Figure 5. It can be observed that the circular particle inside the shear flow eventually reaches a steady state. The obtained stable angular velocity is −0.4948 with h = 1/40 and −0.4945 with h = 1/60 by Algorithm 1, and −0.4946 with h = 1/40 and −0.4947 with h = 1/60 by Algorithm 2, respectively. They are all in good agreement with Jeffery's solution (57). Table 2 and Figure 6 present angular speeds obtained by two moving mesh algorithms with different radii of circular particles when h = 1/60. Observations indicate that as the radius of circular particles decreases, the angular velocity tends to approach −0.5. The inverse process of this behavior can be ascribed to the increasingly subtle influence of particles on the surrounding fluid due to the increasing particle size.

Elliptical Particle
Next, we consider an elliptical particle with semi-major axes r a = 0.2 and r b = 0.1. The particle's mass center is initially located at (1, 1) with the inclination angle θ = π/4 toward the positive x 1 -axis direction. The other parameters are the same as those used in the case of circular particles.
According to (57), while aligning along the x 1 axis, the elliptical particle reaches the minimum angular velocity of −0.2 since it is almost parallel to the direction of fluid velocity, resulting in the minimum hydrodynamic torque experienced. On the other hand, while aligning with the x 2 axis, the particle attains the maximum angular velocity of −0.8 since it is nearly perpendicular to the direction of fluid velocity, leading to the maximum undergoing hydrodynamic torque. Figure 7 shows snapshots of the velocity field and the pressure field surrounding elliptical particles obtained by moving-mesh Algorithms 1 and 2 with h = 1/40 and h = 1/60, respectively. The same phenomenon, i.e., the pressure near the fluid-particle interface being higher on the upstream side and lower on the downstream side, can also be observed in this case.
The angular speeds obtained by the proposed numerical method are also compared with Jeffery's solution in Figure 8, illustrating that as the angle changes from kπ to kπ + π/2, the major axis of the elliptical particle progressively aligns perpendicular to the direction of fluid velocity. This alignment leads to an augmentation in hydrodynamic torque, concurrently increasing the angular velocity. Conversely, as the angle changes from kπ + π/2 to kπ + π, the major axis of the elliptical particle gradually aligns parallel to the direction of fluid velocity. As a consequence, there is a reduction in hydrodynamic torque, causing a gradual decrease in the angular velocity.

Conclusions and Future Work
In this paper, we propose an interface-fitted fictitious domain finite element method for fluid-rigid particle interaction problems. A main ingredient of the approach is to introduce a unified velocity defined in the entire domain based on the fictitious domain method to eliminate the requirement of a fixed number of nodes in the fluid subdomain that is usually required by the ALE method. Another main ingredient is to use an interface-fitted mesh to capture the interface of fluid and rigid particles to improve the accuracy of the numerical simulation. Both ingredients are based upon ideas that an interface-fitted mesh may lead to a more accurate numerical solution, and, an interface-fitted mesh generation algorithm in a fixed domain is easier to develop in contrast to that in a time-dependent domain.
For interface-fitted mesh algorithms, the classical ALE method usually makes the fluid mesh move by solving a PDE-type of ALE mapping that subjects to the Dirichlet boundary condition given by the structural displacement on the interface, which does not support structural motions with large displacement/deformation, as they require a fixed number of mesh nodes/elements in the fluid domain. Otherwise, if the fluid mesh quality becomes too poor due to a large structural motion, then the fluid domain has to be re-meshed at each time step to resume the ALE-like computations on the newly generated mesh which, however, no longer holds a fixed number of mesh nodes in the fluid domain, resulting in additional interpolation errors. But there is no such obstacle in our developed numerical method. In the future, more attempts at simulating realistic FSI problems with large structural displacement or deformation will be studied by following similar ideas presented in this paper.  Data Availability Statement: Data sharing not applicable. No new data were created or analyzed in this study. Data sharing is not applicable to this article.

Conflicts of Interest:
The authors declare no conflict of interest.