# Plot frequency responses based on the ODE results

17 views (last 30 days)

Show older comments

Yu
on 23 May 2024

Commented: Star Strider
on 28 May 2024 at 13:22

Hello all,

I would like to plot the frequency responses based on the ode results to see the behaviors of the system. I have already established the state-space model for my ode equations. Additionally, I use ss2tf to get the transfer fucntion. Now I am to plot the frequency responses of all degrees of freedom on one figure. However, 'freqz' does not work in my case.

I checked the details from turtorial. unfortunately, I failed. I hope you can help me with this. My codes can be seen in the following. I appreciate your help.

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')

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

Best wishes,

Yu

##### 0 Comments

### Accepted Answer

Star Strider
on 23 May 2024

Edited: Star Strider
on 23 May 2024 at 22:02

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’.

.

##### 14 Comments

Star Strider
on 28 May 2024 at 13:22

### More Answers (0)

### See Also

### Categories

### Community Treasure Hunt

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

Start Hunting!