# Flow #1: A Smooth Motivation to Fluid Dynamics

**Pedro Ferro Pereira**Co-author

**Fábio Cruz**Co-author | Technical Reviewer

**David Carvalho**Co-author | Editor

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

exactlymean by afluidand itsdynamics?

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:

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.:

### 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:

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:

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**:

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

## Recent posts from our blog

Hugo Penedones

Luís Sarmento

The Inductiva API v0.4 release brings MPI clusters, the latest Google Cloud CPUs, two new simulators, a lighter Python package, a CLI interface, a template engine and totally revamped documentation. Get started in minutes!

Sofia Guerreiro

Cristiana Carpinteiro

In this series of blog posts we will explore a specific case of the use of AI in the pharmaceutical industry - using Graph Neural Networks for predicting binding affinity. But for now, let’s start by understanding the problem of drug discovery and some fundamental concepts like binding affinity.