jns

3D Navier-stokes solver for incompressible fluids using Java 8 for regions including obstacles and surface.

Status: pre-alpha

Interested in this generally having seen water movement models in use at Australian Maritime Safety Authority to predict the location of people that have fallen in the ocean for example.

Would like to document in java 8 (enjoying the conciseness of lambdas) a basic implementation of a Navier-Stokes equations solver for an incompressible fluid in 3D using the pressure-correction method. See Architecture below for more specific aims.

In terms of reference material I've found very detailed mathematical discussions of the pressure-correction method or high level descriptions and nothing in between. Ideally I'd like some pseudo code describing the algorithm with enough detail to code from. Still looking and in the meantime am fleshing out my best guesses.

Update: These references look useful:

• Fractional step methods for the Navier Stokes equations on non-staggered grids by Armfield and Street. I might run with this shortly.
• Numerical Fluid Mechanics MIT lectures
• Computational Fluid Dynamics lectures Uni Michigan might have sufficent explanation of the SIMPLER algorithm to run with also.
• The Simple Algorithm lecture from Purdue Uni
• Computational Methods for Fluid Dynamics by Ferziger and Peric
• 2D solver in mathematica
• A finite volume method to solve the 3D Navier-Stokes equations on unstructured colocated meshes
• Validation using Taylor-Green Vortex or Lid Driven Cavity Problem

Aims:

• decouple the algorithm, the mesh, derivative methods, root solver so that the program can be altered with ease.
• the algorithm should be clearly visible in the code
• minimize pollution of code with data structure choices (array iterations for example)
• accept performance degradation arising from the decoupling but seek to later leverage concurrency possibly in a distributed fashion.
• consider using lazy computation

Unit tests will be created for regular grid meshes.

The Solver class is the main entry point.

Navier-Stokes Equations

The momentum equation is:

    ρ( δv/δt + (v ⋅ ∇)v ) = -∇p + μ∇²v + f

The governing equation for pressure computation (derived via conservation of mass) is:

    ∇²p = -∇ ⋅ (v ⋅ ∇)v

where

    ρ = fluid density

    v = velocity vector

    p = pressure

    t = time

    f = other forces vector (for example gravity)

    μ = viscosity (1 / Reynolds number)

Explanatory notes

Note that for the term f above in the case of gravity this is a vector with z coordinate equal to the force exerted by gravity on a cubic metre of fluid. This is thus

    f_z = ma = density * volume * accelerationGravity
        = 1025kg/m³ * 1m³ * -9.81m/s² using an approximation for the density of salt water
        = -10055.25kgm/s²

So in this case f = (0,0,-10055.25) assuming gravity is the only extraneous force.

∇ is the Gradient operator and is defined by:

    ∇f = ( δf/δx, δf/δy, δf/δz)

∇² is the Laplacian operator and is defined by:

    ∇²f = δ²f/δx² + δ²f/δy² + δ²f/δz²

For a vector v,

    ∇²v = (∇²v_x,∇²v_y,∇²v_z)

v ⋅ ∇ is the Jacobian matrix of all first order partial derivatives of velocity:

    δv_x/δx  δv_x/δy  δv_x/δz
    δv_y/δx  δv_y/δy  δv_y/δz
    δv_z/δx  δv_z/δy  δv_z/δz

The Divergence operator is defined by:

    ∇ ⋅ v = δv_x/δx + δv_y/δy + δv_z/δz

Computation of derivatives

Numerical approximations for the derivatives (first and second) of a function f are given by:

  Suppose the point x has two close neighbours a,b (a<x<b)

  Then

    f '(x) ≈ (f(b) - f(a)) / (b - a)

    f ''(x) ≈ (f(b) + f(a) - 2f(x)) / (b-a)²

A second order approximation to the first derivative is:

    f '(x) ≈ (h₂² - h₁²)f(b) + h₁²f(c) - h₂²f(a))/(h₁²h₂ + h₁h₂²)

where

    h₁ = b - a

    h₂ = c - b

Lazy computation

What would be nice is the ability to only calculate what we need for the output. For instance in a 10x10x10 grid if I want to know the value of velocity at (5,5,5) after 2 time steps then clearly the whole grid does not have to be computed for the two time steps.

When the full grid is computed, rather than map-reduce (which might be the best bet for distributed processing) seek to enable Rx to improve performance.