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.
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 variables , are where , , and are the density, pressure, and viscosity of the fluid, respectively; is the Eulerian Cartesian coordinates. The source term characterizes the interaction between the fluid and the complex particle (labeled as ) located at for a total of particles: in which 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 -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 where is the squared difference between the components of and ; is the squared distance from a point in the domain to the axis of symmetry of the upper half sphere labeled with the superscript with the radius ; the parameter 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 where , , and are semiaxes. The kernel function takes the form where , , and are the relative distance to its geometrical center in each direction, respectively.
The term in (2) is the vector sum of forces from the particle 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 , and external forces on the particle such as interparticle collision. The term consists of where is the mass of the particle ; is the mass of displaced fluid; is the nominal velocity of the particle ; and is the specific force due to gravity. The counteracting force , elaborated in the next section of the paper, is from the particle to the fluid. This force balances the pressure difference, the tangential friction from the fluid to the particle.
The term 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 , 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 reduce the integral averaged strain rate tensor inside the particle to zero. If the particle 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: where is the strain rate tensor; 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: where is the vector of moment of inertia of the complex particle with respect to the moment axes through its center of mass; and are zero matrices of dimensions and , respectively. 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: where is the counteract force as in (6); is the component of the counteracting torque in the - plane; and are the actual velocity and the actual angular velocity of the particle; is the permutation symbol for the spatial indices (the superscript is the particle label); and 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 equals the actual velocity of the particle , and the unknown counteracting force becomes independent of , and the unknown is found. Similarly, the component of the counteracting torque acting in the - 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 particle is computed from the nonlinear filtering, that is, a volume integral of the fluid velocity with the kernel function Similarly, the nominal angular velocity of the particle is computed with another volume integral using the fluid vorticity and the kernel The translational position vector and rotational angular displacement of the particle are computed explicitly from the Lagrangian constitutive description: The problems studied in this paper are resistance 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/ 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 semiaxes , , in dimensionless units (same as below), is fixed at the origin of a Cartesian coordinate system in which a channel of size , , and is located. A parabolic fluid velocity of profile is specified at the inlet plane . The semiaxes of this spheroid are aligned with all there Eulerian coordinate axes and the long semiaxis is along , the main flow direction. Since this spheroid is off the centerline in the channel, the undisturbed approach velocity towards its center is . The density and dynamic viscosity of the fluid are scaled to 1. Based on the undisturbed approach velocity , the length scale , and the kinematic viscosity , the nominal Reynolds number is 2, although the actual peak Reynolds number is a litter larger. Periodic boundary conditions were imposed at . Nonslip wall boundary conditions were specified at and . A plane view of this flow is given in Figure 1(a).
(a)
(b)
Since this spheroid is off-centered in the channel, it receives a drag force in and a lateral lift force in . 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 a resistance 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 Table 1. The relatively larger error in the lift force is because it is proportional to the shear stress which is sensitive to the spatial resolution.
To compare the overall flow field characteristics between the VIP run and DNS run, Figure 1(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.
Figure 1(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. Figure 2(a) shows the distribution of the mainstream fluid velocity along the line parallel to through the center of mass of this spheroid (). 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 direction. Pressure distribution curves from the VIP and DNS runs match very well except near the surface of the spheroid ( 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.
(a)
(b)
Figure 3(a) presents the mainstream fluid velocity profiles, normalized by , versus the lateral distance between two side walls, in the plane . The fluid velocity profiles were taken at varied locations along , such as the center of the spheroid , tangential planes at two tips upstream and downstream and further away from its center at . 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 () versus , 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 in Figure 3(b) stop at .
(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 is , , and , in dimensionless units (same as below). This half sphere of unit radius is fixed stationary and its geometric center is at in the middle between two parallel planar walls. However, its center of mass is at . A Dirichlet boundary condition was imposed at the inlet plane so that a parabolic profile of fluid velocity, which has a peak at the center and at both walls and 6, was established. Periodic fluid boundary conditions were set on the planes . Nonslip wall boundary conditions were specified on the surface of the half sphere and on both planar walls at 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 -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 Table 2. The negative sign for lift force means that it is against the positive direction of -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.
Although 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. Figure 4(a) compares contour lines of the mainstream fluid velocity () between DNS (upper) and VIP (lower) runs in the plane through the spherical center of the half sphere. Figure 4(b) is the similar comparison but in the plane 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.
(a)
(b)
Figure 5 compares results from the DNS run and the VIP run about the field details in the plane through the center of the half sphere perpendicular to the main flow. Due to imposed nonslip velocity boundary conditions on the top () and bottom () 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.
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 in the middle between two parallel planar walls. However, its center of mass is at . 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 , on the top () and bottom () walls, the spherical surface, and two sides at planes . Therefore, the same parabolic inlet fluid velocity, a peak at the center and at both walls 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 -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 Table 2. The negative sign for torque values means that it is along the negative direction of -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.
Figure 6(a) compares contour lines of the mainstream fluid velocity () between DNS (upper) and VIP (lower) runs in the plane through the spherical center of the lower half sphere. Figure 6(b) is the similar comparison but in the plane 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.
(a)
(b)
Figure 7 compares results from the DNS and the VIP runs about the field details in the plane through the center of the half sphere perpendicular to the main flow. Due to imposed nonslip velocity boundary conditions on the top () and bottom () 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.
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 Table 2. 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.
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
Acknowledgments
This research was supported by the National Science Foundation under Grants DMS-1115546 and DMS-1318988 and Tongji University, Shanghai, China. The computational resources were provided by the Extreme Science and Engineering Discovery Environment (XSEDE), Louisiana Optical Network Initiative (LONI), and Louisiana State University.