Using the force output calculated in one script as an input to mass/spring damper

Hi, I'm very much new to Matlab so would appreciate some guidance with a problem I'm having.

Background: I'm trying to build a very simple model of an offshore wind turbine using a spring and damper model. I have a script whose output is a force that varies with time. I also have a function called 'myspring' that solves the position and velocity of a mass on a spring, considering damping and driving force.

I want to use the force I have calculated using my script as the driving force input to 'myspring', but I can't figure out how I can make it work.

Here is what I have. Firstly the 'Morison' script to calculate the force. The content is not important - it could be any force.

% Morison equation force on body

% Define the variables

H = 3.69;           % Wave height (m)
A = H/2;            % Wave amplitude (m)
Tw = 9.87;           % Wave period (s)
omega = (2*pi)/Tw;  % Angular frequency (rad/s)
lambda = 128.02;    % Wavelength (m)
k = (2*pi)/lambda;  % Wavenumber (1/m)
dw = 25;             % Water depth (m)
Cm = 2;             % Mass coefficient (-)
Cd = 0.7;           % Drag coefficient (-)
rho = 1025;         % Density of water (kg/m^3)
D = 3;              % Diameter of structure (m)

x = 0;              % Fix the position at x = 0 for simplicity
n = 100;            % Number of time steps

t = linspace(0,6*pi,n);   % Define the vector t with n time steps
eta = A*cos(k*x - omega*t); % Create the surface elevation vector with size equal to t

F_M = zeros(1,n); % Initiate an inertia force vector with same dimension as t
F_D = zeros(1,n); % Initiate a drag force vector with same dimension as t
F = zeros(1,n); % Initiate a total force vector with same dimension as t

fun_inertia = @(z)cosh(k*(z+dw));  % Define the inertia function to be integrated
fun_drag = @(z)(cosh(k*(z+dw)).*abs(cosh(k*(z+dw))));   % Define the drag function to be integrated

for i = 1:n
    F_D(i) = abs(((H*pi)/Tw) * (1/sinh(k*dw)) * cos(k*x - omega*t(i))) * ((H*pi)/Tw) * (1/sinh(k*dw)) * cos(k*x - omega*t(i)) * integral(fun_drag,-dw,eta(i));
    F_M(i) = (Cm*rho*pi*(D^2)/4) * ((2*H*pi^2)/(Tw^2)) * (1/(sin(k*dw))) * sin(k*x - omega*t(i)) * integral(fun_inertia,-dw,eta(i));
    F(i) = F_D(i) + F_M(i);
end

and now 'myspring.m':

function pdot = myspring(t,p,c,k,m)

w = sqrt(k/m);

g = sin(t);  % This is the forcing function

pdot = zeros(size(p));
pdot(1) = p(2);
pdot(2) = g - c*p(2) - (w^2)*p(1);

end

where k = spring stiffness, c = damping, m = mass. These values are specified before running ode45 with 'myspring'.

So in summary, is there a way I can replace my 'g = sin(t)' driving function which I have written into 'myspring' with the force calculated in my 'Morison' script? Any support is appreciated.

Answer 1

[Thorsten]

Turn your script into a function by adding these lines on top of your script

function F = morison(t)
if nargin == 0 % define default argument
  t = linspace(0,6*pi,n);   % Define the vector t with n time steps
end

and remove from the body of the function

t = linspace(0,6*pi,n);

You can then use the function "morison" as follows:

     g = morison(t);

Just like a build-in function like sin.

[Richard Gilchrist]

Thanks for the reply. Could you please elaborate? I don't understand why this would allow me to use the output of Morison as the input forcing function for 'myspring'.

[Richard Gilchrist]

Hi Thorsten,

My lack of understanding is more fundamental. The Morison thing will come later, but what I'm trying to do is make the function 'myspring' accept any arbitrary driving force that I specify when running ode45. I've managed to make it accept k, c and m values, but when I try to give it a driving force it gives me this

In an assignment  A(:) = B, the number of elements in A and B must be the same.

Error in myspring (line 9)
pdot(2) = g - c*p(2) - (w^2)*p(1);

Error in @(t,p)myspring(t,p,c,k,m,g)

Error in odearguments (line 87)
f0 = feval(ode,t0,y0,args{:});   % ODE15I sets args{1} to yp0.

Error in ode45 (line 115)
    odearguments(FcnHandlesUsed, solver_name, ode, tspan, y0, options, varargin);

If I give it a 'g' that doesn't depend on time (for troubleshooting purposes) it works e.g. g=2 allows it to run.

To get this error message, I had created a variable g = sin(t), which I intended to use as the driving force.

>> k = 2; m = 5; c = 1.5; t = linspace(0,10,100);
>> x0 = [1 0];
>> g = sin(t);
>> [t,x] = ode45(@(t,p)myspring(t,p,c,k,m,g),t,x0);

Here's what 'myspring' looks like now:

function pdot = myspring(t,p,c,k,m,g)

%g = sin(t); % Driving force

w = sqrt(k/m);

pdot = zeros(size(p));
pdot(1) = p(2);
pdot(2) = g - c*p(2) - (w^2)*p(1);

end

I realise this is probably really simple, but I can't figure out what I'm doing wrong.

[Thorsten]

When g is a vector, g - c*p(2) - (w^2)*p(1); is a vector. You cannot assign a vector to a single variable

pdot(2)

[Richard Gilchrist]

Thanks for your help. I have it working now.

Answer 2

[Ondrej Piroh]

I would like to ask if this equation is applicable for deep waves and if the ratio between D / lamba> 0.2 ?