How to extract the velocity from the differential equation and use the velocity to calculate the external excitations
3 visualizaciones (últimos 30 días)
Mostrar comentarios más antiguos
Yu
el 5 de Abr. de 2024
Comentada: Yu
el 29 de Abr. de 2024
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
0 comentarios
Respuesta aceptada
Abhishek Chakram
el 23 de Abr. de 2024
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
Más respuestas (0)
Ver también
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!