mass spring inpulse ode45

Good evening everyone!

I am trying to model a straight blowback bolt opening (that can be seen as a simple mass and spring).

From another software I have calculated the pressure of gases respect to time and via a spline interpolation I am now trying to feed them into matlab.

for now the code I wrote is the following:

%%%units
%g	mm	ms	N	MPa	N-mm	
function dydt = apertura(t,y)

load('9x19workspace.mat')

Press = fit(tempo,Pcalcolata,'smoothingspline','normalize','on')
Pressione = @(t) Press(t)


A=(0.5*diametrointerno)^2*pi
circ=2*pi*(0.5*diametrointerno)
SL=altezza*circ
mu=0.44;
k=3.000;
M=950;

Eb=100000
Ecc=210000
ni_b=0.35
ni_cc=0.3
a=0.5*9.65
b=0.5*9.68
c=8
alfa=6.06
C=(1/(Eb*(b^2-a^2)))*[(1+ni_b)*a^2+(1-ni_b)*b^2];
B=(1/(Ecc*(c^2-b^2)))*[(1+ni_cc)*c^2+(1-ni_cc)*b^2];
AA=2*a^2/(Eb*(b^2-a^2));
C_delta=1/(b*(B+C));
C_p=AA/(B+C);
pcont= @(t) Pressione(t)*C_p % contact pressure between case and chamber

ub=@(t) (1/Eb).*[-(1+ni_b).*((a^2*b^2.*(pcont(t)-Pressione(t)))./(b*(b^2-a^2)))+(1-ni_b).*(a^2.*Pressione(t)-b.^2.*pcont(t)).*b./(b.^2-a.^2)]
ucc=@(t) (1/Ecc)*[-(1+ni_cc).*((c^2*b^2.*(-pcont(t))/(b*(c^2-b^2))))+(1-ni_cc)*(b^2.*pcont(t))*b/(c^2-b^2)];

interferenza=@(t) (ucc(t)-ub(t))
festr=@(t) pcont(t)*2*pi*altezza*b*mu %extraction force
Forza_spingente=@(t)Pressione(t).*A %bolt thrust
pl=C_delta*tand(alfa)*altezza/(C_p)

if Forza_spingente(t)>festr(t)
    
 dydt=@(t) [y(2);@(t)((-2*pi*(altezza-y(1))*pcont(t)*b*mu*heaviside(altezza-y(1)) +Forza_spingente(t)*heaviside(tcanna-t))-k*y(1))/M];
else
y(2)=0
y(1)=0

end

I have added an if/else cycle in order to compute the fact that, if the bolt thrust is bigger than the extracting force the bolt will move backwards, otherwise it will be held in place by the friction between case and chamber.

Concerning the solver i wrote:

  clear all
  close all
  clc

tin=zeros(1,1)

[t,y] = ode45(@(t,y) apertura(t,y),[0 70],tin,1e-5);


% Plot the solutions for $y_1$ and $y_2$ against |t|.
plot(t,y(:,1),'-o',t,y(:,2),'-o')



title('bolt displacement cal 9x19 Parabellum');
xlabel('Time t [ms]');
ylabel('displacement[mm]  /  velocity[mm/s]');
grid on
legend('displacement','velocity')

when I execute the code i receive several errors such as:

Error using odearguments (line 113)
Inputs must be floats, namely single or double.

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

Error in risolutore (line 10)
[t,y] = ode45(@(t,y) apertura(t,y),[0 70],tin,1e-5);

where am I getting wrong?

How can I fix this problem?

Thank you in advance for your help and advices, they will be really appreciated.

Regards.

G.

Answer 1

[Jan]

Your function to be integrated replies a function handle instead of a numerical vector in some cases:

if Forza_spingente(t)>festr(t)
   dydt=@(t) [y(2);@(t)((-2*pi*(altezza-y(1))*pcont(t)*b*mu*heaviside(altezza-y(1)) +Forza_spingente(t)*heaviside(tcanna-t))-k*y(1))/M];

This cannot work.

If the IF condition is FALSE, you create a [1 x 2] vector y, but this is not used anywhere. y and dydt are both vectors, but you provide a scalar tin as initial value. In addition a colum vector is expected.

Loading parameters in each iteration will slow down the code dramatically.

[Giulio Rossi]

thank you for your answer Jan.

How do you suggest me to write those lines of code?

I am really stuck with them.

Thanks again.

G.

[Jan]

Sorry, the code is rather confusing. Therefore I cannot reliably guess its intention. You are mixing arrays and anonymous functions too frequently.

Move this part:

load('9x19workspace.mat')
Press = fit(tempo,Pcalcolata,'smoothingspline','normalize','on')
Pressione = @(t) Press(t)
A=(0.5*diametrointerno)^2*pi
circ=2*pi*(0.5*diametrointerno)
SL=altezza*circ

out of the integration and perform it once only. Store all needed parameters in a struct and provide the details by a anonymous function:

ode45(@(t,y) apertura(t,y, Param), [0 70], tin, 1e-5);

But even this must fail due to the scalar tin.

I'd suggest to start writing an integration of a simple function, e.g. y*sin(t). If this works, expand the code to include parameters and to operate on vectors also. This allows to learn step by step how to use ODE45.

[Giulio Rossi]

The intent of the code is compute velocity and displacement of a mass spring system under a time variable force, that the code has to compute if it is bigger than the friction force (in that case the system won't move) or smaller and in that case calculate y(1) and y(2).

is corret how I have fitted the datas with the spline and function calling?

So If I have understood correctly, you are suggesting to create a script with all the variables, apart from the function code and the solver? Is that correct?

[Jan]

"is corret how I have fitted the datas with the spline and function calling?" - there is no chance for me to assess this reliably.

Your description of the purpose of the code does not clarify e.g. what these lines of code should do:

ub=@(t) (1/Eb).*[-(1+ni_b).*((a^2*b^2.*(pcont(t)-Pressione(t)))./(b*(b^2-a^2)))+(1-ni_b).*(a^2.*Pressione(t)-b.^2.*pcont(t)).*b./(b.^2-a.^2)]
ucc=@(t) (1/Ecc)*[-(1+ni_cc).*((c^2*b^2.*(-pcont(t))/(b*(c^2-b^2))))+(1-ni_cc)*(b^2.*pcont(t))*b/(c^2-b^2)];

I cannot guess, why your code calculates dydt as anonymous function and y as a vector in some cases.

I suggest to restart from scratch, because the current code contains so many severe problems, that fixing this is more work than rewriting it using clean methods. Start with learning, how an integration of a tiny system works in Matlab. For the implementation of your project move the load and the processing of the constants out of the function to be integrated. Otherwise this part is called repeatedly and this wastes a lot of ressources.

[Giulio Rossi]

these 2 lines of codes are for calculating the strain both of case and camber under pressure..and so calculating friction force.

Thanks for your advice..I'll try to restart from zero adding then small lines of code and see if it works.

I'll let you know!

Thanks again

G.