Computational Fluid Dynamics Fluids Fluid materials Smoke like

Computational Fluid Dynamics

Fluids �Fluid materials �Smoke like phenomena �Clouds �Dyeing

State of the fluid �Velocity (vector field) x = (x, y) position u = (u, v) velocity T time u(x, t) = (u(x, t), v(x, t)) Cartesian grid

Navier-Stokes equations (1822) �Claude Navier and George Gabriel Stokes Movements of fluids �Basic assumptions Two components of stress in fluids ▪ Diffusion proportional to the gradient of the velocity ▪ Pressure

Navier-Stokes equations (1822) �Can describe several physical phenomena Weather Flows in ducts with non-circular cross-section Flows around the wings of aircrafts Movements of solid bodies through fluid like materials (movements in stars in galaxies) Can be connected to the Maxwell equations (Magnetohydrodynamics)

Navier-Stokes equations (1822) �It is also important from theoretical point of view Existence and smoothness of Navier-Stokes solutions in 3 D One of the seven Millennium Problems by Clay Mathematics Institute ▪ One million dollar for the solution

Navier-Stokes equations (1822) density Viscosity (kinematic) External forces Incompressible, homogeneous fluids

Parts of the equations I. �Advection �Forward movement, transport �Any quantity �Own vector field as well

Parts of the equations II. �Pressure �Force transmission is not instaneous in fluids �Molecule collision creates pressure �It causes acceleration

Parts of the equations III. �Diffusion �Different fluids move differently: there are more dense and more viscous ones �Viscosity: the resistance of the fluid against the flow �This resistance causes velocity diffusion

Parts of the equations IV. �External forces �Can be local or global ones (e. g. gravity)

Operators Operator Gradient Divergence Laplace Definition Finite difference form

Solving the equations �Analytic solution only in special (simple) cases �Numeric method, incremental solution �If we would like to animate the flow the time integration could be useful �The problem is divided into smaller parts (Stam, J. 1999. "Stable Fluids. " In Proceedings of SIGGRAPH 1999)

Helmholtz-Hodge decomposition (Projection step) �(A vector can be written as a weighted sum of basis vectors) �A vector field can be written as a sum of a divergence-free vector field and gradient of a scalar field

Helmholtz-Hodge decomposition (Projection step) W

How to calculate pressure? Poisson equation

What to do in a step of the simulation? Projection External Force Diffusion Advection

Advection �Euler method, forward step �Not stable, hard to implement in Open. CL �A solution is the backward step :

Advection

Diffusion �Explicit solution: �Not stable! Implicit solution: Poisson equation

Projection Poisson equation

Solution of Poisson equations �Iterative solution, starting from an initial state and progressive refinement �Form of the equation: �Where is the Laplace operator �The most simple solution is the Jacobi method

Jacobi method Diffusion Pressure x velocity (u) pressure (p) b velocity (u) divergence of the velocity field -1. 0 0. 25

Boundary conditions �Calculations are performed only in a finite volume, boundary conditions are needed �If the fluid is surrounded with solid walls, the conditions for the velocity and pressure are: Velocity: it is zero at the boundaries (no-slip condition) Pressure: its gradient is zero at the boundaries (Neumann condition)

Extension: vorticity �The simulation and the error of the numeric discretization makes some details of the flow vanish �Trying to get it back:

Implementation �Quantities are stored in 2 D arrays �During the calculations sometime information is needed form the neighborhood, therefore some quantities are double buffered �Arrays are updated by Open. CL kernels �Separate kernels are needed for different steps of the calculation �The display is a simple drawing to the screen �Kernel functions. . .

Advection __kernel void advection(const int grid. Resolution, __global float 2* input. Velocity. Buffer, __global float 2* output. Velocity. Buffer) { int 2 id = (int 2)(get_global_id(0), get_global_id(1)); if(id. x > 0 && id. x < grid. Resolution - 1 && id. y > 0 && id. y < grid. Resolution – 1){ float 2 velocity = input. Velocity. Buffer[id. x + id. y * grid. Resolution]; float 2 p = (float 2)((float)id. x - dt * velocity. x, (float)id. y - dt * velocity. y); output. Velocity. Buffer[id. x + id. y * grid. Resolution] = get. Bil(p, grid. Resolution, input. Velocity. Buffer); } else{ //határfeltételek if(id. x == 0) output. Velocity. Buffer[id. x + id. y * grid. Resolution] = - input. Velocity. Buffer[id. x + 1 + id. y * grid. Resolution]; . . . } }

Divergence __kernel void divergence(const int grid. Resolution, __global float 2* velocity. Buffer, __global float* divergence. Buffer) { int 2 id = (int 2)(get_global_id(0), get_global_id(1)); if(id. x > 0 && id. x < grid. Resolution - 1 && id. y > 0 && id. y < grid. Resolution - 1){ float 2 v. L v. R v. B v. T = = velocity. Buffer[id. x + + + 1 + id. y * (id. y - 1) (id. y + 1) grid. Resolution]; * grid. Resolution]; divergence. Buffer[id. x + id. y * grid. Resolution] = 0. 5 f * ((v. R. x - v. L. x) + (v. T. y - v. B. y)); } else{ divergence. Buffer[id. x + id. y * grid. Resolution] = 0. 0 f; } }

Calculating pressure, Jacobi method __kernelvoid pressure. Jacobi(const int grid. Resolution, __global float* input. Pressure. Buffer, __global float* output. Pressure. Buffer, __global float* divergence. Buffer) { int 2 id = (int 2)(get_global_id(0), get_global_id(1)); if(id. x > 0 && id. x < grid. Resolution - 1 && id. y > 0 && id. y < grid. Resolution - 1){ float float alpha = -1. 0 f; beta = 0. 25 f; v. L = input. Pressure. Buffer[id. x - 1 + id. y * grid. Resolution]; v. R = input. Pressure. Buffer[id. x + 1 + id. y * grid. Resolution]; v. B = input. Pressure. Buffer[id. x + (id. y - 1) * grid. Resolution]; v. T = input. Pressure. Buffer[id. x + (id. y + 1) * grid. Resolution]; divergence = divergence. Buffer[id. x + id. y * grid. Resolution]; output. Pressure. Buffer[id. x + id. y * grid. Resolution] = (v. L + v. R + v. B + v. T + alpha * divergence) * beta; }else{ //határfeltételek if(id. x == 0) output. Pressure. Buffer[id. x + id. y * grid. Resolution] = input. Pressure. Buffer[id. x + 1 + id. y * grid. Resolution]; . . . } }

Projection __kernel void projection(const int grid. Resolution, __global float 2* input. Velocity. Buffer, __global float* pressure. Buffer, __global float 2* output. Velocity. Buffer) { int 2 id = (int 2)(get_global_id(0), get_global_id(1)); if(id. x > 0 && id. x < grid. Resolution - 1 && id. y > 0 && id. y < grid. Resolution - 1){ float p. L = pressure. Buffer[id. x - 1 + id. y * grid. Resolution]; float p. R = pressure. Buffer[id. x + 1 + id. y * grid. Resolution]; float p. B = pressure. Buffer[id. x + (id. y - 1) * grid. Resolution]; float p. T = pressure. Buffer[id. x + (id. y + 1) * grid. Resolution]; float 2 velocity = input. Velocity. Buffer[id. x + id. y * grid. Resolution]; output. Velocity. Buffer[id. x + id. y * grid. Resolution] = velocity – (float 2)(p. R - p. L, p. T - p. B); } else {//határfeltételek if(id. x == 0) output. Velocity. Buffer[id. x + id. y * grid. Resolution] = -input. Velocity. Buffer[id. x + 1 + id. y * grid. Resolution]; . . . } }

Diffusion __kernel void diffusion(const int grid. Resolution, __global float 2* input. Velocity. Buffer, __global float 2* output. Velocity. Buffer) { int 2 id = (int 2)(get_global_id(0), get_global_id(1)); float viscousity = 0. 01 f; float alpha = 1. 0 f / (viscousity * dt); float beta = 1. 0 f / (4. 0 f + alpha); if(id. x > 0 && id. x < grid. Resolution - 1 && id. y > 0 && id. y < grid. Resolution - 1){ float 2 float 2 v. L = input. Velocity. Buffer[id. x - 1 + id. y * v. R = input. Velocity. Buffer[id. x + 1 + id. y * v. B = input. Velocity. Buffer[id. x + (id. y - 1) v. T = input. Velocity. Buffer[id. x + (id. y + 1) velocity = input. Velocity. Buffer[id. x + id. y grid. Resolution]; * grid. Resolution]; output. Velocity. Buffer[id. x + id. y * grid. Resolution] = (v. L + v. R + v. B + v. T + alpha * velocity) * beta; } else { output. Velocity. Buffer[id. x + id. y * grid. Resolution] = input. Velocity. Buffer[id. x + id. y * grid. Resolution]; }}

Vorticity __kernel void vorticity(const int grid. Resolution, __global float 2* velocity. Buffer, __global float* vorticity. Buffer) { int 2 id = (int 2)(get_global_id(0), get_global_id(1)); if(id. x > 0 && id. x < grid. Resolution - 1 && id. y > 0 && id. y < grid. Resolution – 1){ float 2 v. L v. R v. B v. T = = velocity. Buffer[id. x + + + 1 + id. y * (id. y - 1) (id. y + 1) grid. Resolution]; * grid. Resolution]; vorticity. Buffer[id. x + id. y * grid. Resolution] = (v. R. y - v. L. y) - (v. T. x - v. B. x); } else{ vorticity. Buffer[id. x + id. y * grid. Resolution] = 0. 0 f; } }

Velocity from vorticity __kernel void add. Vorticity(const int grid. Resolution, __global float* vorticity. Buffer, __global float 2* velocity. Buffer) { int 2 id = (int 2)(get_global_id(0), get_global_id(1)); const float scale = 0. 2 f; if(id. x > 0 && id. x < grid. Resolution - 1 && id. y > 0 && id. y < grid. Resolution - 1){ float v. L = vorticity. Buffer[id. x - 1 + id. y * grid. Resolution]; float v. R = vorticity. Buffer[id. x + 1 + id. y * grid. Resolution]; float v. B = vorticity. Buffer[id. x + (id. y - 1) * grid. Resolution]; float v. T = vorticity. Buffer[id. x + (id. y + 1) * grid. Resolution]; float 4 grad. V = (float 4)(v. R - v. L, v. T - v. B, 0. 0 f); float 4 z = (float 4)(0. 0 f, 1. 0 f, 0. 0 f); if(dot(grad. V, grad. V)){ float 4 vorticity. Force = scale * cross(grad. V, z); velocity. Buffer[id. x + id. y * grid. Resolution] += vorticity. Force. xy * dt; } } }

External forces __kernel void add. Force(const float x, const float y, const float 2 force, const int grid. Resolution, __global float 2* velocity. Buffer, const float 4 density, __global float 4* density. Buffer) { int 2 id = (int 2)(get_global_id(0), get_global_id(1)); float dx = ((float)id. x / (float)grid. Resolution) - x; float dy = ((float)id. y / (float)grid. Resolution) - y; float radius = 0. 001 f; float c = exp( - (dx * dx + dy * dy) / radius ) * dt; velocity. Buffer[id. x + id. y * grid. Resolution] += c * force; density. Buffer[id. x + id. y * grid. Resolution] += c * density; }

Possible extensions �Buoyancy and gravity �Thermodynamic simulation � 3 D �Other grids: FCC �Interaction with solid bodies (voxelization, boundary conditions)