Time- series inputs for ODE Function

2 次查看（过去 30 天）
Yu
Yu 2024-5-1
回答： UDAYA PEDDIRAJU 2024-5-6
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

请先登录，再回答此问题。

采纳的回答

UDAYA PEDDIRAJU
UDAYA PEDDIRAJU 2024-5-6
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 个）

类别

在 Help Center 和 File Exchange 中查找有关 Ordinary Differential Equations 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by