problem at 50Hz!

3 次查看（过去 30 天）
Paul Rogers
Paul Rogers 2020-11-2
评论： Cris LaPierre 2020-11-4
Hello everyone,
I have the 3 files in attached, first run mappa_jerzak to load the parameters, then main_jerzak to call the function.
In the function file greitzer_Jerzak I set up the frequency I want to investigate, line 10:
hertz = 0.01;
and it all goes fine. The prooblem is for higher frequencies, like 20-30-50hz becausee the time step of the ode (45 oor 113)
To overcamee this problem I was told to call the function by setting the time step in the call, like this:
tspan = linspace(0, 100 , 60000);
[t,y]=ode113(@greitzer_Jerzak,[tspan],[0,0]);
Unfortunatley it's not working and I still can't read the 50Hz in output.
At this point i have no idea on how too go on and solve this problem.
  4 个评论
显示 2更早的评论
Mathieu NOE
Mathieu NOE 2020-11-2
hello
maybe you should have fixed time steps , with sampling freq proportionnal to your input frequency (so number of samples per period is constant)
Paul Rogers
Paul Rogers 2020-11-2
yeah, but the ode in matlab have their own adaptive step, so I can't change much.
I tough that by doing this:
tspan = linspace(0, 100 , 60000);
[t,y]=ode113(@greitzer_Jerzak,[tspan],[0,0]);
I would increase the time step, but all I am doing is to evaluete the function on more points, while the integration points remain the same.
Do you know any tricks?
I don't know anything else but ode to solve a differential system.

请先登录，再进行评论。

请先登录，再回答此问题。

采纳的回答

Paul Rogers
Paul Rogers 2020-11-3
this is the solution I was looking for:
init=[0 0]';
options= odeset('MaxStep',0.001); %maximum time-step size
[t,y]=ode45(@greitzer_Jerzak,[0,20],init,options);
so I can make sure the time step the ode chose is small enough to see the wawe I expect.

更多回答（1 个）

Cris LaPierre
Cris LaPierre 2020-11-2
编辑：Cris LaPierre 2020-11-2
Ran in:
As was described a couple times in your previous post, changing tspan does not change the time step interval MATLAB uses to solve the ODE. MATLAB still solves it the same way under the hood, and then interpolates the results to the time points defined in tspan. If the underlying sovler is not getting the correct reseult, passing in a finer and finer tspan is not going to fix that. You need to first choose the correct solver.
I have no idea what the results should look like, but when I use non-stiff solvers (ode34, ode113), The plots all look the same when hertz > 5 (i didn't exhaustively check for where the cutoff point is). When I switched to a stiff solver, I get results that appear to change with Hz.
I've combined everything into a single script below. I did change you equation for gammaT by removing (v), since v is not used in your equations. This meant updating your use of gammaT in defining y(2). I also just define tspan as [0 100] for the reason stated previously.
% Define constants
U = 68; %m/s
Vp = 0.1; %m^3
Lc = 0.41; %m
H = 0.18;
a0 = 340; %m/s
Ac = 0.0038; %m^2
W = 0.25;
gamma_T = 0.61;
psi_c0 = 0.352;
m = [-0.2:0.0001:1];
indice_m0 = find(m==0);
wh = a0*(Ac/Vp*Lc)^0.5;
B = U/(2*wh*Lc);
p01 = 1e5; %Compressor inlet pressure(Pa)
ro1 = 1.15; %Gas density(kg/m^3)
%% Adimensional Parameters
phi = m./(ro1.*U.*Ac);
phi0 = find(phi==0);%trovo l'indice di phi dove la massa è 0
pc = psi_c0+H.*(1+1.5.*((m./W)-1)-0.5.*((m./W)-1).^3);
psi_t = (1/gamma_T.^2).*m(indice_m0:end).^2;
pc_adim = psi_c0+H.*(1+1.5.*((phi./W)-1)-0.5.*((phi./W)-1).^3);
psi_t_adim = (1/gamma_T.^2).*phi(indice_m0:end).^2;
psi = psi_t./(0.5.*ro1.*U.^2);
PHI_T = gamma_T.*(psi).^0.5;
PSI_c = pc./(0.5.*ro1.*U.^2);
%% Equilibrium
a=pc_adim(indice_m0:end); %prendo i valori della caratteristica del compressore solo per m>0
b=psi_t_adim; %prendo i valori della caratteristica della Valvola a farfalla solo per m>0
c = abs(bsxfun(@minus, a, b)); %faccio a-b e ne faccio il valore assoluto della differenza per trovare la disttanza minima
[d,e]=min(c,[],2); %trovo il minimo di c e l'indice corrisponendente
f=e+indice_m0; %riaggiungo indice_m_0 così da trovarmi allineato con gli indici della massa
phi_0=phi(f); %estrae solo i valori di m corrispondenti
psi_0= pc_adim(f); %estrae solo i valori di m corrispondenti
save parametri_Jerzak_main
% Solve ODE using ode15s (stiff solver)
tspan = [0 100];
[t,y]=ode15s(@greitzer_Jerzak,[tspan],[0,0]);
% Plot results
m=[-0.2:0.01:0.8];
figure
subplot(2,2,1)
plot(t,y(:,1))
xlabel('t')
ylabel('\Phi')
grid on
subplot(2,2,2) %plottare la mappa del compressore in questo quadrante
plot(y(:,1),y(:,2))
xlabel('\phi_c')
ylabel('\psi')
grid on
hold on
%plot(y(:,1),((psi_c0+H*(1+(3/2)*((y(:,1)/W)-1)-(1/2)*(((y(:,1)/W)-1).^3)))))
plot(m,((psi_c0+H*(1+(3/2)*((m/W)-1)-(1/2)*(((m/W)-1).^3)))))
subplot(2,2,3)
plot(t,y(:,2))
xlabel('t')
ylabel('\Psi')
grid on
axis([0 100 0 1]);
●
% Define differential equations
function [ dy ] = greitzer_Jerzak( t,y )
load parametri_Jerzak_main.mat
gammaT_max = 0.9;
gammaT_min = 0.61;
A = (gammaT_max - gammaT_min)/2; %amplitude
b = gammaT_max - ((gammaT_max - gammaT_min)/2);
hertz = 50;
w = hertz*2*pi;
gamma_T = A*sin(w*t)+b;
dy = [ B*((psi_c0+H*(1+(3/2)*((y(1)/W)-1)-(1/2)*(((y(1)/W)-1).^3)))-y(2));
(1/B.*(y(1)-gamma_T.*(y(2).^0.5))) ];
end
Let me again state I have no idea if this is correct or not.

类别

在 Help Center 和 File Exchange 中查找有关 Ordinary Differential Equations 的更多信息

标签

产品


版本

R2014b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by