Monte Carlo Vortical Smoothed Particle Hydrodynamics for Simulating Turbulent Flows

For vortex particle methods relying on SPH‐based simulations, the direct approach of iterating all fluid particles to capture velocity from vorticity can lead to a significant computational overhead during the Biot‐Savart summation process. To address this challenge, we present a Monte Carlo vortical smoothed particle hydrodynamics (MCVSPH) method for efficiently simulating turbulent flows within an SPH framework. Our approach harnesses a Monte Carlo estimator and operates exclusively within a pre‐sampled particle subset, thus eliminating the need for costly global iterations over all fluid particles. Our algorithm is decoupled from various projection loops which enforce incompressibility, independently handles the recovery of turbulent details, and seamlessly integrates with state‐of‐the‐art SPH‐based incompressibility solvers. Our approach rectifies the velocity of all fluid particles based on vorticity loss to respect the evolution of vorticity, effectively enforcing vortex motions. We demonstrate, by several experiments, that our MCVSPH method effectively preserves vorticity and creates visually prominent vortical motions.


Introduction
Fluid animation constitutes a captivating area of research within the computer graphics community.Among various fluid simulation techniques, the Navier-Stokes equations [CF88] serve as the cornerstone for modern fluid solvers.In the conventional operator splitting scheme [Bri15], the particle velocity field u n at time step n is updated in two steps.First, non-pressure effects (such as gravity and viscous forces) are applied to calculate an intermediate velocity projection for three-dimensional fluid simulation; and vortex sheets [BKB12] for boundaries.While such methods can capture vorticity features and preserve fine turbulent visual details, the computational cost of iterating over the vortex elements is significant [NWRC22].For instance, the direct approach to integrating the vortex particle method into SPH-based fluid simulation frameworks involves directly considering all fluid particles as vortex particles to track vorticity evolution and compute velocity through the Biot-Savart summation.Typical simulations use many fluid particles, so this direct method imposes high computational costs.
Inspired by the work of Rioux-Lavoie et al. [RSÖ*22] using a Monte Carlo estimator to simulate fluids in an Eulerian framework, we propose the Monte Carlo Vortical Smoothed Particle Hydrodynamics (MCVSPH) method to simulate turbulent fluids in an SPH framework.Our method operates independently of pressure projection loops, allowing the use of existing state-of-the-art incompressibility solvers [ICS*14; BK17] and simplifying the process of implementing a complete Lagrangian fluid solver.We incorporate our vortex particle method into an SPH-based framework using the Monte Carlo stochastic method.This adaption allows us to derive a corrective velocity from vorticity using the Biot-Savart summation within a pre-sampled small particle subset, and makes it feasible to implement the vortex particle method efficiently within a pure particle framework, all without the need for expensive global iterations.To organise vortex particles, we treat vortex and SPH particles as separate entities.Following Monte Carlo method principles, we randomly sample a subset of SPH particles and pair each of them with a corresponding vortex particle.This enables us to adopt probabilistic methods to estimate velocities that would otherwise necessitate iterating through all SPH fluid particles.
Our main contributions are summarised as follows: • A Lagrangian approach enhancing turbulent details; • Improved Biot-Savart summation in SPH with a Monte Carlo method; • An organisation scheme for vortex particles for turbulent detail preservation.
The remainder of this paper is organised as follows.We start by reviewing related work (Sec.2) and also briefly introduce SPH and the vortex particle method (Sec.3).Next, we introduce our MCVSPH approach (Sec.4).We illustrate our approach by a series of three-dimensional simulation experiments and comparisons with other SPH-based methods for turbulent flow (Sec.5).Finally, we conclude the paper and provide directions for future work (Sec.6).

Related work
A substantial body of research exists within the computer graphics community focusing on fluid simulation and SPH methods.We next provide a concise review of relevant work focusing on SPH and techniques for enhancing turbulent details.
SPH-based Incompressible Fluids.Simulating visually authentic fluids often requires incompressibility.Many previous studies focus on achieving incompressible fluid simulations through SPHbased methods.Becker and Teschner [BT07] developed an explicit approach to compute pressure from densities using the Tait equation.Yet, the simulated fluids remain weakly compressible.To overcome this, Solenthaler and Pajarola [SP09] proposed a predictivecorrective scheme (PCISPH), iteratively solving particle pressures from density change using a global stiffness constant.Ihmsen et al. [ICS*14] enhanced the pressure-solving iterations with the implicit incompressible SPH (IISPH) using density deviation computed from velocities and computing stiffness for each particle per time step.However, PCISPH and IISPH solve the Pressure Poisson Equation (PPE) without explicitly enforcing a divergence-free velocity field.
Bender and Koschier [BK17] proposed divergence-free SPH (DF-SPH), where position updates correct densities and velocity updates mitigate velocity divergence.In our work, we use DFSPH for its ability to enforce invariant densities and divergence-free velocities.Yet, typical SPH-based incompressible fluid solvers exhibit noticeable numerical damping from coarse discretisations [IOS*14], leading to significant loss of turbulent details and the need for additional approaches to preserve vorticity and turbulent motions.Our MCVSPH method enhances turbulent details for general fluid simulations independently of the pressure projection process and can be integrated with most existing incompressibility solvers.
Numerical Dissipation Mitigation.The semi-Lagrangian method is commonly employed to achieve accurate advection.Yet, it can inadvertently blur high-frequency fluid motion information due to the use of low-order accuracy interpolation computations.Eulerian-Lagrangian hybrid methods, such as PIC [Har62] and FLIP [ZB05], mitigate numerical dissipation, yet frequent interpolations during particle-grid transfers can cause momentum loss.Jiang et al. [JSS*15] proposed the affine particle-in-cell (APIC) scheme, which considers angular momentum and shearing deformation during grid-particle transfers.Fu et al. [FGG*17] expand upon the APIC method by generalising the local function on each particle.Qu et al. [QLDJ22] compute weights for particle-grid transfers using the power particle method, which ensures uniform particle distributions and volume preservation.
In the realm of pure Eulerian simulation, Qu et al. [QZG*19] introduce the BiMocq scheme which capitalises on dual mesh characteristics and multi-level mapping to effectively compensate error in a back-and-forth manner.Zehnder et al. [ZNT18] propose an advection-reflection scheme applying an energy-preserving reflection operator halfway through the advection step to reduce the dissi- pation at the projection step.Nabizadeh et al. [NWRC22] enhance the advection-reflection scheme by incorporating covectors to preserve circulation.Previous efforts to mitigate numerical dissipation primarily address advection dissipation present in Eulerian methods.In contrast, our work focuses on maintaining intricate turbulent details within a Lagrangian framework.
Vorticity Confinement.The vorticity confinement method, introduced to the computer graphics community by Fedkiw et al. [FSJ01], amplifies existing vortices using corrective forces.Letine et al. [LAF11] enhance the computation of these forces for improved energy and momentum conservation.While initially proposed for grid-based scenarios [JKB*10; WM19], vorticity confinement was applied to SPH-based methods fluids by Macklin and Müller [MM13] to enhance particle vortical motions.Although the vorticity confinement method is relatively simple to implement, its sensitivity to parameters can result in the generation of excessive energy and unstable simulations.Moreover, since the amplification forces are contingent on existing vortices, this approach cannot generate new vortices, which restricts its use in simulating fluids with rich turbulent features.
Micropolar Fluids.The micropolar fluid method [Eri66; Luk99] describes fluid micro-structures with non-symmetric stress tensors.This approach finds applications in areas where fluid micro-motion is significant, such as heat transfer [PMZ22], blood flows [KSPS20], and magnetised fluids [ANK22].Bender et al. [BKKW19] introduced the micropolar method to computer graphics and applied it to SPH-based methods (MPSPH) to simulate turbulent flows.However, as pointed out by Chen et al. [CLL10], micro-rotations represent nonsolid-body-like rotation (gyration).Given that this method does not explicitly model solid-body-like rotation (vorticity), the resulting turbulent features often manifest themselves as small-scale rotational motions.Our method instead employs the vorticity equation to explicitly model vortical motions.
Non-Navier-Stokes Methods.To circumvent the challenges faced by general Navier-Stokes solvers, such as vorticity loss, some methods have turned to non-Navier-Stokes methods.Chern [ZBG15] propose the IVOCK scheme to restore lost vorticity and expedite the Biot-Savart summation using the fast multipole algorithm [Dar00].Liu et al. [LWB*21] adopt the vorticity equation and the stream function to reconstruct turbulent details in SPH fluid simulations, where the Biot-Savart summation traverses vortex particles within the support radius.However, the computational overhead of traversing all vortex particles in the Biot-Savart summation is high.Previous studies have introduced solutions to expedite this process within Eulerian frameworks, such as PPPM [ZB14] and the Monte Carlo method [RSÖ*22].Yet, the application of vortex methods to SPH frameworks remains relatively unexplored.
Given the effectiveness of the vortex particle method in modelling vortex motion, we harness its potential to generate turbulent motion within SPH-based fluids.Our approach leverages the vortex particle method within an SPH-based framework using a Monte Carlo estimator, conducting the Biot-Savart summation via a pre-sampled subset of particles.

Fundamentals of SPH and the vortex particle method
As we use SPH-based approaches in our work, we next revisit the governing Navier-Stokes equations and general SPH-based numerical computations for various differential operators needed in fluid simulations (Sec.3.1).We also cover the basics of the vortex particle method and explain the bottleneck of its use in SPH (Sec.3.2).
The Navier-Stokes equations [CF88] read where Du Dt = ∂u ∂t + (u • ∇)u is the material derivative; p is the pressure; υ is the kinematic viscosity; and g models external forces.The equations govern fluid dynamics and form the basis of fluid solvers.However, for computational purposes, it needs to be discretised.

SPH discretisation
SPH is used as a spatial discretisation method to numerically approximate the differential operators within the governing equations [KBST19].For general approximation without any differential operators, values on fluid particle i are evaluated by where j denotes a neighbor particle of particle i; m j and ρ j are the mass and the density, respectively, of particle j; and W i j = W (x i − x j , h), where x denote particle positions, is a kernel function with support radius h.In this paper, we use for W the cubic spline function [KBST19] given by where a = Vorticity.Vorticity ω is defined as the curl ∇ × u of the velocity.
In the context of discretisation, two SPH formulations exist for representing this curl.The symmetric curl form is given by Alternatively, with A ji = A j − A i , the difference curl is given by As Eqn. (6) exhibits lower error in boundary areas [BKKW19], we use this form for vorticity computation in our method.
The difference form of SPH formulations can also be applied to the discretisation of the gradient operator where m ⊗ n = mn T denotes the dyadic product.
Viscosity.In scenarios where the simulated fluid is not entirely inviscid, the computations of the Laplacian of velocity ∆u and the Laplacian of vorticity ∆ω become necessary.We use the following SPH approximation to discretise the Laplacian operator: where x i j = x i − x j , and d is the spatial dimension of the simulated scenario.This approach reduces the order of the differential equation, enhances numerical stability, and helps conserve linear and angular momenta [Mon92].
Pressure.The pressure force enforces incompressibility of the fluid field.In our work, we opt for the divergence-free SPH solver (DF-SPH) [BK17], which maintains constant density and a divergencefree velocity field.
Boundary handling.We address interactions between fluid and surrounding solids by the approach of Akinci et al. [AIA*12].This involves representing solid boundaries using particles, treating them similarly to fluid particles.The smoothed mass of a solid boundary particle i with density ρ 0 and neighbouring boundary particles k is computed as where m i is the non-smoothed mass of particle i.

The vortex particle method
The vortex particle method has two key components: the vortex strength vector and the Biot-Savart summation, as follows.
The vortex strength vector β represents the circulation at a specific position x within the fluid field.In fluid dynamics, the general circulation is computed as where ∂V is a closed curve in 2D or a closed surface in 3D surrounding position x (see Fig. 3 left).Using Stokes' theorem, we can transform Eqn.(10) in 3D to where V is the volume enclosed by the surface ∂V (see Fig. 3

middle).
Ideally, when computing β at a given position x, the closed surface ∂V and the enclosed volume V should be infinitesimal.In practice, V should be small enough so that ω can be considered constant when evaluating β.In this case, the vortex strength vector β of vortex particle k can be approximated as where ω k is the vorticity of vortex particle k and v k is its (small) volume, respectively (see Fig. 3

right).
The Biot-Savart summation [WL93] recovers the velocity field from vorticity as where K× represents the Biot-Savart kernel, * denotes the convolution operator, and ω is the approximated vorticity.In threedimensional scenarios, the Biot-Savart kernel is given by Expanding the convolution operator * , Eqn. (13) can be written as where χ denotes the domain enclosed by ∂V .
However, Eqn. (15) involves an integration that is not conducive to numerical computation.To discretise this equation, we use the vortex particle method, wherein vortices are abstracted as particles that store position and vortex strength information.The approximated vorticity ω is represented as where x k denotes the position of vortex particle k; δ is the Dirac delta function; and Np is the total number of vortex particles used to cover the whole fluid domain.
By substituting Eqn. ( 16) into the Biot-Savart law (Eqn.( 15)), the velocity computation becomes This shows that the velocity field is recovered from vorticity by iterating over all vortex particles.In SPH, physical attributes (such as position, velocity, and pressure) of a fluid are tracked using a collection of SPH fluid particles -which could also serve as vortex particles that hold vorticity information.Yet, due to the high number of fluid particles, iterating over all these particles becomes impractical when applying the vortex particle method to SPH.An optimised approach, which we present next, solves this issue.

MCVSPH Scheme
We apply the vortex particle method in SPH-based fluid simulations to recover turbulent features.For this, we use the vorticity equation for vorticity evolution to compute vorticity loss in the velocity field (Sec.4.1).We optimise computing velocity corrections from vorticity loss by Monte Carlo sampling, as described next (Sec.4.2).

Vorticity loss
To recover turbulent details, we start from the vorticity equation representing the time evolution of vorticity ω = ∇ × u where υ is the same kinematic viscosity as the one in Eqn.(1).The first term (ω • ∇)u gives the vortex stretching governing the evolution of vortex filaments.
In an ideal scenario, the curl ∇ × u of velocity is identical to the vorticity ω.This is not the case in practice, due to distinct equations governing the evolution of velocity and vorticity and the numerical dissipation during the splitting scheme for updating velocity.Zhang [ZBG15] proposed a method that evaluates lost vorticity by the difference between updated vorticity and the curl of updated velocity.Inspired by this, in our SPH-based approach, we 'save' the advection step and approximate differential operators by SPH formulations to evaluate vorticity loss.Specifically, we quantify vorticity loss for each particle by the following steps: • Compute vorticity ω n at current time step n using the curl of current velocity field ω n = ∇ × u n ; • Evolve current vorticity ω n with the vorticity equation ( 18) to obtain the vorticity ω n+1 at the next time step; • Apply non-pressure forces fnon-pressure to the current velocity u n to acquire the intermediate velocity ūn by ūn = u n + ∆t fnon-pressure; • Quantify the vorticity loss for each fluid particle as It is worth noting that we have the flexibility to choose either ∇ × ūn or ∇ × u n+1 as the second term in Eqn.(20), where u n+1 is the pressure projection outcome of ūn , since the curl of the pressure gradient is zero.We adopt the difference curl (Eqn.( 6)) to discretise the velocity curl as Computing the vortex-stretching term needs the gradient of velocity ∇u, which is given by For viscid fluids, we compute viscosity with the Laplacian of velocity and vorticity using Eqn.(8).Putting it all together, Alg. 1 summarises the computation of vorticity loss.

Algorithm 1 Vorticity loss evaluation
Input: ∆t, u n Output: ω n+1 loss 1: for all SPH fluid particles i do 2: where y i are all sampled N positions measured in the abovementioned Eulerian grid -that is, much fewer than the total number of globally distributed vortex particles Np used in Eqn.(16); and P denotes the probability of the corresponding position x.
We adapt the Monte Carlo estimator to particle methods, making it feasible to correct the velocity field with the vortex particle method.In our method, rather than iterating over all N f SPH fluid particles, we stochastically sample N ≪ N f fluid particles to form a small subset Γ of the entire particle-set.Each sampled SPH fluid particle is paired with a single vortex particle as a moving unit (Fig. 4).Since implementing an infinitely small vortex particle is numerically impractical, we assign a small fixed volume v k to each vortex particle k in Γ, given by where V k is the volume of the corresponding SPH fluid particle, and c ∈ (0, 1), which we next call the volume coefficient, controls the volume of the vortex particle.
By substituting Eqns.( 12) and ( 24) into the Monte Carlo-sampled Biot-Savart formulation (Eqn.( 23)), we obtain Similarly, by converting the vorticity ω k to the smoothed vorticity loss ωloss, k , we derive the velocity correction u corr,i for SPH particle i as where N is the number of vortex particles and P k is the probability of vortex particle k, given by a stochastic distribution.In practice, we use the uniform distribution where N f is the number of fluid particles from which Γ is sampled.
By replacing A in Eqn.
(2) with the vorticity loss, we find the smoothed vorticity loss ωloss,k of vortex particle k as where j are particle k's neighbouring SPH fluid particles (see Fig. 4).
For solid particles (used to model boundaries, see Sec. 3.1), we set the vorticity loss to 0.
We leverage the velocity correction ucorr to compensate the original fluid velocity field and enhance turbulent details by setting Algorithm 2 outlines our MCVSPH approach.We start by stochastically sampling the fluid particles to yield the subset Γ (line 1).Each sampled fluid particle is linked to a vortex particle; these paired particles next move together as cohesive units.During the simulation loop, we perform five steps, making our approach work independently of pressure-projection loops: First, after applying non-pressure forces (line 4), we enforce incompressibility by methods such as DFSPH [BK17] (line 6).Next, we compute vorticity loss for each fluid particle using the vorticity equation on the current velocity field (line 8).Smoothed vorticity loss values for each vortex particle are computed by an SPH approximation (line 10).Corrective velocities for each fluid particle are computed using the Biot-Savart summation method, enhanced by the Monte Carlo technique (line 12).Finally, these corrective velocities are added to the original velocities to obtain the updated velocity field (line 13).

Results
We next present several experiments to validate MCVSPH, which we implemented using the Taichi graphics programming language [HLA*19].All computations are done on a desktop PC with a 3GHz 8-Core Intel Core i7-9700 processor, an NVIDIA Quadro RTX 4000 GPU, and 32GB memory.Rendering uses the Karma renderer of Houdini 19.5.All experiments are visualised in the accompanying video.

Effect of sample size N
We next present experiments to evaluate the impact of the sample size N (Sec.5.1.1)and the volume coefficient c (Sec.5.1.2) on the performance of turbulence generation of our method.
We evaluate the impact of different sample sizes N in our MCVSPH method by simulations of a dam break scenario with 1.2 M fluid particles and a rotating propeller at 120 r/min (see Fig. 5).
In these experiments, we set the volume coefficient c = 0.05 and vary N ∈ {100, 500, 1200, 5000}; note that N = 1200 corresponds to 1‰ of fluid particles.We visualise the results by rendering particles colour-coded by vorticity (white for high vorticity, blue for low vorticity, see Fig. 5).We see that the fluid simulation becomes noisy and unrealistic for smaller sample sizes (N = 100 and N = 500).
A visually plausible and moderately turbulent flow is obtained for N = 1200.Increasing the sample size to N = 5000 results in smoothing of turbulent details due to excessive vortex particles with lower vorticity being distributed throughout the fluid field.To further assess the effect of the sample size N, and compare our method with DFSPH, we show in Fig. 6(top) the dimensionless vorticity ω = ∥ω∥/∥ω DFSPH ∥, the dimensionless time overhead per frame t 0 = t 0 /t DFSPH 0 , and the dimensionless energy variation dE = dE/dE DFSPH .Smaller sample sizes (below 1200) lead to significantly higher vorticity values ω; larger sample sizes converge to stable values.Regarding dE, smaller sample sizes result in a much noisier and unstable fluid energy variation.Figure 6(bottom) further refines this observation by illustrating the energy variation dE evolving with time for different sample sizes.Analysing t 0 (Fig. 6 top), we see a non-linear behaviour when N < 1200, with linear evolution for N > 1200.Interestingly, the t 0 curve for N < 1200 says that the time overhead of MCVSPH is even lower than that of DFSPH.This is due to the sparser distribution of fluid particles with smaller sample sizes, leading to fewer density correction iterations in DFSPH.Once the sample size exceeds 1200 and the noise is suppressed, the number of density correction iterations aligns with typical DFSPH, resulting in linear increase in the t 0 curve with respect to N. This also demonstrates the impracticability of treating all fluid particles as vortex particles due to the heavy computational burden (approximately 63 times slower than DFSPH per frame).By exclusively iterating over a pre-sampled subset N ≪ N f of particles, we substantially reduce the computational overhead, making our method comparable with the baseline DFSPH.
From these results, we see that setting the sample size N to roughly 1‰ of the total fluid particle count yields a well-optimised fluid animation with visually realistic vorticity, reduced noise, and reasonable computation time.We use this setting of N for all further results in this paper, and indicate it next by the notation N ‰ .

Effect of volume coefficient c
To assess the effect of c, we show in Fig. 7 a dynamic scenario featuring fast streaming water with 2.3 M fluid particles, a stationary  ball, and open boundaries on the left and right sides.Fluid particles are initialised with leftward velocities of 3 m/s.We show results for different volume coefficients c ∈ {0.04, 0.05, 0.06, 0.07} for a constant sample size N ‰ = 2300.
We see that, as the volume coefficient c increases, the water flow becomes progressively more turbulent -simply put, this shows how one can control the intensity of turbulent flow via tuning c.

Validation of generating turbulent flows
We further confirm the effectiveness of our MCVSPH approach in capturing turbulent details by five simulation scenarios.We also compare our method with the micropolar SPH (MPSPH) [BKKW19], the vorticity refinement SPH (VRSPH) [LWB*21], and the divergencefree SPH (DFSPH) [BK17] methods.
Falling water column.Our first experiment simulates a falling water column with 1.1 M fluid particles (Fig. 8).The fluid, initially shaped as a box column with a height of 6.0 m and a cross-section  area of 4 m 2 , experiences free fall under gravity.We compare our MCVSPH approach (c = 0.05, N ‰ = 1100) with DFSPH, MPSPH (with the transfer coefficient νt = 0.08), and VRSPH (with the gain α = 1.0).From Fig. 8, we see that both DFSPH and VRSPH fail to capture turbulent details effectively, leading to a quick loss of such features.Also, the left image for VRSPH shows noticeable artefacts in the water column, where fluid particles clump together unnaturally in local regions.We attribute this behaviour to VRSPH's usage of the Biot-Savart summation within the support radius of each SPH particle, which causes significant local movements.We also see that both MPSPH and MCVSPH can effectively preserve turbulent details, but with differing patterns: MPSPH creates smallscale ripples, whereas MCVSPH produces visible vortical motions resulting from our vorticity-explicit method.
Dam breaking with a rotating propeller.We further show the effectiveness of MCVSPH by a dynamic scenario involving the  The graphs in Figure 10 illustrate the temporal progression of the total energy E per particle (top) and vorticity ω (bottom) over a 15 second period.The bottom graph shows that all compared methods generate a higher vorticity than the baseline DFSPH method.Notably, both our MCVSPH approach with c = 0.06 and the MPSPH method with νt = 0.15 exhibit similarly elevated vorticity levels, the highest in this comparison.However, the top graph reveals that MCVSPH yields a greater total energy, suggesting that the increased vorticity is primarily attributed to intensified vortical motions in the MCVSPH method.At a higher level, we see that our method does not score significantly negatively with respect to estimation of E or ω with respect to the compared methods.The E and ω differences are small and, as the other images in this paper show, our method produces more visually-realistic effects than the compared methods.
Fast streaming milk with a static chocolate ball.In this simulation (Fig. 11), the setup is identical to the experiment in Fig. 7.We compare DFSPH, MPSPH with νt = 0.1, VRSPH with α = 1.2, and our MCVSPH with c = 0.07 and N ‰ = 2300 to demonstrate the effectiveness of our method.Observing the left images, we see that both MPSPH and MCVSPH produce rich turbulent details, with MCVSPH showing a more prolonged persistence of turbulent flows compared to MPSPH.As the fluid gradually slows down and becomes shallower (images on the right in Fig. 11), our MCVSPH provides a more credible simulation of a tranquil shallow fluid, more similar to DFSPH as compared with the other studied methods.Table 1: Average kinetic energy E k , average total energy E, average kinetic energy percentage ε k = E k /E × 100% per second per particle, and average computational time per frame t 0 (three methods, four fluid particle counts n, first 450 frames).While it is slightly slower than the other two, our approach can produce a larger percentage of kinetic energy, which means our approach has a higher efficiency in generating turbulent motions.
1.1 (Fig. 8 Forward moving boat on tranquil water with a propeller.In this experiment, we simulate a scenario involving a boat moving forward with a propeller at 60 r/min and 2.5 M fluid particles (Fig. 12).We Dam breaking with a static Stanford dragon and two static Stanford bunnies.In our final experiment, we simulate a dam breaking scenario with 1.5 M fluid particles alongside a static Stanford dragon and two static Stanford bunnies (Fig. 13).We compare DFSPH, MPSPH with νt = 0.08, and our MCVSPH with c = 0.05 and N ‰ = 1500.Differences are primarily noticeable in the middle and right figures.The figures illustrate that DFSPH shows almost no turbulent features, while MPSPH and our MCVSPH preserve clear turbulent motions (see accompanying video).
Performance comparison.Table 1 summarises the average computational times t 0 per frame, the average kinetic energy E k , the average total energy E, and the percentage of kinetic energy ε k per second per particle of the first 450 frames for the compared methods.Our approach (MCVSPH, bottom row) exhibits a slightly lower speed than MPSPH and DFSPH, as traversals over particle samples are still required.This somewhat limits the scalability of our approach, indicating the need for research into further enhancements.However, the comparison of the percentage of kinetic energy ε k demonstrates that our approach can generate a higher proportion of kinetic energy per second.Thus, even if it is slightly slower than the other two methods, our approach can produce turbulent motions with a higher efficiency in general simulation scales.

Conclusion
We have presented a Monte Carlo vortical SPH approach for simulating turbulent flows with SPH-based frameworks.Our method works independently of pressure projection loops and flexibly integrates with existing efficient pressure solvers.We leverage the vortex particle method within an SPH-based framework to model turbulent motions.To enhance the efficiency of computing corrective velocity from vorticity loss, we use a Monte Carlo estimator to compute the Biot-Savart summation on a pre-sampled particle subset.We pair vortex particles with a randomly sampled subset of fluid particles making them move together as cohesive units.Vortex particles, which are different in volume from SPH particles, carry smoothed vorticity loss approximated using SPH.
The effectiveness of our method is mainly influenced by the sample size N and the volume c of vortex particles.Setting N to roughly 1‰ of fluid particles yields an optimal balance among turbulent details, computational efficiency, and visual quality.Additionally, increasing the volume c of vortex particles enhances turbulent effects.Our experiments show that our MCVSPH method can effectively generate and preserve visually prominent vortical motions.
While our method enhances turbulent details in SPH-based fluid simulations, several unresolved issues remain.Firstly, in terms of incompressibility and volume preservation, our method does not provide performance gains compared to DFSPH.Additionally, traversals over particle samples are still indispensable in our approach, which limits the scalability to large-scale simulations.Furthermore, as the distribution of particle samples is randomly initialized, this induces slightly stochastic fluid motions in each simulation, even if all parameter settings are the same.Apart from that, a fully robust setting of our method's parameters is still a topic to be exploredagain, similar with comparable methods.We aim to develop viable solutions for these challenges in our future research endeavors.

Figure 2 :
Figure 2: Left: Velocity field u n at time step n.Middle: Non-pressure forces push the original velocity u n into a divergent intermediate velocity ūn .Right: The pressure force mitigates the divergent (normal) component of ūn , causing turbulent detail loss.

Figure 3 :
Figure3: Left: Circulation loop.The vortex strength vector β is the circulation around position x.Middle: Stokes' theorem converts a surface integral into a volume integral.Right: Convergence to a vortex particle.The loop size is small enough so that volume integration can be approximated by the product of a constant vorticity and the volume.

Figure 4 :
Figure4: Vortex particle organisation.Every sampled SPH fluid particle (white, with volume V k ) is paired with a vortex particle (cyan, with volume cV k ) as a cohesive moving unit.The smoothed vorticity loss ωloss of the vortex particle is computed using SPH (28).

Figure 5 :Figure 6 :
Figure 5: Dam breaking scenario featuring 1.2 M SPH fluid particles and a rotating propeller, with particles colour-coded by vorticity (white for high, blue for low vorticity).The volume coefficient c is set to 0.05.Results are shown for various sample sizes N: 100 (a), 500 (b), 1200 (c), and 5000 (d).Results for N = 100 and N = 500 exhibit pronounced noise.In contrast, results for N = 1200 and N = 5000 show reduced noise levels and enhanced realism.

Figure 7 :
Figure 7: Illustration of a fast streaming water scenario with 2.3 M fluid particles.The simulation employs MCVSPH with a constant sample size N ‰ = 2300 and varying volume coefficients c ∈ {0.04, 0.05, 0.06, 0.07}.As c increases, turbulent details become more pronounced.

Figure 8 :
Figure 8: Falling water column with 1.1 M fluid particles coloured by vorticity.Comparison of DFSPH, MPSPH with νt = 0.08, VR-SPH with α = 1.0, and our MCVSPH with volume coefficient c = 0.05 and sample size N ‰ = 1100.We see that VRSPH creates unrealistic local aggregation patterns resulting from its Biot-Savart summation within local neighbourhoods.Turbulent details get lost quickly in DFSPH and VRSPH.In contrast, MPSPH generates fine ripples on the fluid, and our MCVSPH preserves vortices.

Figure 9 :Figure 10 :
Figure 9: Dam breaking with 1.2 M SPH fluid particles and a rotating propeller at 120 r/min.Comparison of DFSPH, MPSPH with νt = 0.15, VRSPH with α = 1.0, and our MCVSPH with c = 0.06 and N ‰ = 1200.MCVSPH yields more visible turbulent motions than the other three methods.

Figure 11 :
Figure 11: Fast streaming milk with 2.3 M fluid particles and a static chocolate ball (same scenario settings and layout as Fig. 7).Comparison of DFSPH, MPSPH with νt = 0.1, VRSPH with α = 1.2, and our MCVSPH with c = 0.07 and N ‰ = 2300.On the left, our MCVSPH approach generates prominent vortical details on fluid surfaces.As the fluid slows down and becomes shallower, as shown on the right, our MCVSPH yields a realistic simulation of tranquil shallow liquid similar to DFSPH.
Figure 12: A small boat moving forward with 2.5 M vorticity colour-coded fluid particles and a propeller at 60 r/min.Comparison of DFSPH, MPSPH with νt = 0.05, and our MCVSPH with c = 0.057 and N ‰ = 2500.DFSPH shows basic wake flow from the boat rear.MPSPH yields a more detailed flow.Our MCVSPH yields explicit visible transportation, merging, and splitting of vortices.

Figure 13 :
Figure 13: A dam breaking scenario with 1.5 M vorticity colour-coded fluid particles, a static Stanford dragon, and two static Stanford bunnies.Comparison of DFSPH, MPSPH with νt = 0.08, and our MCVSPH with c = 0.05 and N ‰ = 1500.Our MCVSPH approach produces evident vortical movements within the fluid domain, which is visible in the middle and right images.
compare three methods: DFSPH, MPSPH with νt = 0.05, and our MCVSPH with c = 0.057 and N ‰ = 2500.DFSPH generates basic wake flows from the rear of the boat.MPSPH adds more details to the flow.Our MCVSPH approach clearly shows the transportation, merging, and splitting of vortices driven by the propeller at the rear of the boat (see the accompanying video).