Indeed, the 13th step now exists! This case is called the case of the backward facing step (BFS)! ⬜ This is an unofficial continuation to this. If I get it approved by Dr. Barba, then it will be official. The original series is in Python but I coded this in MATLAB without using many MATLAB specific functions so the code can be translated to other programing languages 🖧 quite easily.
I write about Propulsion, Aerodynamics and Renewable Energy (Wind/Hydro Turbines).
Sunday 31 March 2024
13th Step of the 12 steps to Navier-Stokes 😑
Saturday 26 August 2023
Open - Source coded 3D Navier–Stokes equations in C++
Here are 3D Navier–Stokes equations configured for lid-driven cavity flow. The syntax is C++. The results from this code are shown in Fig. 1. From the pressure-Poisson equation, I removed mixed derivative terms to improve solution stability. 😆
p[0][j][k] = p[1][j][k];
p[num_i - 1][j][k] = p[num_i - 2][j][k];
p[i][0][k] = p[i][1][k];
p[i][num_j - 1][k] = p[i][num_j - 2][k];
p[i][j][0] = p[i][j][1];
p[i][j][num_k - 1] = 0.0;
u[i][j][k] = u[i][j][k] - u[i][j][k] * dt / (2 * h) * (u[i+1][j][k] - u[i-1][j][k]) - v[i][j][k] * dt / (2 * h) * (u[i][j+1][k] - u[i][j-1][k]) - w[i][j][k] * dt / (2 * h) * (u[i][j][k+1] - u[i][j][k-1]) - dt / (2 * rho * h) * (p[i+1][j][k] - p[i-1][j][k]) + nu * (dt / (h * h) * (u[i+1][j][k] - 2 * u[i][j][k] + u[i-1][j][k]) + dt / (h * h) * (u[i][j+1] [k] - 2 * u[i][j][k] + u[i][j-1][k]) + dt / (h * h) * (u[i][j][k+1] - 2 * u[i][j][k] + u[i][j][k-1]));
u[0][j][k] = 0.0;
u[num_i - 1][j][k] = 0.0;
u[i][0][k] = 0.0;
u[i][num_j - 1][k] = 0.0;
u[i][j][0] = 0.0;
u[i][j][num_k - 1] = -Uinf;
v[i][j][k] = v[i][j][k] - u[i][j][k] * dt / (2 * h) * (v[i+1][j][k] - v[i-1][j][k]) - v[i][j][k] * dt / (2 * h) * (v[i][j+1][k] - v[i][j-1][k]) - w[i][j][k] * dt / (2 * h) * (v[i][j][k+1] - v[i][j][k-1]) - dt / (2 * rho * h) * (p[i][j+1][k] - p[i][j-1][k]) + nu * (dt / (h * h) * (v[i+1][j][k] - 2 * v[i][j][k] + v[i-1][j][k]) + dt / (h * h) * (v[i][j+1][k] - 2 * v[i][j][k] + v[i][j-1][k]) + dt / (h * h) * (v[i][j][k+1] - 2 * v[i][j][k] + v[i][j][k-1]));
v[0][j][k] = 0.0;
v[num_i - 1][j][k] = 0.0;
v[i][0][k] = 0.0;
v[i][num_j - 1][k] = 0.0;
v[i][j][0] = 0.0;
v[i][j][num_k - 1] = 0.0;
w[i][j][k] = w[i][j][k] - u[i][j][k] * dt / (2 * h) * (w[i+1][j][k] - w[i-1][j][k]) - v[i][j][k] * dt / (2 * h) * (w[i][j+1][k] - w[i][j-1][k]) - w[i][j][k] * dt / (2 * h) * (w[i][j][k+1] - w[i][j][k-1]) - dt / (2 * rho * h) * (p[i][j][k+1] - p[i][j][k-1]) + nu * (dt / (h * h) * (w[i+1][j][k] - 2 * w[i][j][k] + w[i-1][j][k]) + dt / (h * h) * (w[i][j+1][k] - 2 * w[i][j][k] + w[i][j-1][k]) + dt / (h * h) * (w[i][j][k+1] - 2 * w[i][j][k] + w[i][j][k-1]));
w[0][j][k] = 0.0;
w[num_i - 1][j][k] = 0.0;
w[i][0][k] = 0.0;
w[i][num_j - 1][k] = 0.0;
w[i][j][0] = 0.0;
w[i][j][num_k - 1] = 0.0;
Fig. 1, Velocity and pressure iso-surfaces |
Of course constants need to be defined, such as grid spacing in space and time, density, kinematic viscosity. These equations have been validated, as you might have read and here already! Happy coding!
If you want to hire me as your PhD student in the research projects related to turbo-machinery, aerodynamics, renewable energy, please reach out. Thank you very much for reading.
Monday 26 March 2018
1D Transient Diffusion (MATLAB code)
Friday 23 March 2018
1D non-Linear Convection (MATLAB code)
1D Linear Convection (MATLAB code) (Updated with CFL condition)
New Code, with the CFL condition satisfied.
Friday 9 March 2018
Solve PDEs using MATLAB (Code)
The example equation is mentioned below.
d2y/dx2 = x2+1 and y(0)=10, y(1)=20
The result of the code is a plot between the dependent and the independent variables, as shown in Fig. 1. The accuracy of the solutions lie around 1.5 % of the analytical solution.
%solution to 1D Poisson's equations with dirichlet boundary conditions, 2nd order accurate, uniform grid
clear; clc;
Nx=101; %Enter the number of grid points
clc;
A=zeros(Nx,Nx); %Stiffness matrix
B=zeros(Nx,1); %Dependent variable vector, also used for boundary conditions and other known terms
x=zeros(Nx,1); %Independent variable vector
x(Nx)=1; %Domain size
dx=x(Nx)/(Nx-1); %Grid spacing
B(1)=(((x(1).^2)+1)*(dx.^2))-10; %Enter the function the PDE equals to multiplied by the square of the grid spacing minus the first boundary condition
B(Nx)=(((x(Nx).^2)+1)*(dx.^2))-20; %Enter the function the PDE equals to multiplied by the square of the grid spacing minus the second boundary condition
for i=2:Nx-1 %Off boundary values for the dependent variable vector
B(i)=((x(i).^2)+1)*(dx.^2); %Enter the function the PDE equals to multiplied by the the square of the grid spacing
end
for i=1:Nx-1 %Loop for independent variable
x(i+1)=(i)*dx;
end
for i=1:Nx %Diagonal terms for the stiffness matrix
A(i,i)=-2;
end
for i=1:Nx-1 %Off diagonal terms for the stiffness matrix
A(i,i+1)=1;
A(i+1,i)=1;
end
X=A\B;
hold on
set(gca,'FontSize',28)
xlabel('Independent Variable','FontSize',38.5)
ylabel('Dependent Variable','FontSize',38.5)
title('Dependent VS Independent Variable (d^2y/dx^2 = x^2+1)','FontSize',38.5);
plot(x,X)
Wednesday 21 February 2018
Solve ODEs using MATLAB (Code)
The equation in question is mentioned below.
dy/dx = x+4 and y(0)=3
A uniform grid was employed, as I have not yet figured out how to make non uniform grid. When I figure it out, I will update it here, which will be soon.
The code, mentioned below, starts with a variable called Nx, which represents the total number of grid points in which the domain is divided in to. A rule of thumb, the more points in the domain, more accurate the solution becomes at the cost of computation resources, mainly CPU time and RAM and storage.
x and y are the vectors in which the known and the unknown values are stored. x has the known values, like the length of a pipe through which a fluid flow or a rod through which heat flows etc. y stores the unknown values like temperature of the rod or velocity of the fluid etc. The values are stored at nodes.
y(1) means first value in the vector y. which also happens to be our boundary condition. y is the value of the function at x=0. x=0 is the first value in in the x vector, written as x(1).
x(Nx-1) represents the domain size. Minimum domain size is usually zero e.g. starting point of a rod or a pipe and maximum can be any real number.
To calculate the grid spacing, represented by dx, the domain size was divided by the number of elements. In this example, we have four nodes, two at corners and two in-between, forming a total of three elements.
First order finite difference method was used to discretize the given equation. And a for loop was used for the solution. Using 4 grid points, the solution has an accuracy of 2.22% against the analytical solution of the equation which can be verified by hand calculations. Thank you for reading.
%solution to dy/dx=x+4 y(0)=3
%uniform grid
clear; clc;
Nx=4; %number of grid points
x=zeros(Nx,1); %x vector, known, goes from some known value to maximum value of the domain
y=zeros(Nx,1); %y vector, one value known due to boundary condition rest unknown
y(1)=3; %boundary condition
x(Nx-1)=1; %domain size
dx=x(Nx-1)/(Nx-1); %grid spacing
for i=1:Nx-1
x(i+1)=(i)*dx;
y(i+1)=(dx.*(x(i)+4))+y(i); %discretized differential equation
end
Friday 23 October 2015
LU Decomposition MATLAB
%Please don't mess around with my code
clc;
R=input('How many variables are there in your system of linear equations? ');
clc;
disp ('Please enter the elements of the Coefficient matrix "A" and the Constant/Known matrix "B" starting from first equation');
A=zeros(R,R);
B=zeros(R,1);
for H=1:R
for i=1:R
fprintf ('Coefficient Matrix Column \b'), disp(i), fprintf ('\b\b Row \b'), disp(H);
I=input(' ');
A(H,i)=A(H,i)+I;
end
fprintf ('Known Matrix Row \b'), disp(H);
J=input(' ');
B(H,1)=B(H,1)+J;
end
clc;
[L,U] = lu(A);
Y = L\B;
X = U\Y;
fprintf ('The Coefficient Matrix "A" is \n'), disp(A);
fprintf ('The Constant/Known Matrix "B" is \n'), disp(B);
fprintf ('The Solution Matrix "X" is \n'), disp(X);