Introduction to CFD (computational fluid dynamics) using Python and Matlab

Computational Fluid Dynamics :

The term computational refers to one of the pathways to solve the fluid flow equations. the other possible pathways are experimental fluid dynamics and theoretical fluid dynamics. In experimental fluid dynamics, various experimental techniques are employed to determine the fluid flow phenomena. Theoretical fluid dynamics also refers to analytical fluid dynamics where we mathematically solve the governing equations to achieve fluid flow solution. In computational fluid dynamics, various numerical techniques are employed to simulate the flow of fluid using computers.

Fluid: Any material that cannot resist the flow of shear stress is called fluid. All Liquids and Gases are referred to as Fluids.

Dynamics: to determine both the kinematics parameters such as velocity, acceleration and the effect these parameters causes such as Force (Lift, Drag, Pressure forces, etc)

Computational fluid dynamics allows estimation of all flow-related phenomena using computers. This helps solve complex fluid flow phenomena and results can be obtained in no time with less cost than experimental technique.

Governing Equations :

2D incompressible flow:

Conservation of Mass :

`(delu)/(delx) + (delv)/(dely) = 0`

Conservation of Momentum :

`(delu)/(delt) + u(delu)/(delx) + v(delu)/(dely) = -1/rho(delp)/(delx)+vartheta((del^2u)/(delx^2)+(del^2u)/(dely^2))+F_x`

`(delv)/(delt) + u(delv)/(delx) + v(delv)/(dely) = -1/rho(delp)/(dely)+vartheta((del^2v)/(delx^2)+(del^2v)/(dely^2))+F_y`

In CFD, these differential equations are converted into Algebraic equations so that they can easily be manipulated. Finite differencing methods are one of the ways to convert them.

Taylor series:

Taylor series helps in determining the value of a function f(x) at multiple grid points.

if the value of a function is known at any given point `x_0`, then by using taylor series, the value of the function at a distance `h` from `x_0` can be approximated.

`f(x_0+h) = f(x_0) to`Zeroth approximation

`f(x_0+h) = f(x_0)+hfprime(x_0) to` First order approximation (Linear)

`f(x_0+h) = f(x_0)+hfprime(x_0)+h^2/(2!)fprimeprime(x_0) to` Second order approximation (curvature)

`f(x_0+h) = f(x_0)+hfprime(x_0)+h^2/(2!)fprimeprime(x_0)+h^3/(3!)fprimeprimeprime(x_0)+...h^n/(n!)f^n(x_0)+... to` General Formulation - (1)

`f(x_0-h) = f(x_0)-hfprime(x_0)+h^2/(2!)fprimeprime(x_0)-h^3/(3!)fprimeprimeprime(x_0)+...(-1)^nh^n/(n!)f^n(x_0) to`General Formulation

Forward differencing :

After rearranging the terms

`fprime(x_0) = (f(x_0+h)- f(x_0))/h -(h/(2!)fprimeprime(x_0)-h^2/(3!)fprimeprimeprime(x_0)+...h^(n-1)/(n!)f^n(x_0)+...) `

in most of the cases h is always less than 1

therefore `h>h^2>h^3>h^4....`

`fprime(x_0) = (f(x_0+h)- f(x_0))/h` - O(h)

`(f(x_0+h)- f(x_0))/h to` Approximation

`O(h) to` Truncation error

Backward differencing :

After rearranging the terms

`fprime(x_0) = (f(x_0)- f(x_0-h))/h +(h/(2!)fprimeprime(x_0)-h^2/(3!)fprimeprimeprime(x_0)+...(-1)^nh^(n-1)/(n!)f^n(x_0)+...) `

`(f(x_0)- f(x_0-h))/h to`Approximation

`O(h) to` Truncation error

Central Differencing :

`fprime(x_0) = (f(x_0+h)- f(x_0-h))/(2h) -(h^2/(3!)fprimeprimeprime(x_0)+...h^(n-1)/(n!)f^n(x_0)+...)`

`(f(x_0+h)- f(x_0-h))/(2h)``to` Approximation

`h^2/(3!)fprimeprimeprime(x_0)+...h^(n-1)/(n!)f^n(x_0)+...``to` O(h^2) `to` Truncation error

* Central differencing method gives better approximation results compared to Forward and Backward differencing methods

Example :

`f(x) = -4x^3 + 7x^2-3x+9`

`fprime(x) = -12x^2+14x-3`

`fprime(0) = -3`

The same can be approximated using finite difference methods. The following program has been executed in python

# Getting the numpy module

import math
import matplotlib.pyplot as plt 
import numpy as np


# polynomial definition

poly_p = np.array([-4, 7, -3, 9])

# Polynomial derivative (numpy)

p_der = np.polyder(poly_p)

# confirmation of correct derivative polynomial

print('Derivative polynomial =' , p_der)

# Analytical result at x = 0

p_der_eval = np.polyval(p_der,0.0)

print('Theoritical Derivative =' , p_der_eval)

#Numerical calculations using FDM (forward)

x_0 = 0.0
h = np.float64(0.25)
forward_difference = (np.polyval(poly_p, (x_0 + h)) - np.polyval(poly_p, x_0))/h
print('Forward derivative =', forward_difference)

#Numerical calculations using FDM (backward)

x_0 = 0.0
h = np.float64(0.25)
backward_difference = (np.polyval(poly_p, (x_0)) - np.polyval(poly_p, x_0 - h))/h
print('Backward derivative =', backward_difference)

#Numerical calculations using FDM (central)

x_0 = 0.0
h = np.float64(0.25)
central_difference = (np.polyval(poly_p, (x_0 + h)) - np.polyval(poly_p, x_0 - h))/(2*h)
print('Central derivative =', central_difference)

Results

Derivative polynomial = [-12  14  -3]
Theoritical Derivative = -3.0
Forward derivative = -1.5
Backward derivative = -5.0
Central derivative = -3.25

Second-order derivative :

From the above forward difference and backward difference equations, the second derivative can be determined as follows

Central differencing :

`fprimeprime(x_0) = (f(x_0+h)-2f(x_0)+f(x_0-h))/h-((2h^2)/(4!)f^(iv)(x_o)+....)`

`fprimeprime(x_0) = (f(x_0+h)-2f(x_0)+f(x_0-h))/(h^2)-O(h^2)`

Forward differencing :

`2f(x_0+h) = 2f(x_0)+2hfprime(x_0)+(2h^2)/(2!)fprimeprime(x_0)+(2h^3)/(3!)fprimeprimeprime(x_0)+...(2h^n)/(n!)f^n(x_0)+...`

`f(x_0+2h) = f(x_0)+2hfprime(x_0)+(4h^2)/(2!)fprimeprime(x_0)+(8h^3)/(3!)fprimeprimeprime(x_0)+...((2h)^n)/(n!)f^n(x_0)+...`

`fprimeprime(x_0) = (f(x_0+2h)-2f(x_0+h)+f(x_0))/(h^2)-O(h^2)`

similarly for

backward differencing :

`fprimeprime(x_0) = (f(x_0)-2f(x_0-h)+f(x_0 - 2h))/(h^2)-O(h^2)`

Introduction to Finite Volume method (FVM)

In CFD, the physical domain is discretized into a computational mesh to solve the algebraic (converted partial differential) equations. In these series of points some are defined at the boundary and the other at the interior points. The function is evaluated at the node (grid) points

The building block of FVM

In a finite volume, the grid point is the centroid of the volume.

There are two methods in which the domain can be discretized in the Finite Volume method

1D mesh FVM

In the first method, the domain is put over a finite-difference mesh. A finite volume is substituted around each grid point. Even the boundary points here have a finite volume. But the volume will be halved.

The whole domain will have an equal number of finite volumes or equal size of finite volumes in the second method. This will reduce the total number of grid points. Also, there will be no grid points at the boundary.

If boundary conditions and boundary values are defined, then the first method is employed. If rather than boundary value boundary flux is defined then the second method is suitable.

Methodology of FVM

The governing equation is integrated over the finite volume in FVM. The volume is generated between adjacent grid points. The faces where this finite volume is generated are east face(e) and west face(w). The fundamental principle in using FVM is to take the governing equations and simply integrate them over the finite volume.

`int(d^2T)/dx^2 dV = 0`

`int_e^w (d^2T)/dx^2dx = 0`

`|quad(dT)/dx|_e - |quad(dT)/dx|_w = 0 `

Depending on the variation of function between two grid points, the derivatives are numerically approximated.

Basic Algorithm used in Matlab and Python to solve equations in FVM

1D heat conduction :

The physics of the simulation is governed by the differential equation of pure diffusion

`(d^2T)/(dx^2) = 0`

Initially,

T(t =0 , x < 1) = 0 & T(t = 0, x = 1) = 1 ; (Initial conditions)

at any other instant T (t=t , x =0) =0 and T(t = t, x = 1) = 1 (Boundary conditions)

Objective :

To calculate the temperature distribution T(x), at equilibrium using the finite difference method.

The general idea of tackling a CFD problem

• Convert physical geometry into a computational mesh
• Discritize the governing equation on this mesh using a numerical scheme

`(d^2T)/(dx^2) = 0`

`darr` central difference scheme

`Tprimeprime(x_o) = (T(x_0+h) - 2T(x_0) + T(x_0-h))/h^2 = 0`

• Obtain an iterative formula from the discretized equation

`T(x_0) = (1/2)*(T(x_0+h) + T(x_0 - h))`

• Implement the iterations on every mesh point until converge

Use the above `T(x_0)` at the second, third, fourth,.... (n-1) point. The edges are fixed in this case.

Converge is achieved until the error reaches the required error `epsilon`

`sum_(x=x_0) abs(T_n - T_o) < epsilon`