Computational Fluid Dynamics

Verified

Added on  2023/03/17

|19
|3639
|50
AI Summary
This document discusses the concept of Computational Fluid Dynamics and its application in power generation. It includes equations, code, and graphs to calculate power generated. The program calculates power for varying values of KC.

Contribute Materials

Your contribution can guide someone’s learning journey. Share your documents today.
Document Page
Running HEAD: COMPUTATIONAL FLUID DYNAMICS 1
Computational fluid dynamics.
Name
Institute of affiliation
Date

Secure Best Marks with AI Grader

Need help grading? Try our AI Grader for instant feedback on your assignments.
Document Page
COMPUTATIONAL FLUID DYNAMICS 2
KC = T Um
D
U= Water’s velocity in x direction
ρ: Water’s density
Cd: Drag force coefficient
Cm: Mass coefficient
A: = Area
Fx ( t )=D ( t )
¿ Cm+0.5 Cd AU |U |……………………….…………………………….Equation 1
velocity is defined as follows
u=U0 +sin ( 2 πt
T )Um ……………….………………………………………….. …Equation 2
U0= Steady flow velocity.
Um= flow oscillatory velocity
T= time period
Given that
Document Page
COMPUTATIONAL FLUID DYNAMICS 3
P= 1
T
0
T
c V 2 t ……………….…………….…………………….……………Equation 3
where P= power stored, and V is the Vibration of the system.
And
KC =U m T
D ……………….……………………………………………………..……Equation 4
It then follows that the power can be derived as follows
P=N1

n=1
N
C ( V n )
2……………….……………………………………………Equation 5
Power from the water is as follows
Fw=md CA
d V r
dt +0.5 ρ A projected CD V rV r¿……………………………………………..Equation
6
Total force Ftacting is sum of the force acting on the spring Fs, damper Fd and the water FW
Ft=FS + Fd +Fw ……………….………………………………..…………………Equation 7
Ft=kxc dx
dt +0.5 ρ Cd A projected|V r|V r +C A md
d2 x
d t2 ………….……………Equation 8
Since awater=¿ d2 x
d t2 equation above becomes
Ft=kxc dx
dt +C A md awater + 1
2 ρ Cd A p |V wV a|(V ¿¿ wV a )¿…………………….Equation 9
Differentiating equation 2 with respect t gives
aw=ω U m cos (ωt )……………………………………………………………………..Equation 10
awater withbe added into the f t equation¿ substitute other terms
Developing the 4th Runge kurta term
dx
dt =v………………………………………………………………………………Equation 11
dv
dt =f ( x , v )…………………………………………………………………………Equation 12
Document Page
COMPUTATIONAL FLUID DYNAMICS 4
Where,
f ( x , v , t )=
kx + (c + 1
2 ρ Cd Ap |v|)v
C A md + m
………………………………………………Equation 13
x(n+ 1) t =xn t + t
6 ( kx 1 +2 kx 2 +2 k x3 +k x 4 ) ………………………..Equation 14
v(n+1) t= ( kv 1 +2 kv 2 +2 kv 3 +k v 4 ) t
6 + vn t …………………………Equation 15
Where,
k xa
=vn t……………………………………………………………………………Equation 16
k xb
=vn t + t
2 k v1…………………………………………………………………Equation 17
k xc
=vn t + t
2 k v2…………………………………………………………..……Equation 18
k xd
=vn t + t k v3…………………………………………………………………Equation 19
k va
=f (xn t , vn t ) …………………………………………….………………….Equation 20
k vb
=f ( xn t + t
2 k x 1 , vn t + t
2 kv 1 ) ……………………………………..…..Equation 21
k vc
=f (xn t + t
2 k x2 , vn t + t
2 kv 2) ………………..………………………Equation 22
k vd
=f (xn t + t k x 3 , vn t + t kv 3 ) ………………………………………………Equation 23
In order to calculate the power from the equations above, the code below
was developed.
clear all
% Runge-Kutta method
j=7; % last digit of ID
U_m = 1+ (j/10); % Flow velocity Amplitude
U_0= 0.1 * U_m; %Velocity of steady flow
T = .2; % oscillatory flow period
D = (10+j)/100; % Diameter of cylinder
% Cylinder length
Length_cyl = 1;
mass_c = 50; % cylinder mass
rho_water = 1000; % water density

Secure Best Marks with AI Grader

Need help grading? Try our AI Grader for instant feedback on your assignments.
Document Page
COMPUTATIONAL FLUID DYNAMICS 5
K_s = 200; % spring constant
c_d = 100; % damping coefficient
C_inertia = 1; % inertial constant
C_drag = 1.8; % Drag coefficient
om_ega = 2*pi/T; % Angular frequency of the flow
md = rho_water*pi*(D*D)/4*Length_cyl; % coefficient due to the added mass
A_projected = D*Length_cyl; % Area
delta_t = T/40; % step in time
n_delta_t = T/delta_t*50; % the total time
X_length(1:n_delta_t+1) = 0; % displacement from a step
response in time
V(1:n_delta_t+1) = 0; % velocity from a step response in
time
time(1:n_delta_t+1) = (0:n_delta_t)*delta_t;
for k=1:n_delta_t
% finding the kx1 and kv1
t_a = time(k); % defining the accelearation
a_water = om_ega*U_m*cos(om_ega*t_a); % water acceleration
V_w = U_0+U_m*sin(om_ega*t_a); % water velocity
V_a = V(k);
X_a = X_length(k);
t_1 = C_inertia*md*a_water;
t_2 = 0.5*rho_water*C_drag*A_projected*abs(V_w-V_a)*(V_w-V_a);
t_3 = -K_s*X_a-c_d*V_a;
k_x_a = V(k);
k_v_a = (t_1+t_2+t_3)/(mass_c+C_inertia*md);
% finding the k_x_b and k_v_b
t_a = time(k)+delta_t/2;
a_water = om_ega*U_m*cos(om_ega*t_a); % acceleration of the
water
V_w = U_0+U_m*sin(om_ega*t_a); % velocity of the water
V_a = V(k)+0.5*delta_t*k_v_a;
X_a = X_length(k)+0.5*delta_t*k_x_a;
t_1 = C_inertia*md*a_water;
t_2 = 0.5*rho_water*C_drag*A_projected*abs(V_w-V_a)*(V_w-V_a);
t_3 = -K_s*X_a-c_d*V_a;
k_x_b = V(k)+0.5*delta_t*k_v_a;
k_v_b = (t_1+t_2+t_3)/(mass_c+C_inertia*md);
% to calculate k_x_c and k_v_c
t_a = time(k)+delta_t/2;
a_water = om_ega*U_m*cos(om_ega*t_a); % acceleration of the
water
V_w = U_0+U_m*sin(om_ega*t_a); % velocity of the water
V_a = V(k)+0.5*delta_t*k_v_b;
X_a = X_length(k)+0.5*delta_t*k_x_b;
t_1 = C_inertia*md*a_water;
t_2 = 0.5*rho_water*C_drag*A_projected*abs(V_w-V_a)*(V_w-V_a);
t_3 = -K_s*X_a-c_d*V_a;
k_x_c = V(k)+0.5*delta_t*k_v_b;
k_v_c = (t_1+t_2+t_3)/(mass_c+C_inertia*md);
% to calculate k_x_d and
k_v_d
t_a = time(k)+delta_t;
a_water = om_ega*U_m*cos(om_ega*t_a); % acceleration of the
water
V_w = U_0+U_m*sin(om_ega*t_a); % velocity of the
water
V_a = V(k)+0.5*delta_t*k_v_c;
X_a = X_length(k)+0.5*delta_t*k_x_c;
t_1 = C_inertia*md*a_water;
t_2 = 0.5*rho_water*C_drag*A_projected*abs(V_w-V_a)*(V_w-V_a);
Document Page
COMPUTATIONAL FLUID DYNAMICS 6
t_3 = -K_s*X_a-c_d*V_a;
k_x_d = V(k)+delta_t*k_v_c;
k_v_d = (t_1+t_2+t_3)/(mass_c+C_inertia*md);
% values from next step
X_length(k+1) = X_length(k)+(delta_t/6)*(k_x_a+2*k_x_b+2*k_x_c+k_x_d);
V(k+1) = V(k)+(delta_t/6)*(k_v_a+2*k_v_b+2*k_v_c+k_v_d);
end
fileID=fopen('output.txt','w');
for k= 1:n_delta_t+1
fprintf(fileID,'%12.6 %12.6 %12.6\r\n', time(k),X_length(k),V(k));
end
fclose(fileID);
%---------------------------------------------------------------------
add=0;
for m = 1:n_delta_t+1
add = add + 1/n_delta_t*c_d*(V(m)*V(m));
end
value= add;
% developing the power graph
subplot(2,2,1)
plot(time(1:n_delta_t),X_length(1:n_delta_t))
title('displacement against time graph')
xlabel('time.....(s)');
ylabel('X length.....(m)');
subplot(2,2,2);
plot(time(1:n_delta_t),V(1:n_delta_t))
xlabel('time.............seconds');
ylabel('V..............(m/s)');
title('velocity against time graph')
19.918218 is the value of the power generated
Document Page
COMPUTATIONAL FLUID DYNAMICS 7
Published with MATLAB® R2018b
For constant value of C and varying value of kc
KC =U m T
D
D=0.18
T= [0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 1.8 2.0]
Um = 1.8
Making T the subject of the formula
T = KC D
U m
Substituting the values of D, Um the formula becomes as follows
T = KC0.18
1.8
Below is the program to calculate the power developed from values of KC from 0 to 20
clear all
% Runge-Kutta method
j=7; % last digit of ID
U_m = 1+ (j/10); % Flow velocity Amplitude
U_0= 0.1 * U_m; %Steady flow velocity
T = .2; % oscillatory flow period
D = (10+j)/100; % cylinder diameter
% length of the cylinder
Length_cyl = 1;
mass_c = 50; % cylinder mass
rho_water = 1000; % water density
K_s = 200; % spring constant
c_d = 100; % damping coefficient
C_inertia = 1; % inertial constant
C_drag = 1.8; % coefficient of drag
om_ega = 2*pi/T; % flow`s angular frequency
md = rho_water*pi*(D*D)/4*Length_cyl; % coefficient due to the added mass
A_projected = D*Length_cyl; % projected area
delta_t = T/40; % step in time
n_delta_t = T/delta_t*50; % the total time
X_length(1:n_delta_t+1) = 0; % displacement from a step
response in time
V(1:n_delta_t+1) = 0; % velocity from a step response in
time

Paraphrase This Document

Need a fresh take? Get an instant paraphrase of this document with our AI Paraphraser
Document Page
COMPUTATIONAL FLUID DYNAMICS 8
time(1:n_delta_t+1) = (0:n_delta_t)*delta_t;
for k=1:n_delta_t
% finding the kx1 and kv1
t_a = time(k); % defining the acceleration
a_water = om_ega*U_m*cos(om_ega*t_a); % water acceleration
V_w = U_0+U_m*sin(om_ega*t_a); % water velocity
V_a = V(k);
X_a = X_length(k);
t_1 = C_inertia*md*a_water;
t_2 = 0.5*rho_water*C_drag*A_projected*abs(V_w-V_a)*(V_w-V_a);
t_3 = -K_s*X_a-c_d*V_a;
k_x_a = V(k);
k_v_a = (t_1+t_2+t_3)/(mass_c+C_inertia*md);
% finding the k_x_b and k_v_b
t_a = time(k)+delta_t/2;
a_water = om_ega*U_m*cos(om_ega*t_a); % acceleration of the
water
V_w = U_0+U_m*sin(om_ega*t_a); % velocity of the water
V_a = V(k)+0.5*delta_t*k_v_a;
X_a = X_length(k)+0.5*delta_t*k_x_a;
t_1 = C_inertia*md*a_water;
t_2 = 0.5*rho_water*C_drag*A_projected*abs(V_w-V_a)*(V_w-V_a);
t_3 = -K_s*X_a-c_d*V_a;
k_x_b = V(k)+0.5*delta_t*k_v_a;
k_v_b = (t_1+t_2+t_3)/(mass_c+C_inertia*md);
% to calculate k_x_c and k_v_c
t_a = time(k)+delta_t/2;
a_water = om_ega*U_m*cos(om_ega*t_a); % acceleration of the
water
V_w = U_0+U_m*sin(om_ega*t_a); % velocity of the water
V_a = V(k)+0.5*delta_t*k_v_b;
X_a = X_length(k)+0.5*delta_t*k_x_b;
t_1 = C_inertia*md*a_water;
t_2 = 0.5*rho_water*C_drag*A_projected*abs(V_w-V_a)*(V_w-V_a);
t_3 = -K_s*X_a-c_d*V_a;
k_x_c = V(k)+0.5*delta_t*k_v_b;
k_v_c = (t_1+t_2+t_3)/(mass_c+C_inertia*md);
% to calculate k_x_d and
k_v_d
t_a = time(k)+delta_t;
a_water = om_ega*U_m*cos(om_ega*t_a); % acceleration of the
water
V_w = U_0+U_m*sin(om_ega*t_a); % velocity of the
water
V_a = V(k)+0.5*delta_t*k_v_c;
X_a = X_length(k)+0.5*delta_t*k_x_c;
t_1 = C_inertia*md*a_water;
t_2 = 0.5*rho_water*C_drag*A_projected*abs(V_w-V_a)*(V_w-V_a);
t_3 = -K_s*X_a-c_d*V_a;
k_x_d = V(k)+delta_t*k_v_c;
k_v_d = (t_1+t_2+t_3)/(mass_c+C_inertia*md);
% values from next step
X_length(k+1) = X_length(k)+(delta_t/6)*(k_x_a+2*k_x_b+2*k_x_c+k_x_d);
V(k+1) = V(k)+(delta_t/6)*(k_v_a+2*k_v_b+2*k_v_c+k_v_d);
end
fileID=fopen('output.txt','w');
for k= 1:n_delta_t+1
fprintf(fileID,'%12.6 %12.6 %12.6\r\n', time(k),X_length(k),V(k));
end
fclose(fileID);
Document Page
COMPUTATIONAL FLUID DYNAMICS 9
%---------------------------------------------------------------------
add=0;
for m = 1:n_delta_t+1
add = add + 1/n_delta_t*c_d*(V(m)*V(m));
end
value= add;
% developing the power graph
subplot(2,2,1)
plot(time(1:n_delta_t),X_length(1:n_delta_t))
title('displacement against time graph')
xlabel('time.....(s)');
ylabel('X length.....(m)');
subplot(2,2,2);
plot(time(1:n_delta_t),V(1:n_delta_t))
xlabel('time.............seconds');
ylabel('V..............(m/s)');
title('velocity against time graph')
24.832199 is the value of power generated
Published with MATLAB® R2014b
The table below gives the power generated for Different values of KC is given in table below:
KC T P
2.00 0.200 19.9
4.00 0.400 21.9
6.00 0.600 24.8
8.00 0.800 28.2
10.00 1.000 31.8
12.00 1.200 35.2
14.00 1.400 38.4
16.00 1.600 41.3
18.00 1.800 43.9
20.00 2.000 46.3
Document Page
COMPUTATIONAL FLUID DYNAMICS
10
KC=2 and kc =4 figures

Secure Best Marks with AI Grader

Need help grading? Try our AI Grader for instant feedback on your assignments.
Document Page
COMPUTATIONAL FLUID DYNAMICS
11
Figure for KC=6 and for KC=8
0 5 10 15 20 25
0
5
10
15
20
25
30
35
40
45
50
KC values against value of Power
KC
Power
Document Page
COMPUTATIONAL FLUID DYNAMICS
12
The graph above shows power increases gradually from 19.91822 to 46.31243.
Holding KC constant at a value of 10, and conducting simulations for c=0 to 2000 N·m/s with at
an interval of 50 power variation becomes as in the table below
KC T P
2.00 0.200 19.9
4.00 0.400 21.9
6.00 0.600 24.8
8.00 0.800 28.2
10.00 1.000 31.8
12.00 1.200 35.2
14.00 1.400 38.4
16.00 1.600 41.3
18.00 1.800 43.9
20.00 2.000 46.3
Power for given k values
c power
0.0 0.0
50.0 17.4
100.0 31.8
150.0 43.4
200.0 52.7
250.0 59.9
300.0 65.4
350.0 69.4
400.0 72.2
450.0 74.2
500.0 75.4
550.0 76.0
600.0 76.2
650.0 76.1
700.0 75.7
750.0 75.1
800.0 74.3
850.0 73.5
900.0 72.5
950.0 71.5
1000.0 70.5
1050.0 69.4
Document Page
COMPUTATIONAL FLUID DYNAMICS
13
1100.0 68.3
1150.0 67.2
1200.0 66.1
1250.0 65.1
1300.0 64.0
1350.0 62.9
1400.0 61.9
1450.0 60.9
1500.0 60.9
1550.0 58.9
1600.0 58.0
1650.0 57.0
1700.0 56.1
1750.0 55.2
1800.0 54.4
1850.0 53.5
1900.0 52.7
1950.0 51.9
2000.0 51.1

Paraphrase This Document

Need a fresh take? Get an instant paraphrase of this document with our AI Paraphraser
Document Page
COMPUTATIONAL FLUID DYNAMICS
14
For c=0 and c= 100
For c= 500 and c= 1200
Document Page
COMPUTATIONAL FLUID DYNAMICS
15
Graph of C against power
0
100
200
300
400
500
600
700
800
900
1000
1100
1200
1300
1400
1500
1600
1700
1800
1900
2000
0
10
20
30
40
50
60
70
80
90
Graph of C against power
KC
C
Document Page
COMPUTATIONAL FLUID DYNAMICS
16
Discussion of results
C has increased up to 70.47486 at 1000 value then it started decreasing.
QUESTION NO>2
Part 2
1-D equation can be wriiten as shown below
dT
dt = Dx ( d2 c )
d x2 ………………………………………………………………………Equation 24
When having the initial conditions shown below
T(x,0) = ¿ …………………………………………………………Equation 25
T(0,t) = T 0……………………………………………………………………………Equation 26
T(L,t) = T 0……………………………………………..……………………………Equation 27
For better numerical stability the following is used
x =x / L
t'= t
T
Where L is the length of the aquifer and T is the time of simulation
dT
dt = dT
dt
d t'
dt = 1
T
dT
d t' ………………………………………………………………………Equation
28
Now equation 4 can be rewritten as
1
T
dT
d t' =Dx
1
l2
d2 T
d x'2 ………………………………………………………………………Equation 29
Since T = L2
Dx
equation 9 now becomes
dT
d t' = d2 T
d x' 2 ……………………………………………………………………Equation 30
Discretizing the equation of diffusion

Secure Best Marks with AI Grader

Need help grading? Try our AI Grader for instant feedback on your assignments.
Document Page
COMPUTATIONAL FLUID DYNAMICS
17
dT
d t' =T i
n+1 Ci
n
Δt …………………………………………………………Equation 31
d2 T
d x' 2 Ti +1
n 2Ci
n +T i1
n
( Δx ) 2 …………………………………………………………Equation 32
Combining the equation 11 and 12 we have
ci
n+ 1=Ci
n + Δ T
( Δ x ) 2 (Ti +1
n 2Ci
n +T i1
n )
2nd instance
Temperature diffusion in soil
T =Soil temperature
t = time
u = Velocity of water
α =Coefficient of heat diffusion
T=30°C at the top boundary
T
x =0 at the bottom boundary (at x= 2 m).
% part 2
j=8; % student id
U = 0.002*(1+j); %flwo velocitiy
temperature = 10+j; %initial condition
depth=2; % depth
m= 200;
delta_x= depth/m;
x = (0:m)*delta_x;
alpha = 0.0002;
time = 0;
delta_time = delta_x/U/100; % time flow for one cell
N = depth/2/U/delta_time;
w = U*delta_time/(2*delta_x);
p = (alpha*delta_time)/(delta_x*delta_x);
T(1:m+1) = temperature;
T(1)= 30; %initial condition
for n=1:N
time = time+delta_time;
Document Page
COMPUTATIONAL FLUID DYNAMICS
18
T1 = T;
for k=2:m
T(k) = (w+p)*T1(k-1)+(1-2*p)*T1(k)+(-w+p)*T1(k+1);
end
T(m+1) = T(m);
end
plot(x,T,'o');
xlabel('length x (m)');
ylabel('Temperature (degree Celsius)');
s = find(T > 25 & T <25.5);
Time = delta_time*min(s);
fprintf('all temperatures at 1.5 m depth for the time = %f sec\n',Time)
At 1.5 deep on time = 0.544444 sec, temperature got to 25 degrees
Soil temperature increases with time due to diffusion
Document Page
COMPUTATIONAL FLUID DYNAMICS
19
1 out of 19
circle_padding
hide_on_mobile
zoom_out_icon
[object Object]

Your All-in-One AI-Powered Toolkit for Academic Success.

Available 24*7 on WhatsApp / Email

[object Object]