Showing posts with label MATLAB. Show all posts
Showing posts with label MATLAB. Show all posts

Wednesday, 3 May 2017

Lid Driven Cavity - Post 3

Code Files$^1$

1. Main.m (The file to be run)
2. streamfunc.m (Solves the stream vorticity function)
3. velocity.m (Solves the velocity function)
4. omega.m (Iterates the vorticity transport equation)
5. BC.m (Implements the Boundary conditions)
6. Upwind.m (Implements the upwind scheme)
7. rmse psi.m (Computes Errors)

Terminologies used in code$^2$

1. w = $\omega$(Vorticity)
2. psi = $\psi$(StreamFunction)
3. u = (U −V elocity)
4. v = (V −V elocity)
5. u0 = $u_0$(Lid Velocity)
6. psi_1 = $c_1$
7. gamma = $\gamma$ (Kinematic viscosity)
8. dx = ∆x
9. dy = ∆y
10. alpha = Relaxation parameter for stream vorticity function
11. alpha1 = Relaxation parameter for vorticity transport equation
12. x = (X −length) of cavity
13. y = (Y −length) of cavity
14. Nx = $N_x$(No of grid points along x-direction
15. Ny = $N_y$(No of grid points along y-direction
16. err - Percentage error in $\psi$
17. err1 - Percentage error in $\omega$
18. ERR1 - Error in stream vorticity equation
19. ERR2 - Error in vorticity transport equation

-------------------------------------------------------------------------------------------------------------------------
$^1$ Link to code given under section Matlab Code URL
$^2$ Variable Names used in the Matlab code implementing the Stream Vorticity approach for the Modified Lid Driven Cavity Problem

Sunday, 19 February 2017

CFD Tutorial - 2 (Laasonen scheme)

                           Computational Fluid Dynamics 

Problem statement:
Given a rod of Length  L ,with boundary conditions,initial conditions as follows:

Boundary condition:
\begin{eqnarray}
T(0,t) = 0^{\circ}C \\
T(L,t)= 1^{\circ}C
\end{eqnarray}

Initial condition:

\begin{eqnarray}
T(x,0) = 0^{\circ}C
\end{eqnarray}

Compute the temperature for t = 0s to 20s for various values of $\Delta T = 0.1s,0.01s,0.001s$ .

Governing Equations:
PDE:
\begin{eqnarray}
\frac{\delta T}{\delta t} &=& \alpha \frac{\delta ^2 T}{\delta x^2}
\end{eqnarray}

Finite difference formulation using BTCS scheme (Implicit)
\begin{equation}
-\gamma T^{n+1}_{i-1}+T^{n+1}_{i}(2 \gamma +1) -\gamma T^{n+1}_{i+1} = T^n_{i}
\end{equation}

Pseudo Code


  • Initialize the variables $\alpha ,\Delta t,T,\Delta x, N_x , L$ .
     (Note here T is a matrix with $N_x$ columns ,and $20/(\Delta t) +1= N_y$ rows)
  •  For n = 2 to $N_y$ .
     Solve the equation below for $T^{n+1}$ using the TDMA algorithm.
     $$AT^{n+1} = T^{n}$$
    Where,
    $A =\begin{bmatrix}
    1 & 0 &0& \cdots &\cdots& 0 \\
    -\gamma &2\gamma + 1 &-\gamma &0 &\cdots &0&\\
    0 & -\gamma &2\gamma + 1 &-\gamma & \ddots &0\\
    0 & 0&\ddots &\ddots &\ddots & \ddots \\
    0 & \cdots&\cdots &\cdots &0 & 1 \\
    \end{bmatrix}$

The TDMA consists of converting the Tridiagonal matrix to an upper triangular matrix , then solving the system of equations by back substitution.

TDMA Algorithm:


  •  For i = 2:$(N_x-1)$
     $A[i,:] = A[i,:] - A[i-1,:]\frac{A[i,i-1]}{A[i,i]}$
     $T[n-1,i] = T[n-1,i] - T[n-1,i-1]\frac{A[i,i-1]}{A[i,i]}$
  •  Back substitution.
  •  End

Results :

Below are sample Graphs showing the Unconditional stability of the Implicit schemes.




Matlab Code:

https://github.com/RAAKASH/Intro-to-CFD-/tree/master/Assignment%202.

Note : 
  • Assignment 1.m is the file of the previous assignment/tutorial in the form of a function.
  • Assignment 2.m is the file of the present assignment/tutorial in the form of a function.
  • BTCS.m  is the same Assignment 2.m but has been optimized in time to get faster results.

Friday, 17 February 2017

FTCS Tutorial -1

Computational Fluid Dynamics


Problem Statement

 Given a rod of Length  L ,with boundary conditions,boundary conditions as follows :  
   $T(0,t) = 0^{\circ}C $
   $T(1,t) = 0^{\circ}C $
Also the initial conditions are as follows:
  $T(x,0) = 0^{\circ}C $
Find the Temperature distribution for the complete length of the rod for a time period of 300s.

Governing Equations:

The physical equation is as follows from the below heat conduction equation.
$ \frac{\delta T}{\delta t} = \alpha \frac{\delta ^2 T}{\delta x^2}$ 

Now we are applying FTCS scheme to the above equation we get: 

$T_{i}^{n+1} = T_{i}^{n} + \alpha \Delta t \frac{(T_{i+1}^n - 2T_{i}^n + T_{i-1}^n)}{(\Delta x)^2}$

Pseudo code :


  1.  Initialize the variables $\alpha ,\Delta t,T,\Delta x,N_x,L$.
    (Note here T is a matrix with $N_x$ columns ,and $20/(\Delta t) +1= N_y$ rows)
  2.  For n = 2 to $N_y$ execute the statements 3 and 4.
  3.  For i = 2 to $N_x-1$ execute statement 4.
  4.  $T[i][n+1] = T[i][n] + \frac{\alpha \Delta t}{(\Delta x)^2} (T[i+1][n] - 2T[i][n] + T[i-1][n] )$
  5.  end


Stability of scheme : 

Let $\gamma = \alpha \frac{\Delta t}{ \Delta x^2}$.
Scheme is stable if $\gamma$ < $\frac{1}{2} $
Derivation : 
$$e(n,x) = A_n  e^{ikx}$$
Here $k$ is the wave number , $i = (-1)^{0.5}$
Therefore substituting this in the finite difference scheme  we get
$$ \frac{A_{n+1}}{A_{n}} = 1 + \frac{\alpha \Delta t}{\Delta x^2}\left ( e^{i k \Delta x} -2 + e^{-i k \Delta x}   \right )$$
$$ \frac{A_{n+1}}{A_{n}} = 1 + \frac{\alpha \Delta t}{\Delta x^2}\left (-4 sin^2 (k \frac{\Delta x}{2} )\right )$$


$$ \left \vert \frac{A_{n+1}}{A_{n}} \right \vert   \leq 1$$
$$  -1 \leq 1 + \frac{\alpha \Delta t}{\Delta x^2} \left (-4 sin ^2  (k \frac{\Delta x}{2} ) \right ) $$
Finally we get the required criterion. 

MATLAB Run : 





Web App Run:




MATLAB /Web App: