%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% EXAMPLE 12-2: Transient Pressure Drop Model
% UnitSystem SI
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Inputs
function [t,P] = flow
%
% FUELCELL Fuel Cell Stack pressure and velocity model
% Best viewed with a monospaced font with 4 char tabs.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Constants
const.N = 6; % Number of layers
const.tfinal = 50; % Simulation time (s)
const.R = 8.314472; % Ideal gas constant (J/K-mol)
const.v_H2_in = 1.7e-8; % Volumetric fl ow rate of wet hydrogen (m^3/s)
const.v_air_in = 1.e-8; % Volumetric fl ow rate of air (m^3/s)
const.P_tot = 2; % Total pressure (bar)
const.Tf_in = 353.2; % Initial fl uid temperature (K)
const.Tf_air = 273.5; % Initial air temperature (K)
const.mu_H2 = 8.6e-6; % Viscosity of wet hydrogen (Pa-s)
const.mu_air = 8.6e-6; % Viscosity of air (Pa-s)
const.be = 0; % # of bends in channel
const.rho = 0.08988; % Density of hydrogen (kg/m^3)
%const.damp = 0.6; % ODE solver damping factor (to avoid ringing of thesolution)
%Convert volumetric fl ow rate to molar fl ow rate using idealgas law
const.n_air_in = const.v_air_in * (const.P_tot./const.Tf_air) * (1/0.0831) * 1000
% mol/s
% Convert volumetric fl ow rate to molar fl ow rate using idealgas law
const.n_H2_in = const.v_H2_in * (const.P_tot./const.Tf_in) * (1/0.0831) * 1000
% mol/s
% Layers
% 1 − Left end plate
% 2 − Gasket
% 3 − Contact (Copper)
% 4 − Contact
% 5 − Gasket
% 6 − End plate
% Parameters defi ned at the layer boundaries
% 1 2 3 4 5 6 7
% Gas temperature
param.T_f = [353.2, 353.2 353.2 353.2, 353.2, 353.2, 353.2]
% Parmeters defi ned at the layer centers
% 1 2 3 4 5 6
% Area (m^2)
param.A = [0.0367, 0.0367, 0.0367, 0.036 7, 0.0367, 0.0367]
% Channel Area (m^2)
param.Ach = [1.5708e-006, 1.5708e-006, 1.5708e-006, 1.5708e-006, 1.5708e006, 1.5708e-006]
% Thickness (m)
param.thick = [0.025, 0.025, 0.002, 0.002, 0.025, 0.025]
% Number of slices within layer
param.M = [1, 1, 1, 1, 1, 1];
% Channel radius
param.r = [0.000 625, 0.000625, 0.000625, 0.000625, 0.000625, 0.000625]
% Channel width (rectangular channels)
param.wc = [0.000625, 0.000625, 0.000625, 0.000625, 0.000625, 0.000625]
% Channel depth (rectangular channels)
param.dc = [0.000625, 0.000625, 0.000625, 0.000625, 0.000625, 0.000625]
% Channel length
param.L = [0.000625, 0.000625, 0.000625, 0.000625, 0.000625, 0.000625]
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Grid defi nition
% x − interslice coordinates
% n − molar fl ow rates at slice boundary
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Grid up mass fl ows. Assume each layer abuts the
% next one. The mass fl ow rate is at the boundary of each slice. x is at the
% edge of each slice (like a stair plot).
x = 0;
layer = [];
for i = 1:const.N
x = [x, x(end) + (1:param.M(i)) * param.thick(i)/param.M(i)] %Boundary points
layer = [layer, i * ones(1,param.M(i))]
end
% Slice thicknesses
dx = diff(x) %gives approximate derivatives between x’s
% Initial pressure (defi ned at layer boundaries)
P = zeros(size(x)) % Pressure
u_m = zeros(size(x)) % Velocity
% Intial pressure equals the outside pressure
P(1) = const.P_tot % mol/s
P(end) = const.P_tot % mol/s
% Mean velocity in a channel with area A and molar fl ow rateof MOL. P is the pressure and T is the temperature
% Use ideal gas law to calculate initial velocity of hydrogen coming in
u_m(1) = (const.n_H2_in ./ (const.P_tot ./ const.Tf_in * (1/0.0831) * 1000))./ param.A(1) % m^3/s
u_m(end) = (const.n_air_in ./ (const.P_tot ./ const.Tf_air * (1/0.0831) * 1000))./param.A(end) % m^3/s
f = @(t,P) pressure(t,P,u_m,x,layer,param,const,dx')
options = odeset('OutputFcn',@(t,P,opt) pressureplot(t,P,opt,x))
[t,P] = ode45(f, [linspace(0,const.tfinal,100)], P, options)
end % of function
%
function dPdt = pressure(t,P,u_m,x,layer,param,const,dx)
% Pressure calculations for fuel cell
% Make a convenient place to set a breakpoint
if (t > 30)
s = 1
end
% Preallocate output
dPdt = zeros(size(P))
% Fluid fl ows from left to right for layers 1-3 and
% right to left for layers 4-6
inlet = [find(layer < 4) fi nd(layer > 3)+1]
outlet = [find(layer < 4)+1 fi nd(layer > 3)]
% Treat P as a row vector
P = P(:)';
% Calculate outlet pressure drop in fl ow channel from each layer
% Calculate velocity (m/s) % Calculation is on the fl ows into and out of each
% block
u_m(outlet) = 6 * u_m(inlet) * (x(outlet)/param.dc(layer)-(x(outlet)/param.dc(layer))^2)
% velocity (m/s)
%Hydraulic diameter % Calculation is dependent upon the number of layers
Dh_outlet = (4 * param.Ach(layer))./(2 * pi * param.r(layer)); %circular fl ow channel(m)
%Reynold’s number at the channel exit
Re_outlet = (const.rho .* u_m(outlet) .* Dh_outlet)./ const.mu_H2
%Friction Factor
f_outlet = 64./Re_outlet; %for circular channels
%Pressure Drop
Kl_outlet = const.be * 30 * f_outlet(layer)
P_outlet = (f_outlet.* param.thick(layer)./ Dh_outlet).* const.rho.* ((u_m(outlet).^2)./2)+(Kl_outlet.* const.rho.* ((u_m(outlet).^2)./ 2)) %bar
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Combine into rate of change of P
%
% Do this in a loop since more than one layer outlet may contribute to
% a given element of dPdt.
for i = 1:length(outlet)
dPdt(outlet(i)) = dPdt(outlet(i)) + P_outlet(i) - P(outlet(i))
end
% Use damping factor to help numerical convergence
%dPdt = const.damp ∗ dPdt;
end % of function
%
function status = pressureplot(t,P,opt,x)
if isempty(opt)
plot(t,P), title(['t = ',num2str(t)]), hold on
status = 0;
drawnow
end
end % of function
이렇게 되어 있는 코드인데요 자꾸"function [t,P] = flow
↑
오류: 이 컨텍스트에서는 함수 정의가 지원되지 않습니다. 함수는 코드 파일에서 로컬 함수 또는 중첩 함수로만 생성할 수 있습니다."
이런 답이 오는데 어떻게 커쳐야 결과를 얻을 수 있을까요?
