ASSIGNMENT ON COMPUTATIONAL FLUID DYNAMICS

SUBMITTED BY SUBMITTED TO
ANKIT KU. SINGH( 2022TH11) DR. ASHWINI KUMAR YADAV

VISHAKH P.NAIR(2022TH10)

G.V. SRIKANT(2022THO8)

BHANWAR SINGH(2022TH07)
 ASSIGNMENT 1
1. Solve the one dimension steady state conduction equation
∂
( )
k
∂T
∂x ∂ x
=0

Subjected to following boundary condition

X = 0, T = 500
X = L, T= 100
Using:
Finite volume method
Finite difference method
System of linear equation should be solved by
 TDMA
 Gauss-seidal iteration method
 Jacobin method
 SOR method
Compare the numerical solution with exact solution of above equation
for different number of node Compare numerical solution with results
of commercial code-FLUENT
 Solution
Assumptions:-

1. Steady state.
2. One dimensional heat flow.
3. Thermal conductivity of material is independent of temperature.
4. Material is homogeneous and isotropic.
5. No heat generation.
6. Negligible radiation exchange with surrounding.

Using Finite Volume Method

Tb Ta
1 2 n-1 n
B A
∂ x /2∂ x ∂ x /2

∂
( )
k
∂T
∂x ∂ x
=0

∫ ∂∂x (k ∂T
∂x )
dV =0
CV

∭ ∂∂x ( k ∂∂ Tx ) dxdydz=0
w

A=dydz
e

∫ ∂∂x ( k ∂∂ Tx ) A dx=0
w

[( kA
∂T
∂x e) (
− kA
∂T
∂x ) ]=0
w

First, we develop a formula valid for nodal points 2, 3,4,5 …n-1

byintroducing the usuallinear approximations for the temperature gradient.
Subsequent division by cross-sectional area A gives

[ (
ke A e
T E −T P
∂ xe )
T −T W
−k w Aw P
∂ xw
=0
( )]
By rearranging all terms
 ( ke A e kw A w
+
∂ xe ∂ xw
T P=
)
k e Ae
∂ xe
T E+
k w Aw
∂ xw W
T

Because Ae=Aw=A, ke= kw=k and ∂ x e =∂ x w =∂ x

( ∂1x + ∂1x ) T = ∂1x T + ∂1x T

P E W

For nodes 2, 3, 4, 5, 6…….n-1

a P T P=aE T E + aW T W + Su

aE aW aP Su

1 1 a P=aW + aE 0
∂x ∂x

Next, we apply the boundary conditions at nodal points 1 and n.

At node 1the west control volume boundary is kept at a specified
temperature.

( ) ( )
T E −T P T −T B
kA −kA P =0
∂x ∂x
2

By rearranging all terms

( ∂1x + ∂2x ) T = ∂1x T + ∂2x T

P E B

a P T P=aE T E + aW T W + Su

aE aW aP Su SP

1 0 a P=aW + aE −S P 2 −2
T
∂x ∂x B ∂x

Similarly for nodes n

 kA ( T ∂−Tx )−kA ( T ∂−Tx )=0
P W A P

By rearranging all terms

( ∂1x + ∂2x ) T = ∂1x T

P W +
2
T
∂x A

a P T P=aE T E + aW T W + Su

aE aW aP Su SP

0 1 a P=aW + aE −S P 2 −2
T
∂x ∂x A ∂x

We can write in a tabular form for different nodes

Nodes aE aW aP Su SP

1 1 0 a P=aW + aE −S P 2 −2
T
∂x ∂x B ∂x

2 1 1 a P=aW + aE −S P 0 0
∂x ∂x

3 1 1 a P=aW + aE −S P 0 0
∂x ∂x
1 1 a P=aW + aE −S P 0 0
∂x ∂x

n-1 1 1 a P=aW + aE −S P 0 0
∂x ∂x

n 0 1 a P=aW + aE −S P 2 −2
T
∂x ∂x A ∂x

Matrix form of equation

( )
a P −ae 0 0 0 0
−aw a P −ae ⋯ 0 0 0
0 −aw aP 0 0 0
⋮ ⋱ ⋮ =Tridiagonal Matrix
0 0 0 aP −ae 0
0 0 0 ⋯ −a w aP −ae
0 0 0 0 −a w aP

[] []
T1 Su
T2 Su
T3 Su
T= ⋮ =Unknown matrix, B= ⋮ =Known matrix
⋮ ⋮
T n−1 Su
Tn Su

( )[ ][ ]
a P −ae 0 0 0 0 T1 Su
−aw a P −ae ⋯ 0 0 0 T2 Su
0 −aw aP 0 0 0 T3 Su
⋮ ⋱ ⋮ * ⋮ = ⋮
0 0 0 aP −ae 0 ⋮ ⋮
0 0 0 ⋯ −a w aP −ae T n−1 Su
0 0 0 0 −a w aP Tn Su

Now above system of linear equation can be solved by any direct or iterative
method.

Using Finite Difference Method

T −1 T 0 T 1 T 2 T m+1

∆x∆x∆x

X=0 ∂
k ( )
∂T
∂x ∂ x
=0 X=L

T i−1−2 T i+T i+1

=0
( ∆ x )2

When we rearrange this equation

 2T i=T i+1 +T i −1
Boundary Conditions
At x=0, T=500
At x=L, T=100
We put i=1, 2, 3, 4… m
2T 1=T 2+T 0
2T 2=T 3+ T 1
2T 3=T 4 +T 2
2T m−2=T m−1+ T m−3
2T m−1=T m +T m −1
2T m=T m+1 +T m−1
[ A ] [ T ] =[ G ]

( )
−2 1 0 0 0 0
1 −2 0 ⋯ 0 0 0
0 1 −2 0 0 0
A= ⋮ ⋱ ⋮ =Tridiagonal matrix
0 0 0 −2 1 0
0 0 0 ⋯ 1 −2 1
0 0 0 0 1 −2

[] []
T1 −T 0
T2 0
T3 0
T= ⋮ =Unknown Vector, G= ⋮ =known Vector
⋮ ⋮
T m −1 0
Tm −T m+1

( )[ ][ ]
−2 1 0 0 0 0 T1 −T 0
1 −2 0 ⋯ 0 0 0 T2 0
0 1 −2 0 0 0 T3 0
⋮ ⋱ ⋮ ∗ ⋮ = ⋮
0 0 0 −2 1 0 ⋮ ⋮
0 0 0 ⋯ 1 −2 1 T m−1 0
0 0 0 0 1 −2 Tm −T m+1
 Now above system of linear equation can be solved by any direct or
iterativemethod.

Analytical Solution
d2T
=0 ………. (1)
d x2

By Integration of equation (1) we get

T=C 1 x +C2 …… (2)
Apply boundary Conditions
x=0, T=T1
x=L, T=T2
Put x=0, T=T1 in equation (2) we get C 2=T1
T 2−T 1
Put x=L, T=T2 in equation (2) we get C 1=
L

By putting values of C 1 , C2 we get the solution of equation

T −T 1 x
=
T 2−T 1 L
Flow Chart
 Sci-lab Code
clear
clc
option=input('select 1. FVM, 2. FDM:')
if(option==1)
// FINITE VOLUME METHOD
l=input('Enter length in (m):')
Tb=input('Enter the left boundary temperature (in Kelvin):')
n=input('Input no of nodes:')
Ta=input('Enter the right boundary Temperature(in Kelvin):')
a(n,n)=0
T(n,1)=0
B(n,1)=0;
dx=l/n;
aw(n)=0
ae(n)=0
su(n)=0
sp(n)=0
ap(n)=0
// COEFFICIENTS OF EQUATION FOR NODE-1
fori=1
aw(i)=0
ae(i)=(1/dx);
su(i)=(2*Tb/dx);
sp(i)=-(2/dx);
ap(i)=aw(i)+ae(i)-sp(i);
end
// COEFFICIENTS OF EQUATIONS FOR MIDDLE INTERIOR NODES
fori=2:n-1
aw(i)=(1/dx);
ae(i)=(1/dx);
su(i)=0;
sp(i)=0;
ap(i)=aw(i)+ae(i)-sp(i);
end
// COEFFICIENTS OF EQUATION FOR LAST NODE
fori=n
aw(i)=(1/dx);
ae(i)=0
su(i)=(2*Ta/dx);
sp(i)=-(2/dx);
ap(i)=aw(i)+ae(i)-sp(i);
end
//CREATING TRIDIAGONAL MATRIX
fori=2:n-1
a(1,1)=ap(1,1);a(1,2)=-ae(1,1);
a(i,i)=ap(i,1);
a(i,i+1)=-ae(i,1);
a(i,i-1)=-aw(i,1);
a(n,n-1)=-aw(n,1);a(n,n)=ap(n,1);
end
for(i=1:n)
B(1)=su(1);
B(i)=0;
B(n)=su(n);
end
printf("1.TDMA,2.Jacobian,3.Seidel,4.SOR\n")
ch=input("Select:")
// TDMA METHOD
if(ch==1)
fori=2:n
forj=1:n
if(i==j)
a(i,j)=a(i,j)-(a(i,j-1)/a(i-1,j-1)*a(i-1,j));
end
end
end
fori=2:n
B(i)=B(i)-a(i,i-1)/a(i-1,i-1)*B(i-1);
end
T(n)=B(n)/a(n,n);
fori=n-1:-1:1
T(i)=(B(i)-T(i+1)*a(i,i+1))/a(i,i);
end
else//if ch is not 1
app=input("Enter initial approximation:")
tol=input("Enter tolerance:")
iter=input('Enter no. of iteration:')
if(ch==4)
omega=input("Enter the value of omega:")
end
TN(n,1)=0;
for(i=1:n)
T(i)=app;
end
for(m=1:iter)
for(i=2:n-1)
if(ch==2)
TN(1)=(1/ap(1))*(B(1)+ae(1)*T(2));//Temperature based on boundary
condition
TN(i)=(1/ap(i))*(aw(i)*T(i-1)+ae(i)*T(i+1)+B(i));
TN(n)=(1/ap(n))*(aw(n)*T(n-1)+B(n));// JACOBIAN METHOD
elseif(ch==3)
TN(1)=(1/ap(1))*(B(1)+ae(1)*T(2));//Temperature based on boundary
condition
TN(i)=(1/ap(i))*(aw(i)*TN(i-1)+ae(i)*T(i+1)+B(i));
TN(n)=(1/ap(n))*(aw(n)*TN(n-1)+B(n));// SEIDEL METHOD
elseif(ch==4)
TN(1)=((1-omega)*T(1))+(omega*(1/ap(1))*(B(1)+ae(1)*T(2)));//SOR
METHOD
TN(i)=((1-omega)*T(i))+(omega*(1/ap(i))*(aw(i)*TN(i-
1)+ae(i)*T(i+1)+B(i)));
TN(n)=((1-omega)*T(n))+(omega*(1/ap(n))*(aw(n)*TN(n-1)+B(n)));
end
end
err=abs(TN-T)/abs(TN)//Error Calculation
T=TN;
if(err<tol)// To eject from loop if error<tolerence
break
end
end
end
printf("The temperature at %d node are:\n",n)
for(i=1:n)
printf("Temperature%d:%f\n",i,T(i));
end
//ANALYTICAL SOLUTION
T(i)=0;
t1=Tb;
t2=Ta;
printf("The analytical temperature at %d nodes are:\n",n)
for(i=1:n)
T(i)=(((i-1/2)*dx)/l)*(t2-t1)+t1;
printf("%f\n",T(i));
end
//INPUTS
//Finite Difference Method
else
l=input("Enter the value of length(m):")
n=input("Number of nodes:")
Tb=input('enter the left side temperature in (kelvin):')
Ta=input('enter the right side temperature in (kelvin):')
dx=l/(n+1);
a=rand(n,n)//intialising arbitrary matrix "A"
fori=1:n//to create a tridiagonal matrix "A"
forj=1:n
if(i==j)thena(i,j)=-2;
elseif(j==i+1)thena(i,j)=1;
elseif(j==i-1)thena(i,j)=1;
elsea(i,j)=0;
end
end
end
T=rand(n,1);//temperature matrix
G=rand(n,1);//RHS of equation [a]*[T]=[G]
G(1)=-Tb,G(n)=-Ta;//boundary conditions
for(i=2:n-1)
G(i)=0;
end
printf("1.TDMA 2.Gauss-Seidel 3.Jacobian 4.Successive over relaxation\n")
S=input("SELECT:")
if(S==1)
fori=2:n//calculation of new diagonal elements
forj=1:n
if(i==j)
a(i,j)=a(i,j)-(a(i,j-1)/a(i-1,j-1)*a(i-1,j));//Based on thomas algorithm
end
end
end
fori=2:n//matrix operation on matrix[D]
G(i)=G(i)-a(i,i-1)/a(i-1,i-1)*G(i-1);//Based on thomas algorithm
end
T(n)=G(n)/a(n,n);//End node temperature
fori=n-1:-1:1//Calculating temperataur of other nodes in reverse order
T(i)=(G(i)-T(i+1)*a(i,i+1))/a(i,i);
end
else//if S is not 1
I=input("Enter initial approximation:")
tol=input("Enter tolerance:")
if(S==4)
w=input("Enter the value of relaxation parameter:")
end
8
TN=rand(n,1);
for(i=1:n)
T(i)=I;
end
for(m=1:1000)
TN(1)=(-0.5)*(G(1)-T(2));//temperature based on boundary conditio
for(i=2:n-1)
if(S==2)
TN(i)=(-0.5)*(G(i)-TN(i-1)-T(i+1));//gauss seidel
TN(n)=(-0.5)*(G(n)-TN(n-1));
elseif(S==3)
TN(i)=(-0.5)*(G(i)-T(i-1)-T(i+1));//jacobian method
TN(n)=(-0.5)*(G(n)-T(n-1));
elseif(S==4)
TN(1)=((1-w)*T(1))+(w*(-0.5)*(G(1)-T(2)));//
Succesiveoverrelaxation(SOR)
TN(i)=((1-w)*TN(i))+(w*(-0.5)*(G(i)-TN(i-1)-T(i+1)));
TN(n)=((1-w)*T(n))+(w*(-0.5)*(G(n)-TN(n-1)));
end
end
err=abs(TN-T)/abs(TN)//error calculation
T=TN;
if(err<tol)//to eject from loop
break
end
end
end
printf("the temperature at %d nodes are:\n",n)
for(i=1:n)
printf("T%d:%f\n",i,T(i));//calculated temperatures
end
//for analytical solution
m(n)=0;
t1=Tb
t2=Ta
printf("The analytical temperature at %d nodes are\n",n)
fori=1:n
m(i)=(i*dx/l)*(t2-t1)+t1;
printf("T%d:%f\n",i,m(i));//calculated analytical temperature
end
for(i=1:n)
 E=rand(n,1)
E(i)=(m(i)-T(i))
printf("Error at each node%d:%f\n",i,E(i))
end
end

Finite Volume Method

Inputs
1 Length 2 (m)
2 Number of Nodes 15
3 Initial Approximation 1
4 Tolerance 0.0001

Outputs
Nod TDM JACOBI SOR(W=1. ANALYT ERR. ERR
G-S ERR. GS
es A AN 5) ICAL JACOBI. SOR
1 486.67 486.55044 486.54 486.621986 486.66667 0.116224 0.128744 0.044681
2 460 459.65133 459.63 459.878174 460 0.348671 0.36956 0.121826
3 433.33 432.78232 432.75 433.150461 433.33333 0.551018 0.584029 0.182872
4 406.67 405.9133 405.9 406.439427 406.66667 0.753366 0.764033 0.22724
5 380 379.10933 379.1 379.744813 380 0.890675 0.903287 0.255187
6 353.33 352.30535 352.34 353.065649 353.33333 1.027985 0.99754 0.267684
7 326.67 325.59012 325.62 326.400392 326.66667 1.07655 1.044667 0.266275
8 300 298.87489 298.96 299.747076 300 1.125115 1.044667 0.252924
9 273.33 272.25678 272.33 273.103461 273.33333 1.07655 0.999574 0.229872
10 246.67 245.63868 245.75 246.46717 246.66667 1.027985 0.913282 0.199497
11 220 219.10933 219.21 219.835817 220 0.890675 0.791293 0.164183
12 193.33 192.57997 192.69 193.207118 193.33333 0.753365 0.640414 0.126215
13 166.67 166.11565 166.2 166.57898 166.66667 0.551019 0.468405 0.087687
14 140 139.65133 139.72 139.949571 140 0.348671 0.283602 0.050429
15 113.33 113.21711 113.24 113.317367 113.33333 0.116223 0.094534 0.015966
 600

500

400
Temperature

300

200

100

0
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15

No of Nodes

Nodes TDMA JACOBIAN

G-S SOR(W=1.5) ANALYTICAL

Figure 1: Temperature vs Number of nodes

1.2

0.8
Absolute Error

0.6

0.4

0.2

0
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15

No of Nodes

ERR. JACOBI. ERR. GS ERR SOR

Figure 2: Absolute Error vs No of Nodes

 Finite Difference Method

Inputs
S.No.    
1 Length 2 (m)
2 Number of Nodes 15
3 Initial Approximation 1
4 Tolerance 0.0001

Outputs

TDM JACOBIA SOR(W=1. ANALY ERR. ERR

Nodes G-S ERR. GS
A N 5) TICAL JACOBI. SOR

1 475 474.769514 474.738044 474.90785 475 0.230486 0.261956 0.092151

2 450 449.539027 449.496029 449.83039 450 0.460973 0.503971 0.169607
3 425 424.343631 424.282406 424.76896 425 0.656369 0.717594 0.231038
4 400 399.148234 399.104225 399.72409 400 0.851766 0.895775 0.275912
5 375 374.017674 373.966919 374.69558 375 0.982326 1.033081 0.30442
6 350 348.887114 348.874159 349.68262 350 1.112886 1.125841 0.317383
7 325 323.841267 323.827778 324.68386 325 1.158733 1.172222 0.316143
8 300 298.795421 298.827778 299.69755 300 1.204579 1.172222 0.302448
9 275 273.841267 273.872393 274.72167 275 1.158733 1.127607 0.278333
10 250 248.887114 248.958227 249.75399 250 1.112886 1.041773 0.246007
11 225 224.017674 224.080445 224.79226 225 0.982326 0.919555 0.20774
12 200 199.148234 199.233009 199.83423 200 0.851766 0.766991 0.165767
13 175 174.343631 174.408958 174.87779 175 0.656369 0.591042 0.122206
14 150 149.539027 149.600706 149.92102 150 0.460973 0.399294 0.078983
15 125 124.769514 124.800353 124.96222 125 0.230486 0.199647 0.037781
 500
450
Temperature 400
350
300
250
200
150
100
50
0
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15

No. of Nodes

Nodes TDMA JACOBIAN

G-S SOR(W=1.5) ANALYTICAL

Figure 3: Temperature vs Nodes

1.4

1.2
Absolute Error

0.8

0.6

0.4

0.2

0
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15

No. of Nodes

ERR. JACOBI. ERR. GS ERR SOR

Figure 4: Absolute Error vs No of Nodes