Trouble solving a 3-DOF equation of motion with MATLABs ODE45

Dear all,

I'm struggeling with solving a 3-DOF equation of motion of a floating wind turbine with Matlab's ODE45. More specifically, after running the simulation I imediately get results in the range of 10^290 and I have not been able to figure out what goes wrong. Could anyone assist me and help me identify where it all goes south? All tips are welcome, also general tips about solving this problem are highly appreciated since documentation on a 3-DOF equation of motion is quite limited..

Thank you in advance!

The equation that I'm trying to solve is:

x(t) is the 3 DOF motion vector containing the motion in surge, heave and pitch direction.

M =

A =

B =

And C =

The external forcing vector is a time dependent force which are all the forces and moments acting on the turbine.

I have reproduced the following Matlab tutorial for multiple EOM's, but I don't seem to get any valid results: https://nl.mathworks.com/help/matlab/math/solve-equations-of-motion-for-baton.html.

Below I will display how I have formed the code at the moment. First, I have defined the matrices as part of a structure P:

clear all
close all
clc

%M + A matrix
P.Mas = zeros(3,3);
P.Mas(1,1) = (8.07+7.98)*10^6;
P.Mas(1,2) = 0;
P.Mas(1,3) = (-6.29-4.94)*10^8;
P.Mas(2,1) = 0;
P.Mas(2,2) = 8.07*10^6+2.23*10^4;
P.Mas(2,3) = 1.12*10^5;
P.Mas(3,1) = (-6.29-4.94)*10^8;
P.Mas(3,2) = 1.12*10^5;
P.Mas(3,3) = (6.8+3.97)*10^10;

%B matrix
P.B = zeros(3,3);
P.B(1,1) = 1*10^5;
P.B(2,2) = 1.3*10^5;

%C matrix 
P.C = zeros(3,3);
P.C(2,2) = 3.34*10^5;
P.C(3,3) = -5.01*10^9;

The ODE45 part with the special mass matrix function (from the https://nl.mathworks.com/help/matlab/math/solve-equations-of-motion-for-baton.html. tutorial) is specified below. The time dependent forces are F_krylov, F_morison and Moment. I will post the code for the force calculation below, but I can imagine that is not something where an outsider can see an easy mistake in. To give some indication:

F_krylov is in the range e^6 [N]

F_morison is in the range e^6 [N]

Moment is in the range e^7 [N/m]

y0 = [0; 0; 0; 0; 1; 0];                    %Initial conditions
opts = odeset('Mass',@(t,q) mass(t,q,P));   
[t,q] = ode45(@(t,q) f(t,q,P,F_krylov,time,F_morison,Moment),time,y0,opts);

figure
plot(t, q(:,5))
xlabel('time [s]')
ylabel('Pitch [deg]')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function dydt = f(t,q,P,time,F_krylov,F_morison,Moment); % Equations to solve

F_V = interp1(time,F_krylov,t)
F_x = interp1(time,F_morison,t)
M_t = interp1(time,Moment,t)

% RHS of system of equations
dydt = [q(2)
        F_x-P.B(1,1)*q(2)-P.B(1,2)*q(4)-P.B(1,3)*q(6)-P.C(1,1)*q(1)-P.C(1,2)*q(3)-P.C(1,3)*q(5)
        q(4)
        F_V-P.B(2,1)*q(2)-P.B(2,2)*q(4)-P.B(2,3)*q(6)-P.C(2,2)*q(3)-P.C(2,1)*q(1)-P.C(2,3)*q(5)
        q(6)
        M_t-P.B(3,1)*q(2)-P.B(3,2)*q(4)-P.B(3,3)*q(6)-P.C(3,3)*q(5)-P.C(3,1)*q(1)-P.C(3,2)*q(3)];
end
%------------------------------------------------
function M = mass(t,q,P,Mas); % Mass matrix function

% Set nonzero elements in mass matrix
M = zeros(6,6);
M(1,1) = 1;
M(2,2) = P.Mas(1,1);
M(2,4) = P.Mas(1,2)
M(2,6) = P.Mas(1,3);
M(3,3) = 1;
M(4,2) = P.Mas(2,1);
M(4,2) = P.Mas(2,2)
M(4,6) = P.Mas(2,3);
M(5,5) = 1;
M(6,2) = P.Mas(3,1);
M(6,4) = P.Mas(3,2);
M(6,6) = P.Mas(3,3);
end

Time dependent force vector calculation:

%%Parameters
time = [0:1:99];    %Seconds

r =4.6;             %Spar radius above taper        [m]
D = 2*r             %Spar diameter                  [m]
g = 9.81;           %Gravity constant               [m/s^2]
rho= 1021;          %Seawater density               [kg/m^3]
CoB = 20            %Distance to center of buoyancy [m]
Cm = 2;             %Inertia coefficient            [-]
Cd = 1.2;           %Drag coefficient               [-]

%Wave force
omega = 0.15;            %wave period in rad/s
a = 6;                   %wave amplitude in m
d = 1500;                %Water depth, deep water conditions
wave = a*sin(omega*time) %Wave profile
k = omega^2/g            %Wave number
z = [0:-1:-1000];        %Depth coordinate

u_x = 0                  %Initial wave velocity
a_x = 0                  %Initial wave acceleration


%Here the wave forces are calculated from the wave velocity and
%acceleration by means of the Morison's equation. The Froude-krylov force
%is calculated as function of the wetted area and surface elevation. The
%induced moment is simplified as the total force per time multiplied by
%arm CoB
for i = 1:1:100
    for j = 1:1:100
        u_x(i,j) = omega*a*exp(k*z(j))*sin(omega*time(i));
        a_x(i,j) = omega^2*a*exp(k*z(j))*cos(omega*time(i));
        F(i,j) = rho*Cm*pi/4*D^2*a_x(i,j)+Cd*0.5*rho*D*u_x(i,j)^2;
        F_krylov(i) = rho*g*(pi*D^2)/4*wave(i);
        
    end
    F_morison(i) = sum(F(i,:))
    Moment(i) = CoB*F_morison(i)
end

Answer 1

[Alan Stevens]

Seems to work ok the following way

r =4.6;             %Spar radius above taper        [m]
D = 2*r;             %Spar diameter                  [m]
g = 9.81;           %Gravity constant               [m/s^2]
rho= 1021;          %Seawater density               [kg/m^3]
CoB = 20;            %Distance to center of buoyancy [m]
Cm = 2;             %Inertia coefficient            [-]
Cd = 1.2;           %Drag coefficient               [-]
%Wave force
omega = 0.15;            %wave period in rad/s
a = 6;                   %wave amplitude in m
d = 1500;                %Water depth, deep water conditions

tspan = 0:99;

wave = a*sin(omega*tspan); %Wave profile
k = omega^2/g ;           %Wave number
z = [0:-1:-1000];        %Depth coordinate
u_x = 0 ;                 %Initial wave velocity
a_x = 0 ;                 %Initial wave acceleration
%Here the wave forces are calculated from the wave velocity and
%acceleration by means of the Morison's equation. The Froude-krylov force
%is calculated as function of the wetted area and surface elevation. The
%induced moment is simplified as the total force per time multiplied by
%arm CoB
for i = 1:1:100
    for j = 1:1:100
        u_x(i,j) = omega*a*exp(k*z(j))*sin(omega*tspan(i));
        a_x(i,j) = omega^2*a*exp(k*z(j))*cos(omega*tspan(i));
        F(i,j) = rho*Cm*pi/4*D^2*a_x(i,j)+Cd*0.5*rho*D*u_x(i,j)^2;
        F_krylov(i) = rho*g*(pi*D^2)/4*wave(i);
        
    end
    F_morison(i) = sum(F(i,:));
    Moment(i) = CoB*F_morison(i);
end

X0 = [0;0;1;0;0;0]; % [x10; x20; x30; dx1dt0; dx2dt0; dx3dt0]
[t, X] = ode45(@(t,X) odefn(t,X,tspan,F_krylov,F_morison,Moment), tspan, X0);
subplot(3,1,1)
plot(t,X(:,1)),grid
xlabel('t'),ylabel('surge')
subplot(3,1,2)
plot(t,X(:,2)),grid
xlabel('t'),ylabel('heave')
subplot(3,1,3)
plot(t,X(:,3)),grid
xlabel('t'),ylabel('pitch')


function dXdt = odefn(t, X,tspan,F_krylov,F_morison,Moment)

        x = X(1:3);
        xdot = X(4:6);
        M = [8.07e6 0 -6.29e8; 0 8.07e6 1.12e5; -6.29e8 1.12e5 6.80e10];
        A = [7.98e6 0 -4.94e8; 0 2.23e4 0; -4.94e8 0 3.97e10];
        B = [1e5 0 0; 0 1.3e5 0; 0 0 0];
        C = [0 0 0; 0 3.34e5 0;0 0 0];
        
        F_V = interp1(tspan,F_krylov,t);
        F_x = interp1(tspan,F_morison,t);
        M_t = interp1(tspan,Moment,t);

        F = [F_V; F_x; M_t];
        xddot = (M+A)\(F - C*x - B*xdot);
        
        dXdt = [xdot; xddot];
       
end

[Tijmen]

Thanks Alan! You got the EOM working!