problem at 50Hz!

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.

[Stephan]

What about parametri_Jerzak_main.mat file?

[Paul Rogers]

added

[Mathieu NOE]

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]

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.

Answer 1

[Paul Rogers]

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.

[Cris LaPierre]

Nice work!

Answer 2

[Cris LaPierre]

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.