# problem at 50Hz!

3 ビュー (過去 30 日間)
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 件のコメント表示 3 件の古いコメント非表示 3 件の古いコメント
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 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 2020 年 11 月 4 日
Nice work!

サインインしてコメントする。

### その他の回答 (1 件)

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)
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;
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 )
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 ExchangeOrdinary Differential Equations についてさらに検索

R2014b

### Community Treasure Hunt

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

Start Hunting!

Translated by