problem at 50Hz!
2 visualizzazioni (ultimi 30 giorni)
Mostra commenti meno recenti
Paul Rogers
il 2 Nov 2020
Commentato: Cris LaPierre
il 4 Nov 2020
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 Commenti
Mathieu NOE
il 2 Nov 2020
hello
maybe you should have fixed time steps , with sampling freq proportionnal to your input frequency (so number of samples per period is constant)
Risposta accettata
Più risposte (1)
Cris LaPierre
il 2 Nov 2020
Modificato: Cris LaPierre
il 2 Nov 2020
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.
0 Commenti
Vedere anche
Categorie
Scopri di più su Ordinary Differential Equations in Help Center e File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!