```
function fo = gaCAES1(x)
% Caution: 杈撳叆鐨勫帇鍔涙槸kPa锛屽昂瀵告槸mm
%% 绯荤粺鎬诲拰
%鍒濆鍖栧弬鏁?
ambient_pressure=101.32;
ambient_temperature=293.15;
P_start=x(6);
P_end=7500;
compression_ratio=(P_end/ambient_pressure)^(1/3);
[h_ambient s_ambient]=refpropm('HS','T',ambient_temperature,'P',ambient_pressure,'air.mix'); %#ok<*NCOMMA>
delta_m=370000;
% vol=4527.23;% m^3
D_end=refpropm('D','T',ambient_temperature,'P',P_end,'air.mix');
% m_end=D_end*vol;
D_start=refpropm('D','T',ambient_temperature,'P',P_start,'air.mix');
% m_start=D_start*vol;
vol=delta_m/(D_end-D_start);
%% 涓?绾у帇缂╂満 COMP1
%浣庡帇鍘嬬缉鏈?
adiabatic_efficiency = 0.85; % 缁濈儹鏁堢巼
% compression_ratio = 7.87; % 鍘嬬缉姣?(宸茶鎬讳綋鎺у埗)
m_rate=20; % 鍏呰兘娈电殑娴侀噺
s_low_comp_in=refpropm('S','T',ambient_temperature,'P',ambient_pressure,'air.mix');
inlet_entropy=s_low_comp_in;% 杈撳叆鍏ュ彛鐔靛?硷紝鍗曚綅锛欽/(kg*K)
P_low_comp_in=ambient_pressure;% 杈撳叆鍏ュ彛鍘嬪己锛屽崟浣嶏細kPa
inlet_pressure=P_low_comp_in;
outlet_pressure=P_low_comp_in*compression_ratio;
% 璁$畻鍑哄彛鐒撳??
h2_adiabatic = refpropm('H','P',outlet_pressure,'S',inlet_entropy,'air.mix'); % 缁濈儹杩囩▼涓嬬殑鍑哄彛鐒撳??
h1 = refpropm('H','P',inlet_pressure,'S',inlet_entropy,'air.mix'); % 鍏ュ彛鐒撳??
h2_actual = h1 + 1/adiabatic_efficiency*(h2_adiabatic - h1); % 瀹為檯鍑哄彛鐒撳?硷紝鑰冭檻缁濈儹鏁堢巼
% 璁$畻鍘嬬缉鏈哄姛鑰?
h_low_comp_work=h2_actual - h1;
low_comp_work = h_low_comp_work*m_rate; % 鍘嬬缉鏈烘墍鍋氱殑鍔熺巼
% 璁$畻鍑哄彛娓╁害
low_comp_Temp_out=refpropm('T','H',h2_actual,'P',outlet_pressure,'air.mix');
s_low_comp_out=refpropm('S','T',low_comp_Temp_out,'P',outlet_pressure,'air.mix');
% fprintf('浣庡帇鍘嬬缉鏈虹殑杈撳叆鍔熶负%.0fkW\n',low_comp_work/1000);
ex_low_comp_in=(h1-h_ambient)-ambient_temperature*(s_low_comp_in-s_ambient);
ex_low_comp_out=(h2_actual-h_ambient)-ambient_temperature*(s_low_comp_out-s_ambient);
ix_low_comp=ex_low_comp_in+low_comp_work/m_rate-ex_low_comp_out;
%% 浣庡帇鎹㈢儹鍣?
mc_low = 0.4*m_rate; % 鍐锋祦浣撹川閲忔祦閲忥紝鍗曚綅 kg/s
T_hi = low_comp_Temp_out; % 鐑祦浣撳叆鍙f俯搴︼紝鍗曚綅 K
cp_h = refpropm('C','T',T_hi,'P',outlet_pressure,'air.mix'); % 鐑祦浣撳畾鍘嬫瘮鐑锛屽崟浣? J/(kg.K)
while 1
%浣庡帇鎹㈢儹鍣?
cond=0; % 鍒ゆ柇椤洪?嗘祦
cp_c = 2100; % 鍐锋祦浣撳畾鍘嬫瘮鐑锛屽崟浣? J/(kg.K)
T_ci = 293.15; % 鍐锋祦浣撳叆鍙f俯搴︼紝鍗曚綅 K
mh = m_rate; % 鐑祦浣撹川閲忔祦閲忥紝鍗曚綅 kg/s
pp=5; % 鑺傜偣娓╁樊
P=outlet_pressure;
omega=0; % 婕忔皵鐜?
% 瑙f柟绋嬪緱鍒板喎娴佷綋鍑哄彛娓╁害
syms T_co T_ho %#ok<*NASGU> % 瀹氫箟绗﹀彿鍙橀噺 T_co,T_ho
if cond==0
T_ho=T_ci+pp; %閫嗘祦
eqn=mc_low*cp_c*(T_ci-T_co)-mh*cp_h*(T_ho-T_hi)== 0;
T_co_sol=solve(eqn,T_co);
T_co_low=double(T_co_sol);
P_out=P*(1-omega);
else
T_ho=T_co+pp; %椤烘祦
H_s=refpropm('H','T',T_ho,'P',P,'air.mix');
H_a=refpropm('H','T',T_ho,'P',P*(1-omega),'air.mix');
de_H=H_s-H_a;
eqn = mc_low * cp_c * (T_ci - T_co) - mh * cp_h * (T_ho - T_hi) + de_H == 0; % 鐑钩琛℃柟绋?
T_co_sol = solve(eqn, T_co); % 姹傝В鏂圭▼
T_co_low = double(T_co_sol); % 灏嗚В杞崲涓烘暟鍊?
P_out=P*(1-omega);
end
if T_co_low<=T_hi-15
break
else
mc_low=mc_low+0.05;
end
end
Q_ch1=mh*cp_h*(T_hi-T_ho);
T_mean_ch1=(T_hi-T_ho+T_ci-T_co_low)/log((T_hi-T_co_low)/(T_ho-T_ci));
A_ch1=Q_ch1/0.8/1100/T_mean_ch1;
h_air_charging_low_in=refpropm('H','T',T_hi,'P',P,'air.mix');
h_air_charging_low_out=refpropm('H','T',T_ho,'P',P_out,'air.mix');
s_air_charging_low_in=refpropm('S','T',T_hi,'P',P,'air.mix');
s_air_charging_low_out=refpropm('S','T',T_ho,'P',P_out,'air.mix');
% fprintf('浣庡帇鎹㈢儹鍣ㄧ┖姘旂殑鍑哄彛娓╁害涓?%fK\n',T_ho);```显示完整的修改后的代码
最新发布