Taking too long to get the output (Wilson theta method)

1 view (last 30 days)
Could I know where I might be wrong in the following code as it is taking too long get the result.
WilsonMethod.m
function [depl,vel,accl,t] = WilsonMethod(M,K,C,R)
clc ;
sdof = length(K) ;
% Time step and time duration
ti = 0. ;
dt = 0.1 ;
tf = 500 ;
t = ti:dt:tf ;
nt = length(t) ;
% Initialize the displacement,velocity and acceleration matrices
depl = zeros(sdof,nt) ;
vel = zeros(sdof,nt) ;
accl = zeros(sdof,nt) ;
Reff = zeros(sdof,nt) ;
% Initial conditions
depl(:,1) = zeros ;
vel(:,1) = zeros ;
accl(:,1) =zeros;% M\(R-K*depl(:,1)-C*vel(:,1)) ;
% Integration constants
tita = 1.4 ; % Can be changed
a0 = 6/(tita*dt)^2 ; a1 = 3/(tita*dt) ; a2 = 2*a1 ;
a3 = tita*dt/2 ; a4 = a0/tita ; a5 = -a2/tita ;
a6 = 1-3/tita ; a7 = dt/2 ; a8 = dt^2/6 ;
% Form Effective Stiffness Matrix
Keff = K+a0*M+a1*C ;
%Time step starts
for it = 1:nt-1
% Calculating Effective Load
Reff(:,it) = R(:,it)+tita*(R(:,it)-R(:,it))+M*(a0*depl(:,it)+a2*vel(:,it)+2*accl(:,it))+....
C*(a1*depl(:,it)+2*vel(:,it)+a3*accl(:,it)) ;
% Solving for displacements at time (t+dt)
depl(:,it+1) = Keff\Reff(:,it) ;
% Calculating displacements, velocities and accelerations at time t+dt
accl(:,it+1) = a4*(depl(:,it+1)-depl(:,it))+a5*(vel(:,it))+a6*accl(:,it) ;
vel(:,it+1) = vel(:,it)+a7*(accl(:,it+1)+accl(:,it)) ;
depl(:,it+1) = depl(:,it)+dt*vel(:,it)+a8*(accl(:,it+1)+2*accl(:,it)) ;
end
Ex.m
MA = [8070000,0,-629468070;0,8070000,112980;-629468070,112980,6.800000000000000e+10];
Ad = [8.053095218400001e+06,0,-4.831857131040000e+08;0,2.167940435676214e+05,0;-4.831857131040000e+08,0,3.865485704832000e+10];
Ca = [0,0,0;0,3.241885080000000e+05,0;0,0,1.301151158327999e+09];
Cm = [4.12e+04,0,-2.82e+06;0,1.19e+04,0;-2.82e6,0,3.11e+08];
M = MA+Ad;
K = Ca+Cm;
C = zeros(size(K)) ; % Damping Matrix
Fg = -79086000; %Gravitational force
Fbuoy = 7.844814740000000e+07; %Buoyancy force
Fp = 2.712318560000001e+06; %Heave force
profile on
t = 0:0.1:500
for i = 1:length(t)
Fh = hydro(t(i))
FhT = transpose(Fh)
R(:,i) = [-334731.8545+27939.6+6.5*10^5+FhT(i);-3517000+Fg+Fbuoy+Fp;-112510430.2+3.44*10^6+266.5*10^5+(FhT(i)*18)];
end
profile off
profile viewer
[depl,vel,accl,t] = WilsonMethod(M,K,C,R) ;
depl'
figure(1), clf
plot(t,depl(1,:)), xlabel('time(s)'), ylabel('surge(m)')
title ('Surge vs Time')
figure(2), clf
plot(t,depl(2,:)), xlabel('time(s)'), ylabel('heave(m)')
title ('heave vs Time')
figure(3), clf
plot(t,depl(3,:)), xlabel('time(s)'), ylabel('Pitch(deg)')
title ('Pitch vs Time')
  3 Comments
Satish Jawalageri
Satish Jawalageri on 20 Jul 2020
Thanks for your reply. Here is the code:
hydro.m
function [Fh] = hydro(t)
Cd = 0.6;
Ca = 0.9699;
R = 4.7;
t=0:0.1:500;
for i=1:length(t)
vel = compute_wavevelocity(t(i));
acc = compute_waveacceleration(t(i));
fun = sym((0.5*997*Cd*2*R*abs(vel)*vel)+(Ca*997*pi*(R^2)*acc)+(pi*(R^2)*997*acc));
Fh(1,i) = double(int(fun,0,-120));
end
plot(t,Fh), xlabel('time(s)'), ylabel('Fh')
title ('Fh vs Time')
end
compute_wavevelocity.m:
function velocity = compute_wavevelocity(t)
T = 10;
WN = (9.8*(T^2))/(2*pi());
k = (2*pi())/WN;
omega = (2*pi())/T;
y = [320;310;300;290;280;270;260;250;240;230;220;210;200];
d = 320;
H = 4;
vel = (omega*H/2)* ((cosh(k*y))/(sinh(k*d)))* cos((k*0)-(omega*t));
velocity = mean(vel,1);
end
compute_waveacceleration.m:
function acceleration = compute_waveacceleration(t)
T = 10;
WN = (9.8*(T^2))/(2*pi());
k = (2*pi())/WN;
omega = (2*pi())/T;
y = [320;310;300;290;280;270;260;250;240;230;220;210;200];
d = 320;
H = 4;
acc = ((omega^2)*H/2)* ((cosh(k*y))/(sinh(k*d)))* sin((k*0)-(omega*t));
acceleration = mean(acc,1);
end

Sign in to comment.

Accepted Answer

Serhii Tetora
Serhii Tetora on 22 Jul 2020
Without loops it runs faster. Please check it carefully
clc;close all;clear;
MA = [8070000,0,-629468070;0,8070000,112980;-629468070,112980,6.800000000000000e+10];
Ad = [8.053095218400001e+06,0,-4.831857131040000e+08;0,2.167940435676214e+05,0;-4.831857131040000e+08,0,3.865485704832000e+10];
Ca = [0,0,0;0,3.241885080000000e+05,0;0,0,1.301151158327999e+09];
Cm = [4.12e+04,0,-2.82e+06;0,1.19e+04,0;-2.82e6,0,3.11e+08];
M = MA+Ad;
K = Ca+Cm;
C = zeros(size(K)); % Damping Matrix
Fg = -79086000; %Gravitational force
Fbuoy = 7.844814740000000e+07; %Buoyancy force
Fp = 2.712318560000001e+06; %Heave force
t = 0:0.1:500;
Fh = hydro(t);
FhT = transpose(Fh);
R(1,:) = -334731.8545 + 27939.6 + 6.5e5 + FhT;
R(2,:) = -3517000 + Fg + Fbuoy + Fp;
R(3,:) = -112510430.2 + 3.44e6 + 266.5e5 + 18*FhT;
[depl,vel,accl,t] = WilsonMethod(M,K,C,R) ;
% depl';
figure(1), clf
plot(t,depl(1,:)), xlabel('time(s)'), ylabel('surge(m)'), grid on
title ('Surge vs Time')
figure(2), clf
plot(t,depl(2,:)), xlabel('time(s)'), ylabel('heave(m)'), grid on
title ('heave vs Time')
figure(3), clf
plot(t,depl(3,:)), xlabel('time(s)'), ylabel('Pitch(deg)'), grid on
title ('Pitch vs Time')
function [depl,vel,accl,t] = WilsonMethod(M,K,C,R)
sdof = length(K) ;
% Time step and time duration
ti = 0. ;
dt = 0.1 ;
tf = 500 ;
t = ti:dt:tf ;
nt = length(t) ;
% Initialize the displacement,velocity and acceleration matrices
depl = zeros(sdof,nt) ;
vel = zeros(sdof,nt) ;
accl = zeros(sdof,nt) ;
Reff = zeros(sdof,nt) ;
% Initial conditions
depl(:,1) = zeros ;
vel(:,1) = zeros ;
accl(:,1) =zeros;% M\(R-K*depl(:,1)-C*vel(:,1)) ;
% Integration constants
tita = 1.4 ; % Can be changed
a0 = 6/(tita*dt)^2 ; a1 = 3/(tita*dt) ; a2 = 2*a1 ;
a3 = tita*dt/2 ; a4 = a0/tita ; a5 = -a2/tita ;
a6 = 1-3/tita ; a7 = dt/2 ; a8 = dt^2/6 ;
% Form Effective Stiffness Matrix
Keff = K+a0*M+a1*C ;
%Time step starts
for it = 1:nt-1
% Calculating Effective Load
Reff(:,it) = R(:,it)+tita*(R(:,it)-R(:,it))+M*(a0*depl(:,it)+a2*vel(:,it)+2*accl(:,it))+....
C*(a1*depl(:,it)+2*vel(:,it)+a3*accl(:,it)) ;
% Solving for displacements at time (t+dt)
depl(:,it+1) = Keff\Reff(:,it) ;
% Calculating displacements, velocities and accelerations at time t+dt
accl(:,it+1) = a4*(depl(:,it+1)-depl(:,it))+a5*(vel(:,it))+a6*accl(:,it) ;
vel(:,it+1) = vel(:,it)+a7*(accl(:,it+1)+accl(:,it)) ;
depl(:,it+1) = depl(:,it)+dt*vel(:,it)+a8*(accl(:,it+1)+2*accl(:,it)) ;
end
end
function [Fh] = hydro(t)
Cd = 0.6;
Ca = 0.9699;
R = 4.7;
vel = compute_wavevelocity(t);
acc = compute_waveacceleration(t);
fun = sym((0.5*997*Cd*2*R*abs(vel).*vel)+(Ca*997*pi*(R^2).*acc)+(pi*(R^2)*997.*acc));
Fh = double(int(fun,0,-120));
end
function velocity = compute_wavevelocity(t)
T = 10;
WN = (9.8*(T^2))/(2*pi);
k = 2*pi/WN;
omega = 2*pi/T;
y = [320;310;300;290;280;270;260;250;240;230;220;210;200];
d = 320;
H = 4;
vel = (omega*H/2)* ((cosh(k*y))/(sinh(k*d)))* cos((k*0)-(omega.*t));
velocity = mean(vel,1);
end
function acceleration = compute_waveacceleration(t)
T = 10;
WN = (9.8*(T^2))/(2*pi);
k = 2*pi/WN;
omega = 2*pi/T;
y = [320;310;300;290;280;270;260;250;240;230;220;210;200];
d = 320;
H = 4;
acc = ((omega^2)*H/2)* ((cosh(k*y))/(sinh(k*d)))* sin((k*0)-(omega.*t));
acceleration = mean(acc,1);
end

More Answers (0)

Categories

Find more on Programming in Help Center and File Exchange

Products


Release

R2019b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!