How to extract the velocity from the differential equation and use the velocity to calculate the external excitations

Dear all,

I hope you are all doing well. I am using ode45 to sovle the differential equations without any external forces. It works well. However, I get problems in employing the external forces. Specifically, the forces should be relavant to the absolute value of the relative velocity of the wave and structure from the differential equations. I have no idea about how to using the velocity of each time step to calculate the force then apply the force to the structure (equations). My codes are following:

clear all; clc, close all
syms z_0 z_1
% ============ definition of the parameters ==============
% q1 = q_T tower... q2 = q_s surge... q3 = q_p platform pitch
rho = 1025; % kg/m3
g = 9.81;

H_T = 77.6; % m checked
h_R = 90; % m checked
h_T = 10; % m checked
h_G = 89.9155; % m checked
h_b = 30.0845; % m checked
h_BG = 62.69; % m
m_N = 110000+240000; % kg checked
m_T = 249718; % kg checked
m_p = 7446330; % kg checked
J_p = 4.22923e9; % kgm^2 checked
Z_1 = 10; % location of tower bottom
xi_TTFA = 0.01; % damping ratio of the tower
C_m = 1.969954; % checked
C_d = 0.6; % drag coefficient1
k1 = 41180; % N/m mooring stiffness in surge checked
k2 = -2.821*10^6; % N/rad mooring stiffness of coupled effect between surge and pitch
k3 = 3.111e8 -4.99918e9; % Nm/rad mooring stiffness in pitch motion3.111*10^8 
c0 = 10^5; % N/(m/s) additional damping in surge checked

% q_s spar surge; q_p spar pitch; q_T TTFA deformation; q_TT TMD disp
% m_N the mass of RNA; m_TT the mass of TMD; m_p mass of spar
% m_T mass of the tower
% m_TTFA modal mass
% J_p inertia moment of spar
% h_R the undisturbed height of nacelle with reference to the MSL
% h_T the undisturbed distance between tower bottom and the MSL
% h_G the undisturbed distance from the MSL to the spar COG
% h_BG undisturbed vertical distance from the spar COG to the center of buoyancy (COB)
% h_b the distance between the platform bottom and the spar COG
% H_T the tower length
% z_1 the tower segment height in the local coordinate system
% z_1 the location of tower bottom
% z_0 the horizontal force acting at the spar location
% mu(z_1) the tower mass density per length
% Phi_TTFA(z_1) the fundamental mode shape of tower in the fore-aft
% k_TTFA bending stiffness of fore-aft displacement of tower
% C_m hydrodynamic inertia coefficients

tspan = 0:0.1:200;

X0 = [0 0 10*pi/180 0 0 0];

% ======================== definition of the tower =======================

mu = 0.093*z_1^2-42.18*z_1+4667; % mass per unit length ?

EI = -2.235e5*z_1^3+8.971e7*z_1^2-1.345e10*z_1+7.296e11;

 % bending stiffness of the tower segment ?

Phi_TTFA = 0.8689*(z_1/H_T)^2+0.2205*(z_1/H_T)^3-0.0908*(z_1/H_T)^4+0.1167*(z_1/H_T)^5-0.1154*(z_1/H_T)^6;

% mass components

fun1 = mu*Phi_TTFA^2; m_TTFA = double(int(fun1,z_1,0,H_T)); % checked

fun2 = mu*Phi_TTFA; m_a = double(int(fun2,z_1,0,H_T)); % checked

fun3 = mu*(z_1+h_T)*Phi_TTFA; m_b = double(int(fun3,z_1,0,H_T)); % checked

fun4 = mu*(z_1+h_T); m_c = double(int(fun4,z_1,0,H_T)); % checked

fun5 = mu*(z_1+h_T).^2; m_d = double(int(fun5,z_1,0,H_T)); % checked

% bending stiffnes of the tower

D2y = diff(Phi_TTFA,z_1,2); Dy = diff(Phi_TTFA,z_1,1);

fun6 = EI*D2y^2; f1 = int(fun6,0,H_T);
fun7 = mu; f2 = int(fun7,z_1,H_T);
fun8 = g*(m_N+f2)*Dy^2; f3 = int(fun8,0,H_T);
k_TTFA = double(f1-f3);

% ============== definition of hydrodynamic properties ============


fun9 = @(z_0) diameter(z_0).^2;
fun10 = @(z_0) diameter(z_0).^2.*z_0;
fun11 = @(z_0) diameter(z_0).^2.*z_0.^2;

m_as = rho*(pi/4*(C_m-1)*integral(fun9,-h_b-h_G,0)); % added mass on surge motion
m_asp = rho*(pi/4*(C_m-1)*integral(fun10,-h_b-h_G,0)); % coupled added mass pitch, surge
m_ap = rho*(pi/4*(C_m-1)*integral(fun11,-h_b-h_G,0)); % added mass pitch 



% ============== definition of matrices ==================
C_add = [0 0 0;
         0 c0 0;
         0 0 0];
M = [m_N+m_TTFA  m_N+m_a             (m_N*h_R+m_b);
     m_N+m_a     m_N+m_T+m_p         m_N*h_R+m_c-m_p*h_G;
     (m_N*h_R+m_b) m_N*h_R+m_c-m_p*h_G m_N*h_R^2+m_d+m_p*h_G^2+J_p];

M_hydro = [0 0 0;
           0 m_as m_asp;
           0 m_asp m_ap];
C = [2*xi_TTFA*sqrt(m_TTFA*k_TTFA) 0 0; 
     0 0 0; 
     0 0 0];
K = [k_TTFA 0 -(m_N+m_T)*g; 
     0 0 0;
    -(m_N+m_T)*g 0 -(m_N*h_R+m_c-m_p*h_G)*g];
K_moor = [0 0 0;
          0 k1 k2;
          0 k2 k3];


% =============== solve the equations and plot the results =============

options = odeset('RelTol',1e-10,'AbsTol',1e-10);

[t,X] = ode45(@(t,X) odefn(tspan,t,X,M,M_hydro,C,C_add,K,K_moor),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')

% ================== definition of the equations =========================

function dXdt = odefn(t,tspan,X,M,M_hydro,C,C_add,K,K_moor)

syms z_0

rho = 1025; C_d = 0.6; h_G = 89.9155; % m checked
h_b = 30.0845; % m checked

dis = X(1:3)';

vel = X(4:6)'; % get the velocity

% ================== definition of the drag force ========================

v_rel = - (vel(:,2) + vel(:,3) * z_0);

fun12 = @(z_0) diameter(z_0) * abs(v_rel);
fun13 = @(z_0) diameter(z_0) * z_0 * abs(v_rel);

c1 = double(1/2*rho*C_d*int(fun12,z_0,-h_b-h_G,0)); % Morison force1
c2 = double(1/2*rho*C_d*int(fun13,z_0,-h_b-h_G,0)); % Morison force2 force2


F1 = interp1(tspan,c1,t);
F2 = interp1(tspan,c2,t);

F = [0;F1;F2];

x = X(1:3);

xdot = X(4:6);

xddot = (M+M_hydro)\(F-(K+K_moor)*x-(C+C_add)*xdot);

dXdt = [xdot; xddot];


end

% ============== definition of D(z) ===============================
function D=diameter(z_0)
    D = (9.5).*(z_0<-12 & z_0>=-120)+((-0.3825)*z_0+4.59).*(z_0>=-12 & z_0<-4)+(6.5).*(z_0>=-4 & z_0<=0);
end

I appreciate your support.

Thank you.

Best wishes,

Yu

Answer 1

[Abhishek Chakram]

Hi Yu,

To incorporate the external forces of the relative velocity between the wave and structure into the differential equations, you need to adjust the part of the code that calculates the drag force. This involves using the velocity at each time step to compute the force and then applying this force to your structure (equations).

Your current approach of calculating the drag force seems to be on the right path, but there is a misunderstanding in the use of “interp1” for “F1” and “F2”. Since “F1” and “F2” are calculated based on the current state “X”, there is no need to interpolate these forces over “tspan”. Instead, these forces should be directly calculated and applied at each time step by the ODE solver, “ode45”. Here is the corrected code for the same:

clear all; clc; close all;
syms z_0 z_1

% ============ definition of the parameters ==============
% q1 = q_T tower... q2 = q_s surge... q3 = q_p platform pitch
rho = 1025; % kg/m3
g = 9.81;
H_T = 77.6; % m checked
h_R = 90; % m checked
h_T = 10; % m checked
h_G = 89.9155; % m checked
h_b = 30.0845; % m checked
h_BG = 62.69; % m
m_N = 110000+240000; % kg checked
m_T = 249718; % kg checked
m_p = 7446330; % kg checked
J_p = 4.22923e9; % kgm^2 checked
Z_1 = 10; % location of tower bottom
xi_TTFA = 0.01; % damping ratio of the tower
C_m = 1.969954; % checked
C_d = 0.6; % drag coefficient1
k1 = 41180; % N/m mooring stiffness in surge checked
k2 = -2.821*10^6; % N/rad mooring stiffness of coupled effect between surge and pitch
k3 = 3.111e8 -4.99918e9; % Nm/rad mooring stiffness in pitch motion3.111*10^8 
c0 = 10^5; % N/(m/s) additional damping in surge checked

tspan = 0:0.1:200;
X0 = [0 0 10*pi/180 0 0 0];

% ======================== definition of the tower =======================
mu = 0.093*z_1^2-42.18*z_1+4667; % mass per unit length ?
EI = -2.235e5*z_1^3+8.971e7*z_1^2-1.345e10*z_1+7.296e11;
 % bending stiffness of the tower segment ?
Phi_TTFA = 0.8689*(z_1/H_T)^2+0.2205*(z_1/H_T)^3-0.0908*(z_1/H_T)^4+0.1167*(z_1/H_T)^5-0.1154*(z_1/H_T)^6;
% mass components
fun1 = mu*Phi_TTFA^2; m_TTFA = double(int(fun1,z_1,0,H_T)); % checked
fun2 = mu*Phi_TTFA; m_a = double(int(fun2,z_1,0,H_T)); % checked
fun3 = mu*(z_1+h_T)*Phi_TTFA; m_b = double(int(fun3,z_1,0,H_T)); % checked
fun4 = mu*(z_1+h_T); m_c = double(int(fun4,z_1,0,H_T)); % checked
fun5 = mu*(z_1+h_T).^2; m_d = double(int(fun5,z_1,0,H_T)); % checked
% bending stiffnes of the tower
D2y = diff(Phi_TTFA,z_1,2); Dy = diff(Phi_TTFA,z_1,1);
fun6 = EI*D2y^2; f1 = int(fun6,0,H_T);
fun7 = mu; f2 = int(fun7,z_1,H_T);
fun8 = g*(m_N+f2)*Dy^2; f3 = int(fun8,0,H_T);
k_TTFA = double(f1-f3);
% ============== definition of hydrodynamic properties ============
fun9 = @(z_0) diameter(z_0).^2;
fun10 = @(z_0) diameter(z_0).^2.*z_0;
fun11 = @(z_0) diameter(z_0).^2.*z_0.^2;
m_as = rho*(pi/4*(C_m-1)*integral(fun9,-h_b-h_G,0)); % added mass on surge motion
m_asp = rho*(pi/4*(C_m-1)*integral(fun10,-h_b-h_G,0)); % coupled added mass pitch, surge
m_ap = rho*(pi/4*(C_m-1)*integral(fun11,-h_b-h_G,0)); % added mass pitch 
% ============== definition of matrices ==================
C_add = [0 0 0;
         0 c0 0;
         0 0 0];
M = [m_N+m_TTFA  m_N+m_a             (m_N*h_R+m_b);
     m_N+m_a     m_N+m_T+m_p         m_N*h_R+m_c-m_p*h_G;
     (m_N*h_R+m_b) m_N*h_R+m_c-m_p*h_G m_N*h_R^2+m_d+m_p*h_G^2+J_p];
M_hydro = [0 0 0;
           0 m_as m_asp;
           0 m_asp m_ap];
C = [2*xi_TTFA*sqrt(m_TTFA*k_TTFA) 0 0; 
     0 0 0; 
     0 0 0];
K = [k_TTFA 0 -(m_N+m_T)*g; 
     0 0 0;
    -(m_N+m_T)*g 0 -(m_N*h_R+m_c-m_p*h_G)*g];
K_moor = [0 0 0;
          0 k1 k2;
          0 k2 k3];
% =============== solve the equations and plot the results =============
options = odeset('RelTol',1e-10,'AbsTol',1e-10);
[t,X] = ode45(@(t,X) odefn(t,X,M,M_hydro,C,C_add,K,K_moor),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')

% ================== definition of the equations =========================
function dXdt = odefn(t,X,M,M_hydro,C,C_add,K,K_moor)
    rho = 1025; C_d = 0.6; h_G = 89.9155; % m checked
    h_b = 30.0845; % m checked
    dis = X(1:3)';
    vel = X(4:6)'; % get the velocity

    % Corrected approach to calculate the drag force
    % Define the functions fun12 and fun13 to include the velocity calculations
    fun12 = @(z_0) diameter(z_0) .* abs(-vel(2) - vel(3) .* z_0);
    fun13 = @(z_0) diameter(z_0) .* z_0 .* abs(-vel(2) - vel(3) .* z_0);
    
    c1 = 1/2*rho*C_d*integral(fun12,-h_b-h_G,0, 'ArrayValued', true); % Morison force1
    c2 = 1/2*rho*C_d*integral(fun13,-h_b-h_G,0, 'ArrayValued', true); % Morison force2
    
    F = [0; c1; c2]; % Directly use calculated forces

    x = X(1:3);
    xdot = X(4:6);
    xddot = (M+M_hydro)\(F-(K+K_moor)*x-(C+C_add)*xdot);
    dXdt = [xdot; xddot];
end

% ============== definition of D(z) ===============================
function D=diameter(z_0)
    D = (9.5).*(z_0<-12 & z_0>=-120)+((-0.3825)*z_0+4.59).*(z_0>=-12 & z_0<-4)+(6.5).*(z_0>=-4 & z_0<=0);
end

Best Regards,

Abhishek Chakram

[Yu]

Hello Abhishek,

Many thanks for your answer and I apologize for my late reply as I am busy recently. Your modification helps me a lot. I just need to adjust the velocity component.

I have another question but I am not sure if you can give me some suggestions on it. That is:

I am calculating the natural frequency and natural period of the structure, however, the results are quite different when I use 'eig'. It is obvious different as you can compare with the period from the plots (you also have generated).

Is the 'eig' function not suitable for this situation?

Again, thank you for your help.

Best wishes,

Yu