Time- series inputs for ODE Function
5 views (last 30 days)
Show older comments
Dear all,
I hope you are all doing well.
I am trying to introduce the external inputs (forces) data which is time series. I have refered to a number of Matlab discussions and used the interp1 method to calculate in the ode function.
However, I found the results of the interplotion are different with the original data which I used. Furthermore, I have no idea about how to include the external inputs in ode. It does not accept the time-series data when I use it directly. I have attached the file and my code is following:
In the script, I just try to find the problem so, only 1 input (force) which is F_wave is enabled.
Thank you for your kind help.
Best wishes,
Yu
clear; clc;
global all_F
global ex_F
all_F = [];
ex_F = [];
syms z_T
filename = 'Inputs_Force.xlsx';
data = readtable(filename);
use_data=table2cell(data(:,1:6));
use_data=cell2mat(use_data);
Waveforce = use_data(:,2);
Wavemoment =use_data(:,3);
Windforce = use_data(:,5);
Windmoment = use_data(:,6);
tspan = 0:0.025:200;
figure,
subplot(3,1,2), plot(tspan, Windforce), xlabel('time/s'),ylabel('WindForce/N');
subplot(3,1,1), plot(tspan, Waveforce), xlabel('time/s'),ylabel('WaveForce/N');
subplot(3,1,3), plot(tspan, Wavemoment), xlabel('time/s'),ylabel('WaveMoment/N.m')
h_R = 144.582; % m the height from the MSL to tower top;
H_T = 129.582; % m tower height from the tower bottom;
h_T = 15; % m height from tower base to tower bottom;
h = 29.94;
h_1 = 164.94;
z = 14.94; % m the distance from rotational centre to mooring line
h_t = 92.61; % m the distance from the centre of gravity of tower to rotational centre
h_p = 14.94; % m the distance from the centre of gravity of platform to rotational centre
g = 9.81;
m = 20093000; % kg /total mass
m_T = 1.263e6; % 8.6e5; % kg / tower
m_N = 1.017e6; % kg / nacelle
m_p = 1.7838e7; % kg /platform mass
xi_TA = 0.01;
I_p = 1.2507*10^10;% kg m2 /mass moment of inertia of platform
m_as = 9.4*10^6; % kg / Added mass for platform surge
I_a = 1.13*10^10; % kg m2 / Added mass for platform pitch
m_asp = -1.01*10^8; % kg m / Coupling effects of added mass bewteen surge and pitch
c_s = 1e5; % N s2/m2 / damping in surge motion (x-axis translation)
c_sp = -2e5; % coupled damping value between surge and platform pitch
c_p = 6e8; % viscous damping in pitch motionS
k_p = 2.190e9; % rotational stiffness of platform
k_mooringS = 7.965e4;
k_mooringSP = 1.162e6;
k_mooringP = 2.65e8;
% definition of TMD parameters
m_TMD = 1.2e5; % kg
k_TMD = 6.064e3; % kg/m
c_TMD = 1.2678e4; % kg/(m s)
X0 = [0 0 0 0 0 0]; % initial pitch motion
% ==================== definition of the tower properties===============
mu = 0.0084*z_T^3-1.077*z_T^2-171.5*z_T+2.96e04; % mass per length
EI = 1.905e06*z_T^3-2.47e08*z_T^2-5.208e10*z_T+6.851e12; % tower bending stiffness
Phi_TA = 0.9414*(z_T/H_T)^2+0.3468*(z_T/H_T)^3-1.073*(z_T/H_T)^4+1.3139*(z_T/H_T)^5-0.5289*(z_T/H_T)^6; % tower fore-aft first mode shape
% mass component
fun1 = mu*Phi_TA^2; m_TA = double(int(fun1,z_T,0,H_T));
fun2 = mu*Phi_TA; m_1 = double(int(fun2,z_T,0,H_T));
fun3 = mu*(z_T)*Phi_TA; m_2 = double(int(fun3,z_T,0,H_T));
fun4 = mu*(z_T); m_3 = double(int(fun4,z_T,0,H_T));
fun5 = mu*(z_T)^2; m_4 = double(int(fun5,z_T,0,H_T));
% fun6 = mu; m_T = double(int(fun6,z_T,0,H_T));
% Stiffness of the tower
D2y = diff(Phi_TA,z_T,2); Dy = diff(Phi_TA,z_T,1);
fun6 = EI*D2y^2; f1 = int(fun6,0,H_T);
fun7 = mu; f2 = int(fun7,z_T,H_T);
fun8 = g*(m_N+f2)*Dy^2; f3 = int(fun8,0,H_T);
k_TA = double(f1-f3);
% ================= definition of matrices ======================
M = [m_N+m_TA m_N+m_1 m_N*h_R+m_2;
m_N+m_1 m_N+m_p+m_T (m_N*h_R-m_p*h_p+m_3+m_T*h_T);
m_N*h_R+m_2 (m_N*h_R-m_p*h_p+m_3) m_N*h_R^2+I_p+m_4];
C = [2*xi_TA*sqrt(m_TA*k_TA) 0 0;
0 0 0;
0 0 0];
K = [k_TA 0 -(m_N+m_T)*g;
0 0 0;
-(m_N+m_T)*g 0 -(m_N*h_R+m_3-m_p*h_p)*g];
M_add = [0 0 0;
0 m_as m_asp;
0 m_asp I_a];
C_add = [0 0 0; % problem_
0 c_s c_sp;
0 c_sp c_p];
K_add = [0 0 0;
0 0 0;
0 0 k_p];
K_mooring = [0 0 0;
0 k_mooringS k_mooringSP;
0 k_mooringSP k_mooringP];
M_1 = M+M_add;
K_1 = K+K_mooring+K_add;
% the forces and moments are extracted from Orcaflex
F_wave = Waveforce;
M_wave = Wavemoment;
F_aero = Windforce;
M_aero = Windmoment;
% [V,D,W] = eig(K_1,M_1);
%
% w = diag(D).^0.5;
%
% T = (2*pi./w);
options = odeset('RelTol',1e-10,'AbsTol',1e-10);
[t,X] = ode45(@(t,X) reducedmodel(z,h_R,tspan,t,X,M,M_add,C,C_add,K,K_add,K_mooring,F_wave,M_wave,F_aero,M_aero),tspan,X0,options);
PtfmPitch_deg = X(:,3)*180/pi;
figure,
subplot(3,1,1), plot(t,X(:,1)),grid, xlabel('time/ s'), ylabel('TTDspFA/ m')
subplot(3,1,2), plot(t,X(:,2)),grid, xlabel('time/ s'), ylabel('surge/ m')
subplot(3,1,3), plot(t,PtfmPitch_deg),grid, xlabel('time/ s'), ylabel('platform pitch/ deg')
function dXdt = reducedmodel(z,h_R,tspan,t,X,M,M_add,C,C_add,K,K_add,K_mooring,F_wave,M_wave,F_aero,M_aero)
global all_F
global ex_F
x = X(1:3);
xdot = X(4:6);
coe = 9.225e5;
coe1 = 1.16e10;
F1 = -coe*xdot(2)*abs(xdot(2));
F2 = -coe1*xdot(3)*abs(xdot(3));
F3 = interp1(tspan, F_wave, t,'spline');
F4 = interp1(tspan, M_wave, 'spline');
F5 = interp1(tspan, F_aero, 'spline');
F6 = interp1(tspan, F_aero*h_R, 'spline');
F_external = [0;F3;0];
F=[0;F1;F2];
all_F = [all_F,F]; % drag force
ex_F = [ex_F, F_external]; % wave and wind forces
xddot = (M+M_add)\(F+F_external-(K+K_add+K_mooring)*x-(C+C_add)*xdot);
dXdt = [xdot; xddot];
end
% function [u,wn]=eigsort(K_1,M_1)
% Omega=sqrt(eig(K_1,M_1));
% [vtem,~]=eig(K_1,M_1);
% [wn,isort]=sort(Omega);
% il=length(wn);
% for i=1:il
% v(:,i)=vtem(:,isort(i));
% end
% disp("The natural frequencies are (rad/sec)")
% wn
% end
0 Comments
Accepted Answer
UDAYA PEDDIRAJU
on 6 May 2024
Hi Yu,,
You've used 'spline' interpolation in some calls to interp1 but not provided the t argument in some cases, which would cause an error. Ensure you always pass the current time t as the third argument to interp1 for it to work correctly.
function dXdt = reducedmodel(z,h_R,tspan,t,X,M,M_add,C,C_add,K,K_add,K_mooring,F_wave,M_wave,F_aero,M_aero)
global all_F
global ex_F
x = X(1:3);
xdot = X(4:6);
% Interpolate the time-series inputs at the current time step
F3 = interp1(tspan, F_wave, t,'spline');
F4 = interp1(tspan, M_wave, t,'spline');
F5 = interp1(tspan, F_aero, t,'spline');
F6 = interp1(tspan, F_aero*h_R, t,'spline');
F_external = [0; F3; 0];
coe = 9.225e5;
coe1 = 1.16e10;
F1 = -coe*xdot(2)*abs(xdot(2));
F2 = -coe1*xdot(3)*abs(xdot(3));
F = [0; F1; F2];
all_F = [all_F, F]; % drag force
ex_F = [ex_F, F_external]; % wave and wind forces
xddot = (M+M_add) \ (F + F_external - (K+K_add+K_mooring)*x - (C+C_add)*xdot);
dXdt = [xdot; xddot];
end
Try this it should work, I hope this helps!
0 Comments
More Answers (0)
See Also
Categories
Find more on Ordinary Differential Equations in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!