# Modal Spectral Element Solutions to Incompressible Flows over Particles of Complex Shape

Authors: Don Liu; Yonglai Zheng
Journal: Journal of Computational Engineering
 (2014)

---

## Abstract

This paper develops the virtual identity particles (VIP) model to simulate two-phase flows involving complex-shaped particles. VIP assimilates the high efficiency of the Eulerian method and the convenience of the Lagrangian approach in tracking particles. It uses one fixed Eulerian mesh to compute the fluid field and the Lagrangian description to handle constitutive properties of particles. The interaction between the fluid and complex particles is characterized with source terms in the fluid momentum equations, while the same source terms are computed iteratively from the particulate Lagrangian equations. The advantage of VIP is its economy in modeling a two-phase flow problem almost at the cost of solving only the fluid phase with added source terms. This high efficiency in computational cost makes VIP viable for simulating particulate flows with numerous particles. Owing to the spectral convergence and high resolvability of the modal spectral element method, VIP provides acceptable resolution comparable to DNS but at much reduced computational cost. Simulation results indicate that VIP is promising for investigating flows with complex-shaped particles, especially abundant complex particles.


---

## Body

## 1. Introduction

Particles of complex shape, especially in particulate flows, occur widely in scientific research and engineering applications. The most recent advances in the preparation and self-assembly of micron-sized colloids with well-defined anisotropic shapes were reviewed by [1]. Complex and precisely shaped polymeric particles have been made in laboratories to achieve certain functions in advanced materials and microfluids technology [2]. Janus particles have unique anisotropic characteristics, tunable chemistry, and physical properties and have been used in colloidal physics and chemistry for applications including optoelectronics, e-ink, drug delivery, and bioimaging [3]. Complex-shaped metal nanoparticles with controlled morphology and architecture have potential application in material science, chemistry, physics, and medicine [4]. Complex-shaped microgel particles were fabricated via selective polymerization of aqueous two-phase systems [5]. Micro- and nanoparticles of complex shape have been used in various drug delivery applications. The complexity of particle shape does influence biological processes such as circulation, vascular adhesion, and phagocytosis [6]. A systematic study of particles with shape anisotropy and directed assembly of particles into ordered structures via capillary interactions was presented by [7]. Complex particles of any shape, size, and distribution were numerically recognized and counted in stereology [8].Limited research has been reported for complex-shaped particles especially in particulate flows. A simulation model for studying earthquake microphysics and rods mechanics [9] used an arbitrary number of piecewise spherical patches to describe the particle shape and variable-shaped particles were simulated to investigate how particle shape affects frictional properties of fault gouge. A numerical scheme for freely moving rigid particles in complex flows [10] treated particles as a part of the fluid and constrained particles which were filled with the fluid to have rigid body motion using mollification kernels. This scheme was effective in modeling multiple particle problems. A study of complex shaped particles settling under gravity [11] showed that the dispersion behavior depends on the initial orientation even in the absence of thermal agitation. Collisions between convex or nonconvex three-dimensional polyhedral particles were simulated with joint methods, molecular dynamics and multibody systems, to determine their different motion [12]. The simulation model [13–15] used tailored force fields on the fluid to model particulate flows. In the limit of Stokes flow, this method produces the exact result for an isolated particle. A three-dimensional molecular dynamics (MD) simulation of a hopper flow showed that a strong influence of nonconvex particles, represented by Minkowski operators, on the jamming process in granular flows [16]. Another MD simulation using Voronoi construction [17] studied how the specific shape of particle grains affects the macroscopic mechanical behavior of granular matter, for example, noncohesive soil. A simulation using discrete element method [18] studied how complex shape could affect mechanical behavior of powder particles in order to optimize production.The research on half sphere is few and far between. Visual observations were made on the experimental research on flow around a sphere which was half-submerged in still water and forced to oscillate up and down [19]. Numerical simulations for the laminar flow past a hemisphere up to Reynolds number 300 [20] showed that the onsets of axisymmetry breaking and unsteadiness take place at lower Reynolds numbers than those for flow past a sphere. Besides, zero mean lift force was found instead of a nonzero mean for flow past a whole sphere. In addition, a time-periodic asymmetric flow appeared but not for a whole sphere case in [20]. An experimental and numerical study of scouring sediment transport flow around a half-buried sphere exposed to a steady turbulent current showed that as the external turbulence generated, the scour depth increases but time scale decreases dramatically [21].This paper develops a computational model to simulate the interaction of complex-shaped particles with a fluid. This model aims at a balance between the computational economy as well as the accuracy, instead of being accurate but expensive for numerous particles. Since Lagrangian descriptions provide the ease in tracking free surfaces and interfaces between different materials and facilitate the treatment of materials with history-dependent constitutive relations [22], the developed model uses combined Eulerian and Lagrangian approaches.


## 2. Mathematical Description

The challenge in simulating two-phase particulate flows involving complex-shaped solid particles lies in handling the hydrodynamic interactions between the fluid and particles especially when they are abundant. One way is to fully resolve the fluid boundary layer around each complex-shaped particle to obtain hydrodynamic interaction such as drag, lift, and torque from the fluid to the particle. Practically this demands prohibitively complex discretization with a massive number of body-fitting elements or grid points. This paper provides an alternative and develops an economical and accurate model to meet this challenge. This model takes most of the computational time for the fluid in an Eulerian reference and computes the hydrodynamic interaction with a penalized Lagrangian description for each particle. Instead of using tiny elements around each particle, only a virtual particle is modeled. Hence, this model is called virtual identity particles (VIP). Instead of resolving the boundary layer outside a complex-shaped particle, VIP imposes two tailored force fields at the current time instant for this particle and resolves the interior volume and the vicinity of this particle. Then a set of Lagrangian kinetic equations are integrated to compute its velocity and position and so forth, using a penalty method. The VIP model is similar to [23, 24] in theory. Due to the fact that VIP avoids resolving boundary layers around a particle, VIP significantly reduces the computational cost by one to two orders of magnitude compared with a direct numerical simulation.


### 2.1. Efficient Eulerian Description of the Fluid

Taking advantage of the efficiency of the Eulerian description for the fluid, the incompressible Navier-Stokes equations with added source terms are used to describe the fluid motion. The momentum and mass conservation for the fluid in terms of the primitive variablesu, p are
(1)ρ∂u∂t+ρu·∇u=-∇p+μ∇2u+∑k=1N‍fk(x,t),∇·u=0,
where ρ, p, and μ are the density, pressure, and viscosity of the fluid, respectively; x=(x,y,z) is the Eulerian Cartesian coordinates. The source term fk(x,t) characterizes the interaction between the fluid and the complex particle (labeled as k) located at Yk(t) for a total of N particles:
(2)fk(x,t)=(Fk+G→k∇)Φ[x,Yk(t),σ,t],
in which Φ[x,Yk(t),σ,t] is a kernel function that specifies where and how the interaction between the fluid and the complex-shaped particle happens in space and σ is the vector of parameters. The specific form of this kernel function depends on the actual shape of a complex particle. For example, for an upper half sphere with its axis of symmetry parallel to x-axis, the parametric vector is σ=(σ,σ,σ). At zero to moderately high Reynolds numbers (defined with the diameter of the particle), the kernel function is a factorized exponential function of the form
(3)Φ[x,Yk(t),σ,t]=2(2π)3/2σ3Exp[-12σ2ξ(x,Yk(t))]×Exp[-12σ2ψ(x,Yk(t))],
where ξ(x,Yk(t)) is the squared difference between the x components of x and Yk(t); ψ(x,Yk(t)) is the squared distance from a point in the domain to the axis of symmetry of the upper half sphere labeled with the superscript k with the radius r; the parameter σ=rπ-(1/2) is a scaling factor for the effective range of this kernel function. It is similar to the setting in [25–28]. For an ellipsoidal particle
(4)σ=(σ1,σ2,σ3),σi=riπ-(1/2),i=1,2,3,
where r1, r2, and r3 are semiaxes. The kernel function takes the form
(5)Φ[x,Yk(t),σ,t]=1(2π)3/2σ1σ2σ3Exp[-x^22σ12]Exp[-y^22σ22]Exp[-z^22σ32],
where x^, y^, and z^ are the relative distance to its geometrical center in each direction, respectively.The termFk in (2) is the vector sum of forces from the particle k to the fluid, such as the net inertial effect due to the added mass and buoyancy and the forces counteracting the drag (called balancing drag) and lift (called balancing lift), denoted by FCk, and external forces FEk on the particle k such as interparticle collision. The term Fk consists of
(6)Fk=(mp-mf)(g-dVkdt)+FCk+FEk,
where mp is the mass of the particle k; mf is the mass of displaced fluid; Vk is the nominal velocity of the particle k; and g is the specific force due to gravity. The counteracting force FCk, elaborated in the next section of the paper, is from the particle k to the fluid. This force balances the pressure difference, the tangential friction from the fluid to the particle.The termG→k in (2) contains the force couples, equivalent to the sum of the counteracting torques from the particle to the fluid. It is a second order tensor mathematically which could be decomposed into the symmetric and antisymmetric parts. The symmetric part neutralizes the pressure difference around the particle and prevents the virtual particle from deviating away from its actual position. The antisymmetric part is equivalent to torques, denoted by T→k, in different planes exerting on the fluid in response to shear stress. Both parts are computed iteratively, similar to the dipole iteration scheme [27, 29], to correct the field of strain rate tensor to help compute the unknown drag, lift, and torque. Specifically for solid particles, components of G→k reduce the integral averaged strain rate tensor inside the particle to zero. If the particle k has an impermeable surface, the integral averaged strain rate tensor inside the volume of the corresponding virtual particle should be a zero tensor of second order:
(7)∭Vol‍ϵ→Φ[x,Yk(t),σ,t]dVol≈0,
where ϵ→ is the strain rate tensor; Vol is the volume of the computational domain.


### 2.2. Counteracting Drag, Lift, and Torque

The total force and torque from the fluid to the particle satisfy the differential form of Newton’s second law of motion:(8)[mp-mf0102I]ddt[VkΩk]=[FkT→k],
where I is the vector of moment of inertia of the complex particle with respect to the moment axes through its center of mass; 01 and 02 are zero matrices of dimensions 1×3 and 3×1, respectively. Ωk is the vector of the angular velocity of this particle.The process of computing the unknown counteracting drag, lift, and torque from the particle to fluid is within the Lagrangian reference. Because it is unnecessary to resolve the nonslip surface boundary layer around each particle in the VIP formulation, tiny curvature-fitting elements are avoided in VIP runs. In order to obtain the unknown drag, lift, and torque on each particle, the nominal velocity and angular velocity of the particle are penalized against the actual values to compute the unknowns which counteract the drag, lift, and torque on from the fluid to the particle. To implement this penalty algorithm, a pseudotimeτ is introduced between two consecutive time steps and a set of pseudotime-dependent problems about the counteracting force and torque is solved:
(9)∂FCk(t,τ)∂τ=λ1[V-k(t)-Vk(t,τ)],(10)∂Tij(t,τ)∂τ=λ2ϵijk[Ω-k(t)-Ωk(t,τ)],
where FCk(t,τ) is the counteract force as in (6); Tij(t,τ) is the component of the counteracting torque T→k in the i-j plane; V-k(t) and Ω-k(t) are the actual velocity and the actual angular velocity of the particle; ϵijk is the permutation symbol for the spatial indices (the superscript k is the particle label); λ1 and λ2 are the associated penalty parameters. Their magnitude controls the speed of convergence. The larger the penalty parameter, the faster the convergence rate. However, due to their explicit nature, the magnitude is subject to a stability condition. Upon convergence, the nominal velocity Vk(t,τ) equals the actual velocity of the particle k, and the unknown counteracting force becomes independent of τ, and the unknown FCk(t,τ) is found. Similarly, the component Tij of the counteracting torque T→k acting in the i-j plane is computed from the difference in angular velocity. Once these forces and torque on each particle are computed, they enter the Navier-Stokes equations to solve for the flow field at the next time step. This penalty algorithm applies to both resistance and mobility problems.


### 2.3. Accurate Lagrangian Description of Particles

The nominal velocity of each particleVk(t) is computed from the nonlinear filtering, that is, a volume integral of the fluid velocity with the kernel function
(11)Vk(t)=∭Vol‍u(x,t)Φ[x,Yk(t),σ,t]dVol.
Similarly, the nominal angular velocity Ωk(t) of the particle is computed with another volume integral using the fluid vorticity ∇×u(x,t) and the kernel
(12)Ωk(t)=12∭Vol‍∇×u(x,t)Φ[x,Yk(t),σ,t]dVol.
The translational position vector and rotational angular displacement of the particle are computed explicitly from the Lagrangian constitutive description:
(13)dYk(t)=Vk(t)dt,dAk(t)=Ωk(t)dt.
The problems studied in this paper areresistance problems [30] in multiphase flows. Therefore, the translational and rotational velocities are zero as particles are fixed. However, during computations with the penalty algorithm, they were penalized to be nonzero until reaching convergence and return to zero.


### 2.4. Discretization Method and Integration Schemes

The modal spectral/hp element method developed by [31, 32] was implemented in space for all simulations and a high order splitting algorithm [33] was used in time to obtain the numerical solution to the incompressible Navier-Stokes equations. Using hierarchical orthogonal functions as the expansion and test bases, this modal spectral element method has high order accuracy with errors decreasing at the exponential rate. Gauss-Lobatto-Legendre quadrature was extensively used to numerically evaluate volume integrals for properties of particles. Adam-Bashforth schemes of third order accuracy were used to integrate Lagrangian equations for particles.


## 2.1. Efficient Eulerian Description of the Fluid

Taking advantage of the efficiency of the Eulerian description for the fluid, the incompressible Navier-Stokes equations with added source terms are used to describe the fluid motion. The momentum and mass conservation for the fluid in terms of the primitive variablesu, p are
(1)ρ∂u∂t+ρu·∇u=-∇p+μ∇2u+∑k=1N‍fk(x,t),∇·u=0,
where ρ, p, and μ are the density, pressure, and viscosity of the fluid, respectively; x=(x,y,z) is the Eulerian Cartesian coordinates. The source term fk(x,t) characterizes the interaction between the fluid and the complex particle (labeled as k) located at Yk(t) for a total of N particles:
(2)fk(x,t)=(Fk+G→k∇)Φ[x,Yk(t),σ,t],
in which Φ[x,Yk(t),σ,t] is a kernel function that specifies where and how the interaction between the fluid and the complex-shaped particle happens in space and σ is the vector of parameters. The specific form of this kernel function depends on the actual shape of a complex particle. For example, for an upper half sphere with its axis of symmetry parallel to x-axis, the parametric vector is σ=(σ,σ,σ). At zero to moderately high Reynolds numbers (defined with the diameter of the particle), the kernel function is a factorized exponential function of the form
(3)Φ[x,Yk(t),σ,t]=2(2π)3/2σ3Exp[-12σ2ξ(x,Yk(t))]×Exp[-12σ2ψ(x,Yk(t))],
where ξ(x,Yk(t)) is the squared difference between the x components of x and Yk(t); ψ(x,Yk(t)) is the squared distance from a point in the domain to the axis of symmetry of the upper half sphere labeled with the superscript k with the radius r; the parameter σ=rπ-(1/2) is a scaling factor for the effective range of this kernel function. It is similar to the setting in [25–28]. For an ellipsoidal particle
(4)σ=(σ1,σ2,σ3),σi=riπ-(1/2),i=1,2,3,
where r1, r2, and r3 are semiaxes. The kernel function takes the form
(5)Φ[x,Yk(t),σ,t]=1(2π)3/2σ1σ2σ3Exp[-x^22σ12]Exp[-y^22σ22]Exp[-z^22σ32],
where x^, y^, and z^ are the relative distance to its geometrical center in each direction, respectively.The termFk in (2) is the vector sum of forces from the particle k to the fluid, such as the net inertial effect due to the added mass and buoyancy and the forces counteracting the drag (called balancing drag) and lift (called balancing lift), denoted by FCk, and external forces FEk on the particle k such as interparticle collision. The term Fk consists of
(6)Fk=(mp-mf)(g-dVkdt)+FCk+FEk,
where mp is the mass of the particle k; mf is the mass of displaced fluid; Vk is the nominal velocity of the particle k; and g is the specific force due to gravity. The counteracting force FCk, elaborated in the next section of the paper, is from the particle k to the fluid. This force balances the pressure difference, the tangential friction from the fluid to the particle.The termG→k in (2) contains the force couples, equivalent to the sum of the counteracting torques from the particle to the fluid. It is a second order tensor mathematically which could be decomposed into the symmetric and antisymmetric parts. The symmetric part neutralizes the pressure difference around the particle and prevents the virtual particle from deviating away from its actual position. The antisymmetric part is equivalent to torques, denoted by T→k, in different planes exerting on the fluid in response to shear stress. Both parts are computed iteratively, similar to the dipole iteration scheme [27, 29], to correct the field of strain rate tensor to help compute the unknown drag, lift, and torque. Specifically for solid particles, components of G→k reduce the integral averaged strain rate tensor inside the particle to zero. If the particle k has an impermeable surface, the integral averaged strain rate tensor inside the volume of the corresponding virtual particle should be a zero tensor of second order:
(7)∭Vol‍ϵ→Φ[x,Yk(t),σ,t]dVol≈0,
where ϵ→ is the strain rate tensor; Vol is the volume of the computational domain.


## 2.2. Counteracting Drag, Lift, and Torque

The total force and torque from the fluid to the particle satisfy the differential form of Newton’s second law of motion:(8)[mp-mf0102I]ddt[VkΩk]=[FkT→k],
where I is the vector of moment of inertia of the complex particle with respect to the moment axes through its center of mass; 01 and 02 are zero matrices of dimensions 1×3 and 3×1, respectively. Ωk is the vector of the angular velocity of this particle.The process of computing the unknown counteracting drag, lift, and torque from the particle to fluid is within the Lagrangian reference. Because it is unnecessary to resolve the nonslip surface boundary layer around each particle in the VIP formulation, tiny curvature-fitting elements are avoided in VIP runs. In order to obtain the unknown drag, lift, and torque on each particle, the nominal velocity and angular velocity of the particle are penalized against the actual values to compute the unknowns which counteract the drag, lift, and torque on from the fluid to the particle. To implement this penalty algorithm, a pseudotimeτ is introduced between two consecutive time steps and a set of pseudotime-dependent problems about the counteracting force and torque is solved:
(9)∂FCk(t,τ)∂τ=λ1[V-k(t)-Vk(t,τ)],(10)∂Tij(t,τ)∂τ=λ2ϵijk[Ω-k(t)-Ωk(t,τ)],
where FCk(t,τ) is the counteract force as in (6); Tij(t,τ) is the component of the counteracting torque T→k in the i-j plane; V-k(t) and Ω-k(t) are the actual velocity and the actual angular velocity of the particle; ϵijk is the permutation symbol for the spatial indices (the superscript k is the particle label); λ1 and λ2 are the associated penalty parameters. Their magnitude controls the speed of convergence. The larger the penalty parameter, the faster the convergence rate. However, due to their explicit nature, the magnitude is subject to a stability condition. Upon convergence, the nominal velocity Vk(t,τ) equals the actual velocity of the particle k, and the unknown counteracting force becomes independent of τ, and the unknown FCk(t,τ) is found. Similarly, the component Tij of the counteracting torque T→k acting in the i-j plane is computed from the difference in angular velocity. Once these forces and torque on each particle are computed, they enter the Navier-Stokes equations to solve for the flow field at the next time step. This penalty algorithm applies to both resistance and mobility problems.


## 2.3. Accurate Lagrangian Description of Particles

The nominal velocity of each particleVk(t) is computed from the nonlinear filtering, that is, a volume integral of the fluid velocity with the kernel function
(11)Vk(t)=∭Vol‍u(x,t)Φ[x,Yk(t),σ,t]dVol.
Similarly, the nominal angular velocity Ωk(t) of the particle is computed with another volume integral using the fluid vorticity ∇×u(x,t) and the kernel
(12)Ωk(t)=12∭Vol‍∇×u(x,t)Φ[x,Yk(t),σ,t]dVol.
The translational position vector and rotational angular displacement of the particle are computed explicitly from the Lagrangian constitutive description:
(13)dYk(t)=Vk(t)dt,dAk(t)=Ωk(t)dt.
The problems studied in this paper areresistance problems [30] in multiphase flows. Therefore, the translational and rotational velocities are zero as particles are fixed. However, during computations with the penalty algorithm, they were penalized to be nonzero until reaching convergence and return to zero.


## 2.4. Discretization Method and Integration Schemes

The modal spectral/hp element method developed by [31, 32] was implemented in space for all simulations and a high order splitting algorithm [33] was used in time to obtain the numerical solution to the incompressible Navier-Stokes equations. Using hierarchical orthogonal functions as the expansion and test bases, this modal spectral element method has high order accuracy with errors decreasing at the exponential rate. Gauss-Lobatto-Legendre quadrature was extensively used to numerically evaluate volume integrals for properties of particles. Adam-Bashforth schemes of third order accuracy were used to integrate Lagrangian equations for particles.


## 3. Simulation Results and Verifications

### 3.1. Flow over Spheroid at Finite Reynolds Number

A spheroid with semiaxesa1=2, a2=a3=1, in dimensionless units (same as below), is fixed at the origin of a Cartesian coordinate system in which a channel of size -14≤x1≤6, -10/3≤x2≤10, and -5≤x3≤5 is located. A parabolic fluid velocity of profile u1=(1+3x2/10)(1-x2/10) is specified at the inlet plane x1=-14. The semiaxes of this spheroid are aligned with all there Eulerian coordinate axes and the long semiaxis is along x1, the main flow direction. Since this spheroid is off the centerline in the channel, the undisturbed approach velocity towards its center is uo=1. The density and dynamic viscosity of the fluid are scaled to 1. Based on the undisturbed approach velocity uo=1, the length scale 2a2=2, and the kinematic viscosity ν=1, the nominal Reynolds number is 2, although the actual peak Reynolds number is a litter larger. Periodic boundary conditions were imposed at x3=±5. Nonslip wall boundary conditions were specified at x2=-10/3 and x2=10. A plane view of this flow is given in Figure 1(a).Comparisons between the DNS run and the VIP run for a spheroid atRe=2 in the plane x3=0 through its center of mass: (a) contour lines of the streamwise fluid velocity; (b) spectral elements used in the VIP run (upper) and the DNS run (lower).
(a)
(b)Since this spheroid is off-centered in the channel, it receives a drag force inx1 and a lateral lift force in x2. Due to the asymmetrical shear stress, this spheroid also receives a torque from the fluid. These forces and torque are unknowns prior to all simulations. Because this spheroid is fixed without any motion, this problem is aresistance problem in multiphase flows [30]. The penalty method described earlier in (9) and (10) was implemented to compute the unknown balancing forces and torque, which counteract the drag, and shear induced lift and torque from the fluid.Because no analytic solution is readily available for comparison, numerical solutions using the VIP model were produced and verified with a direct numerical simulation (DNS) for accuracy. The computational domain for the VIP run is the entire channel including the volume occupied by the spheroid since the VIP run simulates a virtual particle. But for the DNS run, the domain is the volume inside the channel and outside the spheroid. Conforming and non-uniform computational meshes were used for both runs, with VIP runs on 1,920 hexahedral elements and DNS runs on 2,220 hexahedral elements. The scaled drag, lift, torque, and their associated relative errors are tabulated in Table1. The relatively larger error in the lift force is because it is proportional to the shear stress which is sensitive to the spatial resolution.Table 1
Comparison of the scaled drag, lift, and torque from DNS run and VIP run for a spheroid in a channel flow at Reynolds number 2. The errors shown in percentage are relative to values from DNS run.
Source
D
/
(
μ
a
2
u
o
)
L
/
(
μ
a
2
u
o
)
τ
/
(
μ
a
2
2
u
o
)
V
1
/
u
o
Ω
3
a
2
/
u
o
VIP run
34.371
3.078
4.102
4.61
×
1
0
-
4
1.47
×
1
0
-
3
DNS run
34.231
2.898
4.12
0
0
Relative error
0.4
%
6.2
%
0.5
%
0.05
%
0.15
%To compare the overall flow field characteristics between the VIP run and DNS run, Figure1(a) presents the contour lines of the mainstream fluid velocity field. In the VIP (upper) plot, a virtual spheroid was simulated and the fluid fills the volume nominally occupied by the spheroid. In the DNS (lower) plot, the nonslip boundary condition was imposed on the surface of the spheroid, which is impermeable to the fluid. Therefore contour lines are outside the spheroid. It is shown in this figure that outside the spheroid, the flow field in VIP run is similar to that in DNS run. However, in terms of the computational cost, the VIP run is more economical than the DNS run.Figure1(b) shows the spectral elements used in the VIP (upper) and DNS (lower) runs. To enhance the resolution and to maintain less elements, finer elements are arranged around the surface of this spheroid in the DNS run. The advantage of VIP is that it is unnecessary to resolve the surface boundary layer in a VIP run. Therefore, much larger and less elements were used in the VIP run as seen in the upper plot. This is why the computational cost is much less than that of the DNS run.Next, the details of the flow field are compared between the VIP and DNS runs. Figure2(a) shows the distribution of the mainstream fluid velocity along the line parallel to x1 through the center of mass of this spheroid (x2=x3=0). Up to the surface of this spheroid, the VIP predicted a velocity profile almost the same as from a DNS. Figure 2(b) demonstrates the pressure distribution along x1 direction. Pressure distribution curves from the VIP and DNS runs match very well except near the surface of the spheroid (x1/a2=±2 in the plot). This was because the VIP does not impose nonslip boundary condition and does not resolve the thin boundary layer on the surface.Comparisons between the DNS and VIP runs for a spheroid atRe=2: (a) mainstream fluid velocity profile versus x1; (b) normalized pressure distribution versus x1. The ellipse is represented in dash-dot-dot.
(a)
(b)Figure3(a) presents the mainstream fluid velocity profiles, u1 normalized by uo, versus the lateral distance x2 between two side walls, in the plane x3=0. The fluid velocity profiles were taken at varied locations along x1, such as the center of the spheroid x1=0, tangential planes at two tips x1=±2 upstream and downstream and further away from its center at x1=±4. Although lines are clustered together, it is still observed that velocity profiles match well and all the way up to the surface of the spheroid. Figure 3(b) plots the lateral velocity (u2) versus x2, at the same locations as in Figure 3(a). Except some small difference at the surface of the spheroid, those lateral velocity profiles from the VIP run match with those from the DNS run. Since the DNS run does not resolve the fluid inside the spheroid, the diamond symbols at x1=0 in Figure 3(b) stop at -1≤x2≤1.Comparisons between the DNS run and VIP run for a spheroid atRe=2. Fluid velocity profile: (a) u1 versus x2; (b) u2 versus x2.
(a)
(b)In summary, the main characteristics of a particulate flow could be simulated by the VIP model with good accuracy, although on the surface of a particle, the flow field is less accurate due to the fact that nonslip boundary conditions are not imposed at all. However, the gain in computational time outweighs the reduced accuracy on the surface of a particle. This indicates that VIP is ideal for modeling flows with numerous particles involved, especially for engineering application.


### 3.2. Plane Poiseuille Flow over Upper Half Sphere

A plane Poiseuille flow over an upper half sphere is numerically simulated with VIP and verified with DNS. The channel is0≤x≤6, -3≤y≤3, and -4≤z≤6, in dimensionless units (same as below). This half sphere of unit radius is fixed stationary and its geometric center is at (3,0,0) in the middle between two parallel planar walls. However, its center of mass is at (3.375,0,0). A Dirichlet boundary condition was imposed at the inlet plane z=-4 so that a parabolic profile of fluid velocity, which has a peak w=10 at the center x=3 and w=0 at both walls x=0 and 6, was established. Periodic fluid boundary conditions were set on the planes y=±3. Nonslip wall boundary conditions were specified on the surface of the half sphere and on both planar walls at x=0 and 6. The density and viscosity of the fluid were scaled to 1. The nominal Reynolds number based on the incoming peak fluid velocity and the diameter is 20. The incurred drag and torque about the y-axis were unknown prior to both VIP and DNS runs. Regarding the spatial resolution, VIP runs used 800 hexahedral elements at 3rd, 4th, and 5th order basis polynomial expansion. DNS runs were on 1,744 hexahedral elements at 4th, 5th, and 10th order basis expansion.The computed dimensionless drag force, lift force, torque, and so forth from 10th order DNS and 5th order VIP runs are tabulated in the first and second rows in Table2. The negative sign for lift force means that it is against the positive direction of y-axis. Comparing the corresponding values between these two rows, the relative differences in drag and lift are within 3%. The relative difference in torque is less than 4.7%. DNS data were at high resolution settings with the highest hierarchical basis function being 10th order orthogonal polynomial. Therefore, a parallel algorithm was implemented for this run using 48 processors in a Linux clusters for almost 30 hours. Results from the VIP run were from a single i5 processor for a little over 2 hours. The computational time in DNS is about 500 times of that in VIP. This indicates that VIP runs are far more economical than DNS runs.Table 2
Summary of dimensionless drag, lift, and torque from DNS and VIP runs at different orders of basis expansion, computational time, and processor information for different cases. [1] Upper half sphere with the torque center set at the center of mass (3.375, 0, 0). [2] Lower half sphere with the torque center set at the center of mass (2.625, 0, 0). [3] Upper half sphere with the torque center set at the spherical center (3, 0, 0). [4] Lower half sphere with the torque center set at the spherical center (3, 0, 0).
Runs
Order
Drag
Lift
Torque
Computational time
CPU number/type
DNS [1]
10
348.09
−10.71
40.55
29:45:14
48 Intel 2.33 GHz
VIP [1]
5
346.66
−10.35
42.44
2:10:07
One i5 3.3 GHz
DNS [2]
5
339.60
9.89
−39.23
3:40:01
One i5 3.3 GHz
VIP [2]
5
345.98
10.05
−42.21
1:46:09
One i5 3.3 GHz
DNS [3]
5
342.27
−9.44
−90.14
4:07:28
One i5 3.3 GHz
DNS [4]
5
339.60
9.89
88.12
4:14:54
One i5 3.3 GHzAlthough the spatial resolution (degrees of freedom) used in DNS is about 14 times higher than in VIP, the latter is still able to predict the main feature of the interaction between the particle and the fluid. Figure4(a) compares contour lines of the mainstream fluid velocity (z) between DNS (upper) and VIP (lower) runs in the plane y=0 through the spherical center of the half sphere. Figure 4(b) is the similar comparison but in the plane x=3 through its spherical center. The computational elements are shown in the background in both plots. Both the DNS (upper) and VIP (lower) plots in Figure 4(a) indicate that the fluid detours mostly from below the half sphere due to the asymmetry and less resistance. Similar patterns of velocity distribution especially in the wake region are demonstrated. The upper plot in Figure 4(b) shows an almost circular stagnation zone. This is exactly the result of an impermeable nonslip boundary condition in the DNS run. The corresponding plot below it from the VIP run has a noticeably smaller stagnation zone. But outside the circle, both plots are in good agreement. There are two reasons. First, the impermeable nonslip boundary condition was never imposed in the VIP run; instead, the volume within the half sphere was completely resolved. Second, the geometric specifics of this half sphere were set to the fluid through a kernel function which distributes the force and torque coupled between the fluid and the half sphere. The lower plot illustrated salient features of the field structure in DNS at a much reduced cost in both time and spatial resolution.Comparisons between the DNS and VIP runs for the top half sphere atRe=20. Fluid velocity distribution in (a) y=0 and (b) x=3 planes. Elements are shown in black lines.
(a)
(b)Figure5 compares results from the DNS run and the VIP run about the field details in the plane z=0 through the center of the half sphere perpendicular to the main flow. Due to imposed nonslip velocity boundary conditions on the top (x=6) and bottom (x=0) walls, the fluid detours through two sides of the half sphere in both the DNS and VIP runs, although the VIP run did not have impermeable nonslip wall boundary condition imposed and the fluid has free access to the interior of the half sphere. The difference between two plots in Figure 5 is partially due to the difference in the actual computational domain. In the DNS run, it is the volume outside the surface of the half sphere but within the channel; in the VIP run, it is the entire volume inside the channel. Considering at only one fourteenth of the resolution of DNS and one five-hundredth of the wall clock time of DNS, the VIP run has provided acceptable results.Figure 5
Comparisons between the DNS and VIP runs for the top half sphere atRe=20. Fluid velocity distribution in z=0 plane. Elements are shown in black lines.


### 3.3. Plane Poiseuille Flow over Lower Half Sphere

Now we study the similar flow but with a lower half sphere in exactly the same channel. This lower half sphere still has unit radius and is fixed stationary so that its geometric center is at(3,0,0) in the middle between two parallel planar walls. However, its center of mass is at (2.625,0,0). The orientation of this half sphere is just rotated 180° from the previous case. The same boundary conditions as the previous case were imposed at the inlet plane z=-4, on the top (x=0) and bottom (x=6) walls, the spherical surface, and two sides at planes y=±3. Therefore, the same parabolic inlet fluid velocity, a peak w=10 at the center x=3 and w=0 at both walls x=0 and 6, was established. The density and viscosity of the fluid were scaled to 1. The nominal Reynolds number based on the incoming peak fluid velocity and the diameter remains 20. The induced drag and torque about the y-axis were unknown prior to both VIP and DNS runs. Similar spatial resolution settings were used as before: VIP runs on 800 hexahedral elements at 3rd, 4th, and 5th order basis polynomial expansion and DNS runs on 1,744 hexahedral elements at 3rd, 4th, and 5th order basis expansion.The computed dimensionless drag force, lift force, torque, and so forth from both DNS and VIP runs are tabulated in the third and fourth rows in Table2. The negative sign for torque values means that it is along the negative direction of y-axis. DNS data were obtained from a run using a relatively low resolution setting at 5th order basis expansion. Therefore, a serial algorithm was used to solve the incompressible Navier-Stokes equations. Results from both DNS and VIP runs were generated from an i5 processor. Comparing the corresponding values between the third and fourth rows, the relative differences in drag and lift are within 2%. The relative difference in torque is less than 7.6%. Next we only look at results from VIP runs and just compare corresponding values between the second rows and the fourth row. It is observed that drag values are very close with only 0.2% relative difference; lift values have opposite signs but within 3% relative difference in magnitude; and torque value have opposite signs as well but within 0.6% relative difference. The opposite signs are just as expected. Both at the 5th order basis expansion, the spatial resolution in the DNS run is about twice of that of the VIP run due to the number of elements used. Therefore, the computational time for the former is about twice of that for the latter. This indicates that VIP handles a fluid-particle interactive flow as if a DNS run for a single fluid flow.Figure6(a) compares contour lines of the mainstream fluid velocity (z) between DNS (upper) and VIP (lower) runs in the plane y=0 through the spherical center of the lower half sphere. Figure 6(b) is the similar comparison but in the plane x=3 through its spherical center. The computational elements are shown in the background in both plots. Both the DNS (upper) and VIP (lower) plots in Figure 6(a) indicate that the fluid detours mostly from below the half sphere due to the asymmetry and less resistance. Similar patterns of velocity distribution especially in the wake region are demonstrated. The upper plot in Figure 6(b) shows an almost circular stagnation zone. This is exactly the result of an impermeable nonslip boundary condition in the DNS run. The corresponding plot below it from the VIP run has a noticeably smaller stagnation zone. But outside the circle, both plots are in good agreement. There are two reasons. First, the impermeable nonslip boundary condition was never imposed in the VIP run; instead, the volume within the half sphere was completely resolved. Second, the geometric specifics of this half sphere were set to the fluid through a kernel function which distributes the force and torque coupled between the fluid and the half sphere. The lower plot illustrated salient features of the field structure in DNS at a much reduced cost in both time and spatial resolution.Comparisons between the DNS and VIP runs for the bottom half sphere atRe=20. Fluid velocity contours at (a) y=0 plane and (b) x=3 plane. The actual elements are also shown.
(a)
(b)Figure7 compares results from the DNS and the VIP runs about the field details in the plane z=0 through the center of the half sphere perpendicular to the main flow. Due to imposed nonslip velocity boundary conditions on the top (x=6) and bottom (x=0) walls, the fluid detours through two sides of the half sphere in both the DNS and VIP runs, although the VIP run did not have impermeable nonslip wall boundary condition imposed and the fluid has free access to the interior of the half sphere. The difference between two plots in Figure 7 is partially due to the difference in the actual computational domain. In the DNS run, it is the volume outside the surface of the half sphere but within the channel; in the VIP run, it is the entire volume inside the channel. Considering at only one-fourteenth of the resolution of DNS and one five-hundredth of the wall clock time of DNS, the VIP run has provided acceptable results.Figure 7
Comparisons between the DNS and VIP runs for the bottom half sphere atRe=20. Fluid velocity contour lines at z=0 plane. The actual elements are also shown.Regarding the torque on the half sphere in DNS runs, it is important to have the torque center set at the center of mass. This is because the mesh for DNS runs only provides coordinates and curvature information, no specifics about the torque center. To see that different settings of the torque center do influence the results in DNS results, two DNS runs for the upper and lower half sphere, respectively, were performed with the torque center intentionally set at the spherical center instead of the mass center. Results are also tabulated in the fifth and sixth rows of Table2. These results are compared with results from the corresponding DNS runs with the torque center set at the center of mass. For the upper half sphere, the first row versus the fifth row in Table 2, results show that the difference in drag force is within 1.7% and the difference in lift force is less than 12%. There are differences due to different expansion orders, 10th versus 5th and parallel versus serial algorithms. However, the differences in torque value are significant. Not only the opposite sign appears, but also the magnitude increases to 222% of the correct value in row one. For the lower half sphere case, the third row versus the sixth row in Table 2, the values for both drag and lift forces are identical because both runs used the serial algorithm at exactly the same order of expansion. However, there are significant differences in the torque value as well. The magnitude increases to 225% with an opposite sign. Therefore, it is crucial to set the torque center at the center of mass. Comparing those two runs with torque center set at the spherical center, that is, values in the fifth and sixth rows in Table 2, it is observed that the drag force is within 1% difference; the lift force has opposite signs as anticipated and is within 4.8% difference; and the torque has opposite signs as well and is about 2% in difference. For the cases with the upper and lower half spheres, the lift force and torque should have opposite signs. This is reasonable from physics point of view.However, the torque on the half sphere in VIP runs is easy to handle—simply specifying nothing about the torque center. It is unnecessary to set the torque center in VIP runs because the torque is generated from the volume averaged rotational rate (VARR), that is, the nominal angular velocity in (12), via the penalty method as described in (10). The torque value was computed from penalizing VARR instead of evaluating the surface integrals around the half sphere involving the normal and shear stresses on the spherical surface. Simulation results indicate that VIP is capable of computing not only the accurate magnitude of the torque, but also the correct direction.


## 3.1. Flow over Spheroid at Finite Reynolds Number

A spheroid with semiaxesa1=2, a2=a3=1, in dimensionless units (same as below), is fixed at the origin of a Cartesian coordinate system in which a channel of size -14≤x1≤6, -10/3≤x2≤10, and -5≤x3≤5 is located. A parabolic fluid velocity of profile u1=(1+3x2/10)(1-x2/10) is specified at the inlet plane x1=-14. The semiaxes of this spheroid are aligned with all there Eulerian coordinate axes and the long semiaxis is along x1, the main flow direction. Since this spheroid is off the centerline in the channel, the undisturbed approach velocity towards its center is uo=1. The density and dynamic viscosity of the fluid are scaled to 1. Based on the undisturbed approach velocity uo=1, the length scale 2a2=2, and the kinematic viscosity ν=1, the nominal Reynolds number is 2, although the actual peak Reynolds number is a litter larger. Periodic boundary conditions were imposed at x3=±5. Nonslip wall boundary conditions were specified at x2=-10/3 and x2=10. A plane view of this flow is given in Figure 1(a).Comparisons between the DNS run and the VIP run for a spheroid atRe=2 in the plane x3=0 through its center of mass: (a) contour lines of the streamwise fluid velocity; (b) spectral elements used in the VIP run (upper) and the DNS run (lower).
(a)
(b)Since this spheroid is off-centered in the channel, it receives a drag force inx1 and a lateral lift force in x2. Due to the asymmetrical shear stress, this spheroid also receives a torque from the fluid. These forces and torque are unknowns prior to all simulations. Because this spheroid is fixed without any motion, this problem is aresistance problem in multiphase flows [30]. The penalty method described earlier in (9) and (10) was implemented to compute the unknown balancing forces and torque, which counteract the drag, and shear induced lift and torque from the fluid.Because no analytic solution is readily available for comparison, numerical solutions using the VIP model were produced and verified with a direct numerical simulation (DNS) for accuracy. The computational domain for the VIP run is the entire channel including the volume occupied by the spheroid since the VIP run simulates a virtual particle. But for the DNS run, the domain is the volume inside the channel and outside the spheroid. Conforming and non-uniform computational meshes were used for both runs, with VIP runs on 1,920 hexahedral elements and DNS runs on 2,220 hexahedral elements. The scaled drag, lift, torque, and their associated relative errors are tabulated in Table1. The relatively larger error in the lift force is because it is proportional to the shear stress which is sensitive to the spatial resolution.Table 1
Comparison of the scaled drag, lift, and torque from DNS run and VIP run for a spheroid in a channel flow at Reynolds number 2. The errors shown in percentage are relative to values from DNS run.
Source
D
/
(
μ
a
2
u
o
)
L
/
(
μ
a
2
u
o
)
τ
/
(
μ
a
2
2
u
o
)
V
1
/
u
o
Ω
3
a
2
/
u
o
VIP run
34.371
3.078
4.102
4.61
×
1
0
-
4
1.47
×
1
0
-
3
DNS run
34.231
2.898
4.12
0
0
Relative error
0.4
%
6.2
%
0.5
%
0.05
%
0.15
%To compare the overall flow field characteristics between the VIP run and DNS run, Figure1(a) presents the contour lines of the mainstream fluid velocity field. In the VIP (upper) plot, a virtual spheroid was simulated and the fluid fills the volume nominally occupied by the spheroid. In the DNS (lower) plot, the nonslip boundary condition was imposed on the surface of the spheroid, which is impermeable to the fluid. Therefore contour lines are outside the spheroid. It is shown in this figure that outside the spheroid, the flow field in VIP run is similar to that in DNS run. However, in terms of the computational cost, the VIP run is more economical than the DNS run.Figure1(b) shows the spectral elements used in the VIP (upper) and DNS (lower) runs. To enhance the resolution and to maintain less elements, finer elements are arranged around the surface of this spheroid in the DNS run. The advantage of VIP is that it is unnecessary to resolve the surface boundary layer in a VIP run. Therefore, much larger and less elements were used in the VIP run as seen in the upper plot. This is why the computational cost is much less than that of the DNS run.Next, the details of the flow field are compared between the VIP and DNS runs. Figure2(a) shows the distribution of the mainstream fluid velocity along the line parallel to x1 through the center of mass of this spheroid (x2=x3=0). Up to the surface of this spheroid, the VIP predicted a velocity profile almost the same as from a DNS. Figure 2(b) demonstrates the pressure distribution along x1 direction. Pressure distribution curves from the VIP and DNS runs match very well except near the surface of the spheroid (x1/a2=±2 in the plot). This was because the VIP does not impose nonslip boundary condition and does not resolve the thin boundary layer on the surface.Comparisons between the DNS and VIP runs for a spheroid atRe=2: (a) mainstream fluid velocity profile versus x1; (b) normalized pressure distribution versus x1. The ellipse is represented in dash-dot-dot.
(a)
(b)Figure3(a) presents the mainstream fluid velocity profiles, u1 normalized by uo, versus the lateral distance x2 between two side walls, in the plane x3=0. The fluid velocity profiles were taken at varied locations along x1, such as the center of the spheroid x1=0, tangential planes at two tips x1=±2 upstream and downstream and further away from its center at x1=±4. Although lines are clustered together, it is still observed that velocity profiles match well and all the way up to the surface of the spheroid. Figure 3(b) plots the lateral velocity (u2) versus x2, at the same locations as in Figure 3(a). Except some small difference at the surface of the spheroid, those lateral velocity profiles from the VIP run match with those from the DNS run. Since the DNS run does not resolve the fluid inside the spheroid, the diamond symbols at x1=0 in Figure 3(b) stop at -1≤x2≤1.Comparisons between the DNS run and VIP run for a spheroid atRe=2. Fluid velocity profile: (a) u1 versus x2; (b) u2 versus x2.
(a)
(b)In summary, the main characteristics of a particulate flow could be simulated by the VIP model with good accuracy, although on the surface of a particle, the flow field is less accurate due to the fact that nonslip boundary conditions are not imposed at all. However, the gain in computational time outweighs the reduced accuracy on the surface of a particle. This indicates that VIP is ideal for modeling flows with numerous particles involved, especially for engineering application.


## 3.2. Plane Poiseuille Flow over Upper Half Sphere

A plane Poiseuille flow over an upper half sphere is numerically simulated with VIP and verified with DNS. The channel is0≤x≤6, -3≤y≤3, and -4≤z≤6, in dimensionless units (same as below). This half sphere of unit radius is fixed stationary and its geometric center is at (3,0,0) in the middle between two parallel planar walls. However, its center of mass is at (3.375,0,0). A Dirichlet boundary condition was imposed at the inlet plane z=-4 so that a parabolic profile of fluid velocity, which has a peak w=10 at the center x=3 and w=0 at both walls x=0 and 6, was established. Periodic fluid boundary conditions were set on the planes y=±3. Nonslip wall boundary conditions were specified on the surface of the half sphere and on both planar walls at x=0 and 6. The density and viscosity of the fluid were scaled to 1. The nominal Reynolds number based on the incoming peak fluid velocity and the diameter is 20. The incurred drag and torque about the y-axis were unknown prior to both VIP and DNS runs. Regarding the spatial resolution, VIP runs used 800 hexahedral elements at 3rd, 4th, and 5th order basis polynomial expansion. DNS runs were on 1,744 hexahedral elements at 4th, 5th, and 10th order basis expansion.The computed dimensionless drag force, lift force, torque, and so forth from 10th order DNS and 5th order VIP runs are tabulated in the first and second rows in Table2. The negative sign for lift force means that it is against the positive direction of y-axis. Comparing the corresponding values between these two rows, the relative differences in drag and lift are within 3%. The relative difference in torque is less than 4.7%. DNS data were at high resolution settings with the highest hierarchical basis function being 10th order orthogonal polynomial. Therefore, a parallel algorithm was implemented for this run using 48 processors in a Linux clusters for almost 30 hours. Results from the VIP run were from a single i5 processor for a little over 2 hours. The computational time in DNS is about 500 times of that in VIP. This indicates that VIP runs are far more economical than DNS runs.Table 2
Summary of dimensionless drag, lift, and torque from DNS and VIP runs at different orders of basis expansion, computational time, and processor information for different cases. [1] Upper half sphere with the torque center set at the center of mass (3.375, 0, 0). [2] Lower half sphere with the torque center set at the center of mass (2.625, 0, 0). [3] Upper half sphere with the torque center set at the spherical center (3, 0, 0). [4] Lower half sphere with the torque center set at the spherical center (3, 0, 0).
Runs
Order
Drag
Lift
Torque
Computational time
CPU number/type
DNS [1]
10
348.09
−10.71
40.55
29:45:14
48 Intel 2.33 GHz
VIP [1]
5
346.66
−10.35
42.44
2:10:07
One i5 3.3 GHz
DNS [2]
5
339.60
9.89
−39.23
3:40:01
One i5 3.3 GHz
VIP [2]
5
345.98
10.05
−42.21
1:46:09
One i5 3.3 GHz
DNS [3]
5
342.27
−9.44
−90.14
4:07:28
One i5 3.3 GHz
DNS [4]
5
339.60
9.89
88.12
4:14:54
One i5 3.3 GHzAlthough the spatial resolution (degrees of freedom) used in DNS is about 14 times higher than in VIP, the latter is still able to predict the main feature of the interaction between the particle and the fluid. Figure4(a) compares contour lines of the mainstream fluid velocity (z) between DNS (upper) and VIP (lower) runs in the plane y=0 through the spherical center of the half sphere. Figure 4(b) is the similar comparison but in the plane x=3 through its spherical center. The computational elements are shown in the background in both plots. Both the DNS (upper) and VIP (lower) plots in Figure 4(a) indicate that the fluid detours mostly from below the half sphere due to the asymmetry and less resistance. Similar patterns of velocity distribution especially in the wake region are demonstrated. The upper plot in Figure 4(b) shows an almost circular stagnation zone. This is exactly the result of an impermeable nonslip boundary condition in the DNS run. The corresponding plot below it from the VIP run has a noticeably smaller stagnation zone. But outside the circle, both plots are in good agreement. There are two reasons. First, the impermeable nonslip boundary condition was never imposed in the VIP run; instead, the volume within the half sphere was completely resolved. Second, the geometric specifics of this half sphere were set to the fluid through a kernel function which distributes the force and torque coupled between the fluid and the half sphere. The lower plot illustrated salient features of the field structure in DNS at a much reduced cost in both time and spatial resolution.Comparisons between the DNS and VIP runs for the top half sphere atRe=20. Fluid velocity distribution in (a) y=0 and (b) x=3 planes. Elements are shown in black lines.
(a)
(b)Figure5 compares results from the DNS run and the VIP run about the field details in the plane z=0 through the center of the half sphere perpendicular to the main flow. Due to imposed nonslip velocity boundary conditions on the top (x=6) and bottom (x=0) walls, the fluid detours through two sides of the half sphere in both the DNS and VIP runs, although the VIP run did not have impermeable nonslip wall boundary condition imposed and the fluid has free access to the interior of the half sphere. The difference between two plots in Figure 5 is partially due to the difference in the actual computational domain. In the DNS run, it is the volume outside the surface of the half sphere but within the channel; in the VIP run, it is the entire volume inside the channel. Considering at only one fourteenth of the resolution of DNS and one five-hundredth of the wall clock time of DNS, the VIP run has provided acceptable results.Figure 5
Comparisons between the DNS and VIP runs for the top half sphere atRe=20. Fluid velocity distribution in z=0 plane. Elements are shown in black lines.


## 3.3. Plane Poiseuille Flow over Lower Half Sphere

Now we study the similar flow but with a lower half sphere in exactly the same channel. This lower half sphere still has unit radius and is fixed stationary so that its geometric center is at(3,0,0) in the middle between two parallel planar walls. However, its center of mass is at (2.625,0,0). The orientation of this half sphere is just rotated 180° from the previous case. The same boundary conditions as the previous case were imposed at the inlet plane z=-4, on the top (x=0) and bottom (x=6) walls, the spherical surface, and two sides at planes y=±3. Therefore, the same parabolic inlet fluid velocity, a peak w=10 at the center x=3 and w=0 at both walls x=0 and 6, was established. The density and viscosity of the fluid were scaled to 1. The nominal Reynolds number based on the incoming peak fluid velocity and the diameter remains 20. The induced drag and torque about the y-axis were unknown prior to both VIP and DNS runs. Similar spatial resolution settings were used as before: VIP runs on 800 hexahedral elements at 3rd, 4th, and 5th order basis polynomial expansion and DNS runs on 1,744 hexahedral elements at 3rd, 4th, and 5th order basis expansion.The computed dimensionless drag force, lift force, torque, and so forth from both DNS and VIP runs are tabulated in the third and fourth rows in Table2. The negative sign for torque values means that it is along the negative direction of y-axis. DNS data were obtained from a run using a relatively low resolution setting at 5th order basis expansion. Therefore, a serial algorithm was used to solve the incompressible Navier-Stokes equations. Results from both DNS and VIP runs were generated from an i5 processor. Comparing the corresponding values between the third and fourth rows, the relative differences in drag and lift are within 2%. The relative difference in torque is less than 7.6%. Next we only look at results from VIP runs and just compare corresponding values between the second rows and the fourth row. It is observed that drag values are very close with only 0.2% relative difference; lift values have opposite signs but within 3% relative difference in magnitude; and torque value have opposite signs as well but within 0.6% relative difference. The opposite signs are just as expected. Both at the 5th order basis expansion, the spatial resolution in the DNS run is about twice of that of the VIP run due to the number of elements used. Therefore, the computational time for the former is about twice of that for the latter. This indicates that VIP handles a fluid-particle interactive flow as if a DNS run for a single fluid flow.Figure6(a) compares contour lines of the mainstream fluid velocity (z) between DNS (upper) and VIP (lower) runs in the plane y=0 through the spherical center of the lower half sphere. Figure 6(b) is the similar comparison but in the plane x=3 through its spherical center. The computational elements are shown in the background in both plots. Both the DNS (upper) and VIP (lower) plots in Figure 6(a) indicate that the fluid detours mostly from below the half sphere due to the asymmetry and less resistance. Similar patterns of velocity distribution especially in the wake region are demonstrated. The upper plot in Figure 6(b) shows an almost circular stagnation zone. This is exactly the result of an impermeable nonslip boundary condition in the DNS run. The corresponding plot below it from the VIP run has a noticeably smaller stagnation zone. But outside the circle, both plots are in good agreement. There are two reasons. First, the impermeable nonslip boundary condition was never imposed in the VIP run; instead, the volume within the half sphere was completely resolved. Second, the geometric specifics of this half sphere were set to the fluid through a kernel function which distributes the force and torque coupled between the fluid and the half sphere. The lower plot illustrated salient features of the field structure in DNS at a much reduced cost in both time and spatial resolution.Comparisons between the DNS and VIP runs for the bottom half sphere atRe=20. Fluid velocity contours at (a) y=0 plane and (b) x=3 plane. The actual elements are also shown.
(a)
(b)Figure7 compares results from the DNS and the VIP runs about the field details in the plane z=0 through the center of the half sphere perpendicular to the main flow. Due to imposed nonslip velocity boundary conditions on the top (x=6) and bottom (x=0) walls, the fluid detours through two sides of the half sphere in both the DNS and VIP runs, although the VIP run did not have impermeable nonslip wall boundary condition imposed and the fluid has free access to the interior of the half sphere. The difference between two plots in Figure 7 is partially due to the difference in the actual computational domain. In the DNS run, it is the volume outside the surface of the half sphere but within the channel; in the VIP run, it is the entire volume inside the channel. Considering at only one-fourteenth of the resolution of DNS and one five-hundredth of the wall clock time of DNS, the VIP run has provided acceptable results.Figure 7
Comparisons between the DNS and VIP runs for the bottom half sphere atRe=20. Fluid velocity contour lines at z=0 plane. The actual elements are also shown.Regarding the torque on the half sphere in DNS runs, it is important to have the torque center set at the center of mass. This is because the mesh for DNS runs only provides coordinates and curvature information, no specifics about the torque center. To see that different settings of the torque center do influence the results in DNS results, two DNS runs for the upper and lower half sphere, respectively, were performed with the torque center intentionally set at the spherical center instead of the mass center. Results are also tabulated in the fifth and sixth rows of Table2. These results are compared with results from the corresponding DNS runs with the torque center set at the center of mass. For the upper half sphere, the first row versus the fifth row in Table 2, results show that the difference in drag force is within 1.7% and the difference in lift force is less than 12%. There are differences due to different expansion orders, 10th versus 5th and parallel versus serial algorithms. However, the differences in torque value are significant. Not only the opposite sign appears, but also the magnitude increases to 222% of the correct value in row one. For the lower half sphere case, the third row versus the sixth row in Table 2, the values for both drag and lift forces are identical because both runs used the serial algorithm at exactly the same order of expansion. However, there are significant differences in the torque value as well. The magnitude increases to 225% with an opposite sign. Therefore, it is crucial to set the torque center at the center of mass. Comparing those two runs with torque center set at the spherical center, that is, values in the fifth and sixth rows in Table 2, it is observed that the drag force is within 1% difference; the lift force has opposite signs as anticipated and is within 4.8% difference; and the torque has opposite signs as well and is about 2% in difference. For the cases with the upper and lower half spheres, the lift force and torque should have opposite signs. This is reasonable from physics point of view.However, the torque on the half sphere in VIP runs is easy to handle—simply specifying nothing about the torque center. It is unnecessary to set the torque center in VIP runs because the torque is generated from the volume averaged rotational rate (VARR), that is, the nominal angular velocity in (12), via the penalty method as described in (10). The torque value was computed from penalizing VARR instead of evaluating the surface integrals around the half sphere involving the normal and shear stresses on the spherical surface. Simulation results indicate that VIP is capable of computing not only the accurate magnitude of the torque, but also the correct direction.


## 4. Conclusion

VIP integrates an Eulerian description of the fluid dynamics and a Lagrangian description of the kinetics of all particles. The computational domain consists of the volume occupied by fluid and complex-shaped particles in it. VIP treats two phases, particles, and the fluid as one through imposing tailored force fields that constraint the fluid in the way that the interaction with a solid particle is approximated. The complexity of the particle shape is handled in the kernel function which distributes force fields according to its particular shape. This model avoids the demand to completely resolve the thin boundary layer around the complex-shaped particles and avoids the requirement of tiny elements and large numbers of nodes near the surface curvature. Therefore, this model significantly reduces the computational intensity by one to two orders of magnitude as indicated in this paper and [14, 34]. Due to the balance between computational accuracy and economical consideration in formulation, simulation results indicate that VIP is promising for investigating flows with complex-shaped particles, especially for the case with abundant complex particles.



---

*Source: 768538-2014-04-30.xml*