Jeffs Code

% Jeff Wood - Homework 7

clear all
close all
clc


% ============ PART I ============


rho = 1.99; % Density of seawater (slug/ft^3)
g = 32.2; % Force of gravity (ft/s^2)
d = 100; % Water depth (ft)
L = 700; % Wave length (ft)
D = 1.6; % Pile diameter (ft)
DeltaZ = 1; % Unit change in depth (ft)
Cm = 1.8; % Mass Coefficient
Cd = .7; % Drag Coefficient
H = 16; % Wave Height (ft)
z = -25; % Depth being analyzed (ft)
k = 2*pi/L; % Wave Number
w = sqrt(k*tanh(k*d) *g); % Wave frequency (rad/s)
zeta = H/2; % Wave Amplitude (ft)
Knot = 1.688; % Knots in ft/s

% Analyzing linear wave theory for proper water depth:
WD = d/L;
Water_d = sprintf('It is appropriate to use Finite water because d/L is: %1.4f', WD);
disp(Water_d)
disp('Finite Water:(d/L)<0.5')
disp(' ')

% Identifying if Morison is applicable:
Mor = D/L; % YES
Morison = sprintf('Morison is applicable because D/L is: %1.4f', Mor);
disp(Morison)
disp('Requirement: (D/L) < 0.2')
disp(' ')

% ---- PART 1-B ----

% Now the wave period needs to be found. From rearranging the transitional
% wavelength equation:
%T = (2*pi)/w;
T=13.82;

% The horizontal water particle velocity can now be analyzed:
ux = w*zeta*((cosh(k*(z+d))/(sinh(k*d))));%*sin(wt-kx)
vel = sprintf('Horizontal Velocity Component: %3.4f', ux);
disp(vel)
disp(' ')

% The horizontal water particle acceleration:
ax = (w^2)*zeta*((cosh(k*(z+d))/(sinh(k*d))));%*cos(wt-kx)
acc = sprintf('Horizontal Acceleration Component: %3.4f', ax);
disp(acc)
disp(' ')

InertF = (rho*((pi*D^2)/4)*DeltaZ*Cm*ax);
Inertial_Force = sprintf('Inertial Force: %3.4f', InertF);
disp(Inertial_Force)
disp(' ')

% The two knot current is added to the velocity in the drag force.
DragF = (.5*rho*Cd*D*DeltaZ*(ux+2*Knot)*(ux+2*Knot));
Drag_Force = sprintf('Drag Force: %3.4f', DragF);
disp(Drag_Force)
disp(' ')

disp('DeltaF = 14.3789*sin(wt-kx) + 67.2357*cos(wt-kx)')
disp(' ')

disp('As you can see, if the phase is not included, the hand calculations')
disp('match up with the matlab evaluation of the forces.')

% ---- Part 1-C ----

x = 1;
delz = 0:1:100;
Dist = [100:-1:0]'*ones(1,139);
t=0:0.1:T;

%t = linspace(0,T,139);
A = (pi*(D^2))/4; % Cross sectional area of the piling (ft^2)

tmp = w*zeta*((cosh(k*(-delz+d)))/(sinh(k*d))); %tmp variable so we don't redo work on the next to calculations
a = w*tmp'*cos(w*t-k*x);
u = tmp'*sin(w*t-k*x);

InertF == rho*A*DeltaZ*Cm*a;

tmp = u+2*Knot*ones(size(u));  %tmp variable so we don't redo work on the next to calculations
DragF = (.5)*rho*Cd*D*DeltaZ*abs(tmp).*tmp;

Inert = sum(InertF);
Drag = sum(DragF);
MF = (Inert + Drag);

Inert_M = InertF.*Dist;
Drag_M = DragF.*Dist;
Inert_Mom = sum(Inert_M);
Drag_Mom = sum(Drag_M);
MF_Mom = Inert_M + Drag_M;
Total_Mom = sum(MF_Mom);

% ---- PART 1-A ----

figure(1)
plot(t,MF,'-r',t,Inert,'-g',t,Drag,'-b')
title('Part 1-A Plot')
xlabel('Time (s)')
ylabel('Forces (lbs)')
legend('Total Force','Inertial Force','Drag Force')
grid on

% ---- PART 1-B ----

figure(2)
plot(t,Total_Mom,'-r',t,Inert_Mom,'-g',t,Drag_Mom,'-b')
title('Part 1-B Plot')
xlabel('Time (s)')
ylabel('Moment (lb-ft)')
legend('Total Moment','Inertial Moment','Drag Moment')
grid on

% ---- PART 1-C ----

figure(3)
plot(t,Total_Mom*.02,'-r',t,MF,'-g')
title('Part 1-C Plot')
xlabel('Time (s)')
ylabel('Force (lb)/Moment (lb-ft)')
legend('Total Moment (Scaled)','Total Force Force')
grid on