Time- series inputs for ODE Function

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

Answer 1

[UDAYA PEDDIRAJU]

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!