Cedric THIEULOT

SPHENE

back

SPHene is a CFD code developed in 2007. At its core is the Smoothed Particle Hydrodynamics method developped in the late 70's independently by Monaghan and Lucy. This method relies on a set of discrete points called particles that sample the fluid and carry a given mass, density, velocity, acceleration ... From these values it is possible to compute spatial gradients and higher order derivatives, such as the pressure gradient, the strain-rate or Laplacian of the velocity for every particle. In order to do so, each particle only 'sees' a given subset of its neighbours that are within a cut-off distance and attributes them a given weight by means of a kernel function.

This method is explicit in time and once the acceleration of each particle is known, it is advected and from the new spatial distribution of these particles one can compute new accelerations, and so forth.

In its standard and commonly described form, the SPH method suffers from many drawbacks, such as boundary conditions implementation, near free surface gradient estimations, conditional stability, lack of high order consistency, non unicity of derivatives discretisation, etc ... SPHene has been primarily designed as a modular code for the testing of various implementations of existing algorithms and/or new ones. This code is therefore two-dimensional only: it would be trivial to produce a three-dimensional version but memory requirements and computational time would rise up with yet no substancial benefit for the understanding of the observed phenomena. Moreover it is for the same reason sequential, even though the parallelisation of some 'bottleneck' routines would most certainly improve its performances while making it impossible to run on my old G4 laptop.

I wish to stress out that this code is currently under development and that the simulations presented on this page are by no means 'perfect'. I nonetheless show these results to later compare them to the ones obtained with improved versions of the algorithms. SPHene has been calibrated for very simple numerical experiments (Poiseuille and Couette flows for instance) only, and what follows is simply an illustration of the type of problems a purely Lagrangian code such as SPHene can tackle. It is also worth noticing that besides the specification of material properties and resolution in the input file of each of these experiments, the only necessary coding required to run them is the specification of the boundary particles and the prescribed boundary conditions, like any other CFD code.

The visualisation is based on the PGPLOT library. Convertion from postscript to png is done w ith Imagemagick, and from images to movies by means of Quicktime for the time being.

I wish to thank Mathieu Labbe for his many improvements and suggestions, the coding of the linked-list neighbour finding algorithm, the coding of the slip/no-slip boundary conditions, and for useful discussions and comments. Some figures displayed on this page come from his 3-month internship report (Spring and Fall 2007). You can download

A good starting point in the world of SPH is the book by Liu and Liu, Smoothed Particle Hydrodynamics: a meshfree particle method , at World S cientific. I must admit I do not like this book very much (structure, relevance of some topics, choice of illustrations) even though I cannot deny the wealth of information this book provides. As complementary reading, I would advise: I am a member of the spheric community. The goal of SPHERIC is to foster the spread of the SPH simulation method within Europe and abroad. It forms a framework for closer co-operation between research groups working on the sub ject and serves as a platform for the information exchange from science to industry. One of the most important goals is the assessment of this method for all its possible applications and its development. I present hereafter various typical (and less typical) fluid mechanics experiments:


Couette flow at low Reynolds number

The flow takes place between two vertical walls. The left wall is fixed while the right wall has a constant velocity v=-10-4m/s. The fluid has a density ρ=1000kg/m3, a viscosity μ=0.001Pa.s and there is no gravitational force. No-slip boundary conditions apply on both walls, so that the right wall is the only driving force of the flow. Vertical periodic boundary conditions are implemented.

On the following snapshots a-e) of the system at different times, a row of particles was colored differently than the others from the bulk. It is then possible to observe how this line first curves as the velocity profile hasn't reached its stationary state:

a) b)
c) d)
e) f)

One can also colour the particles with a color scale indicating their velocity, as shown on the following a-d) plots:

a) b) c) d)

The non-stationary and stationary analytical velovity profiles are given by

so it is possible to plot on the same graph these profiles at different times, and the measured ones on the particles:

The error between the measurements and the theoretical curves is found to be less than 0.3%.


Poiseuille flow at low Reynolds number

It is a simulation of the flow of a viscous fluid in a channel. The fluid has a density ρ=1000 and a viscosity μ=0.005. The channel is 1cm wide and 5mm high. Periodic boundary conditions in the y direction are applied. Rows of darker colored particles are placed so as to facilitate the observation of the flow, and especially its onset. There are 20000 particles in the simulation (200 in the x direction, 100 in the y direction). No slip boundary conditions apply on lateral walls. The Earth gravity is the only driving force acting on the fluid.

a) b)
c) d)
e) f)

One can also colour the particles with a color scale indicating their velocity, as shown on the following a-d) plots:

a) b) c) d)

The non-stationary and stationary analytical velovity profiles are given by

so it is possible to plot on the same graph these profiles at different times, and the measured ones on the particles:

The error between the measurements and the theoretical curves is found to be less than 0.8%.


The dam break problem

This is a classical fluid mechanics benchmarks, see for instance spheric page (Validation tests). The box is 20cm long, 10cm high. The initial column of water is 6.5cm wide and 8.3cm high. The density of the fluid is ρ=1000kg.m-3 and its viscosity μ=5 Pa.s . The gravity vector points vertically downwards with g=9.81m.s-2. The initial particle distribution is 400x200, and only the ones in the water column are kept in the system, i.e. 23111 particles.

the collapse of a viscous fluid column

A bloc of viscous fluid (μ=10Pa.s) is placed on a 45o inclined slope (this is achieved by modifying the gravity acceleration vector). Periodic boundary conditions in the y direction are implemented, and after some time, the system evolves towards a fluid flow with a free surface, as shown on figures a-l).

a) b) c) d)
e) f) g) h)
i) j) k) l)

The expected stationary flow velocity profile is given by the following formula:

and one can plot the measured velocity profile on the SPH particle with the theoretical profile:

The dispersion error is approximately 5%.


Shortening of a viscous fluid

A (really) viscous fluid (μ=100 Pa.s, ρ=1000kg.m3) is placed in a 2.5m x 1m tank up to a height of 50cm under normal gravity. The right wall is pushed inwards at a steady speed of 5cm.s-1. No-slip boundary conditions apply on all walls. One sees the unexpected apparition of some buckling at the surface (figures c,d,e), probably due to some (numerical) instability near the free surface, but this is later accomodated by the shortening and the fluid displays an expected horizontal free surface shortly after, and until the fluid level reaches to point of overflowing the tank.

a) b) c)
d) e) f)
g) h) i)


The splash of a liquid drop

We consider a droplet of fluid falling into a fluid-filled basin. Both droplet and bulk are from the same fluid of density ρ=1000kg/m3 and viscosity μ=50Pa.s. The gravity field is vertical (g=9.81m/s2). The droplet falls under its own weight with no initial velocity. The bassin is 10cm large and the height of the water is 5cm. The drop has a radius of 1.25cm and is placed at x=5cm, and y=6.5cm.


The lid-driven cavity

This is a simulation of a driven cavity flow in a square domain. The physical configuration consists of a square container filled with a fluid. The lid of the container moves at a given, constant velocity ν, thereby setting the fluid in motion. No-slip conditions are imposed on all four segments of the boundary. The used parameters are ρ=1000.0, μ=1.0, ν=1, Lx=Ly=1 which give Re=1000.


The rotating drum

A fluid ( ρ=1000 kg/m3, μ=0.1 Pa.s) ) is placed in a two-dimensional rotating drum of diameter D=25cm, with a counter-clockwise angular velocity ω= 1 rad.s-1. No-slip boundary conditions are imposed on the inside walls of the drum. Because of its relatively high viscosity, the fluid tends to stick to the walls and its free surface is higher on the right part of the simulation domain than on the left due to the rotation. On the other hand, one expects the free surface of the fluid to remain flat. It is indeed what we observe: fluid particles that are taken along the rotating drum on the right finally flow inwards at the surface to restore the free surface flatness. Colours are meaningless: they simply help the visualisation of fluid flow.