I am not certain what result you want. The freqz function only works on Signal Ptocessing Toolbox filter results. To use the Control System Toolbox to get the transfer function plots, do this instead:
sys = ss(A,B,E,D);
figure
bodeplot(sys)
grid
I added that at the end of my code here. (It is not going to be straightforward to interpret those results, at least for me.) I also added my one-sided Fourier transform function ‘FFT1’ since it may be convenient to do element-wise division of its outputs (complex) to determine whatever transfer function you want from the results of integrating your differential equations. I used it to calculate the Fourier transforms of your results, and then plotted them.
I am not certain what result you want, however this may get you started —
clear; clc;
global all_F
global ex_F
all_F = [];
ex_F = [];
syms z_T
h_f = 150;
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.263e5; % 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.64*10^6; % kg / Added mass for platform surge
I_a = 1.16*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 = 7.94e4; % N s2/m2 / damping in surge motion (x-axis translation)
c_sp = -2.25e5; % coupled damping value between surge and platform pitch
c_p = 5.60e8; % damping in pitch motionS 5
k_p = 2.190e9; % rotational stiffness of platform
k_mooringS = 7.965e4;
k_mooringSP = 1.162e6;
k_mooringP = 2.65e8;
tspan = 0:0.025:200;
% definition of TMD parameters
X0 = [0 0 10*pi/180 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.94e04; % 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;
C_1 = C+C_add;
% state space model--------------
O = zeros(3,3);
O1 = zeros(3,1);
N = ones(3,1);
I = eye(3);
A = [O I;-inv(M_1)*K_1 -inv(M_1)*C_1];
B = [O O; O inv(M_1)];
E = [I O; O O];
D = zeros(6);
% the forces and moments are extracted from Orcaflex
options = odeset('RelTol',1e-10,'AbsTol',1e-10);
[t,X] = ode45(@(t,X) reducedmodel(t,A,B,X),tspan,X0,options);
[b,a] = ss2tf(A,B,E,D,1);
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')
[FTx1,Fv1] = FFT1(X(:,1),t);
[FTx2,Fv2] = FFT1(X(:,2),t);
[FTPitch,Fv3] = FFT1(PtfmPitch_deg,t);
figure
tiledlayout(3,1)
nexttile
plot(Fv1, abs(FTx1)*2)
grid
set(gca, 'YScale','log')
ylabel('Amplitude')
title('X_1')
nexttile
plot(Fv1, abs(FTx2)*2)
grid
set(gca, 'YScale','log')
ylabel('Amplitude')
title('X_2')
nexttile
plot(Fv1, abs(FTx1)*2)
grid
set(gca, 'YScale','log')
xlabel('Frequency')
ylabel('Amplitude')
title('PtfmPitch\_deg')
sgtitle('Fourier Transforms')
sys = ss(A,B,E,D,1);
figure
bodeplot(sys)
grid
function dXdt = reducedmodel(t,A,B,X)
coes = 9.225e5;
coesp = 8.92e6;
coep = 1.68e10;
F1 = -(coes*X(5)*abs(X(5))-coesp*X(6)*abs(X(6)));
F2 = -(coep*X(6)*abs(X(6))-coesp*X(5)*abs(X(5)));
F=[0;0;0;0;F1;F2];
dXdt = A*X+B*F;
% disp(t)
end
function [FTs1,Fv] = FFT1(s,t)
t = t(:);
L = numel(t);
if size(s,2) == L
s = s.';
end
Fs = 1/mean(diff(t));
Fn = Fs/2;
NFFT = 2^nextpow2(L);
FTs = fft((s - mean(s)) .* hann(L).*ones(1,size(s,2)), NFFT)/sum(hann(L));
Fv = Fs*(0:(NFFT/2))/NFFT;
% Fv = linspace(0, 1, NFFT/2+1)*Fn;
Iv = 1:numel(Fv);
FTs1 = FTs(Iv,:);
end
EDIT — (23 May 2024 at 22:02)
Changed ‘Fourier Transform’ figure xlabel to correctly label it as ‘Frequency’.
.
