mass spring inpulse ode45
1 次查看(过去 30 天)
显示 更早的评论
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.
0 个评论
回答(1 个)
Jan
2021-4-3
编辑:Jan
2021-4-3
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.
5 个评论
Jan
2021-4-5
"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.
另请参阅
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!