CFD: 2D thermal solver using Octave
Introduction - Phase 1:
CFD which stands for Computational fluid dynamics, is a field which emerged in late 20th century with advent of computers and fast computing, processing capabilities. Intially Fluid studies were predominantly experimental based, then theories were constructed to explain the observations from experiments. Sometimes, theories were hypothised and was experimentally proved by experiments. However, in late 20th century, with the advancement in fluid dynamics field, it came with the challenge of setting up experimental apparatus for the test condition such as engine running at very low barometric pressure or wind tunnels that could accomodate high speed and high temperature simulataneously. It also proved challenging in determining the underlying nuances in physics of phenomenon.
Hence, the third field, which is an amalgmation of both theory and experiments took birth. CFD, provides every analyst with personal virtual lab where he could run numerous experiments to understand the behaviour, predict the outcomes.
Unlike, any field of science, fluids is goverened by three conservation laws, namely
1) Conservation of Mass
2) Conservation of momentum
3) Conservation of energy.
All three together form the Navier stokes equation and there is no closed form solution existing till date. This challenge was partially resolved by utilising numerical schemes and techniques of linear algebra. The geometry of the object under test was broken down into numerous small sized elements by the process of meshing. The property values were calculated in each cell, an array of numerical equations are generated either by means of Finite element method or Finite difference or Finite volume method.
These array of numerical equations are solved computationally and results are post processed to analyse various phenomenon. All this, happens virtually without the need to setup lab equipment and results from simulation could then be experimentally validated
Thermal Solver
Conservation of energy equation gives raise to determining temperature of the fluid element. In this work, a 2D thermal solver is built to solve a simple laplacian equation of temperature. The user is provided with the option of entering various boundary condition
1) Dirchlit's boundary condition: Constant Temperature is applied at the boundary
2) Neuman Boundary condition: Temperature gradient is applied at the boundary
3) Mixed boundary condition: Combination of Dirchlit and Neuman boundary condition
The compuatation time can be awfully high and hence various means are employed to reduce the compuatation time. In this work, numerical iteration is used to reduce the computation time.
Two numerical iterations are provided as an option to user, Jacobian iteration or Gauss siedel iteration. The option is also provided to add additional numerical iteration techniques such as successive over relaxation method.
The laplacian equation that is considered is as follows:
Finite difference method is employed for generating numerical equations.
The temperature is calculated at every node based on the boudary condition and geometry, user has chosen and contour and surface plots are generated for the visualisation. Convergence plot is also generated for better understanding of improving the comuputation time.
These results are of vital importance in engineering as they help us understand the thermal behaviour within fluid element, thermal stresses and better means of managing thermal boundaries. Heat flow, insulation and allied comcepts could be calculated from the work. Overall, it determines the energy the system would carry due to changes in temperature and boundary conditions.
Objectives:
1) To caclculate the temperature at each node for the given user geometry, boundary condition and iteration methodology
2) To plot the temperature profile for better means of understanding the temperature distribution in the element
3) To improve computation time by employing various numerical iterations
4) To study the advantage of Gauss siedel over Jacobian iterations.
Motivation:
To better understand thermal behaviour in fluid element for various boundary condition, importance of insulation and by the means, designing of better thermal management systems
Geometry:
User is allowed to enter the length and breadth of the fluid element under consideration, number of discretisation elements are also taken in as the user input
fprintf("\n............................Enter Geometric details for the plate.....................\n")
len = input("\n Enter length of the fluid element in mm: ");
bre = input("Enter breadth of the fluid element in mm: ");
cellsx = input("Enter number of discretization elements required in x direction: ");This generates a 4 sided objects and meshing is done based on user requirement.
Boundary condition and Initial condition:
Boundary conditions are either taken in Dirchlit's format or Neuman format. Insulated boundary layer is an example of Neuman boudary condition which is used in the code. Dirchlit's boundary conditon uses imposed wall temperature.
All node elements are intialised to zero.
fprintf("\n.............................Enter thermal boundary condition..........................\n")
fprintf("This code offers flexibity to enter dirchlit or Neuman or combination of both for boundary condition\n")
fprintf("\nEnter 1 if the side is insulated\n")
in1 = input("Whether left vertical wall is insulated: ");
in2 = input("Whether bottom wall is insulated: ");
in3 = input("Whether top wall is insulated: ");
in4 = input("Whether right vertical wall is insulated: ");
if(in1==1)
tem1 = 100;
else
tem1 = input("Enter temperature of left vertical wall: ");
endif
if(in2==1)
tem2 = 100;
else
tem2 = input("Enter temperature of bottom wall: ");
endif
if(in3==1)
tem3=100;
else
tem3 = input("Enter temperature of top wall: ");
endif
if(in4==1)
tem4 = 0;
else
tem4 = input("Enter temperature of right vertical wall: ");
endif
..
.
.
.
.
.
.
%Setting up thermal boundaries of the plate
temp = zeros(s,t);
%For left and right vertical walls
for i=1:s
if(in1==1)
temp(i,1) = temp(i,2);
else
temp(i,1) = T(1);
endif
if(in4==1)
temp(i,t) = temp(i,t-1);
else
temp(i,t) = T(4);
endif
endforOctave Code:
1) Geometric details, thermal conditions and Numerical iteration preference is read from user input
2) Meshing is done based on inputs
3) All nodes of the fluid element is intialised, wall temperatures are set using the user inputs
4) The temperature is calculated at each node by discretising the laplacian equation shown before
5) Process is iterated till the error is under accepatable zone
6) Switch-Case is used to switch between various iteration technique
7) Plots are constructed based on the outputs.
%This is coded by Sweekar, please reach out at [email protected]
%This code is the intellectual property of Sweekar and please contact before reproduction
%CFD which stands for computational fluid dynamics is a field which emerged in late 1950s with advent of computers
%This has made possible in understanding physics of fluids and heat transfer at a much higher resolution than we could have from experiments
%It has also made possible couple of experimets which could be done in the lab, like simulating engine at low barometric pressure
%High speed wind tunnel with high temperature boundary condition, both imposed simulataneously and safely controlled
%CFD gave a virtual lab for every coder to experiment on computer without risk of damaging lab apparatus
%This is the code to solve energy equation in the fluid element without heat generation and heat source
%This code provides option to choose either dirchlit's boundary condition ot Neuman's or combination of both
%From the solver approach, this follows explicit solution methodology
%This code also provides option of choosing numerical scheme of iteration
%The code provides option to choose between Jacobian iteration and Gauss siedel iteration
%The code also allows easy addition of other numerical iteration schemes such as Successive over relaxation.
%Cleanup
close all
clear all
clc
fprintf("....................................2D CFD Thermal Solver...............................\n")
fprintf(".........................This is based on Finite difference method......................\n")
fprintf(".................Laplacian equation is solved, no heat source is considered.............\n")
fprintf(".........................This code follows Explicit solution methodogy...................n")
fprintf(".......This is an elliptical system, a disturbance is propagated in all directions.......\n")
fprintf("......The code follows finite difference princicples for solution with simple meshing.....\n")
%User inputs
%Geometric inputs
fprintf("\n............................Enter Geometric details for the plate.....................\n")
len = input("\n Enter length of the fluid element in mm: ");
bre = input("Enter breadth of the fluid element in mm: ");
cellsx = input("Enter number of discretization elements required in x direction: ");
%Thermal inputs
fprintf("\n.............................Enter thermal boundary condition..........................\n")
fprintf("This code offers flexibity to enter dirchlit or Neuman or combination of both for boundary condition\n")
fprintf("\nEnter 1 if the side is insulated\n")
in1 = input("Whether left vertical wall is insulated: ");
in2 = input("Whether bottom wall is insulated: ");
in3 = input("Whether top wall is insulated: ");
in4 = input("Whether right vertical wall is insulated: ");
if(in1==1)
tem1 = 100;
else
tem1 = input("Enter temperature of left vertical wall: ");
endif
if(in2==1)
tem2 = 100;
else
tem2 = input("Enter temperature of bottom wall: ");
endif
if(in3==1)
tem3=100;
else
tem3 = input("Enter temperature of top wall: ");
endif
if(in4==1)
tem4 = 0;
else
tem4 = input("Enter temperature of right vertical wall: ");
endif
%Solver approach and Iteration method
fprintf("\n...............................Enter solver method................................\n")
fprintf("This code provides option to solve via Jacobian iteration or Gauss siedel iteration\n")
fprintf("\nPlease contact Sweekar if additional Iteration schemes are to be added\n")
fprintf("Enter 1 for Jacobian iteration\n")
fprintf("Enter 2 for gauss siedel iteration\n")
operator = input("Enter the integer for the iteration method: ");
%Meshing
cellsx = cellsx+1;
x = linspace(0,len,cellsx);
y = linspace(0,bre,cellsx);
s = length(x);
t = length(y);
T = [tem1 tem2 tem3 tem4];
%Setting up thermal boundaries of the plate
temp = zeros(s,t);
%For left and right vertical walls
for i=1:s
if(in1==1)
temp(i,1) = temp(i,2);
else
temp(i,1) = T(1);
endif
if(in4==1)
temp(i,t) = temp(i,t-1);
else
temp(i,t) = T(4);
endif
endfor
%For Bottom and top walls
for j=1:t
if(in2==1)
temp(1,j)=temp(2,j);
else
temp(1,j) = T(2);
endif
if(in3==1)
temp(s,j)= temp(s-1,j);
else
temp(s,j) = T(3);
endif
endfor
%Initialising the corners
if(T(1)>=T(2))
temp(1,1) = T(1);
else
temp(1,1) = T(2);
endif
if(T(2)>=T(4))
temp(1,t) = T(2);
else
temp(1,t) = T(4);
endif
if(T(4)>=T(3))
temp(s,t) = T(4);
else
temp(s,t) = T(3);
endif
if(T(1)>=T(3))
temp(s,1) = T(1);
else
temp(s,1) = T(3);
endif
Te = temp;
e = 2;
%Counter
countj = 0;
countg = 0;
counto = 0;
p=1;
switch(operator)
case 1
%Jacobian iteration
while abs(e)>power(10,-14)
sum = 0;
for i=1:s
if(in1==1)
Te(i,1) = Te(i,2);
endif
endfor
for j=1:t
if(in2==1)
Te(1,j) = Te(2,j);
endif
endfor
for j=1:t
if(in3==1)
Te(s,j) = Te(s-1,j);
endif
endfor
for i=1:s
if(in4==1)
Te(i,t) = Te(i,t-1);
endif
endfor
for i=2:s-1
for j=2:t-1
Te(i,j) = 0.25*(temp(i+1,j)+temp(i-1,j)+temp(i,j+1)+temp(i,j-1)); %Average of temperatures in all surrounding nodes
endfor
endfor
d = Te-temp;
for i=2:s-1
for j=2:t-1
sum = sum+d(i,j);
endfor
endfor
e = sum;%Calculation of error
err(p)=e;
p=p+1;
countj = countj+1;
ite(p) = countj;
fprintf("%f error after %d iteration\n",e,countj)
temp = Te;
endwhile
case 2
%Gauss Siedel iteration
while abs(e)>power(10,-14)
sum = 0;
for i=1:s
if(in1==1)
Te(i,1) = Te(i,2);
endif
endfor
for j=1:t
if(in2==1)
Te(1,j) = Te(2,j);
endif
endfor
for j=1:t
if(in3==1)
Te(s,j) = Te(s-1,j);
endif
endfor
for i=1:s
if(in4==1)
Te(i,t) = Te(i,t-1);
endif
endfor
for i=2:s-1
for j=2:t-1
Te(i,j) = 0.25*(temp(i+1,j)+Te(i-1,j)+temp(i,j+1)+Te(i,j-1));
endfor
endfor
d = Te-temp;
for i=2:s-1
for j=2:t-1
sum = sum+d(i,j);
endfor
endfor
e = sum;
err(p)=e;
p=p+1;
countg = countg+1;
ite(p) = countg;
fprintf("%f error after %d iteration\n",e,countg)
temp = Te;
endwhile
end
%Plotting
mn=length(ite)-1;
figure(1)
contourf(x,y,temp) %Contour plot of temperature
colorbar
xlabel('Length')
ylabel('Breadth')
title('Temperature contour')
figure(2)
plot(ite(1:mn),err) %Convergence plot
xlabel('Iterations')
ylabel('Error')
title('Convergence plot')
figure(3)
surf(x,y,temp) %Surface plot
colorbar
xlabel('Length')
ylabel('Breadth')
title('Temperature surface')Results and Validation
The UI is as follows
Case 1) Jacobian Iteration
Len = 20; Bre = 20; Number of cells = 400;
Temperature of the left vertical wall = 50C;
Temperature of the right vertical wall = Insulated
Temperature of the bottom wall = 30C
Temperature of the top wall = 25C
Number of Iteration to solve using the Jacobian iteration is 4393
Case 2: When same inputs is repeated for the Gauss siedel iteration
The Gauss siedel takes 2320 iteration which is almost half the time taken to solve Jacobian iteration.
Validation:
Check1: Checking for the boundary conditions,
Temperture along the walls should be same as what was entered that is 50, 30, 25C
Check 2: Along the insulated boundary the dT/dx = 0, hence temperature in i-1th column should match
Check 3: None of the temperature in any node can be greater than 50C and lesser than 25C.
Check4: Temperature transition should be smooth
All these could be easily verified with the contour plot
Check 5: Similarly if all surfaces are insulated, temperature would be that of initialised temperature, here temperatures were intialised to 0, It shouldnt even take one iteration to throw result
Interesting plot: When all three surfaces are insualted and left wall is set at 100C, temperature changes in folds.
Discussion on results and report:
1) Gauss siedel iteration is much faster than Jacobian iteration
2) CFD results should be validated as JUNK IN is equal to JUNK OUT in CFD domain
3) CFD is a thrid approach, amalgm of experiment and theory
4) Laplacian equation could be replaced by the poission's equation, Wave equation with no additional changes to the code
Resources:
1) Introduction to CFD by Versateeg and Malalasekara
2) Computational fluid dynamics: The basics with applications by Anderson Junior
3) Fluid mechanics by Kundu Cohen