ODE45 Multiple Degrees of Freedom Problem
6 vues (au cours des 30 derniers jours)
Afficher commentaires plus anciens
Hello,
I am currently trying to model the displacement of a 3DOF system using matlabs ODE45 solver. However as much I try to fix the errors I am getting I do not seem to be making any progress. At the moment I am just trying to get the solver to run without error and will then plot this at a later date. Any advice / help would be greatly appreciated.
These are the errors that I am getting.
Thanks in advance.
2 commentaires
Stephen23
le 10 Nov 2022
Original question by Maximilian retrieved from Google Cache:
"ODE45 Multiple Degrees of Freedom Problem"
Hello,
I am currently trying to model the displacement of a 3DOF system using matlabs ODE45 solver. However as much I try to fix the errors I am getting I do not seem to be making any progress. At the moment I am just trying to get the solver to run without error and will then plot this at a later date. Any advice / help would be greatly appreciated.
% This script produces a transient response of a 3DOF dynamic model of a
% jacket platform
clear all
clc
% Define global variables which will be used in the other functions
% global M k F I
%% Input Parameters
E=204e+9; % Youngs modulus (GN/m^2)
R=1.25; % Outer radius of steel jacket (m)
r=1.23; % Inner radius of steel jacket (m)
M=9000e+3; % Deck mass (kg)
I=1.35e+9; % Inertia (kgm^2)
a= pi*(R^2-r^2); % Cross-sectional area (m^2)
I2nd= (pi/4)*(R^4-r^4); % Second moment of area (m^4)
l= [16 14 17 13] ; % Length to centre of mass (m)
L = 20 ; % Leg Length (m)
Nrun = length(L);
omega = zeros(3, Nrun); % Natural Frequencies
%% Defining forcing
% Wave specs : Height=2m ; Period=4s
fxa= 0.00;
fxb= 0.00;
fxc= 0.00;
fxd= 0.00;
fya= 371.56;
fyb= 378.96;
fyc= 371.56;
fyd= 351.44;
K=(3*E*I2nd) / (L^3) % Jacket stiffness (N/m)
% Define the mass, stiffness, and damping matirces m, k and c
% and amplitude of force F. We are using proportional damping
m=[ M 0 0 ; 0 M 0 ; 0 0 I ];
k= [ 4.*K 0 (2.*K*(l(1)-l(2))) ; 0 4.*K (2.*K*(l(4)-l(3))) ; 2.*K*(l(1)-l(2)) 2.*K*(l(4)-l(3)) 2.*K*(l(1)^2+l(2)^2+l(3)^2+l(4)^2) ];
c= 0.2*k;
F=[ fxa+fxb+fxc+fxd ; fya+fyb+fyc+fyd ; (fxa+fxb)*l(1)-(fxc+fxd)*l(2)-(fya+fyd)*l(3)+(fyb+fyc)*l(4) ];
% ODE Solution
% Initial Conditions
tspan = [0 1000] ;
z0 = [0 0 0 0 0 0] ;
% Solution
[t, z] = ode45(@(t, z) doffunc(t, z, k, M, I), tspan, z0) ;
%% ODE Function
function zdot=doffunc(t, z, k, M, I)
zdot = zeros(6, 1) ;
zdot(1)= z(4);
zdot(2)= z(5);
zdot(3)= z(6);
zdot(4)= ( - 0.2*((k(1)*z(4)) - (k(3)*z(6))) - (k(1)*z(1)) - (k(3)*z(3)) + F(1))/M;
zdot(5)= ( - 0.2*((k(5)*z(5)) - (k(6)*z(6))) - (k(5)*z(2)) - (k(6)*z(3)) + F(2))/M;
zdot(6)= ( - 0.2*((k(7)*z(4)) - (k(8)*z(5)) - (k(9)*z(6))) - (k(7)*z(1)) - (k(8)*z(2)) - (k(9)*z(3)) + F(3))/I;
end
These are the errors that I am getting.
Réponse acceptée
Jan
le 9 Nov 2022
Modifié(e) : Jan
le 9 Nov 2022
The error message tells you, that Matlab fails in the function F. Unfortunately you post the other code, but not the failing function. Please edit the question and add this function also to allow to help you.
Or do you mean this vector?
F=[fxa+fxb+fxc+fxd; ...
fya+fyb+fyc+fyd; ...
(fxa+fxb)*l(1)-(fxc+fxd)*l(2)-(fya+fyd)*l(3)+(fyb+fyc)*l(4) ];
If so, include in in the list of provided parameters:
[t, z] = ode45(@(t, z) doffunc(t, z, k, M, I, F), tspan, z0) ;
...
function zdot=doffunc(t, z, k, M, I, F)
Then inside doffunc the values of the vector F is used, not the function stored in the File F.m .
Plus de réponses (1)
Torsten
le 9 Nov 2022
% This script produces a transient response of a 3DOF dynamic model of a
% jacket platform
clear all
clc
% Define global variables which will be used in the other functions
% global M k F I
%% Input Parameters
E=204e+9; % Youngs modulus (GN/m^2)
R=1.25; % Outer radius of steel jacket (m)
r=1.23; % Inner radius of steel jacket (m)
M=9000e+3; % Deck mass (kg)
I=1.35e+9; % Inertia (kgm^2)
a= pi*(R^2-r^2); % Cross-sectional area (m^2)
I2nd= (pi/4)*(R^4-r^4); % Second moment of area (m^4)
l= [16 14 17 13] ; % Length to centre of mass (m)
L = 20 ; % Leg Length (m)
Nrun = length(L);
omega = zeros(3, Nrun); % Natural Frequencies
%% Defining forcing
% Wave specs : Height=2m ; Period=4s
fxa= 0.00;
fxb= 0.00;
fxc= 0.00;
fxd= 0.00;
fya= 371.56;
fyb= 378.96;
fyc= 371.56;
fyd= 351.44;
K=(3*E*I2nd) / (L^3) % Jacket stiffness (N/m)
% Define the mass, stiffness, and damping matirces m, k and c
% and amplitude of force F. We are using proportional damping
m=[ M 0 0 ; 0 M 0 ; 0 0 I ];
k= [ 4.*K 0 (2.*K*(l(1)-l(2))) ; 0 4.*K (2.*K*(l(4)-l(3))) ; 2.*K*(l(1)-l(2)) 2.*K*(l(4)-l(3)) 2.*K*(l(1)^2+l(2)^2+l(3)^2+l(4)^2) ];
c= 0.2*k;
F=[ fxa+fxb+fxc+fxd ; fya+fyb+fyc+fyd ; (fxa+fxb)*l(1)-(fxc+fxd)*l(2)-(fya+fyd)*l(3)+(fyb+fyc)*l(4) ];
% ODE Solution
% Initial Conditions
tspan = [0 10] ;
z0 = [0 0 0 0 0 0] ;
% Solution
[t, z] = ode45(@(t, z) doffunc(t, z, k, M, I, F), tspan, z0) ;
plot(t,z)
%% ODE Function
function zdot=doffunc(t, z, k, M, I, F)
zdot = zeros(6, 1) ;
zdot(1)= z(4);
zdot(2)= z(5);
zdot(3)= z(6);
zdot(4)= ( - 0.2*((k(1)*z(4)) - (k(3)*z(6))) - (k(1)*z(1)) - (k(3)*z(3)) + F(1))/M;
zdot(5)= ( - 0.2*((k(5)*z(5)) - (k(6)*z(6))) - (k(5)*z(2)) - (k(6)*z(3)) + F(2))/M;
zdot(6)= ( - 0.2*((k(7)*z(4)) - (k(8)*z(5)) - (k(9)*z(6))) - (k(7)*z(1)) - (k(8)*z(2)) - (k(9)*z(3)) + F(3))/I;
end
Voir également
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!