The past few months at Inductiva have been of great development in the computational infrastructure necessary to perform realistic simulations of fluids. A significant part of our focus has been on simulating fluids with Smoothed Particle Hydrodynamics (SPH).

We are particularly interested in unleashing the potential SPH offers when scaling its capacity to generate high-throughput flow data with high resolution. On the computational side this is a tremendous task but, if properly represented and implemented, can be tamed with the power of… Machine Learning.

Granted, a whole lot of words here. Allow us to go step-by-step.

The Flow Series
A Journey on ML for Fluid Dynamics Simulation using SPH

In this series, we will dive deep into how Machine Learning (ML) models can be built to learn Fluid Dynamics. To this end, we will use a simple physical setup of a dam break as our case study and simulate the corresponding water dynamics with SPlisHSPlasH [1] — an open-source state-of-the-art SPH solver.

We will slowly build up all the tools – from SPH fluid simulations to model architecture optimization — to get our very SPH ML model running mighty fine.

• In #1 (this post): First, fluids & thoughts. We lay the foundations of Fluid Dynamics and motivate why it is difficult to solve their governing equations, the Navier-Stokes Equations (NSEs), even if numerically.

Alright. 😏 Let’s get smooth in here! 😏

Flow #1:
A Smooth Motivation to Fluid Dynamics

The Complex World
… of Fluids

The development of technology and infrastructure of much of what we consider mundane these days could not subsist without a thorough understanding of physical phenomena involving fluids.

Whether it is air flowing through a wing that maintains a plane airborne, hydroelectric energy being generated by a dam, fans ventilating an IT server with cool air or natural gas flowing along pipelines, Fluid Mechanics affords us a framework where we can describe such events with the aid of deterministic equations that govern these physical phenomena.

Fig. 1: Given an initial condition of a fluid (on the left), physical laws dictate its evolution. After a short period of time, extremely complex behaviour may be observed (on the right). Credits: 3Blue1Brown

Modelling adequately fluids is, in general, a pretty challenging task. Fluids can have extremely complex structure - think of honey or a chemical dye. Moreover, fluids can show intricately detailed behaviour (even if their structure is not) — shedding across a wing plane or the turbulent flow around a bonfire [2]. In Fig. 1 you can see the resulting smoke dispersion in air — wow!

Before we move on though, we must answer the question:

What do we exactly mean by a fluid and its dynamics?

Time to brief you on your Fluid Dynamics 101.

A Primer
on Fluid Dynamics

Think of a glass of water: it will contain an estimated \(4 \times 10^{24}\) atoms in! As this unimaginable number of particles drift around and interact with each other, each particle will change its motion according to well-established laws.

Our go-to dynamical variable is the (linear) momentum \(\mathbf{p}\). Each atom \(i\), of mass \(m_i\), stationed at position \(\mathbf{x}_i\) and moving with a velocity vector \(\mathbf{v}_i\), carries it as the vector:

\[\mathbf{p}_i = m_i \mathbf{v}_i\]

It simply is not just prohibitively expensive — it is unfeasible to know the momentum of all individual atoms or molecules in large domains.

Anyway, what would we accomplish if such information were known? At face value, not much. We need a notion of medium (and not the individual constituent particles). As a fluid, specific interactions between the particles — chemical, electromagnetic or even thermal — dictate the global properties of the ensemble.

The basic principle behind Fluid Dynamics is that materials are approximated by a continuum. In this framework, the contributions stemming out of these discrete particle interactions are lumped into a field which quantifies the aggregate role of all particles.

In this fashion, discrete information about the particle momenta \(\mathbf{p}_i\) can lead us to the fluid momentum \(\mathbf{p}(\mathbf{x})\). Naturally, as interactions take place, these fields can also change, becoming time-dependent i.e. \(\mathbf{p}(\mathbf{x},t)\).

The greatest advantage? We can now access fluid properties at any location and time and see how they evolve.

But how do they?

Towards the Holy Grail of Fluid Dynamics:
the Navier-Stokes Equations (NSEs)

A tale of two forces in (hydro)static equilibrium

Consider a column of water in hydrostatic equilibrium. All the forces on the water are in balance and the water is motionless.

On any drop of water, two forces are balanced: gravity, a body force which acts directly on each atom and molecule inside and a surface force exerted by the surrounding water.

Consider the fluid element in grey in Fig. 2. It has a mass \(\delta m\), a height \(\delta h\), base area \(A\) and volume \(\delta V = A \delta h\).

Fig. 2: A fluid element inside a column of water in hydrostatic equilibrium. This motionlessness state can be reached whenever the pressure-gradient forces \(\mathbf{F}_{\rm up}\) and \(\mathbf{F}_{\rm down}\) balance the gravitational force \(\mathbf{F}_g\). Credits: David Carvalho / Inductiva

In addition to its gravitational force \(\delta \mathbf{F}_g\), the fluid sets up a combined force due to the water above and below the element, respectively \(\delta \mathbf{F}_{\rm up}\) and \(\delta \mathbf{F}_{\rm down}\). Let’s analyse what has to happen in the vertical direction for this motionlessness state to be reached.

The gravitational force is simply \(F_g = m g\), where \(g\) is the gravitational acceleration. Using the element density \(\rho = \delta m/ \delta V\), the force element is \(\delta F_g = g \rho A \delta h\).

Now, the force above is simply \(\delta F_{\rm up} = P A\), whereas the force below is \(\delta F_{\rm down} = (P + \delta P) A\). So the fluid force in this direction is \(\delta F_{\rm up} - \delta F_{\rm down} = -A \delta P\).

For this force to balance the gravitational force, we find that \(\delta P / \delta h = \rho g\). In the infinitesimal limit, this is a condition for the partial derivative \(\partial P/\partial z\)!

In this case, since our acceleration is entirely vertical, so is this pressure-driven force. Since the pressure does not change in any other directions, the contribution of these coordinates to this force must be zero.

The effect of the gradient of pressure \(P\) throughout the flow is to accelerate particles from high to low pressure, resulting in a pressure-gradient force \(\mathbf{F} = -\nabla P\) per unit volume on the surface. On any volume element, these forces must then be balanced with the total acceleration vector \(\mathbf{a}\) i.e.:

\[-\nabla P + \rho \ \mathbf{a} = 0\]

A dropping drop (out of equilibrium)

If the forces are not balanced, the droplet accelerates.

We know that the acceleration \(\mathbf{a}\) is the rate of change \(\mathbf{a} = d \mathbf{v}/dt\). To compute this, we cannot simply take the partial derivative \(\partial \mathbf{v} / \partial t\) because the flow changes in space as it changes in time! In essence, the material derivative [3] considers these contributions as the differential operator:

\[\frac{D}{Dt} \equiv \frac{\partial }{\partial t} + v \cdot \nabla\]

Intuitively, this derivative is taken following the motion of the particle. The interesting consequence is that a new term — quantifying advection [4] shows up.

What, then, is this nonzero acceleration felt by the droplet? If we apply the material derivative to the fluid momentum, each unit volume experiences the force:

\[\rho \frac{Dv}{Dt} = \rho \left( \frac{\partial v }{\partial t} + (v \cdot \nabla) \ v \right)\]

… we obtain the net force acting on the droplet! The second contribution, which originated from the advection term of the derivative, is the convective acceleration [5].

Gettin’ sticky: adding viscosity

However, other types of forces can still be exerted.

Fluids do not flow entirely freely — viscous forces generated throughout it drag particles, making them resist motion. Fluids become viscous whenever they withstand deformation forces. When deformed, the fluid will be moving faster in one region than in another, developing a velocity gradient — and hence a force.

In the simplest case, the viscous stress response of the fluid \(\mathbf{\tau}\) [6] is linearly proportional to the rate of the deformation. The proportionality constant is the fluid viscosity \(\mu\). For example, if \(v_x\) varies with \(z\), a force in the direction \(x\) develops and is proportional to \(\tau_{zx} = -\mu (\partial v_{x}/\partial z )\).

Whenever this approximation is valid, the fluid is Newtonian. This relation reduces hugely the possible channels between the deformation tensor and the viscous stress tensor! Fortunately for us, most everyday fluids – water, air and even glycerine — are Newtonian.

Okay - we are now ready to assemble our soon-to-be favourite set of equations.

Enter the Navier-Stokes Equations

Proposed in 1822 by Claude-Louis Navier (no one other than the father of modern Structural Analysis) and Sir George Stokes (yes, the one from Stokes’ Theorem), the Navier-Stokes Equations (NSEs) are a set of fundamental equations which describe the flow of everyday fluids.

Video 1: A great motivation to the inherent complexity of turbulent flow (alongside mesmerising vortex shedding in smoke with lasers!) with the Navier-Stokes Equations. Credits: 3Blue1Brown

The NSEs are a fluid analogue of Newton’s 2nd Law of Motion, setting how the fluid momentum changes due to all forces \(\mathbf{F}_i\) exerted on it:

\[\mathbf{F} = m \mathbf{a} = \frac{\partial \mathbf{p}}{\partial t} = \sum_{j} \mathbf{F}_j\]

For a Newtonian fluid, the NSEs look like this:

\[\underbrace{\rho \frac{D \mathbf{v}}{D t}}_{\rm 'total \ force'} = -\underbrace{\nabla P}_{\rm pressure} + \underbrace{\mu \nabla^2 \mathbf{v}}_{\rm viscosity} + \underbrace{\rho \mathbf{g}}_{\rm gravity} + \underbrace{\mathbf{F}^{\rm ext}}_{\rm external \ forces}\]

where we lumped all the intrinsic contributions we discussed before plus any other external force \(\mathbf{F}^{\rm ext}\) which may be impinged on the fluid.

Additionally, the fluid total mass does not change over time. As fluid flows, the rate at which its mass enters a volume element must then be equal to the rate at which it leaves it plus any accumulation of mass within that volume. This notion is expressed with the continuity equation:

\[\frac{D \rho}{D t} + \rho (\nabla \cdot \mathbf{v}) = 0\]

Five variables are present in these two sets of equations: three components of the velocity \(\mathbf{v}\), the pressure \(P\) and the density \(\rho\). With only four equations available, a fifth one must be introduced.

A reasonable condition is to apply a type of equation of state, which can relate the pressure to the density on general principle considerations e.g. a linear relation \(P(\rho) = k \rho\), for a stiffness factor \(k>0\).

All in all, the NSEs comprise a coupled system of nonlinear Partial Differential Equations (PDEs) whose solutions (probably) cannot be found analytically (except for a couple of simplistic scenarios).

Unsurprising. The task of solving the NSEs requires the estimation of five functions of the form \(f(x,y,z,t)\) which must always satisfy these conditions:

\[\small{ \left\{ \begin{aligned} & \rho \frac{\partial v_x}{\partial t} + v_x \frac{\partial v_x}{\partial x} + v_y \frac{\partial v_x}{\partial y} + v_z \frac{\partial v_x}{\partial z} + \frac{\partial P}{\partial x} - \mu \left( \frac{\partial^2 v_x}{\partial x^2} + \frac{\partial^2 v_x}{\partial y^2} + \frac{\partial^2 v_x}{\partial z^2} \right) - F_x^{\rm ext} = 0 \\ & \rho \frac{\partial v_y}{\partial t} + v_x \frac{\partial v_y}{\partial x} + v_y \frac{\partial v_y}{\partial y} + v_z \frac{\partial v_y}{\partial z} + \frac{\partial P}{\partial y} - \mu \left( \frac{\partial^2 v_y}{\partial x^2} + \frac{\partial^2 v_y}{\partial y^2} + \frac{\partial^2 v_y}{\partial z^2} \right) - F_y^{\rm ext} = 0 && \\ & \rho \frac{\partial v_z}{\partial t} + v_z \frac{\partial v_z}{\partial x} + v_z \frac{\partial v_y}{\partial y} + v_z \frac{\partial v_y}{\partial z} + \frac{\partial P}{\partial y} - \mu \left( \frac{\partial^2 v_z}{\partial x^2} + \frac{\partial^2 v_z}{\partial y^2} + \frac{\partial^2 v_z}{\partial z^2} \right) - \rho g - F_y^{\rm ext} = 0 && \\ & \frac{\partial \rho}{\partial t} + v_x \frac{\partial \rho}{\partial x} + v_y \frac{\partial \rho}{\partial y} + v_z \frac{\partial \rho}{\partial z} + \rho \left( \frac{\partial v_x}{\partial x} + \frac{\partial v_y}{\partial y} + \frac{\partial v_z}{\partial z}\right) = 0 \\ & P = P(\rho) \end{aligned} \right. }\]

Ouch - what a nested kerfuffle this is to solve!

Gaining analytical insights of Fluid Dynamics seems notoriously hard. In fact, making substantial progress towards a mathematical theory which will unlock the secrets hidden in the Navier-Stokes equations could earn you a million dollars, as stated by the Millennium Prize Problems proposed by the Clay Mathematics Institute [7]

However, we can live with that.

Numerical methods to the rescue!

The Computational Fluid

The NSEs can be estimated numerically by applying algorithms that converge trial solutions to the true result. Even in the (simpler) Newtonian regime, this beastly system of equations needs to be solved across the entire domain and for a set time window.

For instance, take a look at the velocity field of the fluid once it moves past a circular obstacle that we obtained by simulating the NSEs for incompressible flow (constant density \(\rho\) or, equivalently, \(\nabla \cdot \mathbf{v} = 0\)), using FEniCS, an open-source in Python/C++ computing platform for solving PDEs with efficient Finite Element Methods:

Video 2: 2D Flow past a circular obstacle obtained by simulating the incompressible NSEs, where \(v_x(x,y,t)\) (top), \(v_y(x,y,t)\) (centre) and velocity magnitude \(|\mathbf{v}(x,y,t)|\) (bottom) are shown. Credits: Fábio Cruz / Inductiva

If the numerical routines estimate these functions unsatisfactorily at some points \(x,y,z\) and time instants \(t\), the trial solutions \(\mathbf{v}(x,y,z,t)\), \(\rho(x,y,z,t)\) and \(P(x,y,z,t)\) may look horribly wrong!

Lesson learnt — the simulation route is riddled with potholes and we must be careful. In order to have guarantees the estimate is sound, we must use robust and stable numerical routines which are able to resolve the physics taking place without introducing any spurious feature or artefact.

At large, the field of Computational Fluid Dynamics (CFD) aims at solving these challenges. Most well-established methods in CFD tend to be mesh-based (e.g. ANSYS’s proprietary software or OpenFOAM’s open-source libraries), meaning the spatiotemporal domain is discretized and quantities of interest estimated at those nodes only.

As we have seen in the Heat series, where we solved the Heat Equation in 2D using a Finite Element Method, discretizing the system of PDEs in a mesh can work beautifully or fail spectacularly.

Video 3 & 4: For the same numerical routine to solve a system of PDEs, the output can approximate the “true” solution (left) or fail spectacularly (right) depending on a huge combination of factors. Credits: Manuel Madeira / Inductiva

That was for an admittedly simple case. As we move to 3D, or to include more complex geometries and custom boundaries and interfaces, the seemingly simple task of subdividing the domain appropriately can become a headache once over/undersampling and singularities occur.

Mesh generation is an expensive and time-consuming task which sets the bar on how good the estimates can potentially be. The best mesh-based CFD solver will always fail whenever an inappropriate mesh is used.

Not only that. Think of a pipeline carrying oil — the pumping is accomplished by mixing two fluids in different phases (gaseous and liquid) in a fine-tuned way. As more complex, heterogeneous and multi-phase media need to be simulated, the framework used must be capable of resolving their combined flow.

Video 5: A simulation of a lava lamp – a multi-phase flow using SPH. From Müller et al., Particle-Based Fluid-Fluid Interaction, SCA 2005. Credits: Barbara Solenthaler

In the
next post…

… we will introduce Smoothed Particle Hydrodynamics (SPH) as a mesh-free, particle-based alternative framework to simulate fluid flow (but not only!).

Starting from its core principles, we will state the technical features that allow SPH to compute the NSEs. After that, we will solve the physical scenario for which we will train our ML model and see how well it fares when benchmarked against real data.

🌊🌊🌊 Until then, flow away! 🌊🌊🌊

---

References & Remarks

[1] The official documentation page of the SPlisHSPlasH library.

[2] Wired: The Weird and Beautiful World of Fluid Dynamics: a beautiful gallery containing ten different physical phenomena that happen within the realm of fluids.

[3] A bit more on the mathematical formulation of the material derivative.

[4] The advection operator that enters the material derivative is \(v \cdot \nabla = v_x \frac{\partial}{\partial x} + v_y \frac{\partial}{\partial y} + v_z \frac{\partial}{\partial z}\)

[5] The convective acceleration \((\mathbf{v} \cdot \nabla) \mathbf{v}\) is a vector of operators whose elements are \([(\mathbf{v} \cdot \nabla) \mathbf{v}]_i = v_j \frac{\partial v_i}{\partial x_j}\)

[6] For a general flow velocity \(\mathbf{v}(\mathbf{x})\), the stress tensor \(\tau\) has nine components i.e.

\[\tau = \begin{bmatrix} \tau_{xx} & \tau_{yx} & \tau_{zx} \\ \tau_{yx} & \tau_{yy} & \tau_{yz} \\ \tau_{zx} & \tau_{zy} & \tau_{zz} \end{bmatrix}.\]

Using \(\mathbf{x} = (x_1, x_2, x_3)\), and assuming the fluid is isotropic (meaning the viscosity \(\mu\) is the same along each direction), each component is simply \(\tau_{ij} = - \mu \dfrac {\partial v_{j}}{\partial x_i}\)

[7] Clay Mathematics Institute: Navier–Stokes Equation