Simple Initial Value Problem

Part one: simple initial value problem
Introduction
The presentation of an energy extraction device that is made up of the cylinder and elastically
mounted to an electric generator where the extraction of power is done using a damper of
coefficient damping as shown the figure one are factors to be considered in this part. The flow
direction is normal to the cylinder. the flowing fluid past cylinder is modelled to comprise of an
oscillatory flow as well as a steady flow having amplitude Um. The MATLAB dynamics
simulation is used in the prediction the problem caused by the varying parameters such as
damping coefficient (c) as well as the displacement Keulegen-carpenter KC, vibration velocity
and extraction resultant power. These problems need a critical hydrodynamics acting forces
analysis on the cylinder thus give an elaborate knowledge of forces which act on any water
vessel with respect to water forces.
Figure 1-1
Definition of motion equation.
Deriving procedure of equation of the motion is illustrated below

Draw the cylinder’s free body diagram first as in the figure 1-2 below.
In the diagram 1-2 shows an indication of three forces acting on the cylinder which include:
Fwater= water forces.
Fdamping= force of damping = - cV (Cv is negative since it is against the cylinder motion)
Fspring=force of the spring = - KX
FTotal = Fwater + Fdamping + Fspring
The marision equation illustrate that Fwater
FWATER=C A md
d V r
dt + 1
2 ρC D AP |V r|V r
Where
 CA= Coefficient of inertia force
 CD= coefficient of drag force
 Md= the mass of the displaced fluid
 Ap= The cylinder DL project area
 Vr= The velocity of water as relate to cylinder similar to Vr equal to u –v
o U is the incoming velocity equivalent to u=U0 +Umsin(ωt)
o V is the velocity of the cylinder.
How the newton second law is applied.
∑ F=ma=m dv
dt

FTotal= m dV
dt
Fwater + Fdamping + Fspring = m dV
dt
m dV
dt =C A
md
d V r
dt + 1
2 ρ CD A P|V r|V r−cV −KX
m dV
dt =C A
md ( du
dt −dV
dt )+ 1
2 ρCD AP |u V |(u−V )−cV −KX
m dV
dt +C A md
dV
dt =C A md
du
dt + 1
2 ρCD AP |u V |(u−V )−cV −KX
dV
dt = 1
m+C A md (C A md
du
dt + 1
2 ρ CD AP|uV |(u−V )−cV −KX )
∴ The illustration of the cylinder motion equation
d2 X
d t2 = 1
m+C A md (CA md
du
dt + 1
2 ρC D AP |u dX
dt |(u− dX
dt )−C dX
dt −KX )

Numerical method.
The application of the Rung Kutta method of fourth order is used in this part in the calculation of
the numerical solution of the second order of differential equations required when changing the
first two order equations
dX
dt =V
dV
dt =f ¿V, X, t)
Where f ( V , X , t ) = 1
m+C A md ( CA md
du
dt + 1
2 ρ CD AP |u V |(u−V )−CV −Kx )
The fourth order- Kutta method is
X ( n+1 ) ∆ t= Xn ∆ t + ∆ t
6 (k X 1+2 k X 2+ 2 kX 3 +k X 4 )
V ( n+1 ) ∆ t =V n ∆ t + ∆ t
6 (k V 1+ 2k V 2+ 2 kV 3 +kV 4 )
where

MATLAB program for solving the problem
% A MATLAB program for solving the problem of the vibration of a
% circular cylinder in a water flow
% The code below is use the Fourth Runge-Kutta method
clear all;
pi = acos (-1); % pi = 3.14159
Um = 1.4; % Fluid Velicity amplitude
Uo=0.14; % Uniform flow vwlocity
T =1.2; % oscillatory flow period KC=10
D = 0.14; % cylinder diamter in meter
L = 1; % cylinder length
m = 50; % cylinder mass

rho = 1024; % fluid density
K = 200; % spring stiffness
c =100; % damping constant
CA = 1; % added mass coefficient
CD = 1.8; % drag coefficient
omega = 2*pi/T; % angular frequency of the flow
md = rho*pi*(D*D)/4*L; % added mass coefficient
Ap = D*L; % projected area
dt = T/40; % time step
ndt = (T/dt*5)*12 ; % %total step to be calculated
X(1:ndt+1) = 0; % save displacement X from step 0 to step ndt
V(1:ndt+1) = 0;
PW(1:ndt+1)=0;
time(1:ndt+1) = (0:ndt)*dt;
for n=1:ndt
% to calculate kx1 and kv1
ta = time(n);
aw = omega*Um*cos(omega*ta); % acceleration of the fluid
Vw = Uo+Um*sin(omega*ta); % velocity of the fluid
Va = V(n);
Xa = X(n);
t1 = CA*md*aw;
t2 = 0.5*rho*CD*Ap*abs(Vw-Va)*(Vw-Va);

t3 = -K*Xa-c*Va;
kx1 = V(n);
kv1 = (t1+t2+t3)/(m+CA*md);
% to calculate kx2 and kv2
ta = time(n)+dt/2;
aw = omega*Um*cos(omega*ta); % acceleration of the fluid
Vw =Uo+Um*sin(omega*ta); % velocity of the fluid
Va = V(n)+0.5*dt*kv1;
Xa = X(n)+0.5*dt*kx1;
t1 = CA*md*aw;
t2 = 0.5*rho*CD*Ap*abs(Vw-Va)*(Vw-Va);
t3 = -K*Xa-c*Va;
kx2 = V(n)+0.5*dt*kv1;
kv2 = (t1+t2+t3)/(m+CA*md);
% to calculate kx3 and kv3
ta = time(n)+dt/2;
aw = omega*Um*cos(omega*ta); % acceleration of the fluid
Vw = Uo+Um*sin(omega*ta); % velocity of the fluid
Va = V(n)+0.5*dt*kv2;
Xa = X(n)+0.5*dt*kx2;
t1 = CA*md*aw;
t2 = 0.5*rho*CD*Ap*abs(Vw-Va)*(Vw-Va);

t3 = -K*Xa-c*Va;
kx3 = V(n)+0.5*dt*kv2;
kv3 = (t1+t2+t3)/(m+CA*md);
% to calculate kx4 and kv4
ta = time(n)+dt;
aw = omega*Um*cos(omega*ta); % acceleration of the fluid
Vw = Uo+Um*sin(omega*ta); % velocity of the fluid
Va = V(n)+0.5*dt*kv3;
Xa = X(n)+0.5*dt*kx3;
t1 = CA*md*aw;
t2 = 0.5*rho*CD*Ap*abs(Vw-Va)*(Vw-Va);
t3 = -K*Xa-c*Va;
kx4 = V(n)+dt*kv3;
kv4 = (t1+t2+t3)/(m+CA*md);
% next step values
X(n+1) = X(n)+(dt/6)*(kx1+2*kx2+2*kx3+kx4);
V(n+1) = V(n)+(dt/6)*(kv1+2*kv2+2*kv3+kv4);
%power loop
PW(n+1)=1/ndt*(c*V(n)*V(n));
end
P=0;
for k=1:ndt