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);
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 件のコメント
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);
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
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 月 4 日
Nice work!


その他の回答 (1 件)

Cris LaPierre
Cris LaPierre 2020 年 11 月 2 日
編集済み: Cris LaPierre 2020 年 11 月 2 日
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];
% Plot results
grid on
subplot(2,2,2) %plottare la mappa del compressore in questo quadrante
grid on
hold on
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))) ];
Let me again state I have no idea if this is correct or not.


Help Center および File ExchangeOrdinary Differential Equations についてさらに検索




Community Treasure Hunt

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

Start Hunting!

Translated by