Add full SHIFT project: Matlab/Simulink models, data, planning docs, and team repos
This commit is contained in:
parent
d24b146912
commit
20b62acf6e
354 changed files with 215247 additions and 0 deletions
7
.github/README.md
vendored
Normal file
7
.github/README.md
vendored
Normal file
|
|
@ -0,0 +1,7 @@
|
||||||
|
# About This Repository
|
||||||
|
This Repository specificly designed to edit the main page profile, so main page of the organization can look nice.
|
||||||
|
|
||||||
|
> If you like to change something, please go to ~ .github/profile
|
||||||
|
|
||||||
|
Edit the parts you want to change or to add photo, then push!
|
||||||
|
- You can see the README file on ~.github/README.md
|
||||||
26
.github/profile/README.md
vendored
Normal file
26
.github/profile/README.md
vendored
Normal file
|
|
@ -0,0 +1,26 @@
|
||||||
|
# Developing an AI-driven Hydrogen Microgrid Management System
|
||||||
|
|
||||||
|
### Our Mission
|
||||||
|
SHIFT's mission is to accelerate the transition to resilient, low-carbon energy systems through AI-driven energy management. We develop intelligent control solutions for critical infrastructure that optimize renewables, storage, hydrogen technologies, and grid interaction. By combining engineering and artificial intelligence, we aim to improve reliability, reduce emissions, and increase energy efficiency where uninterrupted power matters most.
|
||||||
|
|
||||||
|
### Our Vision
|
||||||
|
SHIFT envisions a future where clean, intelligent energy systems power the world's most critical infrastructure. Hybrid-Energy technologies, coordinated by AI, create autonomous energy networks that eliminate waste, cut emissions, and guarantee reliability.
|
||||||
|
|
||||||
|
<img src="https://github.com/Team-SHIFT-Space/.github/blob/main/profile/photos/1.png" width="1000">
|
||||||
|
|
||||||
|
### About team SHIFT
|
||||||
|
|
||||||
|
SHIFT has one goal in mind: creating the next generation of intelligent, resilient, and sustainable energy systems powered by AI. Our mission is to tackle some of the most urgent challenges in today's energy landscape: How can critical infrastructure secure reliable power in an unstable grid? How can we reduce dependence on fossil fuels without compromising safety? And how can we enable this transition in a way that is both technologically innovative and economically viable? Our answer is an AI-driven energy management system that coordinates renewables, hydrogen, batteries, and grid interaction to deliver clean, uninterrupted power where it matters most. As a multidisciplinary student team, all our members are working hard to achieve our goals. Founded in 2021, we work from the TU/e Innovation Space and Neuron, the center of expertise for challenge-based learning and student entrepreneurship at TU/e.
|
||||||
|
|
||||||
|
<img src="https://github.com/Team-SHIFT-Space/.github/blob/main/profile/photos/2.png" width="1000">
|
||||||
|
|
||||||
|
### Addressing the Challenge
|
||||||
|
Hospitals face a growing energy challenge: their demand for reliable, uninterrupted power is increasing, while the electricity grid is becoming more congested, less predictable, and more dependent on weather-driven renewable production. Starting in 2025, energy producers will only receive compensation for a portion of the energy they return to the grid, due to the Dutch government phasing out the "balancing system."
|
||||||
|
<img src="https://github.com/Team-SHIFT-Space/.github/blob/main/profile/photos/3.png" width="1000">
|
||||||
|
SHIFT is developing an AI-driven hydrogen microgrid designed specifically for healthcare environments. Our intelligent energy management system forecasts hospital demand, optimizes solar and battery usage, and controls electrolyzers and fuel cells to ensure every component operates at its full potential.
|
||||||
|
|
||||||
|
### Teams
|
||||||
|
- Energy Management Cluster
|
||||||
|
- Simulations Cluster
|
||||||
|
- Business & Economics Cluster
|
||||||
|
- AI & Control Systems Cluster
|
||||||
BIN
.github/profile/photos/1.png
vendored
Normal file
BIN
.github/profile/photos/1.png
vendored
Normal file
Binary file not shown.
|
After Width: | Height: | Size: 110 KiB |
BIN
.github/profile/photos/2.png
vendored
Normal file
BIN
.github/profile/photos/2.png
vendored
Normal file
Binary file not shown.
|
After Width: | Height: | Size: 32 KiB |
BIN
.github/profile/photos/3.png
vendored
Normal file
BIN
.github/profile/photos/3.png
vendored
Normal file
Binary file not shown.
|
After Width: | Height: | Size: 75 KiB |
1
.github/profile/photos/temp.md
vendored
Normal file
1
.github/profile/photos/temp.md
vendored
Normal file
|
|
@ -0,0 +1 @@
|
||||||
|
|
||||||
18
.gitignore
vendored
Normal file
18
.gitignore
vendored
Normal file
|
|
@ -0,0 +1,18 @@
|
||||||
|
# Node / web build artifacts (website/)
|
||||||
|
node_modules/
|
||||||
|
dist/
|
||||||
|
.astro/
|
||||||
|
.wrangler/
|
||||||
|
|
||||||
|
# Claude Code local config / session data
|
||||||
|
.claude/
|
||||||
|
|
||||||
|
# OS junk
|
||||||
|
Thumbs.db
|
||||||
|
ehthumbs.db
|
||||||
|
Desktop.ini
|
||||||
|
.DS_Store
|
||||||
|
|
||||||
|
# MATLAB / Simulink autosaves
|
||||||
|
*.asv
|
||||||
|
*.autosave
|
||||||
BIN
Forecasting requirements.pdf
Normal file
BIN
Forecasting requirements.pdf
Normal file
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
File diff suppressed because it is too large
Load diff
Binary file not shown.
File diff suppressed because it is too large
Load diff
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
File diff suppressed because it is too large
Load diff
File diff suppressed because one or more lines are too long
Binary file not shown.
Binary file not shown.
BIN
Shift Matlab Drive/Shift Matlab Drive/RL_ML/Empty.mlx
Normal file
BIN
Shift Matlab Drive/Shift Matlab Drive/RL_ML/Empty.mlx
Normal file
Binary file not shown.
|
|
@ -0,0 +1,333 @@
|
||||||
|
component (Propagation = blocks) Electrolyzer
|
||||||
|
% Electrolyzer
|
||||||
|
% This block models a stack of membrane electrode assemblies (MEA) for a
|
||||||
|
% proton exchange membrane (PEM) electrolyzer. Water is consumed from the
|
||||||
|
% connected anode flow channels. Oxygen and hydrogen are produced at the
|
||||||
|
% connected anode and cathode flow channels, respectively. Waste heat is
|
||||||
|
% dissipated to the connected thermal mass.
|
||||||
|
|
||||||
|
% Copyright 2021-2024 The MathWorks, Inc.
|
||||||
|
|
||||||
|
nodes
|
||||||
|
% Electrodes
|
||||||
|
p = foundation.electrical.electrical; % +
|
||||||
|
n = foundation.electrical.electrical; % -
|
||||||
|
% Moist air ports to obtain domain parameters and property tables for
|
||||||
|
% the anode network and the cathode network only. There is no flow through ports A and C
|
||||||
|
A = foundation.moist_air.moist_air; % A
|
||||||
|
C = foundation.moist_air.moist_air; % C
|
||||||
|
% Anode side thermal liquid port to remove water from channels
|
||||||
|
H2O = foundation.thermal_liquid.thermal_liquid; % H2O
|
||||||
|
% Anode side moist air source port to add oxygen to channels
|
||||||
|
O2 = foundation.moist_air.moist_air_source; % O2
|
||||||
|
% Cathode side moist air source port to add hydrogen and water vapor to channels
|
||||||
|
H2 = foundation.moist_air.moist_air_source; % H2
|
||||||
|
% Thermal port for heat generation
|
||||||
|
H = foundation.thermal.thermal; % H
|
||||||
|
end
|
||||||
|
|
||||||
|
inputs
|
||||||
|
% Input ports for the internal states of the cathode gas channels
|
||||||
|
FC = [101325; 293.15; 0.5; 0.01; 0.01; 0.01; 0.001; 0.001; 0]; % FC
|
||||||
|
end
|
||||||
|
|
||||||
|
annotations
|
||||||
|
[p, n] : Side = top;
|
||||||
|
[A, H2O, O2] : Side = left;
|
||||||
|
[C, H2, FC] : Side = right;
|
||||||
|
H : Side = bottom;
|
||||||
|
end
|
||||||
|
|
||||||
|
parameters
|
||||||
|
N_cell = {50, '1' }; % Number of cells in stack
|
||||||
|
area_cell = {280, 'cm^2' }; % Cell area
|
||||||
|
t_membrane = {125, 'um' }; % Membrane thickness
|
||||||
|
t_gdl_A = {25, 'um' }; % Anode gas diffusion layer (GDL) thickness
|
||||||
|
t_gdl_C = {250, 'um' }; % Cathode gas diffusion layer (GDL) thickness
|
||||||
|
io = {8e-05, 'A/cm^2'}; % Exchange current density
|
||||||
|
alpha = {0.5, '1' }; % Charge transfer coefficient
|
||||||
|
D_H2O_gdl_A = {0.07, 'cm^2/s'}; % Water diffusivity in anode GDL
|
||||||
|
D_H2O_gdl_C = {0.07, 'cm^2/s'}; % Water diffusivity in cathode GDL
|
||||||
|
rho_membrane = {2000, 'kg/m^3'}; % Density of dry membrane
|
||||||
|
M_membrane = {1.1, 'kg/mol'}; % Equivalent weight of dry membrane
|
||||||
|
end
|
||||||
|
|
||||||
|
% Parameter checks
|
||||||
|
equations
|
||||||
|
assert(N_cell > 0)
|
||||||
|
assert(area_cell > 0)
|
||||||
|
assert(t_membrane > 0)
|
||||||
|
assert(t_gdl_A > 0)
|
||||||
|
assert(t_gdl_C > 0)
|
||||||
|
assert(io > 0)
|
||||||
|
assert(alpha > 0)
|
||||||
|
assert(D_H2O_gdl_A > 0)
|
||||||
|
assert(D_H2O_gdl_C > 0)
|
||||||
|
assert(rho_membrane > 0)
|
||||||
|
assert(M_membrane > 0)
|
||||||
|
end
|
||||||
|
|
||||||
|
% Internal parameters
|
||||||
|
parameters (Access = protected)
|
||||||
|
% Constants
|
||||||
|
R_u = {8.31446261815324, 'J/K/mol'}; % Universal gas constant
|
||||||
|
F = {96485.33212, 'C/mol' }; % Faraday constant
|
||||||
|
G_H2O = {-237.14, 'kJ/mol' }; % Gibbs free energy of water
|
||||||
|
HHV_H2 = {285.8, 'kJ/mol' }; % Higher heating value of hydrogen
|
||||||
|
LHV_H2 = HHV_H2 - MW_H2O * hfg_H2O_std; % Lower heating value of hydrogen
|
||||||
|
T_std = {25, 'degC' }; % Standard temperature
|
||||||
|
p_std = {1, 'atm' }; % Standard pressure
|
||||||
|
|
||||||
|
% Standard cell potential
|
||||||
|
E_cell = G_H2O/(-2*F);
|
||||||
|
|
||||||
|
% Membrane permeability to water
|
||||||
|
K_darcy = {1.58e-14, 'cm^2'};
|
||||||
|
|
||||||
|
% Compute molar masses from specific gas constants
|
||||||
|
MW_H2O = R_u/A.R_w;
|
||||||
|
MW_O2 = R_u/A.R_g;
|
||||||
|
MW_H2 = R_u/C.R_g;
|
||||||
|
|
||||||
|
% Species specific enthalpy at standard temperature
|
||||||
|
h_H2O_std = tablelookup(A.T_TLU, A.h_w_TLU, T_std, interpolation = linear, extrapolation = linear);
|
||||||
|
hfg_H2O_std = tablelookup(A.T_TLU, A.h_w_vap_TLU, T_std, interpolation = linear, extrapolation = linear);
|
||||||
|
h_O2_std = tablelookup(A.T_TLU, A.h_g_TLU, T_std, interpolation = linear, extrapolation = linear);
|
||||||
|
h_H2_std = tablelookup(C.T_TLU, C.h_g_TLU, T_std, interpolation = linear, extrapolation = linear);
|
||||||
|
end
|
||||||
|
|
||||||
|
variables (ExternalAccess = observe)
|
||||||
|
i = {0, 'A' }; % Current (positive in)
|
||||||
|
Q = {0, 'kW'}; % Heat flow rate (positive in)
|
||||||
|
|
||||||
|
a_acl = {0.5, '1'}; % Water vapor activity at anode catalyst layer (ACL)
|
||||||
|
a_ccl = {0.5, '1'}; % Water vapor activity at cathode catalyst layer (CCL)
|
||||||
|
end
|
||||||
|
|
||||||
|
variables (ExternalAccess = none)
|
||||||
|
% Through variables for thermal liquid port
|
||||||
|
mdot_H2O = {0, 'kg/s'}; % TL mass flow rate into port H2O
|
||||||
|
Phi_H2O = {0, 'kW' }; % TL energy flow rate into port H2O
|
||||||
|
end
|
||||||
|
|
||||||
|
branches
|
||||||
|
i : p.i -> n.i;
|
||||||
|
Q : H.Q -> *;
|
||||||
|
|
||||||
|
mdot_H2O : H2O.mdot -> *
|
||||||
|
Phi_H2O : H2O.Phi -> *
|
||||||
|
end
|
||||||
|
|
||||||
|
% For logging
|
||||||
|
intermediates
|
||||||
|
v = p.v - n.v; % Voltage
|
||||||
|
|
||||||
|
i_cell = if ge(i, 0), i/area_cell else 0 end; % Cell current density
|
||||||
|
v_cell = v_nernst + v_act + v_ohm; % Cell voltage
|
||||||
|
R_ohm = t_membrane / sigma; % Membrane resistance
|
||||||
|
|
||||||
|
power_elec = N_cell * v_cell * i_cell * area_cell; % Electrical power consumed
|
||||||
|
power_dissipated = power_elec - power_net; % Power dissipated
|
||||||
|
efficiency_HHV = simscape.function.limit(HHV_H2/(2*F) / v_cell, 0, 1, false); % Thermal efficiency based on HHV
|
||||||
|
efficiency_LHV = simscape.function.limit(LHV_H2/(2*F) / v_cell, 0, 1, false); % Thermal efficiency based on LHV
|
||||||
|
|
||||||
|
H2O_consumed = N_cell * MW_H2O * area_cell * i_cell/(2*F); % Water consumed
|
||||||
|
O2_produced = N_cell * MW_O2 * area_cell * i_cell/(4*F); % Oxygen produced
|
||||||
|
H2_produced = N_cell * MW_H2 * area_cell * i_cell/(2*F); % Hydrogen produced
|
||||||
|
|
||||||
|
H2O_transport = (nflux_H2O_diff + nflux_H2O_drag + nflux_H2O_hydraulic) ...
|
||||||
|
* MW_H2O * area_cell * N_cell; % Net water transport from anode to cathode
|
||||||
|
end
|
||||||
|
|
||||||
|
annotations
|
||||||
|
[v, v_cell] : LoggingUnit = 'V';
|
||||||
|
i_cell : LoggingUnit = 'A/cm^2';
|
||||||
|
R_ohm : LoggingUnit = 'Ohm*cm^2';
|
||||||
|
[power_elec, power_dissipated] : LoggingUnit = 'kW';
|
||||||
|
[efficiency_HHV, efficiency_LHV] : LoggingUnit = '1';
|
||||||
|
[H2O_consumed, O2_produced, H2_produced, H2O_transport] : LoggingUnit = 'g/s';
|
||||||
|
end
|
||||||
|
|
||||||
|
intermediates (ExternalAccess = none)
|
||||||
|
% Extract cathode gas states from input FC
|
||||||
|
% Note: Can also connect Measurement Select (MA) block to choose values
|
||||||
|
p_C = {FC(1), 'Pa'}; % Pressure
|
||||||
|
y_H2O_C = FC(5); % Water vapor mole fraction
|
||||||
|
y_H2_C = FC(8); % Hydrogen mole fraction
|
||||||
|
|
||||||
|
% Anode pressure and temperature
|
||||||
|
p_A = H2O.p;
|
||||||
|
T_A = H2O.T;
|
||||||
|
|
||||||
|
% Mole fractions at anode
|
||||||
|
% Assume saturated water vapor at stack temperature
|
||||||
|
y_H2O_A = if ge(p_ws_ratio_A, 1), 1/p_ws_ratio_A else 1 end;
|
||||||
|
y_O2_A = 1 - y_H2O_A;
|
||||||
|
|
||||||
|
% Stack temperature
|
||||||
|
T_stack = H.T;
|
||||||
|
|
||||||
|
% Ratio of pressure to water vapor saturation pressure
|
||||||
|
p_ws_ratio_A = exp(log(value(p_A, 'Pa')) - tablelookup(A.T_TLU, A.log_p_ws_TLU, T_stack, interpolation = linear, extrapolation = linear));
|
||||||
|
p_ws_ratio_C = exp(log(value(p_C, 'Pa')) - tablelookup(C.T_TLU, C.log_p_ws_TLU, T_stack, interpolation = linear, extrapolation = linear));
|
||||||
|
|
||||||
|
% Activities
|
||||||
|
a_H2O_A_ = y_H2O_A * p_ws_ratio_A;
|
||||||
|
a_O2_A_ = y_O2_A * p_A / p_std;
|
||||||
|
a_H2_C_ = y_H2_C * p_C / p_std;
|
||||||
|
a_H2O_A = if ge(a_H2O_A_, 1e-9), a_H2O_A_ else 1e-6 end;
|
||||||
|
a_O2_A = if ge(a_O2_A_, 1e-9), a_O2_A_ else 1e-6 end;
|
||||||
|
a_H2_C = if ge(a_H2_C_, 1e-9), a_H2_C_ else 1e-6 end;
|
||||||
|
|
||||||
|
% Nernst voltage
|
||||||
|
v_nernst = E_cell + R_u*T_stack/(2*F) * log((a_H2_C * a_O2_A^0.5) / a_H2O_A);
|
||||||
|
|
||||||
|
% Activation losses from Tafel equation
|
||||||
|
b = R_u * T_stack / (2 * alpha * F);
|
||||||
|
v_act = if ge(i_cell, io), b*log(i_cell/io) else 0 end;
|
||||||
|
|
||||||
|
% Water content
|
||||||
|
lambda_acl = membrane_water(a_acl);
|
||||||
|
lambda_ccl = membrane_water(a_ccl);
|
||||||
|
lambda_membrane = (lambda_acl + lambda_ccl)/2;
|
||||||
|
|
||||||
|
% Membrane conductivity
|
||||||
|
sigma_30 = {if ge(lambda_membrane, 1), 0.005139*lambda_membrane - 0.00326 else 0.005139 - 0.00326 end, '1/(Ohm*cm)'};
|
||||||
|
sigma = sigma_30 * exp(1268*(1/303.15 - 1/value(T_stack, 'K')));
|
||||||
|
|
||||||
|
% Resistive voltage loss
|
||||||
|
v_ohm = R_ohm * i_cell;
|
||||||
|
|
||||||
|
% Water diffusion coefficient across membrane
|
||||||
|
D_H2O_membrane = {1.25e-10, 'm^2/s'} * exp(2416*(1/303.15 - 1/value(T_stack, 'K')));
|
||||||
|
|
||||||
|
% Water concentrations at anode and cathode catalyst layers
|
||||||
|
Conc_H2O_acl = rho_membrane / M_membrane * lambda_acl;
|
||||||
|
Conc_H2O_ccl = rho_membrane / M_membrane * lambda_ccl;
|
||||||
|
|
||||||
|
% Molar flux of water across membrane due to diffusion
|
||||||
|
nflux_H2O_diff = D_H2O_membrane*(Conc_H2O_acl - Conc_H2O_ccl)/t_membrane;
|
||||||
|
|
||||||
|
% Water electro-osmotic drag coefficient
|
||||||
|
nd_H2O_membrane = ...
|
||||||
|
if ge(lambda_acl, 0), ...
|
||||||
|
0.0029*lambda_acl^2 + 0.05*lambda_acl - 3.4e-19 ...
|
||||||
|
else ...
|
||||||
|
0.05*lambda_acl - 3.4e-19 ...
|
||||||
|
end;
|
||||||
|
|
||||||
|
% Molar flux of water across membrane due to electro-osmotic drag
|
||||||
|
nflux_H2O_drag = nd_H2O_membrane * i_cell / F;
|
||||||
|
|
||||||
|
% Water vapor mole fraction at anode and cathode catalyst layers
|
||||||
|
y_H2O_acl = a_acl/p_ws_ratio_A;
|
||||||
|
y_H2O_ccl = a_ccl/p_ws_ratio_C
|
||||||
|
|
||||||
|
% Molar flux of water across anode and cathode diffusion layers
|
||||||
|
nflux_H2O_A = p_A*D_H2O_gdl_A/(R_u*T_stack) * (y_H2O_A - y_H2O_acl)/t_gdl_A;
|
||||||
|
nflux_H2O_C = p_C*D_H2O_gdl_C/(R_u*T_stack) * (y_H2O_ccl - y_H2O_C)/t_gdl_C;
|
||||||
|
|
||||||
|
% Molar flux of water across membrane due to hydraulic pressure difference based on Darcy's law
|
||||||
|
mu_H2O = tablelookup(A.T_TLU, A.mu_w_TLU, T_stack, interpolation = linear, extrapolation = nearest);
|
||||||
|
nflux_H2O_hydraulic = ...
|
||||||
|
if gt(p_A, p_C), ...
|
||||||
|
(p_A - p_C) * K_darcy * p_A * y_H2O_acl / (R_u * T_stack * mu_H2O * t_membrane) ...
|
||||||
|
else...
|
||||||
|
(p_A - p_C) * K_darcy * p_C * y_H2O_ccl / (R_u * T_stack * mu_H2O * t_membrane) ...
|
||||||
|
end
|
||||||
|
|
||||||
|
% Specific enthalpy of water at port H2O
|
||||||
|
% based on fluid properties of thermal liquid network
|
||||||
|
u_H2O = tablelookup(H2O.T_TLU, H2O.p_TLU, H2O.u_TLU, T_A, p_A, interpolation = linear, extrapolation = linear);
|
||||||
|
rho_H2O = tablelookup(H2O.T_TLU, H2O.p_TLU, H2O.rho_TLU, T_A, p_A, interpolation = linear, extrapolation = linear);
|
||||||
|
h_H2O = u_H2O + p_A/rho_H2O;
|
||||||
|
|
||||||
|
% Specific enthalpy of water at anode
|
||||||
|
% based on fluid properties of moist air network
|
||||||
|
h_H2O_A = tablelookup(A.T_TLU, A.h_w_TLU, T_A, interpolation = linear, extrapolation = linear);
|
||||||
|
hfg_H2O_A = tablelookup(A.T_TLU, A.h_w_vap_TLU, T_A, interpolation = linear, extrapolation = linear);
|
||||||
|
|
||||||
|
% Energy consumed by reaction at standard temperature
|
||||||
|
power_rxn = HHV_H2 * H2_produced / MW_H2;
|
||||||
|
|
||||||
|
% Energy gain in membrane due to bringing reactants and products to standard temperature
|
||||||
|
power_delta_std = H2O_consumed*(h_H2O_A - hfg_H2O_A) - H2O_consumed*(h_H2O_std - hfg_H2O_std) ...
|
||||||
|
+ (source_O2_A.Phi_S - source_O2_A.mdot_g_S*h_O2_std) ...
|
||||||
|
+ (source_H2_C.Phi_S - source_H2_C.mdot_g_S*h_H2_std);
|
||||||
|
|
||||||
|
% Energy gain in membrane due to water transport
|
||||||
|
power_trans = H2O_transport*(h_H2O_A - hfg_H2O_A) + transport_H2O_C.Phi_S;
|
||||||
|
|
||||||
|
% Net energy consumption in stack
|
||||||
|
power_net = power_rxn - power_delta_std - power_trans;
|
||||||
|
end
|
||||||
|
|
||||||
|
equations
|
||||||
|
% Stack voltage
|
||||||
|
v == N_cell * v_cell;
|
||||||
|
|
||||||
|
% Heat generated
|
||||||
|
-Q == power_dissipated;
|
||||||
|
|
||||||
|
% Equate water vapor mass flow rates at GDL and membrane
|
||||||
|
% to solve for water vapor activity at ACL and CCL
|
||||||
|
nflux_H2O_A * MW_H2O * area_cell * N_cell == H2O_transport + H2O_consumed;
|
||||||
|
nflux_H2O_C * MW_H2O * area_cell * N_cell == H2O_transport;
|
||||||
|
|
||||||
|
% Through variables for removal of water at port H2O
|
||||||
|
mdot_H2O == H2O_consumed + H2O_transport;
|
||||||
|
Phi_H2O == mdot_H2O*h_H2O;
|
||||||
|
|
||||||
|
% Assign mass flow rate to internal trace gas source and moisture source blocks
|
||||||
|
% to model mass consumption/generation due to reaction and water transport
|
||||||
|
source_O2_A.M == O2_produced;
|
||||||
|
source_H2_C.M == H2_produced;
|
||||||
|
transport_H2O_C.M == H2O_transport;
|
||||||
|
|
||||||
|
source_O2_A.T == T_stack;
|
||||||
|
source_H2_C.T == T_stack;
|
||||||
|
transport_H2O_C.T == T_stack;
|
||||||
|
end
|
||||||
|
|
||||||
|
% Internal components to add or remove mass at the anode and cathode
|
||||||
|
components (ExternalAccess = none)
|
||||||
|
% Anode trace gas source block for oxygen production
|
||||||
|
source_O2_A = foundation.moist_air.sources.moisture_trace_gas.trace_gas_source( ...
|
||||||
|
source_type = foundation.enum.constant_controlled.controlled);
|
||||||
|
% Cathode trace gas source and moisture source blocks for hydrogen production and water transport
|
||||||
|
source_H2_C = foundation.moist_air.sources.moisture_trace_gas.trace_gas_source( ...
|
||||||
|
source_type = foundation.enum.constant_controlled.controlled);
|
||||||
|
transport_H2O_C = foundation.moist_air.sources.moisture_trace_gas.moisture_source( ...
|
||||||
|
source_type = foundation.enum.constant_controlled.controlled, ...
|
||||||
|
moisture_source_spec = foundation.enum.MoistureSourceSpec.WaterVapor, ...
|
||||||
|
vapor_from_liquid = true);
|
||||||
|
end
|
||||||
|
|
||||||
|
connections
|
||||||
|
connect(O2, source_O2_A.S)
|
||||||
|
connect(H2, source_H2_C.S, transport_H2O_C.S)
|
||||||
|
end
|
||||||
|
|
||||||
|
end
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
function lambda = membrane_water(a)
|
||||||
|
% Compute the water content of the membrane as a function of the water
|
||||||
|
% activity. This function is based on fits to FEA simulations.
|
||||||
|
%
|
||||||
|
% See: Dutta, et.al., Numerical prediction of mass-exchange between
|
||||||
|
% cathode and anode channels in a PEM fuel cell. (Equation 16)
|
||||||
|
|
||||||
|
definitions
|
||||||
|
lambda = ...
|
||||||
|
if lt(a, 0), ...
|
||||||
|
0.043 + 17.81*a ...
|
||||||
|
elseif le(a, 1), ...
|
||||||
|
0.043 + 17.81*a - 39.85*a^2 + 36*a^3 ...
|
||||||
|
else ...
|
||||||
|
14.003 + 1.4*(a - 1) ...
|
||||||
|
end;
|
||||||
|
end
|
||||||
|
|
||||||
|
end
|
||||||
|
|
@ -0,0 +1,80 @@
|
||||||
|
<?xml-stylesheet href="physmodIcon.css" type="text/css"?>
|
||||||
|
<!-- Generated by Microsoft Visio, SVG Export cell_voltage6.svg Page-1 -->
|
||||||
|
<svg xmlns="http://www.w3.org/2000/svg" xmlns:xlink="http://www.w3.org/1999/xlink"
|
||||||
|
xmlns:ev="http://www.w3.org/2001/xml-events"
|
||||||
|
xmlns:v="http://schemas.microsoft.com/visio/2003/SVGExtensions/"
|
||||||
|
xmlns:d="http://www.mathworks.com/blockgraphics"
|
||||||
|
width="215" height="222"
|
||||||
|
viewBox="136 417.37 215.09 222"
|
||||||
|
xml:space="preserve" color-interpolation-filters="sRGB" class="st5">
|
||||||
|
<v:documentProperties v:langID="1033" v:metric="true" v:viewMarkup="false">
|
||||||
|
<v:userDefs>
|
||||||
|
<v:ud v:nameU="msvNoAutoConnect" v:val="VT0(1):26"/>
|
||||||
|
</v:userDefs>
|
||||||
|
</v:documentProperties>
|
||||||
|
|
||||||
|
<style type="text/css">
|
||||||
|
<![CDATA[
|
||||||
|
.st1 {fill:url(#ptrn12-4);shape-rendering:crispEdges;stroke:#000000;stroke-linecap:round;stroke-linejoin:round;stroke-width:0.24}
|
||||||
|
.st2 {fill:#a5a5a5;stroke:#000000;stroke-linecap:round;stroke-linejoin:round;stroke-width:0.24}
|
||||||
|
.st3 {fill:#bfbfbf;stroke:#000000;stroke-linecap:round;stroke-linejoin:round;stroke-width:0.24}
|
||||||
|
.st4 {fill:url(#grad0-15);stroke:#000000;stroke-linecap:round;stroke-linejoin:round;stroke-width:0.24}
|
||||||
|
.st5 {fill:none;fill-rule:evenodd;font-size:12px;overflow:visible;stroke-linecap:square;stroke-miterlimit:3}
|
||||||
|
]]>
|
||||||
|
</style>
|
||||||
|
|
||||||
|
<defs id="Patterns_And_Gradients">
|
||||||
|
<pattern id="ptrn12-4" v:fillPattern="12" v:foreground="#000000" v:background="#ffffff" patternUnits="userSpaceOnUse"
|
||||||
|
width="6" height="6" viewBox="0 0 64 64">
|
||||||
|
<image x="0" y="0" width="64" height="64" image-rendering="optimizeSpeed"
|
||||||
|
xlink:href="data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAAgAAAAICAYAAADED76LAAAAAXNSR0IArs4c6QAAAARnQU1BAACxjwv8YQUAAAAJcEhZcwAADsMAAA7DAcdvqGQAAAAgSURBVChTY/iPBBgYULhggCmCBvAqAJlImgkD4ob//wHXAPMN73PsmQAAAABJRU5ErkJggg=="/>
|
||||||
|
</pattern>
|
||||||
|
<linearGradient id="grad0-15" x1="0" y1="0" x2="1" y2="0" gradientTransform="rotate(180 0.5 0.5)">
|
||||||
|
<stop offset="0" stop-color="#31859b" stop-opacity="1"/>
|
||||||
|
<stop offset="0.48" stop-color="#50aec7" stop-opacity="1"/>
|
||||||
|
<stop offset="1" stop-color="#92cddc" stop-opacity="1"/>
|
||||||
|
</linearGradient>
|
||||||
|
</defs>
|
||||||
|
<g v:mID="0" v:index="1" v:groupContext="foregroundPage">
|
||||||
|
<title>Page-1</title>
|
||||||
|
<v:pageProperties v:drawingScale="0.0393701" v:pageScale="0.0393701" v:drawingUnits="24" v:shadowOffsetX="8.50394"
|
||||||
|
v:shadowOffsetY="-8.50394"/>
|
||||||
|
<v:layer v:name="Electrical" v:index="0"/>
|
||||||
|
<v:layer v:name="Connector" v:index="1"/>
|
||||||
|
<g d:options="Port:R0" id="shape1-1" v:mID="1" v:groupContext="shape" transform="translate(136.164,-202.677)">
|
||||||
|
<title>Rectangle</title>
|
||||||
|
<v:userDefs>
|
||||||
|
<v:ud v:nameU="visVersion" v:val="VT0(15):26"/>
|
||||||
|
</v:userDefs>
|
||||||
|
<rect d:options="Port:T0" x="0" y="619.37" width="50.9224" height="222.52" class="st1"/>
|
||||||
|
</g>
|
||||||
|
<g id="shape4-5" v:mID="4" v:groupContext="shape" transform="translate(314.889,-202.677)">
|
||||||
|
<title>Rectangle.4</title>
|
||||||
|
<v:userDefs>
|
||||||
|
<v:ud v:nameU="visVersion" v:val="VT0(15):26"/>
|
||||||
|
</v:userDefs>
|
||||||
|
<rect d:options="Port:T0" x="0" y="619.37" width="38.9224" height="222.52" class="st1"/>
|
||||||
|
</g>
|
||||||
|
<g id="shape5-8" v:mID="5" v:groupContext="shape" transform="translate(180.904,-202.677)">
|
||||||
|
<title>Rectangle.5</title>
|
||||||
|
<v:userDefs>
|
||||||
|
<v:ud v:nameU="visVersion" v:val="VT0(15):26"/>
|
||||||
|
</v:userDefs>
|
||||||
|
<rect d:options="Port:T0" x="0" y="619.37" width="32.3465" height="222.52" class="st2"/>
|
||||||
|
</g>
|
||||||
|
<g id="shape6-10" v:mID="6" v:groupContext="shape" transform="translate(282.504,-202.677)">
|
||||||
|
<title>Rectangle.6</title>
|
||||||
|
<v:userDefs>
|
||||||
|
<v:ud v:nameU="visVersion" v:val="VT0(15):26"/>
|
||||||
|
</v:userDefs>
|
||||||
|
<rect d:options="Port:T0" x="0" y="619.37" width="32.3465" height="222.52" class="st3"/>
|
||||||
|
</g>
|
||||||
|
<g id="shape7-12" v:mID="7" v:groupContext="shape" transform="translate(212.425,-202.677)">
|
||||||
|
<title>Rectangle.7</title>
|
||||||
|
<v:userDefs>
|
||||||
|
<v:ud v:nameU="visVersion" v:val="VT0(15):26"/>
|
||||||
|
</v:userDefs>
|
||||||
|
<rect d:options="Port:T0" x="0" y="619.37" width="69.6929" height="222.52" class="st4"/>
|
||||||
|
</g>
|
||||||
|
</g>
|
||||||
|
</svg>
|
||||||
|
After Width: | Height: | Size: 4 KiB |
|
|
@ -0,0 +1,359 @@
|
||||||
|
component (Propagation = blocks) FuelCell
|
||||||
|
% Fuel Cell
|
||||||
|
% This block models a stack of membrane electrode assemblies (MEA) for a
|
||||||
|
% proton exchange membrane (PEM) fuel cell. Hydrogen and oxygen are
|
||||||
|
% consumed from the connected anode and cathode flow channels,
|
||||||
|
% respectively. Water vapor is produced at the cathode flow channel. The
|
||||||
|
% block also models water transport between the anode and cathode. Waste
|
||||||
|
% heat is dissipated to the connected thermal mass.
|
||||||
|
|
||||||
|
% Copyright 2020-2024 The MathWorks, Inc.
|
||||||
|
|
||||||
|
nodes
|
||||||
|
% Electrodes
|
||||||
|
n = foundation.electrical.electrical; % -
|
||||||
|
p = foundation.electrical.electrical; % +
|
||||||
|
% Moist air ports to obtain domain parameters and property tables for
|
||||||
|
% the anode network and the cathode network only. There is no flow through ports A and C
|
||||||
|
A = foundation.moist_air.moist_air; % A
|
||||||
|
C = foundation.moist_air.moist_air; % C
|
||||||
|
% Source ports to add and remove mass from the anode and cathode gas channels
|
||||||
|
% for the reaction in the fuel cell
|
||||||
|
SA = foundation.moist_air.moist_air_source; % SA
|
||||||
|
SC = foundation.moist_air.moist_air_source; % SC
|
||||||
|
end
|
||||||
|
|
||||||
|
inputs
|
||||||
|
% Input ports for the internal states of the anode and cathode gas channels
|
||||||
|
% Inflow states
|
||||||
|
FAi = [101325; 293.15; 0.5; 0.01; 0.01; 0.01; 0.001; 0.001; 0]; % FAi
|
||||||
|
FCi = [101325; 293.15; 0.5; 0.01; 0.01; 0.01; 0.001; 0.001; 0]; % FCi
|
||||||
|
% Outflow states
|
||||||
|
FAo = [101325; 293.15; 0.5; 0.01; 0.01; 0.01; 0.001; 0.001; 0]; % FAo
|
||||||
|
FCo = [101325; 293.15; 0.5; 0.01; 0.01; 0.01; 0.001; 0.001; 0]; % FCo
|
||||||
|
end
|
||||||
|
|
||||||
|
nodes
|
||||||
|
% Thermal port for heat generation
|
||||||
|
H = foundation.thermal.thermal; % H
|
||||||
|
end
|
||||||
|
|
||||||
|
annotations
|
||||||
|
[n, p] : Side = top;
|
||||||
|
[A, FAi, FAo, SA] : Side = left;
|
||||||
|
[C, FCi, FCo, SC] : Side = right;
|
||||||
|
H : Side = bottom;
|
||||||
|
end
|
||||||
|
|
||||||
|
parameters
|
||||||
|
N_cell = {400, '1' }; % Number of cells in stack
|
||||||
|
area_cell = {280, 'cm^2' }; % Cell area
|
||||||
|
t_membrane = {125, 'um' }; % Membrane thickness
|
||||||
|
t_gdl_A = {250, 'um' }; % Anode gas diffusion layer (GDL) thickness
|
||||||
|
t_gdl_C = {250, 'um' }; % Cathode gas diffusion layer (GDL) thickness
|
||||||
|
io = {8e-05, 'A/cm^2'}; % Exchange current density
|
||||||
|
iL = {1.4, 'A/cm^2'}; % Max (limiting) current density
|
||||||
|
alpha = {0.5, '1' }; % Charge transfer coefficient
|
||||||
|
D_H2O_gdl_A = {0.07, 'cm^2/s'}; % Water vapor diffusivity in anode GDL
|
||||||
|
D_H2O_gdl_C = {0.07, 'cm^2/s'}; % Water vapor diffusivity in cathode GDL
|
||||||
|
rho_membrane = {2000, 'kg/m^3'}; % Density of dry membrane
|
||||||
|
M_membrane = {1.1, 'kg/mol'}; % Equivalent weight of dry membrane
|
||||||
|
end
|
||||||
|
|
||||||
|
% Parameter checks
|
||||||
|
equations
|
||||||
|
assert(N_cell > 0)
|
||||||
|
assert(area_cell > 0)
|
||||||
|
assert(t_membrane > 0)
|
||||||
|
assert(t_gdl_A > 0)
|
||||||
|
assert(t_gdl_C > 0)
|
||||||
|
assert(io > 0)
|
||||||
|
assert(iL > io)
|
||||||
|
assert(alpha > 0)
|
||||||
|
assert(D_H2O_gdl_A > 0)
|
||||||
|
assert(D_H2O_gdl_C > 0)
|
||||||
|
assert(rho_membrane > 0)
|
||||||
|
assert(M_membrane > 0)
|
||||||
|
end
|
||||||
|
|
||||||
|
% Internal parameters
|
||||||
|
parameters (Access = protected)
|
||||||
|
% Constants
|
||||||
|
R_u = {8.31446261815324, 'J/K/mol'}; % Universal gas constant
|
||||||
|
F = {96485.33212, 'C/mol' }; % Faraday constant
|
||||||
|
G_H2O = {-237.14, 'kJ/mol' }; % Gibbs free energy of water
|
||||||
|
HHV_H2 = {285.8, 'kJ/mol' }; % Higher heating value of hydrogen
|
||||||
|
LHV_H2 = HHV_H2 - MW_H2O * hfg_H2O_std; % Lower heating value of hydrogen
|
||||||
|
T_std = {25, 'degC' }; % Standard temperature
|
||||||
|
p_std = {1, 'atm' }; % Standard pressure
|
||||||
|
|
||||||
|
% Standard cell potential
|
||||||
|
E_cell = G_H2O/(-2*F);
|
||||||
|
|
||||||
|
% Membrane permeability to water
|
||||||
|
K_darcy = {1.58e-14, 'cm^2'};
|
||||||
|
|
||||||
|
% Compute molar masses from specific gas constants
|
||||||
|
MW_H2O = R_u/C.R_w;
|
||||||
|
MW_O2 = R_u/C.R_g;
|
||||||
|
MW_H2 = R_u/A.R_g;
|
||||||
|
|
||||||
|
% Species specific enthalpy at standard temperature
|
||||||
|
h_H2_std = tablelookup(A.T_TLU, A.h_g_TLU, T_std, interpolation = linear, extrapolation = linear);
|
||||||
|
h_O2_std = tablelookup(C.T_TLU, C.h_g_TLU, T_std, interpolation = linear, extrapolation = linear);
|
||||||
|
h_H2O_std = tablelookup(C.T_TLU, C.h_w_TLU, T_std, interpolation = linear, extrapolation = linear);
|
||||||
|
hfg_H2O_std = tablelookup(C.T_TLU, C.h_w_vap_TLU, T_std, interpolation = linear, extrapolation = linear);
|
||||||
|
end
|
||||||
|
|
||||||
|
variables (ExternalAccess = observe)
|
||||||
|
i = {0, 'A' }; % Current (positive in)
|
||||||
|
Q = {0, 'kW'}; % Heat flow rate (positive in)
|
||||||
|
|
||||||
|
a_H2O_acl = {0.5, '1'}; % Water vapor activity at anode catalyst layer (ACL)
|
||||||
|
a_H2O_ccl = {0.5, '1'}; % Water vapor activity at cathode catalyst layer (CCL)
|
||||||
|
end
|
||||||
|
|
||||||
|
branches
|
||||||
|
i : p.i -> n.i;
|
||||||
|
Q : H.Q -> *;
|
||||||
|
end
|
||||||
|
|
||||||
|
% For logging
|
||||||
|
intermediates
|
||||||
|
v = p.v - n.v; % Voltage
|
||||||
|
|
||||||
|
i_cell = if le(i, 0), -i/area_cell else 0 end; % Cell current density
|
||||||
|
v_cell = v_nernst - v_act - v_ohm - v_conc; % Cell voltage
|
||||||
|
R_ohm = t_membrane / sigma; % Membrane resistance
|
||||||
|
|
||||||
|
power_elec = N_cell * v_cell * i_cell * area_cell; % Electrical power produced
|
||||||
|
power_dissipated = power_net - power_elec; % Power dissipated
|
||||||
|
efficiency_HHV = simscape.function.limit(v_cell / (HHV_H2/(2*F)), 0, 1, false); % Thermal efficiency based on HHV
|
||||||
|
efficiency_LHV = simscape.function.limit(v_cell / (LHV_H2/(2*F)), 0, 1, false); % Thermal efficiency based on LHV
|
||||||
|
|
||||||
|
H2_consumed = N_cell * MW_H2 * area_cell * i_cell/(2*F); % Hydrogen consumed
|
||||||
|
O2_consumed = N_cell * MW_O2 * area_cell * i_cell/(4*F); % Oxygen consumed
|
||||||
|
H2O_produced = N_cell * MW_H2O * area_cell * i_cell/(2*F); % Water produced
|
||||||
|
|
||||||
|
H2O_transport = (nflux_H2O_diff + nflux_H2O_drag + nflux_H2O_hydraulic) ...
|
||||||
|
* MW_H2O * area_cell * N_cell; % Net water transport from anode to cathode
|
||||||
|
end
|
||||||
|
|
||||||
|
annotations
|
||||||
|
[v, v_cell] : LoggingUnit = 'V';
|
||||||
|
i_cell : LoggingUnit = 'A/cm^2';
|
||||||
|
R_ohm : LoggingUnit = 'Ohm*cm^2';
|
||||||
|
[power_elec, power_dissipated] : LoggingUnit = 'kW';
|
||||||
|
[efficiency_HHV, efficiency_LHV] : LoggingUnit = '1';
|
||||||
|
[H2_consumed, O2_consumed, H2O_produced, H2O_transport] : LoggingUnit = 'g/s';
|
||||||
|
end
|
||||||
|
|
||||||
|
intermediates (ExternalAccess = none)
|
||||||
|
% Extract anode gas states from input FAi
|
||||||
|
% Note: Can also connect Measurement Selector (MA) block to choose values
|
||||||
|
p_Ai = {FAi(1), 'Pa'}; % Pressure
|
||||||
|
y_H2O_Ai = FAi(5); % Water vapor mole fraction
|
||||||
|
y_H2_Ai = FAi(8); % Hydrogen mole fraction
|
||||||
|
|
||||||
|
% Extract anode gas states from input FAo
|
||||||
|
% Note: Can also connect Measurement Selector (MA) block to choose values
|
||||||
|
p_Ao = {FAo(1), 'Pa'}; % Pressure
|
||||||
|
y_H2O_Ao = FAo(5); % Water vapor mole fraction
|
||||||
|
y_H2_Ao = FAo(8); % Hydrogen mole fraction
|
||||||
|
|
||||||
|
% Extract cathode gas states from input FCi
|
||||||
|
% Note: Can also connect Measurement Select (MA) block to choose values
|
||||||
|
p_Ci = {FCi(1), 'Pa'}; % Pressure
|
||||||
|
y_H2O_Ci = FCi(5); % Water vapor mole fraction
|
||||||
|
y_O2_Ci = FCi(8); % Oxygen mole fraction
|
||||||
|
|
||||||
|
% Extract cathode gas states from input FCo
|
||||||
|
% Note: Can also connect Measurement Select (MA) block to choose values
|
||||||
|
p_Co = {FCo(1), 'Pa'}; % Pressure
|
||||||
|
y_H2O_Co = FCo(5); % Water vapor mole fraction
|
||||||
|
y_O2_Co = FCo(8); % Oxygen mole fraction
|
||||||
|
|
||||||
|
% Stack temperature
|
||||||
|
T_stack = H.T;
|
||||||
|
|
||||||
|
% Anode and cathode pressure
|
||||||
|
p_A = (p_Ai + p_Ao)/2;
|
||||||
|
p_C = (p_Ci + p_Co)/2;
|
||||||
|
|
||||||
|
% Ratio of pressure to water vapor saturation pressure
|
||||||
|
p_ws_ratio_A = exp(log(value(p_A, 'Pa')) - tablelookup(A.T_TLU, A.log_p_ws_TLU, T_stack, interpolation = linear, extrapolation = linear));
|
||||||
|
p_ws_ratio_C = exp(log(value(p_C, 'Pa')) - tablelookup(C.T_TLU, C.log_p_ws_TLU, T_stack, interpolation = linear, extrapolation = linear));
|
||||||
|
|
||||||
|
% Mole fractions at anode and cathode
|
||||||
|
y_H2_A = (y_H2_Ai + y_H2_Ao)/2;
|
||||||
|
y_O2_C = (y_O2_Ci + y_O2_Co)/2;
|
||||||
|
y_H2O_C = (y_H2O_Ci + y_H2O_Co)/2;
|
||||||
|
|
||||||
|
% Activities
|
||||||
|
a_H2_A_ = y_H2_A * p_A / p_std;
|
||||||
|
a_O2_C_ = y_O2_C * p_C / p_std;
|
||||||
|
a_H2O_C_ = y_H2O_C * p_ws_ratio_C;
|
||||||
|
a_H2_A = if ge(a_H2_A_, 1e-9), a_H2_A_ else 1e-6 end;
|
||||||
|
a_O2_C = if ge(a_O2_C_, 1e-9), a_O2_C_ else 1e-6 end;
|
||||||
|
a_H2O_C = if ge(a_H2O_C_, 1e-9), a_H2O_C_ else 1e-6 end;
|
||||||
|
|
||||||
|
% Nernst voltage
|
||||||
|
v_nernst = E_cell + R_u*T_stack/(2*F) * log((a_H2_A * a_O2_C^0.5) / a_H2O_C);
|
||||||
|
|
||||||
|
% Activation losses from Tafel equation
|
||||||
|
b = R_u * T_stack / (2 * alpha * F);
|
||||||
|
v_act = if ge(i_cell, io), b*log(i_cell/io) else 0 end;
|
||||||
|
|
||||||
|
% Gas transport voltage loss
|
||||||
|
v_conc = -R_u*T_stack/(2*F) * ...
|
||||||
|
if le(i_cell, 0.999*iL), ...
|
||||||
|
log(1 - i_cell/iL) ...
|
||||||
|
else ...
|
||||||
|
log(1 - 0.999) - (i_cell/iL - 0.999)/(1 - 0.999) ...
|
||||||
|
end;
|
||||||
|
|
||||||
|
% Water content
|
||||||
|
lambda_acl = membrane_water(a_H2O_acl);
|
||||||
|
lambda_ccl = membrane_water(a_H2O_ccl);
|
||||||
|
lambda_membrane = (lambda_acl + lambda_ccl)/2;
|
||||||
|
|
||||||
|
% Membrane conductivity
|
||||||
|
sigma_30 = {if ge(lambda_membrane, 1), 0.005139*lambda_membrane - 0.00326 else 0.005139 - 0.00326 end, '1/(Ohm*cm)'};
|
||||||
|
sigma = sigma_30 * exp(1268*(1/303.15 - 1/value(T_stack, 'K')));
|
||||||
|
|
||||||
|
% Resistive voltage loss
|
||||||
|
v_ohm = R_ohm * i_cell;
|
||||||
|
|
||||||
|
% Water diffusion coefficient across membrane
|
||||||
|
D_H2O_membrane = {1.25e-10, 'm^2/s'} * exp(2416*(1/303.15 - 1/value(T_stack, 'K')));
|
||||||
|
|
||||||
|
% Water concentrations at anode and cathode catalyst layers
|
||||||
|
Conc_H2O_acl = rho_membrane / M_membrane * lambda_acl;
|
||||||
|
Conc_H2O_ccl = rho_membrane / M_membrane * lambda_ccl;
|
||||||
|
|
||||||
|
% Molar flux of water across membrane due to diffusion
|
||||||
|
nflux_H2O_diff = D_H2O_membrane*(Conc_H2O_acl - Conc_H2O_ccl)/t_membrane;
|
||||||
|
|
||||||
|
% Water electro-osmotic drag coefficient
|
||||||
|
nd_H2O_membrane = ...
|
||||||
|
if ge(lambda_acl, 0), ...
|
||||||
|
0.0029*lambda_acl^2 + 0.05*lambda_acl ...
|
||||||
|
else ...
|
||||||
|
0.05*lambda_acl ...
|
||||||
|
end;
|
||||||
|
|
||||||
|
% Molar flux of water across membrane due to electro-osmotic drag
|
||||||
|
nflux_H2O_drag = nd_H2O_membrane * i_cell / F;
|
||||||
|
|
||||||
|
% Water vapor mole fraction at anode and cathode catalyst layers
|
||||||
|
y_H2O_acl = a_H2O_acl/p_ws_ratio_A;
|
||||||
|
y_H2O_ccl = a_H2O_ccl/p_ws_ratio_C
|
||||||
|
|
||||||
|
% Molar flux of water across anode and cathode diffusion layers
|
||||||
|
nflux_H2O_A = p_A*D_H2O_gdl_A/(R_u*T_stack) * (y_H2O_Ao - y_H2O_acl)/t_gdl_A;
|
||||||
|
nflux_H2O_C = p_C*D_H2O_gdl_C/(R_u*T_stack) * (y_H2O_ccl - y_H2O_Co)/t_gdl_C;
|
||||||
|
|
||||||
|
% Molar flux of water across membrane due to hydraulic pressure difference based on Darcy's law
|
||||||
|
mu_H2O = tablelookup(A.T_TLU, A.mu_w_TLU, T_stack, interpolation = linear, extrapolation = nearest);
|
||||||
|
nflux_H2O_hydraulic = ...
|
||||||
|
if gt(p_A, p_C), ...
|
||||||
|
(p_A - p_C) * K_darcy * p_A * y_H2O_acl / (R_u * T_stack * mu_H2O * t_membrane) ...
|
||||||
|
else...
|
||||||
|
(p_A - p_C) * K_darcy * p_C * y_H2O_ccl / (R_u * T_stack * mu_H2O * t_membrane) ...
|
||||||
|
end
|
||||||
|
|
||||||
|
% Energy generation due to reaction at standard temperature
|
||||||
|
power_rxn = LHV_H2 * H2_consumed / MW_H2;
|
||||||
|
|
||||||
|
% Energy gain in membrane due to bringing reactants and products to standard temperature
|
||||||
|
power_delta_std = (source_H2_A.Phi_S - source_H2_A.mdot_g_S*h_H2_std) ...
|
||||||
|
+ (source_O2_C.Phi_S - source_O2_C.mdot_g_S*h_O2_std) ...
|
||||||
|
+ (source_H2O_C.Phi_S - source_H2O_C.mdot_w_S*h_H2O_std);
|
||||||
|
|
||||||
|
% Energy gain in membrane due to water transport
|
||||||
|
power_trans = transport_H2O_A.Phi_S + transport_H2O_C.Phi_S;
|
||||||
|
|
||||||
|
% Net energy gain in stack
|
||||||
|
power_net = power_rxn + power_delta_std + power_trans;
|
||||||
|
end
|
||||||
|
|
||||||
|
equations
|
||||||
|
% Stack voltage
|
||||||
|
v == N_cell * v_cell;
|
||||||
|
|
||||||
|
% Heat generated
|
||||||
|
-Q == power_dissipated;
|
||||||
|
|
||||||
|
% Equate water vapor mass flow rates at GDL and membrane
|
||||||
|
% to solve for water vapor activity at ACL and CCL
|
||||||
|
nflux_H2O_A * MW_H2O * area_cell * N_cell == H2O_transport;
|
||||||
|
nflux_H2O_C * MW_H2O * area_cell * N_cell == H2O_transport + H2O_produced;
|
||||||
|
|
||||||
|
% Assign mass flow rate to the internal moisture source blocks
|
||||||
|
% to model transport of water across membrane
|
||||||
|
transport_H2O_A.M == -H2O_transport;
|
||||||
|
transport_H2O_C.M == H2O_transport;
|
||||||
|
transport_H2O_A.T == T_stack;
|
||||||
|
transport_H2O_C.T == T_stack;
|
||||||
|
|
||||||
|
% Assign mass flow rate to the internal moisture source and trace gas source blocks
|
||||||
|
% to model mass consumption/generation due to reaction
|
||||||
|
source_H2_A.M == -H2_consumed;
|
||||||
|
source_O2_C.M == -O2_consumed;
|
||||||
|
source_H2O_C.M == H2O_produced;
|
||||||
|
source_H2_A.T == T_stack;
|
||||||
|
source_O2_C.T == T_stack;
|
||||||
|
source_H2O_C.T == T_stack;
|
||||||
|
end
|
||||||
|
|
||||||
|
% Internal components to add or remove mass at the anode and cathode
|
||||||
|
components (ExternalAccess = none)
|
||||||
|
% Moisture source blocks for water transport
|
||||||
|
transport_H2O_A = foundation.moist_air.sources.moisture_trace_gas.moisture_source( ...
|
||||||
|
source_type = foundation.enum.constant_controlled.controlled, ...
|
||||||
|
moisture_source_spec = foundation.enum.MoistureSourceSpec.WaterVapor, ...
|
||||||
|
vapor_from_liquid = false);
|
||||||
|
transport_H2O_C = foundation.moist_air.sources.moisture_trace_gas.moisture_source( ...
|
||||||
|
source_type = foundation.enum.constant_controlled.controlled, ...
|
||||||
|
moisture_source_spec = foundation.enum.MoistureSourceSpec.WaterVapor, ...
|
||||||
|
vapor_from_liquid = false);
|
||||||
|
|
||||||
|
% Moisture source and trace gas source blocks for reaction
|
||||||
|
source_H2_A = foundation.moist_air.sources.moisture_trace_gas.trace_gas_source( ...
|
||||||
|
source_type = foundation.enum.constant_controlled.controlled);
|
||||||
|
source_O2_C = foundation.moist_air.sources.moisture_trace_gas.trace_gas_source( ...
|
||||||
|
source_type = foundation.enum.constant_controlled.controlled);
|
||||||
|
source_H2O_C = foundation.moist_air.sources.moisture_trace_gas.moisture_source( ...
|
||||||
|
source_type = foundation.enum.constant_controlled.controlled, ...
|
||||||
|
moisture_source_spec = foundation.enum.MoistureSourceSpec.WaterVapor, ...
|
||||||
|
vapor_from_liquid = false);
|
||||||
|
end
|
||||||
|
|
||||||
|
connections
|
||||||
|
connect(SA, transport_H2O_A.S, source_H2_A.S)
|
||||||
|
connect(SC, transport_H2O_C.S, source_O2_C.S, source_H2O_C.S)
|
||||||
|
end
|
||||||
|
|
||||||
|
end
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
function lambda = membrane_water(a)
|
||||||
|
% Compute the water content of the membrane as a function of the water
|
||||||
|
% activity. This function is based on fits to FEA simulations.
|
||||||
|
%
|
||||||
|
% See: Dutta, et.al., Numerical prediction of mass-exchange between
|
||||||
|
% cathode and anode channels in a PEM fuel cell. (Equation 16)
|
||||||
|
|
||||||
|
definitions
|
||||||
|
lambda = ...
|
||||||
|
if lt(a, 0), ...
|
||||||
|
0.043 + 17.81*a ...
|
||||||
|
elseif le(a, 1), ...
|
||||||
|
0.043 + 17.81*a - 39.85*a^2 + 36*a^3 ...
|
||||||
|
else ...
|
||||||
|
14.003 + 1.4*(a - 1) ...
|
||||||
|
end;
|
||||||
|
end
|
||||||
|
|
||||||
|
end
|
||||||
|
|
@ -0,0 +1,80 @@
|
||||||
|
<?xml-stylesheet href="SimscapeIconStyles.css" type="text/css"?>
|
||||||
|
<!-- Generated by Microsoft Visio, SVG Export cell_voltage6.svg Page-1 -->
|
||||||
|
<svg xmlns="http://www.w3.org/2000/svg" xmlns:xlink="http://www.w3.org/1999/xlink"
|
||||||
|
xmlns:ev="http://www.w3.org/2001/xml-events"
|
||||||
|
xmlns:v="http://schemas.microsoft.com/visio/2003/SVGExtensions/"
|
||||||
|
xmlns:d="http://www.mathworks.com/blockgraphics"
|
||||||
|
width="215" height="222"
|
||||||
|
viewBox="136 417.37 215.09 222"
|
||||||
|
xml:space="preserve" color-interpolation-filters="sRGB" class="st5">
|
||||||
|
<v:documentProperties v:langID="1033" v:metric="true" v:viewMarkup="false">
|
||||||
|
<v:userDefs>
|
||||||
|
<v:ud v:nameU="msvNoAutoConnect" v:val="VT0(1):26"/>
|
||||||
|
</v:userDefs>
|
||||||
|
</v:documentProperties>
|
||||||
|
|
||||||
|
<style type="text/css">
|
||||||
|
<![CDATA[
|
||||||
|
.st1 {fill:url(#ptrn12-4);shape-rendering:crispEdges;stroke:#000000;stroke-linecap:round;stroke-linejoin:round;stroke-width:0.24}
|
||||||
|
.st2 {fill:#a5a5a5;stroke:#000000;stroke-linecap:round;stroke-linejoin:round;stroke-width:0.24}
|
||||||
|
.st3 {fill:#bfbfbf;stroke:#000000;stroke-linecap:round;stroke-linejoin:round;stroke-width:0.24}
|
||||||
|
.st4 {fill:url(#grad0-15);stroke:#000000;stroke-linecap:round;stroke-linejoin:round;stroke-width:0.24}
|
||||||
|
.st5 {fill:none;fill-rule:evenodd;font-size:12px;overflow:visible;stroke-linecap:square;stroke-miterlimit:3}
|
||||||
|
]]>
|
||||||
|
</style>
|
||||||
|
|
||||||
|
<defs id="Patterns_And_Gradients">
|
||||||
|
<pattern id="ptrn12-4" v:fillPattern="12" v:foreground="#000000" v:background="#ffffff" patternUnits="userSpaceOnUse"
|
||||||
|
width="6" height="6" viewBox="0 0 64 64">
|
||||||
|
<image x="0" y="0" width="64" height="64" image-rendering="optimizeSpeed"
|
||||||
|
xlink:href="data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAAgAAAAICAYAAADED76LAAAAAXNSR0IArs4c6QAAAARnQU1BAACxjwv8YQUAAAAJcEhZcwAADsMAAA7DAcdvqGQAAAAgSURBVChTY/iPBBgYULhggCmCBvAqAJlImgkD4ob//wHXAPMN73PsmQAAAABJRU5ErkJggg=="/>
|
||||||
|
</pattern>
|
||||||
|
<linearGradient id="grad0-15" x1="0" y1="0" x2="1" y2="0" gradientTransform="rotate(180 0.5 0.5)">
|
||||||
|
<stop offset="0" stop-color="#31859b" stop-opacity="1"/>
|
||||||
|
<stop offset="0.48" stop-color="#50aec7" stop-opacity="1"/>
|
||||||
|
<stop offset="1" stop-color="#92cddc" stop-opacity="1"/>
|
||||||
|
</linearGradient>
|
||||||
|
</defs>
|
||||||
|
<g v:mID="0" v:index="1" v:groupContext="foregroundPage">
|
||||||
|
<title>Page-1</title>
|
||||||
|
<v:pageProperties v:drawingScale="0.0393701" v:pageScale="0.0393701" v:drawingUnits="24" v:shadowOffsetX="8.50394"
|
||||||
|
v:shadowOffsetY="-8.50394"/>
|
||||||
|
<v:layer v:name="Electrical" v:index="0"/>
|
||||||
|
<v:layer v:name="Connector" v:index="1"/>
|
||||||
|
<g d:options="Port:R0" id="shape1-1" v:mID="1" v:groupContext="shape" transform="translate(136.164,-202.677)">
|
||||||
|
<title>Rectangle</title>
|
||||||
|
<v:userDefs>
|
||||||
|
<v:ud v:nameU="visVersion" v:val="VT0(15):26"/>
|
||||||
|
</v:userDefs>
|
||||||
|
<rect d:options="Port:T0" x="0" y="619.37" width="50.9224" height="222.52" class="st1"/>
|
||||||
|
</g>
|
||||||
|
<g id="shape4-5" v:mID="4" v:groupContext="shape" transform="translate(314.889,-202.677)">
|
||||||
|
<title>Rectangle.4</title>
|
||||||
|
<v:userDefs>
|
||||||
|
<v:ud v:nameU="visVersion" v:val="VT0(15):26"/>
|
||||||
|
</v:userDefs>
|
||||||
|
<rect d:options="Port:T0" x="0" y="619.37" width="38.9224" height="222.52" class="st1"/>
|
||||||
|
</g>
|
||||||
|
<g id="shape5-8" v:mID="5" v:groupContext="shape" transform="translate(180.904,-202.677)">
|
||||||
|
<title>Rectangle.5</title>
|
||||||
|
<v:userDefs>
|
||||||
|
<v:ud v:nameU="visVersion" v:val="VT0(15):26"/>
|
||||||
|
</v:userDefs>
|
||||||
|
<rect d:options="Port:T0" x="0" y="619.37" width="32.3465" height="222.52" class="st2"/>
|
||||||
|
</g>
|
||||||
|
<g id="shape6-10" v:mID="6" v:groupContext="shape" transform="translate(282.504,-202.677)">
|
||||||
|
<title>Rectangle.6</title>
|
||||||
|
<v:userDefs>
|
||||||
|
<v:ud v:nameU="visVersion" v:val="VT0(15):26"/>
|
||||||
|
</v:userDefs>
|
||||||
|
<rect d:options="Port:T0" x="0" y="619.37" width="32.3465" height="222.52" class="st3"/>
|
||||||
|
</g>
|
||||||
|
<g id="shape7-12" v:mID="7" v:groupContext="shape" transform="translate(212.425,-202.677)">
|
||||||
|
<title>Rectangle.7</title>
|
||||||
|
<v:userDefs>
|
||||||
|
<v:ud v:nameU="visVersion" v:val="VT0(15):26"/>
|
||||||
|
</v:userDefs>
|
||||||
|
<rect d:options="Port:T0" x="0" y="619.37" width="69.6929" height="222.52" class="st4"/>
|
||||||
|
</g>
|
||||||
|
</g>
|
||||||
|
</svg>
|
||||||
|
After Width: | Height: | Size: 4 KiB |
Binary file not shown.
Binary file not shown.
Binary file not shown.
|
|
@ -0,0 +1,154 @@
|
||||||
|
%% PEM Electrolysis System
|
||||||
|
%
|
||||||
|
% This example shows how to model a proton exchange membrane (PEM) water
|
||||||
|
% electrolyzer with a custom Simscape(TM) block. The PEM electrolyzer
|
||||||
|
% consumes electrical power to split water into hydrogen and oxygen. The
|
||||||
|
% custom block represents the membrane electrode assembly (MEA) and is
|
||||||
|
% connected to a thermal liquid network and two separate moist air
|
||||||
|
% networks: the thermal liquid network models the water supply, the anode
|
||||||
|
% moist air network models the oxygen flow, and the cathode moist air
|
||||||
|
% network models the hydrogen flow.
|
||||||
|
%
|
||||||
|
% The circulation pump provides a continuous flow of water to the anode
|
||||||
|
% side of the electrolyzer. Consumed water is removed from the thermal
|
||||||
|
% liquid network and excess water is recirculated. Oxygen produced at the
|
||||||
|
% anode is carried away by the excess water flow; it is modeled separately
|
||||||
|
% by the anode moist air network. The separator tank models the balance of
|
||||||
|
% water and oxygen in the return flow before the oxygen is vented. The
|
||||||
|
% supply pump replenishes the system with fresh water.
|
||||||
|
%
|
||||||
|
% Hydrogen produced at the cathode side along with any water that has been
|
||||||
|
% transported across the MEA is modeled by the cathode moist air network.
|
||||||
|
% The dehumidifier removes the unwanted water vapor from the hydrogen. A
|
||||||
|
% pressure regulator valve maintains a pressure of 3 MPa at the cathode in
|
||||||
|
% comparison to atmospheric pressure at the anode. The differential
|
||||||
|
% pressure across the MEA results in water transport due to hydraulic
|
||||||
|
% pressure that helps counteract the electro-osmosis drag and reduces the
|
||||||
|
% amount of water at the cathode side.
|
||||||
|
%
|
||||||
|
% Unlike a fuel cell stack, a separate cooling network is not needed. Heat
|
||||||
|
% dissipated by the electrolyzer is carried away by the excess water and
|
||||||
|
% then rejected to the environment via the heat exchanger. The
|
||||||
|
% recirculating water is controlled to maintain a temperature of 80 degC in
|
||||||
|
% the electrolylzer.
|
||||||
|
%
|
||||||
|
% The custom MEA block is implemented in the Simscape code
|
||||||
|
% |<matlab:edit('Electrolyzer.ssc') Electrolyzer.ssc>|. The thermal liquid
|
||||||
|
% port H2O is used to remove water from the thermal liquid network. The
|
||||||
|
% produced H2 and O2 and the transported H2O are added to the two moist air
|
||||||
|
% networks using Controlled Trace Gas Source (MA) and Controlled Moisture
|
||||||
|
% Source (MA) blocks. The excess heat is sent through the thermal port H to
|
||||||
|
% the connected Thermal Mass block. Refer to the comments in the code for
|
||||||
|
% additional details on the implementation.
|
||||||
|
%
|
||||||
|
% See also the <docid:simscape_ug#example-ssc_fuel_cell PEM Fuel Cell
|
||||||
|
% System> example.
|
||||||
|
%
|
||||||
|
% References:
|
||||||
|
%
|
||||||
|
% Liso, Vincenzo, et al. "Modelling and experimental analysis of a polymer
|
||||||
|
% electrolyte membrane water electrolysis cell at different operating
|
||||||
|
% temperatures." _Energies_ 11.12 (2018): 3273.
|
||||||
|
%
|
||||||
|
% Mo, Jingke, et al. "Thin liquid/gas diffusion layers for high-efficiency
|
||||||
|
% hydrogen production from water splitting." _Applied Energy_ 177 (2016):
|
||||||
|
% 817-822.
|
||||||
|
|
||||||
|
% Copyright 2021 The MathWorks, Inc.
|
||||||
|
|
||||||
|
|
||||||
|
%% Model
|
||||||
|
|
||||||
|
open_system('PEMElectrolysisSystem')
|
||||||
|
|
||||||
|
set_param(find_system(bdroot,'FindAll','on','type','annotation','Tag','ModelFeatures'),'Interpreter','off');
|
||||||
|
|
||||||
|
%% Anode Fluid Channels Subsystem
|
||||||
|
|
||||||
|
open_system('PEMElectrolysisSystem/Anode Fluid Channels','force')
|
||||||
|
|
||||||
|
%% Cathode Gas Channels Subsystem
|
||||||
|
|
||||||
|
open_system('PEMElectrolysisSystem/Cathode Gas Channels','force')
|
||||||
|
|
||||||
|
%% Dehumidifier Subsystem
|
||||||
|
|
||||||
|
open_system('PEMElectrolysisSystem/Dehumidifier','force')
|
||||||
|
|
||||||
|
%% Electrical Supply Subsystem
|
||||||
|
|
||||||
|
open_system('PEMElectrolysisSystem/Electrical Supply','force')
|
||||||
|
|
||||||
|
%% Heat Exchanger Subsystem
|
||||||
|
|
||||||
|
open_system('PEMElectrolysisSystem/Heat Exchanger','force')
|
||||||
|
|
||||||
|
%% Hydrogen Output Subsystem
|
||||||
|
|
||||||
|
open_system('PEMElectrolysisSystem/Hydrogen Output','force')
|
||||||
|
|
||||||
|
%% Recirculation Subsystem
|
||||||
|
|
||||||
|
open_system('PEMElectrolysisSystem/Recirculation','force')
|
||||||
|
|
||||||
|
%% Separator Tank Subsystem
|
||||||
|
|
||||||
|
open_system('PEMElectrolysisSystem/Recirculation/Separator Tank','force')
|
||||||
|
|
||||||
|
%% Water Supply Subsystem
|
||||||
|
|
||||||
|
open_system('PEMElectrolysisSystem/Water Supply','force')
|
||||||
|
|
||||||
|
%% Simulation Results from Scopes
|
||||||
|
|
||||||
|
set_param('PEMElectrolysisSystem/Scope','open','on');
|
||||||
|
sim('PEMElectrolysisSystem');
|
||||||
|
|
||||||
|
%%
|
||||||
|
|
||||||
|
set_param('PEMElectrolysisSystem/Scope','open','off');
|
||||||
|
%% Simulation Results from Simscape Logging
|
||||||
|
%%
|
||||||
|
%
|
||||||
|
% This plot shows the current-voltage (i-v) curve and the power consumed by
|
||||||
|
% a cell in the stack. As the current ramps up, an initial rise in voltage
|
||||||
|
% occurs due to electrode activation losses, followed by a gradual increase
|
||||||
|
% in voltage due to Ohmic resistances. The cell voltage is about 1.71 V at
|
||||||
|
% a current density of 2 A/cm^2.
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
PEMElectrolysisSystemPlot1IV;
|
||||||
|
%%
|
||||||
|
%
|
||||||
|
% This plot shows the electrical power consumed by the electrolyzer. The
|
||||||
|
% electrical power is greater than the power needed to produce hydrogen due
|
||||||
|
% to various losses. The difference is the heat dissipated.
|
||||||
|
%
|
||||||
|
% This plot also shows the thermal efficiency of the electrolyzer, which
|
||||||
|
% indicates the fraction of the electrical power used to generate hydrogen
|
||||||
|
% based on hydrogen's heating value. This electrolyzer is about 87%
|
||||||
|
% efficient at a current density of 2 A/cm^2.
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
PEMElectrolysisSystemPlot2Power;
|
||||||
|
%%
|
||||||
|
%
|
||||||
|
% This plot shows the rate of hydrogen produced, the rate of water consumed
|
||||||
|
% at the anode, as well as the rate of water transported to the cathode due
|
||||||
|
% to diffusion, electro-osmosis drag, and hydraulic pressure difference. As
|
||||||
|
% a result, a dehumidification step is necessary to produce hydrogen at the
|
||||||
|
% desired purity.
|
||||||
|
%
|
||||||
|
% This plot also shows the total mass of hydrogen produced and the
|
||||||
|
% equivalent energy based on its higher heating value. This provides an
|
||||||
|
% indication of the amount of energy available if the hydrogen is used to
|
||||||
|
% generate power in a fuel cell.
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
PEMElectrolysisSystemPlot3Hydrogen;
|
||||||
|
|
||||||
|
%%
|
||||||
|
|
||||||
|
|
@ -0,0 +1,72 @@
|
||||||
|
%% Code to define parameters for PEMElectrolysisSystem
|
||||||
|
% Open Model Workspace in the Model Explorer to view and modify parameter
|
||||||
|
% values. Click 'Reinitialize from Source' to reset to the parameter values
|
||||||
|
% in this script.
|
||||||
|
|
||||||
|
% Copyright 2021 The MathWorks, Inc.
|
||||||
|
|
||||||
|
% Environment conditions
|
||||||
|
env_p = 0.101325; % [MPa] Pressure
|
||||||
|
env_T = 20; % [degC] Temperature
|
||||||
|
|
||||||
|
% Example solar power profile over 24 hr
|
||||||
|
solar_profile_time = (0:23)' * 3600; % s
|
||||||
|
solar_profile_power = [0 0 0 0 0 0 0.572 13.1 44.9 74.5 96.4 100 97.5 89.5 65.8 34.0 7.61 0 0 0 0 0 0 0]'; % kW
|
||||||
|
|
||||||
|
% Fuel cell stack
|
||||||
|
stack_num_cells = 50; % [-] Number cells
|
||||||
|
stack_area = 280; % [cm^2] Cell area
|
||||||
|
stack_t_membrane = 125; % [um] Membrane thickness
|
||||||
|
stack_t_gdl_A = 25; % [um] Anode gas diffusion layer thickness
|
||||||
|
stack_t_gdl_C = 250; % [um] Cathode gas diffusion layer thickness
|
||||||
|
stack_w_channels = 1; % [cm] Gas channel width and height
|
||||||
|
stack_num_channels = 8; % [-] Number of gas channels per cell
|
||||||
|
stack_io = 1e-04; % [A/cm^2] Exchange current density
|
||||||
|
stack_alpha = 0.7; % [-] Charge transfer coefficient
|
||||||
|
stack_D_gdl_A = 0.07; % [cm^2/s] Water diffusivity in anode GDL
|
||||||
|
stack_D_gdl_C = 0.07; % [cm^2/s] Water diffusivity in cathode GDL
|
||||||
|
stack_membrane_rho = 2000; % [kg/m^3] Density of dry membrane
|
||||||
|
stack_membrane_MW = 1.1; % [kg/mol] Equivalent weight of dry membrane
|
||||||
|
stack_mea_rho = 1800; % [kg/s] Overall density of membrane electrode assembly
|
||||||
|
stack_mea_cp = 870; % [J/(kg*K)] Overall specific heat of membrane electrode assembly
|
||||||
|
|
||||||
|
water_pipe_D = 0.01; % m
|
||||||
|
gas_pipe_D = 0.01; % m
|
||||||
|
|
||||||
|
% Heat exchanger dimensions
|
||||||
|
exchanger_L = 1; % [m] Overall radiator length
|
||||||
|
exchanger_W = 0.025; % [m] Overall radiator width
|
||||||
|
exchanger_H = 0.5; % [m] Overal radiator height
|
||||||
|
exchanger_N_tubes = 25; % [-] Number of coolant tubes
|
||||||
|
exchanger_tube_H = 0.0015; % [m] Height of each coolant tube
|
||||||
|
exchanger_fin_spacing = 0.002; % [-] Fin spacing
|
||||||
|
exchanger_eta_fin = 0.7; % [-] Fin efficiency
|
||||||
|
|
||||||
|
exchanger_gap_H = (exchanger_H - exchanger_N_tubes*exchanger_tube_H) / (exchanger_N_tubes - 1); % [m] Height between coolant tubes
|
||||||
|
exchanger_air_area_primary = 2 * (exchanger_N_tubes - 1) * exchanger_W * (exchanger_L + exchanger_gap_H); % [m^2] Primary air heat transfer surface area
|
||||||
|
exchanger_N_fins = (exchanger_N_tubes - 1) * exchanger_L / exchanger_fin_spacing; % [-] Total number of fins
|
||||||
|
exchanger_air_area_fins = 2 * exchanger_N_fins * exchanger_W * exchanger_gap_H; % [m^2] Total fin surface area
|
||||||
|
exchanger_tube_Leq = 2*(exchanger_H + 20*exchanger_tube_H*exchanger_N_tubes); % [m] Additional equivalent tube length for losses due to manifold and splits
|
||||||
|
|
||||||
|
% Hydrogen property tables
|
||||||
|
H2_R = 4124.48151675695; % [J/(kg*K)] Specific gas constant
|
||||||
|
H2_T = [-56.55, -50:10:-10, -5:1:5, 10:10:350]; % [degC] Temperature vector
|
||||||
|
H2_h = [2783.66044045879, 2873.94181645932, 3012.62601981068, 3152.22537301871, 3292.62486344326, 3433.72138359680, 3504.50184295996, 3518.67512520997, 3532.85394543712, 3547.03822323045, 3561.22787913325, 3575.42283463482, 3589.62301216223, 3603.82833507219, 3618.03872764290, 3632.25411506595, 3646.47442343830, 3717.64725530348, 3860.32198734890, 4003.38288163040, 4146.77354613486, 4290.44463593559, 4434.35318483605, 4578.46197843019, 4722.73896833090, 4867.15672726297, 5011.69194456718, 5156.32496143709, 5301.03934494342, 5445.82149962224, 5590.66031514150, 5735.54684833147, 5880.47403768083, 6025.43644826542, 6170.43004499140, 6315.45199199528, 6460.50047604533, 6605.57455182556, 6750.67400704912, 6895.79924543541, 7040.95118568927, 7186.13117473630, 7331.34091359044, 7476.58239435534, 7621.85784698693, 7767.16969456766, 7912.52051596217, 8057.91301483815, 8203.34999414341, 8348.83433523064, 8494.36898091475, 8639.95692183309]; % [kJ/kg] Specific enthalpy vector
|
||||||
|
H2_mu = [7.13280985169147, 7.28192581803172, 7.50683010986396, 7.72860158311813, 7.94743907374985, 8.16352107178216, 8.27057966039486, 8.29191522386247, 8.31322574339419, 8.33451136221049, 8.35577222214406, 8.37700846365822, 8.39822022586488, 8.41940764654233, 8.44057086215268, 8.46171000785900, 8.48282521754219, 8.58804682376172, 8.79676834158601, 9.00329345223078, 9.20773222364170, 9.41018554308578, 9.61074613691738, 9.80949945111700, 10.0065244149257, 10.2018941058184, 10.3956763308132, 10.5879341365094, 10.7787262581519, 10.9681075163184, 11.1561291684380, 11.3428392212160, 11.5282827091027, 11.7125019431712, 11.8955367341236, 12.0774245926121, 12.2582009096067, 12.4378991191677, 12.6165508456572, 12.7941860371557, 12.9708330866186, 13.1465189421111, 13.3212692072909, 13.4951082331655, 13.6680592020247, 13.8401442043429, 14.0113843093520, 14.1817996299056, 14.3514093821846, 14.5202319407340, 14.6882848892642, 14.8555850676093]; % [uPa*s] Dynamic viscosity vector
|
||||||
|
H2_k = [142.834673487453, 146.602208511061, 152.236222042957, 157.732550941948, 163.097652939664, 168.338595936873, 170.914793442425, 171.426608923333, 171.937300226766, 172.446874689218, 172.955339634979, 173.462702374715, 173.968970204106, 174.474150402550, 174.978250231913, 175.481276935350, 175.983237736166, 178.477304495221, 183.389581255575, 188.206454189023, 192.934525090876, 197.580042887656, 202.148890500710, 206.646585109232, 211.078287130673, 215.448814629342, 219.762660868078, 224.024013435935, 228.236773894856, 232.404577248187, 236.530810785993, 240.618632037420, 244.670985680918, 248.690619345047, 252.680098287219, 256.641818973249, 260.578021602720, 264.490801638200, 268.382120402950, 272.253814814231, 276.107606318911, 279.945109095949, 283.767837587053, 287.577213413001, 291.374571728986, 295.161167068166, 298.938178718576, 302.706715674597, 306.467821200562, 310.222477040649, 313.971607306067, 317.716082067666]; % [mW/(m*K)] Thermal conductivity vector
|
||||||
|
H2_D = 74; % [mm^2/s] Diffusivity
|
||||||
|
|
||||||
|
% Oxygen property tables
|
||||||
|
O2_R = 259.836612622973; % [J/(kg*K)] Specific gas constant
|
||||||
|
O2_T = [-56.55, -50:10:-10, -5:1:5, 10:10:350]; % [degC] Temperature vector
|
||||||
|
O2_h = [196.314045635230, 202.302438242727, 211.445560194021, 220.590899228792, 229.740231593225, 238.895348952500, 243.475641264001, 244.391947430149, 245.308340402054, 246.224821996439, 247.141394030059, 248.058058319568, 248.974816681386, 249.891670931563, 250.808622885637, 251.725674358497, 252.642827164235, 257.230174605453, 266.413508632715, 275.609852765421, 284.820965750972, 294.048557917965, 303.294277512928, 312.559698673577, 321.846311315623, 331.155513043735, 340.488603075706, 349.846778083808, 359.231129801281, 368.642644208585, 378.082202097877, 387.550580810768, 397.048456949970, 406.576409877154, 416.134925824813, 425.724402467528, 435.345153816409, 444.997415318718, 454.681349062211, 464.397049000008, 474.144546126762, 483.923813550211, 493.734771414086, 503.577291638585, 513.451202453549, 523.356292706982, 533.292315937896, 543.258994207727, 553.256021688858, 563.283068012223, 573.339781378746, 583.425791441447]; % [kJ/kg] Specific enthalpy vectdor
|
||||||
|
O2_mu = [15.7665575307288, 16.1727626078850, 16.7849249433608, 17.3877315055497, 17.9815180942032, 18.5666055809413, 18.8559834763191, 18.9136105724852, 18.9711555251620, 19.0286186248432, 19.0860001606775, 19.1433004204737, 19.2005196907073, 19.2576582565260, 19.3147164017561, 19.3716944089079, 19.4285925591820, 19.7118952084836, 20.2726692662067, 20.8258892364295, 21.3718096794610, 21.9106736606882, 22.4427133237112, 22.9681504449081, 23.4871969666892, 24.0000555077772, 24.5069198496509, 25.0079753988870, 25.5033996255713, 25.9933624782616, 26.4780267762054, 26.9575485796604, 27.4320775392600, 27.9017572254189, 28.3667254387956, 28.8271145028300, 29.2830515393565, 29.7346587282699, 30.1820535521826, 30.6253490269743, 31.0646539190931, 31.5000729504200, 31.9317069914666, 32.3596532436319, 32.7840054111974, 33.2048538637040, 33.6222857893065, 34.0363853396718, 34.4472337669430, 34.8549095532629, 35.2594885333161, 35.6610440103170]; % [uPa*s] Dynamic viscosity vector
|
||||||
|
O2_k = [19.6661951881093, 20.2224765128182, 21.0645207724504, 21.8980296960484, 22.7232456337983, 23.5404090717681, 23.9460452976483, 24.0269406552993, 24.1077592592045, 24.1885013406030, 24.2691671301751, 24.3497568580346, 24.4302707537208, 24.5107090461917, 24.5910719638169, 24.6713597343701, 24.7515725850235, 25.1515209478718, 25.9459275198455, 26.7331962712293, 27.5135401678269, 28.2871691203210, 29.0542797087419, 29.8150645640764, 30.5697094830221, 31.3183935589296, 32.0612893462030, 32.7985630478626, 33.5303747183955, 34.2568784758844, 34.9782227188418, 35.6945503442869, 36.4059989644644, 37.1127011202719, 37.8147844899821, 38.5123720922485, 39.2055824826964, 39.8945299436419, 40.5793246666704, 41.2600729279479, 41.9368772562485, 42.6098365937651, 43.2790464498297, 43.9445990477190, 44.6065834647486, 45.2650857658856, 45.9201891311213, 46.5719739768538, 47.2205180715361, 47.8658966458414, 48.5081824975986, 49.1474460917422]; % [mW/(m*K)] Thermal conductivity vector
|
||||||
|
O2_D = 18; % [mm^2/s] Diffusivity
|
||||||
|
|
||||||
|
% Nitrogen property tables
|
||||||
|
N2_R = 296.802103844292; % [J/(kg*K)] Specific gas constant
|
||||||
|
N2_T = [-56.55, -50:10:-10, -5:1:5, 10:10:350]; % [degC] Temperature vector
|
||||||
|
N2_h = [224.312019475783, 231.141016528663, 241.564116974996, 251.984329606098, 262.402210767198, 272.818266126753, 278.025759131002, 279.067222480148, 280.108675169779, 281.150117667373, 282.191550440608, 283.232973957487, 284.274388686473, 285.315795096618, 286.357193657682, 287.398584840260, 288.439969115894, 293.646803474573, 304.060229225710, 314.473742203777, 324.887858161051, 335.303123493654, 345.720119003427, 356.139461810715, 366.561805672305, 376.987839945795, 387.418287430292, 397.853901301030, 408.295461341278, 418.743769658649, 429.199646054884, 439.663923199077, 450.137441734658, 460.621045430963, 471.115576471389, 481.621870952319, 492.140754650669, 502.673039103050, 513.219518026408, 523.780964098585, 534.358126107445, 544.951726469054, 555.562459108678, 566.190987693008, 576.837944197865, 587.503927792527, 598.189504019663, 608.895204248440, 619.621525377624, 630.368929765337, 641.137845362314, 651.928666026144]; % [kJ/kg] Specific enthalpy vector
|
||||||
|
N2_mu = [13.7948920665312, 14.1367187351510, 14.6513506072199, 15.1575394445859, 15.6556177860831, 16.1459023391225, 16.3882166949119, 16.4364580271216, 16.4846261830610, 16.5327214469937, 16.5807441017456, 16.6286944287114, 16.6765727078615, 16.7243792177482, 16.7721142355126, 16.8197780368914, 16.8673708962236, 17.1042805670559, 17.5729330929829, 18.0349108462106, 18.4904598524977, 18.9398140034241, 19.3831957199835, 19.8208165930053, 20.2528779961256, 20.6795716689317, 21.1010802692360, 21.5175778943589, 21.9292305719283, 22.3361967211248, 22.7386275855675, 23.1366676391963, 23.5304549665897, 23.9201216191847, 24.3057939488592, 24.6875929202991, 25.0656344035241, 25.4400294478834, 25.8108845387614, 26.1783018381670, 26.5423794103038, 26.9032114331518, 27.2608883970194, 27.6154972909606, 27.9671217778900, 28.3158423591670, 28.6617365293689, 29.0048789219164, 29.3453414461711, 29.6831934165754, 30.0185016743671, 30.3513307023604]; % [uPa*s] Dynamic viscosity vector
|
||||||
|
N2_k = [19.6296532607374, 20.1533289725195, 20.9436243615696, 21.7231523445233, 22.4923063139090, 23.2514823266989, 23.6274403302585, 23.7023472871786, 23.7771601872884, 23.8518793809024, 23.9265052168265, 24.0010380423605, 24.0754782033001, 24.1498260439393, 24.2240819070729, 24.2982461339996, 24.3723190645241, 24.7413260601341, 25.4726839899413, 26.1954347636057, 26.9098870600572, 27.6163358040134, 28.3150627088262, 29.0063368423447, 29.6904152006066, 30.3675432785589, 31.0379556302294, 31.7018764131555, 32.3595199136277, 33.0110910505971, 33.6567858570436, 34.2967919382928, 34.9312889072676, 35.5604487970156, 36.1844364511046, 36.8034098926428, 37.4175206727930, 38.0269141997105, 38.6317300488700, 39.2321022557508, 39.8281595918420, 40.4200258249072, 41.0078199644172, 41.5916564930240, 42.1716455849090, 42.7478933117983, 43.3205018373934, 43.8895696009277, 44.4551914905138, 45.0174590069098, 45.5764604182929, 46.1322809065919]; % [mW/(m*K)] Thermal conductivity vector
|
||||||
|
|
@ -0,0 +1,50 @@
|
||||||
|
% Code to plot simulation results from PEMElectrolysisSystem
|
||||||
|
%% Plot Description:
|
||||||
|
%
|
||||||
|
% This plot shows the current-voltage (i-v) curve and the power consumed by
|
||||||
|
% a cell in the stack. As the current ramps up, an initial rise in voltage
|
||||||
|
% occurs due to electrode activation losses, followed by a gradual increase
|
||||||
|
% in voltage due to Ohmic resistances. The cell voltage is about 1.71 V at
|
||||||
|
% a current density of 2 A/cm^2.
|
||||||
|
|
||||||
|
% Copyright 2021 The MathWorks, Inc.
|
||||||
|
|
||||||
|
% Generate simulation results if they don't exist or are outdated
|
||||||
|
if ~exist('simlog_PEMElectrolysisSystem', 'var') || ...
|
||||||
|
(get_param('PEMElectrolysisSystem', 'RTWModifiedTimeStamp') ~= ...
|
||||||
|
simscape.logging.timestamp(simlog_PEMElectrolysisSystem))
|
||||||
|
sim('PEMElectrolysisSystem')
|
||||||
|
end
|
||||||
|
|
||||||
|
% Reuse figure if it exists, else create new figure
|
||||||
|
if ~exist('h1_PEMElectrolysisSystem', 'var') || ...
|
||||||
|
~isgraphics(h1_PEMElectrolysisSystem, 'figure')
|
||||||
|
h1_PEMElectrolysisSystem = figure('Name', 'PEMElectrolysisSystem');
|
||||||
|
end
|
||||||
|
figure(h1_PEMElectrolysisSystem)
|
||||||
|
clf(h1_PEMElectrolysisSystem)
|
||||||
|
|
||||||
|
plotIV(simlog_PEMElectrolysisSystem)
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
% Plot electrolyzer i-v curve
|
||||||
|
function plotIV(simlog)
|
||||||
|
|
||||||
|
% Get simulation results
|
||||||
|
i_cell = simlog.Membrane_Electrode_Assembly.i_cell.series.values('A/cm^2');
|
||||||
|
v_cell = simlog.Membrane_Electrode_Assembly.v_cell.series.values('V');
|
||||||
|
|
||||||
|
% Plot results
|
||||||
|
yyaxis left
|
||||||
|
plot(i_cell, v_cell, 'LineWidth', 1)
|
||||||
|
grid on
|
||||||
|
title('Fuel Cell I-V Curve')
|
||||||
|
ylabel('Cell Voltage (V)')
|
||||||
|
yyaxis right
|
||||||
|
plot(i_cell, i_cell.*v_cell, 'LineWidth', 1)
|
||||||
|
ylabel('Power Density (W/cm^2)')
|
||||||
|
xlabel('Current Density (A/cm^2)')
|
||||||
|
set(gca, 'LineWidth', 1)
|
||||||
|
|
||||||
|
end
|
||||||
|
|
@ -0,0 +1,63 @@
|
||||||
|
% Code to plot simulation results from PEMElectrolysisSystem
|
||||||
|
%% Plot Description:
|
||||||
|
%
|
||||||
|
% This plot shows the electrical power consumed by the electrolyzer. The
|
||||||
|
% electrical power is greater than the power needed to produce hydrogen due
|
||||||
|
% to various losses. The difference is the heat dissipated.
|
||||||
|
%
|
||||||
|
% This plot also shows the thermal efficiency of the electrolyzer, which
|
||||||
|
% indicates the fraction of the electrical power used to generate hydrogen
|
||||||
|
% based on hydrogen's heating value. This electrolyzer is about 87%
|
||||||
|
% efficient at a current density of 2 A/cm^2.
|
||||||
|
|
||||||
|
% Copyright 2021 The MathWorks, Inc.
|
||||||
|
|
||||||
|
% Generate simulation results if they don't exist or are outdated
|
||||||
|
if ~exist('simlog_PEMElectrolysisSystem', 'var') || ...
|
||||||
|
(get_param('PEMElectrolysisSystem', 'RTWModifiedTimeStamp') ~= ...
|
||||||
|
simscape.logging.timestamp(simlog_PEMElectrolysisSystem))
|
||||||
|
sim('PEMElectrolysisSystem')
|
||||||
|
end
|
||||||
|
|
||||||
|
% Reuse figure if it exists, else create new figure
|
||||||
|
if ~exist('h2_PEMElectrolysisSystem', 'var') || ...
|
||||||
|
~isgraphics(h2_PEMElectrolysisSystem, 'figure')
|
||||||
|
h2_PEMElectrolysisSystem = figure('Name', 'PEMElectrolysisSystem');
|
||||||
|
end
|
||||||
|
figure(h2_PEMElectrolysisSystem)
|
||||||
|
clf(h2_PEMElectrolysisSystem)
|
||||||
|
|
||||||
|
plotPower(simlog_PEMElectrolysisSystem)
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
% Plot power and efficiency
|
||||||
|
function plotPower(simlog)
|
||||||
|
|
||||||
|
% Get simulation results
|
||||||
|
t = simlog.Membrane_Electrode_Assembly.power_elec.series.time;
|
||||||
|
power_elec = simlog.Membrane_Electrode_Assembly.power_elec.series.values('kW');
|
||||||
|
power_dissipated = simlog.Membrane_Electrode_Assembly.power_dissipated.series.values('kW');
|
||||||
|
eta_HHV = simlog.Membrane_Electrode_Assembly.efficiency_HHV.series.values('1');
|
||||||
|
eta_LHV = simlog.Membrane_Electrode_Assembly.efficiency_LHV.series.values('1');
|
||||||
|
|
||||||
|
% Plot results
|
||||||
|
handles(1) = subplot(2, 1, 1);
|
||||||
|
plot(t, power_elec, 'LineWidth', 1)
|
||||||
|
hold on
|
||||||
|
plot(t, power_dissipated, 'LineWidth', 1)
|
||||||
|
grid on
|
||||||
|
legend('Electrical Supply', 'Heat Dissipated', 'Location', 'best')
|
||||||
|
title('Power (kW)')
|
||||||
|
handles(2) = subplot(2, 1, 2);
|
||||||
|
plot(t, eta_HHV*100, 'LineWidth', 1)
|
||||||
|
hold on
|
||||||
|
plot(t, eta_LHV*100, 'LineWidth', 1)
|
||||||
|
grid on
|
||||||
|
legend('Based on HHV', 'Based on LHV', 'Location', 'best')
|
||||||
|
title('Thermal Efficiency (%)')
|
||||||
|
xlabel('Time (s)')
|
||||||
|
|
||||||
|
linkaxes(handles, 'x')
|
||||||
|
|
||||||
|
end
|
||||||
|
|
@ -0,0 +1,75 @@
|
||||||
|
% Code to plot simulation results from PEMElectrolysisSystem
|
||||||
|
%% Plot Description:
|
||||||
|
%
|
||||||
|
% This plot shows the rate of hydrogen produced, the rate of water consumed
|
||||||
|
% at the anode, as well as the rate of water transported to the cathode due
|
||||||
|
% to diffusion, electro-osmosis drag, and hydraulic pressure difference. As
|
||||||
|
% a result, a dehumidification step is necessary to produce hydrogen at the
|
||||||
|
% desired purity.
|
||||||
|
%
|
||||||
|
% This plot also shows the total mass of hydrogen produced and the
|
||||||
|
% equivalent energy based on its higher heating value. This provides an
|
||||||
|
% indication of the amount of energy available if the hydrogen is used to
|
||||||
|
% generate power in a fuel cell.
|
||||||
|
|
||||||
|
% Copyright 2021 The MathWorks, Inc.
|
||||||
|
|
||||||
|
% Generate simulation results if they don't exist or are outdated
|
||||||
|
if ~exist('simlog_PEMElectrolysisSystem', 'var') || ...
|
||||||
|
(get_param('PEMElectrolysisSystem', 'RTWModifiedTimeStamp') ~= ...
|
||||||
|
simscape.logging.timestamp(simlog_PEMElectrolysisSystem))
|
||||||
|
sim('PEMElectrolysisSystem')
|
||||||
|
end
|
||||||
|
|
||||||
|
% Reuse figure if it exists, else create new figure
|
||||||
|
if ~exist('h3_PEMElectrolysisSystem', 'var') || ...
|
||||||
|
~isgraphics(h3_PEMElectrolysisSystem, 'figure')
|
||||||
|
h3_PEMElectrolysisSystem = figure('Name', 'PEMElectrolysisSystem');
|
||||||
|
end
|
||||||
|
figure(h3_PEMElectrolysisSystem)
|
||||||
|
clf(h3_PEMElectrolysisSystem)
|
||||||
|
|
||||||
|
plotHydrogen(simlog_PEMElectrolysisSystem)
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
% Plot hydrogen production
|
||||||
|
function plotHydrogen(simlog)
|
||||||
|
|
||||||
|
% Get simulation results
|
||||||
|
t = simlog.Membrane_Electrode_Assembly.H2_produced.series.time;
|
||||||
|
H2_produced = simlog.Membrane_Electrode_Assembly.H2_produced.series.values('g/s');
|
||||||
|
H2O_consumed = simlog.Membrane_Electrode_Assembly.H2O_consumed.series.values('g/s');
|
||||||
|
H2O_transport = simlog.Membrane_Electrode_Assembly.H2O_transport.series.values('g/s');
|
||||||
|
power_elec = simlog.Membrane_Electrode_Assembly.power_elec.series.values('kW');
|
||||||
|
power_dissipated = simlog.Membrane_Electrode_Assembly.power_dissipated.series.values('kW');
|
||||||
|
|
||||||
|
M_H2 = cumtrapz(t, H2_produced) * 1e-3;
|
||||||
|
E_H2 = cumtrapz(t, power_elec - power_dissipated) / 3600;
|
||||||
|
|
||||||
|
% Plot results
|
||||||
|
handles(1) = subplot(2, 1, 1);
|
||||||
|
plot(t, H2_produced, 'LineWidth', 1)
|
||||||
|
hold on
|
||||||
|
plot(t, H2O_consumed, 'LineWidth', 1)
|
||||||
|
plot(t, H2O_transport, 'LineWidth', 1)
|
||||||
|
grid on
|
||||||
|
legend('H2 Produced', 'H2O Consumed', 'H2O Transported', 'Location', 'best')
|
||||||
|
title('Hydrogen and Water in MEA')
|
||||||
|
ylabel('Mass Flow Rate (g/s)')
|
||||||
|
handles(2) = subplot(2, 1, 2);
|
||||||
|
yyaxis left
|
||||||
|
plot(t, M_H2, 'LineWidth', 1)
|
||||||
|
grid on
|
||||||
|
title('Hydrogen Produced')
|
||||||
|
ylabel('Mass (kg)')
|
||||||
|
yyaxis right
|
||||||
|
plot(t, E_H2, 'LineWidth', 1)
|
||||||
|
ylim(max(ylim, 0))
|
||||||
|
ylabel('Energy (kWh)')
|
||||||
|
xlabel('Time (s)')
|
||||||
|
set(gca, 'LineWidth', 1)
|
||||||
|
|
||||||
|
linkaxes(handles, 'x')
|
||||||
|
|
||||||
|
end
|
||||||
Binary file not shown.
Binary file not shown.
Binary file not shown.
|
|
@ -0,0 +1,219 @@
|
||||||
|
%% PEM Fuel Cell System
|
||||||
|
%
|
||||||
|
% This example shows how to model a proton exchange membrane (PEM) fuel
|
||||||
|
% cell stack with a custom Simscape(TM) block. The PEM fuel cell generates
|
||||||
|
% electrical power by consuming hydrogen and oxygen and producing water
|
||||||
|
% vapor. The custom block represents the membrane electrode assembly (MEA)
|
||||||
|
% and is connected two separate moist air networks: one for the anode gas
|
||||||
|
% flow and one for the cathode gas flow.
|
||||||
|
%
|
||||||
|
% The two moist air networks represent different gas mixtures. The anode
|
||||||
|
% network consists of nitrogen (N2), water vapor (H2O), and hydrogen (H2),
|
||||||
|
% representing the fuel. The hydrogen is stored in the fuel tank at 70 MPa.
|
||||||
|
% A pressure-reducing valve releases hydrogen to the fuel cell stack at
|
||||||
|
% around 0.16 MPa. Unconsumed hydrogen is recirculated back to the stack.
|
||||||
|
% The cathode network consists of nitrogen (N2), water vapor (H2O), and
|
||||||
|
% oxygen (O2), representing air from the environment. A compressor brings
|
||||||
|
% air to the fuel cell stack at a controlled rate to ensure that the fuel
|
||||||
|
% cell is not starved of oxygen. A back pressure relief valve maintains a
|
||||||
|
% pressure of around 0.16 MPa in the stack and vents the exhaust to the
|
||||||
|
% environment.
|
||||||
|
%
|
||||||
|
% The temperature and relative humidity in the fuel cell stack must be
|
||||||
|
% maintained at an optimal level to ensure efficient operation under
|
||||||
|
% various loading conditions. Higher temperatures increase thermal
|
||||||
|
% efficiency but reduce relative humidity, which causes higher membrane
|
||||||
|
% resistance. Therefore, in this model, the fuel cell stack temperature is
|
||||||
|
% kept at 80 degC. The cooling system circulates coolant between the cells
|
||||||
|
% to absorb heat and rejects it to the environment via the radiator. The
|
||||||
|
% humidifiers saturate the gas with water vapor to keep the membrane
|
||||||
|
% hydrated and minimize electrical resistance.
|
||||||
|
%
|
||||||
|
% The custom MEA block is implemented in the Simscape code
|
||||||
|
% |<matlab:edit('FuelCell.ssc') FuelCell.ssc>|. The output port F of the
|
||||||
|
% anode and cathode gas channel pipe blocks provide the gas mole fractions
|
||||||
|
% needed to model the fuel cell reaction. The removal of H2 and O2 from the
|
||||||
|
% anode and cathode gas flows are implemented by Controlled Trace Gas
|
||||||
|
% Source (MA) blocks. The production of H2O and the transport of water
|
||||||
|
% vapor across the MEA are implemented by Controlled Moisture Source (MA)
|
||||||
|
% blocks. The heat generated by the reaction is sent through the thermal
|
||||||
|
% port H to the connected Thermal Mass block. Refer to the comments in the
|
||||||
|
% code for additional details on the implementation.
|
||||||
|
%
|
||||||
|
% See also the <docid:simscape_ug#example-ssc_electrolyzer PEM Electrolysis
|
||||||
|
% System System> example.
|
||||||
|
%
|
||||||
|
% References:
|
||||||
|
%
|
||||||
|
% Dutta, Sandip, Sirivatch Shimpalee, and J. W. Van Zee. "Numerical
|
||||||
|
% prediction of mass-exchange between cathode and anode channels in a PEM
|
||||||
|
% fuel cell." _International Journal of Heat and Mass Transfer_ 44.11
|
||||||
|
% (2001): 2029-2042.
|
||||||
|
%
|
||||||
|
% EG&G Technical Services, Inc. _Fuel Cell Handbook (Seventh Edition)._ US
|
||||||
|
% Department of Energy, Office of Fossil Energy, National Energy Technology
|
||||||
|
% Laboratory, 2004.
|
||||||
|
%
|
||||||
|
% Pukrushpan, Jay T., Anna G. Stefanopoulou, and Huei Peng. _Control of
|
||||||
|
% fuel cell power systems: principles, modeling, analysis and feedback
|
||||||
|
% design._ Springer-Verlag London, 2004.
|
||||||
|
%
|
||||||
|
% Spiegel, Colleen. _PEM fuel cell modeling and simulation using MATLAB._
|
||||||
|
% Elsevier, 2008.
|
||||||
|
|
||||||
|
% Copyright 2020-2021 The MathWorks, Inc.
|
||||||
|
|
||||||
|
|
||||||
|
%% Model
|
||||||
|
|
||||||
|
open_system('PEMFuelCellSystem')
|
||||||
|
|
||||||
|
set_param(find_system(bdroot,'FindAll','on','type','annotation','Tag','ModelFeatures'),'Interpreter','off');
|
||||||
|
|
||||||
|
%% Anode Humidifier Subsystem
|
||||||
|
|
||||||
|
open_system('PEMFuelCellSystem/Anode Humidifier','force')
|
||||||
|
|
||||||
|
%% Anode Exhaust Subsystem
|
||||||
|
|
||||||
|
open_system('PEMFuelCellSystem/Anode Exhaust','force')
|
||||||
|
|
||||||
|
%% Anode Gas Channels Subsystem
|
||||||
|
|
||||||
|
open_system('PEMFuelCellSystem/Anode Gas Channels','force')
|
||||||
|
|
||||||
|
%% Cathode Humidifier Subsystem
|
||||||
|
|
||||||
|
open_system('PEMFuelCellSystem/Cathode Humidifier','force')
|
||||||
|
|
||||||
|
%% Cathode Exhaust Subsystem
|
||||||
|
|
||||||
|
open_system('PEMFuelCellSystem/Cathode Exhaust','force')
|
||||||
|
|
||||||
|
%% Pressure Relief Valve Subsystem
|
||||||
|
|
||||||
|
open_system('PEMFuelCellSystem/Cathode Exhaust/Pressure Relief Valve','force')
|
||||||
|
|
||||||
|
%% Cathode Gas Channels Subsystem
|
||||||
|
|
||||||
|
open_system('PEMFuelCellSystem/Cathode Gas Channels','force')
|
||||||
|
|
||||||
|
%% Cooling System Subsystem
|
||||||
|
|
||||||
|
open_system('PEMFuelCellSystem/Cooling System','force')
|
||||||
|
|
||||||
|
%% Coolant Tank Subsystem
|
||||||
|
|
||||||
|
open_system('PEMFuelCellSystem/Cooling System/Coolant Tank','force')
|
||||||
|
|
||||||
|
%% Electrical Load Subsystem
|
||||||
|
|
||||||
|
open_system('PEMFuelCellSystem/Electrical Load','force')
|
||||||
|
|
||||||
|
%% Hydrogen Source Subsystem
|
||||||
|
|
||||||
|
open_system('PEMFuelCellSystem/Hydrogen Source','force')
|
||||||
|
|
||||||
|
%% Pressure-Reducing Valve Subsystem
|
||||||
|
|
||||||
|
open_system('PEMFuelCellSystem/Hydrogen Source/Pressure-Reducing Valve','force')
|
||||||
|
|
||||||
|
%% Oxygen Source Subsystem
|
||||||
|
|
||||||
|
open_system('PEMFuelCellSystem/Oxygen Source','force')
|
||||||
|
|
||||||
|
%% Recirculation Subsystem
|
||||||
|
|
||||||
|
open_system('PEMFuelCellSystem/Recirculation','force')
|
||||||
|
|
||||||
|
%% Simulation Results from Scopes
|
||||||
|
|
||||||
|
set_param('PEMFuelCellSystem/Scope','open','on');
|
||||||
|
sim('PEMFuelCellSystem');
|
||||||
|
|
||||||
|
%%
|
||||||
|
|
||||||
|
set_param('PEMFuelCellSystem/Scope','open','off');
|
||||||
|
%% Simulation Results from Simscape Logging
|
||||||
|
%%
|
||||||
|
%
|
||||||
|
% This plot shows the current-voltage (i-v) curve of a fuel cell in the
|
||||||
|
% stack. As the current ramps up, an initial drop in voltage occurs due to
|
||||||
|
% electrode activation losses, followed by a gradual decrease in voltage
|
||||||
|
% due to Ohmic resistances. Near maximum current, a sharp drop in voltage
|
||||||
|
% occurs due to gas-transport-related losses.
|
||||||
|
%
|
||||||
|
% This plot also shows the power produced by the cell. When the ramp
|
||||||
|
% scenario is selected, the power increases until a maximum power output,
|
||||||
|
% then decreases due to the high losses near maximum current.
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
PEMFuelCellSystemPlot1IV;
|
||||||
|
%%
|
||||||
|
%
|
||||||
|
% This plot shows the electrical power produced by the fuel cell stack as
|
||||||
|
% well as the power consumed by the cathode air compressor and the coolant
|
||||||
|
% pump to maintain stable and efficient system operation. As a result, the
|
||||||
|
% net power produced by the system is a few percent less than the power
|
||||||
|
% produced by the stack. Note that this model assumes an isentropic
|
||||||
|
% compressor. Accounting for compressor efficiency would decrease net power
|
||||||
|
% gain by another couple of percent.
|
||||||
|
%
|
||||||
|
% This plot also shows the excess heat generated by the fuel cell stack,
|
||||||
|
% which must be removed by the cooling system. The maximum power produced
|
||||||
|
% by the fuel cell stack is 110 kW.
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
PEMFuelCellSystemPlot2Power;
|
||||||
|
%%
|
||||||
|
%
|
||||||
|
% This plot show the thermal efficiency of the fuel cell and its reactant
|
||||||
|
% utilization fraction. The thermal efficiency indicates the fraction of
|
||||||
|
% the hydrogen fuel's energy that the fuel cell has converted to useful
|
||||||
|
% electrical work. The theoretical maximum efficiency for a PEM fuel cell
|
||||||
|
% is 83%. However, actual efficiency is around 60% due to internal losses.
|
||||||
|
% Near maximum current, the efficiency drops to around 45%.
|
||||||
|
%
|
||||||
|
% The reactant utilization is the fraction of the reactants, H2 and O2,
|
||||||
|
% flowing into the fuel cell stack that has been consumed by the fuel cell.
|
||||||
|
% While higher utilization makes better use of the gases flowing through
|
||||||
|
% the fuel cell, it decreases the concentration of the reactants and thus
|
||||||
|
% reduces the voltage produced. Unconsumed O2 is vented to the environment,
|
||||||
|
% but unconsumed H2 is recirculated back to the anode to avoid waste.
|
||||||
|
% However, in practice, the H2 is periodically purged to remove
|
||||||
|
% contaminants.
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
PEMFuelCellSystemPlot3Efficiency;
|
||||||
|
%%
|
||||||
|
%
|
||||||
|
% This plot shows the temperatures at various locations in the system. The
|
||||||
|
% fuel cell stack temperature is maintained at a maximum of 80 degC by the
|
||||||
|
% cooling system. Fuel flowing to the anode is warmed by the recirculated
|
||||||
|
% flow. Air flowing to the cathode is warmed by the compressor.
|
||||||
|
%
|
||||||
|
% Maintaining an optimal temperature is critical to the operation of the
|
||||||
|
% fuel cell because higher temperatures lower the relative humidity which
|
||||||
|
% increases the membrane resistance. In this model, the cooling system is
|
||||||
|
% operated by a simple control of the coolant pump flow rate. The plot
|
||||||
|
% shows the temperature of the coolant after it has absorbed heat from the
|
||||||
|
% fuel cell stack and after it has rejected heat in the radiator.
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
PEMFuelCellSystemPlot4T;
|
||||||
|
%%
|
||||||
|
%
|
||||||
|
% This plot shows the mass of hydrogen used during operation and the
|
||||||
|
% corresponding decrease in the hydrogen tank pressure. The energy of the
|
||||||
|
% consumed hydrogen fuel is converted to electrical energy.
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
PEMFuelCellSystemPlot5Energy;
|
||||||
|
|
||||||
|
%%
|
||||||
|
|
||||||
|
|
@ -0,0 +1,95 @@
|
||||||
|
%% Code to define parameters for PEMFuelCellSystem
|
||||||
|
% Open Model Workspace in the Model Explorer to view and modify parameter
|
||||||
|
% values. Click 'Reinitialize from Source' to reset to the parameter values
|
||||||
|
% in this script.
|
||||||
|
|
||||||
|
% Copyright 2020 The MathWorks, Inc.
|
||||||
|
|
||||||
|
load PEMFuelCellSystemDriveCycle.mat
|
||||||
|
|
||||||
|
% Environment conditions
|
||||||
|
env_p = 0.101325; % [MPa] Pressure
|
||||||
|
env_T = 20; % [degC] Temperature
|
||||||
|
env_RH = 0.5; % [-] Relative humidity
|
||||||
|
env_yO2 = 0.21; % [-] Oxygen mole fraction
|
||||||
|
|
||||||
|
% Hydrogen fuel
|
||||||
|
tank_p = 70; % [MPa] Fuel tank pressure
|
||||||
|
tank_yH2 = 1 - 3e-4; % [-] Hydrogen mole fraction
|
||||||
|
tank_V = 120; % [l] Fuel tank volume
|
||||||
|
|
||||||
|
% Fuel cell stack
|
||||||
|
stack_num_cells = 400; % [-] Number cells
|
||||||
|
stack_area = 280; % [cm^2] Cell area
|
||||||
|
stack_t_membrane = 125; % [um] Membrane thickness
|
||||||
|
stack_t_gdl = 250; % [um] Gas diffusion layer thickness
|
||||||
|
stack_w_channels = 1; % [cm] Gas channel width/height
|
||||||
|
stack_num_channels = 8; % [-] Number of gas channels per cell
|
||||||
|
stack_io = 1e-04; % [A/cm^2] Exchange current density
|
||||||
|
stack_iL = 1.4; % [A/cm^2] Limiting current density
|
||||||
|
stack_alpha = 0.7; % [-] Charge transfer coefficient
|
||||||
|
stack_D_gdl = 0.07; % [cm^2/s] Water diffusivity in GDL
|
||||||
|
stack_membrane_rho = 2000; % [kg/m^3] Density of dry membrane
|
||||||
|
stack_membrane_MW = 1.1; % [kg/mol] Equivalent weight of dry membrane
|
||||||
|
stack_mea_rho = 1800; % [kg/s] Overall density of membrane electrode assembly
|
||||||
|
stack_mea_cp = 870; % [J/(kg*K)] Overall specific heat of membrane electrode assembly
|
||||||
|
|
||||||
|
anode_tube_D = 0.01; % [m] Hydrogen tube diameter
|
||||||
|
cathode_tube_D = 0.05; % [m] Air tube diameter
|
||||||
|
|
||||||
|
% Coolant system
|
||||||
|
coolant_w_channels = 1; % [cm] Coolant channel width/height
|
||||||
|
coolant_num_passes = 12; % [-] Number of coolant channel passes per layer
|
||||||
|
coolant_num_layers = 20; % [-] Number of coolant layers in stack
|
||||||
|
coolant_tube_D = 0.05; % [m] Coolant tube diameter
|
||||||
|
|
||||||
|
% Radiator dimensions
|
||||||
|
radiator_L = 1; % [m] Overall radiator length
|
||||||
|
radiator_W = 0.025; % [m] Overall radiator width
|
||||||
|
radiator_H = 0.5; % [m] Overal radiator height
|
||||||
|
radiator_N_tubes = 25; % [-] Number of coolant tubes
|
||||||
|
radiator_tube_H = 0.0015; % [m] Height of each coolant tube
|
||||||
|
radiator_fin_spacing = 0.002; % [-] Fin spacing
|
||||||
|
radiator_eta_fin = 0.7; % [-] Fin efficiency
|
||||||
|
radiator_t_wall = 1e-4; % [m] Material thickness
|
||||||
|
radiator_rho = 2700; % [kg/s] Radiator material density
|
||||||
|
radiator_cp = 910; % [J/(kg*K)] Radiator material specific heat
|
||||||
|
|
||||||
|
radiator_gap_H = (radiator_H - radiator_N_tubes*radiator_tube_H) / (radiator_N_tubes - 1); % [m] Height between coolant tubes
|
||||||
|
radiator_air_area_primary = 2 * (radiator_N_tubes - 1) * radiator_W * (radiator_L + radiator_gap_H); % [m^2] Primary air heat transfer surface area
|
||||||
|
radiator_N_fins = (radiator_N_tubes - 1) * radiator_L / radiator_fin_spacing; % [-] Total number of fins
|
||||||
|
radiator_air_area_fins = 2 * radiator_N_fins * radiator_W * radiator_gap_H; % [m^2] Total fin surface area
|
||||||
|
radiator_tube_Leq = 2*(radiator_H + 20*radiator_tube_H*radiator_N_tubes); % [m] Additional equivalent tube length for losses due to manifold and splits
|
||||||
|
|
||||||
|
% Compressor map
|
||||||
|
comp_p_ratio_TLU = [1; 1.25; 1.5; 1.75; 2]; % [-] Pressure ratio vector
|
||||||
|
comp_rpm_TLU = [0, 1800, 3600]; % [rpm] Shaft speed vector
|
||||||
|
comp_mdot_corr_TLU = [
|
||||||
|
0, 0.05, 0.1;
|
||||||
|
0, 0.0375, 0.075;
|
||||||
|
0, 0.025, 0.05;
|
||||||
|
0, 0.0125, 0.025;
|
||||||
|
0, 0, 0] * 4; % [kg/s] Corrected mass flow rate table
|
||||||
|
|
||||||
|
% Hydrogen property tables
|
||||||
|
H2_R = 4124.48151675695; % [J/(kg*K)] Specific gas constant
|
||||||
|
H2_T = [-56.55, -50:10:-10, -5:1:5, 10:10:350]; % [degC] Temperature vector
|
||||||
|
H2_h = [2783.66044045879, 2873.94181645932, 3012.62601981068, 3152.22537301871, 3292.62486344326, 3433.72138359680, 3504.50184295996, 3518.67512520997, 3532.85394543712, 3547.03822323045, 3561.22787913325, 3575.42283463482, 3589.62301216223, 3603.82833507219, 3618.03872764290, 3632.25411506595, 3646.47442343830, 3717.64725530348, 3860.32198734890, 4003.38288163040, 4146.77354613486, 4290.44463593559, 4434.35318483605, 4578.46197843019, 4722.73896833090, 4867.15672726297, 5011.69194456718, 5156.32496143709, 5301.03934494342, 5445.82149962224, 5590.66031514150, 5735.54684833147, 5880.47403768083, 6025.43644826542, 6170.43004499140, 6315.45199199528, 6460.50047604533, 6605.57455182556, 6750.67400704912, 6895.79924543541, 7040.95118568927, 7186.13117473630, 7331.34091359044, 7476.58239435534, 7621.85784698693, 7767.16969456766, 7912.52051596217, 8057.91301483815, 8203.34999414341, 8348.83433523064, 8494.36898091475, 8639.95692183309]; % [kJ/kg] Specific enthalpy vector
|
||||||
|
H2_mu = [7.13280985169147, 7.28192581803172, 7.50683010986396, 7.72860158311813, 7.94743907374985, 8.16352107178216, 8.27057966039486, 8.29191522386247, 8.31322574339419, 8.33451136221049, 8.35577222214406, 8.37700846365822, 8.39822022586488, 8.41940764654233, 8.44057086215268, 8.46171000785900, 8.48282521754219, 8.58804682376172, 8.79676834158601, 9.00329345223078, 9.20773222364170, 9.41018554308578, 9.61074613691738, 9.80949945111700, 10.0065244149257, 10.2018941058184, 10.3956763308132, 10.5879341365094, 10.7787262581519, 10.9681075163184, 11.1561291684380, 11.3428392212160, 11.5282827091027, 11.7125019431712, 11.8955367341236, 12.0774245926121, 12.2582009096067, 12.4378991191677, 12.6165508456572, 12.7941860371557, 12.9708330866186, 13.1465189421111, 13.3212692072909, 13.4951082331655, 13.6680592020247, 13.8401442043429, 14.0113843093520, 14.1817996299056, 14.3514093821846, 14.5202319407340, 14.6882848892642, 14.8555850676093]; % [uPa*s] Dynamic viscosity vector
|
||||||
|
H2_k = [142.834673487453, 146.602208511061, 152.236222042957, 157.732550941948, 163.097652939664, 168.338595936873, 170.914793442425, 171.426608923333, 171.937300226766, 172.446874689218, 172.955339634979, 173.462702374715, 173.968970204106, 174.474150402550, 174.978250231913, 175.481276935350, 175.983237736166, 178.477304495221, 183.389581255575, 188.206454189023, 192.934525090876, 197.580042887656, 202.148890500710, 206.646585109232, 211.078287130673, 215.448814629342, 219.762660868078, 224.024013435935, 228.236773894856, 232.404577248187, 236.530810785993, 240.618632037420, 244.670985680918, 248.690619345047, 252.680098287219, 256.641818973249, 260.578021602720, 264.490801638200, 268.382120402950, 272.253814814231, 276.107606318911, 279.945109095949, 283.767837587053, 287.577213413001, 291.374571728986, 295.161167068166, 298.938178718576, 302.706715674597, 306.467821200562, 310.222477040649, 313.971607306067, 317.716082067666]; % [mW/(m*K)] Thermal conductivity vector
|
||||||
|
H2_D = 74; % [mm^2/s] Diffusivity
|
||||||
|
|
||||||
|
% Oxygen property tables
|
||||||
|
O2_R = 259.836612622973; % [J/(kg*K)] Specific gas constant
|
||||||
|
O2_T = [-56.55, -50:10:-10, -5:1:5, 10:10:350]; % [degC] Temperature vector
|
||||||
|
O2_h = [196.314045635230, 202.302438242727, 211.445560194021, 220.590899228792, 229.740231593225, 238.895348952500, 243.475641264001, 244.391947430149, 245.308340402054, 246.224821996439, 247.141394030059, 248.058058319568, 248.974816681386, 249.891670931563, 250.808622885637, 251.725674358497, 252.642827164235, 257.230174605453, 266.413508632715, 275.609852765421, 284.820965750972, 294.048557917965, 303.294277512928, 312.559698673577, 321.846311315623, 331.155513043735, 340.488603075706, 349.846778083808, 359.231129801281, 368.642644208585, 378.082202097877, 387.550580810768, 397.048456949970, 406.576409877154, 416.134925824813, 425.724402467528, 435.345153816409, 444.997415318718, 454.681349062211, 464.397049000008, 474.144546126762, 483.923813550211, 493.734771414086, 503.577291638585, 513.451202453549, 523.356292706982, 533.292315937896, 543.258994207727, 553.256021688858, 563.283068012223, 573.339781378746, 583.425791441447]; % [kJ/kg] Specific enthalpy vectdor
|
||||||
|
O2_mu = [15.7665575307288, 16.1727626078850, 16.7849249433608, 17.3877315055497, 17.9815180942032, 18.5666055809413, 18.8559834763191, 18.9136105724852, 18.9711555251620, 19.0286186248432, 19.0860001606775, 19.1433004204737, 19.2005196907073, 19.2576582565260, 19.3147164017561, 19.3716944089079, 19.4285925591820, 19.7118952084836, 20.2726692662067, 20.8258892364295, 21.3718096794610, 21.9106736606882, 22.4427133237112, 22.9681504449081, 23.4871969666892, 24.0000555077772, 24.5069198496509, 25.0079753988870, 25.5033996255713, 25.9933624782616, 26.4780267762054, 26.9575485796604, 27.4320775392600, 27.9017572254189, 28.3667254387956, 28.8271145028300, 29.2830515393565, 29.7346587282699, 30.1820535521826, 30.6253490269743, 31.0646539190931, 31.5000729504200, 31.9317069914666, 32.3596532436319, 32.7840054111974, 33.2048538637040, 33.6222857893065, 34.0363853396718, 34.4472337669430, 34.8549095532629, 35.2594885333161, 35.6610440103170]; % [uPa*s] Dynamic viscosity vector
|
||||||
|
O2_k = [19.6661951881093, 20.2224765128182, 21.0645207724504, 21.8980296960484, 22.7232456337983, 23.5404090717681, 23.9460452976483, 24.0269406552993, 24.1077592592045, 24.1885013406030, 24.2691671301751, 24.3497568580346, 24.4302707537208, 24.5107090461917, 24.5910719638169, 24.6713597343701, 24.7515725850235, 25.1515209478718, 25.9459275198455, 26.7331962712293, 27.5135401678269, 28.2871691203210, 29.0542797087419, 29.8150645640764, 30.5697094830221, 31.3183935589296, 32.0612893462030, 32.7985630478626, 33.5303747183955, 34.2568784758844, 34.9782227188418, 35.6945503442869, 36.4059989644644, 37.1127011202719, 37.8147844899821, 38.5123720922485, 39.2055824826964, 39.8945299436419, 40.5793246666704, 41.2600729279479, 41.9368772562485, 42.6098365937651, 43.2790464498297, 43.9445990477190, 44.6065834647486, 45.2650857658856, 45.9201891311213, 46.5719739768538, 47.2205180715361, 47.8658966458414, 48.5081824975986, 49.1474460917422]; % [mW/(m*K)] Thermal conductivity vector
|
||||||
|
O2_D = 18; % [mm^2/s] Diffusivity
|
||||||
|
|
||||||
|
% Nitrogen property tables
|
||||||
|
N2_R = 296.802103844292; % [J/(kg*K)] Specific gas constant
|
||||||
|
N2_T = [-56.55, -50:10:-10, -5:1:5, 10:10:350]; % [degC] Temperature vector
|
||||||
|
N2_h = [224.312019475783, 231.141016528663, 241.564116974996, 251.984329606098, 262.402210767198, 272.818266126753, 278.025759131002, 279.067222480148, 280.108675169779, 281.150117667373, 282.191550440608, 283.232973957487, 284.274388686473, 285.315795096618, 286.357193657682, 287.398584840260, 288.439969115894, 293.646803474573, 304.060229225710, 314.473742203777, 324.887858161051, 335.303123493654, 345.720119003427, 356.139461810715, 366.561805672305, 376.987839945795, 387.418287430292, 397.853901301030, 408.295461341278, 418.743769658649, 429.199646054884, 439.663923199077, 450.137441734658, 460.621045430963, 471.115576471389, 481.621870952319, 492.140754650669, 502.673039103050, 513.219518026408, 523.780964098585, 534.358126107445, 544.951726469054, 555.562459108678, 566.190987693008, 576.837944197865, 587.503927792527, 598.189504019663, 608.895204248440, 619.621525377624, 630.368929765337, 641.137845362314, 651.928666026144]; % [kJ/kg] Specific enthalpy vector
|
||||||
|
N2_mu = [13.7948920665312, 14.1367187351510, 14.6513506072199, 15.1575394445859, 15.6556177860831, 16.1459023391225, 16.3882166949119, 16.4364580271216, 16.4846261830610, 16.5327214469937, 16.5807441017456, 16.6286944287114, 16.6765727078615, 16.7243792177482, 16.7721142355126, 16.8197780368914, 16.8673708962236, 17.1042805670559, 17.5729330929829, 18.0349108462106, 18.4904598524977, 18.9398140034241, 19.3831957199835, 19.8208165930053, 20.2528779961256, 20.6795716689317, 21.1010802692360, 21.5175778943589, 21.9292305719283, 22.3361967211248, 22.7386275855675, 23.1366676391963, 23.5304549665897, 23.9201216191847, 24.3057939488592, 24.6875929202991, 25.0656344035241, 25.4400294478834, 25.8108845387614, 26.1783018381670, 26.5423794103038, 26.9032114331518, 27.2608883970194, 27.6154972909606, 27.9671217778900, 28.3158423591670, 28.6617365293689, 29.0048789219164, 29.3453414461711, 29.6831934165754, 30.0185016743671, 30.3513307023604]; % [uPa*s] Dynamic viscosity vector
|
||||||
|
N2_k = [19.6296532607374, 20.1533289725195, 20.9436243615696, 21.7231523445233, 22.4923063139090, 23.2514823266989, 23.6274403302585, 23.7023472871786, 23.7771601872884, 23.8518793809024, 23.9265052168265, 24.0010380423605, 24.0754782033001, 24.1498260439393, 24.2240819070729, 24.2982461339996, 24.3723190645241, 24.7413260601341, 25.4726839899413, 26.1954347636057, 26.9098870600572, 27.6163358040134, 28.3150627088262, 29.0063368423447, 29.6904152006066, 30.3675432785589, 31.0379556302294, 31.7018764131555, 32.3595199136277, 33.0110910505971, 33.6567858570436, 34.2967919382928, 34.9312889072676, 35.5604487970156, 36.1844364511046, 36.8034098926428, 37.4175206727930, 38.0269141997105, 38.6317300488700, 39.2321022557508, 39.8281595918420, 40.4200258249072, 41.0078199644172, 41.5916564930240, 42.1716455849090, 42.7478933117983, 43.3205018373934, 43.8895696009277, 44.4551914905138, 45.0174590069098, 45.5764604182929, 46.1322809065919]; % [mW/(m*K)] Thermal conductivity vector
|
||||||
|
|
@ -0,0 +1,52 @@
|
||||||
|
% Code to plot simulation results from PEMFuelCellSystem
|
||||||
|
%% Plot Description:
|
||||||
|
%
|
||||||
|
% This plot shows the current-voltage (i-v) curve of a fuel cell in the
|
||||||
|
% stack. As the current ramps up, an initial drop in voltage occurs due to
|
||||||
|
% electrode activation losses, followed by a gradual decrease in voltage
|
||||||
|
% due to Ohmic resistances. Near maximum current, a sharp drop in voltage
|
||||||
|
% occurs due to gas-transport-related losses.
|
||||||
|
%
|
||||||
|
% This plot also shows the power produced by the cell. When the ramp
|
||||||
|
% scenario is selected, the power increases until a maximum power output,
|
||||||
|
% then decreases due to the high losses near maximum current.
|
||||||
|
|
||||||
|
% Copyright 2020 The MathWorks, Inc.
|
||||||
|
|
||||||
|
% Generate simulation results if they don't exist
|
||||||
|
if ~exist('simlog_PEMFuelCellSystem', 'var')
|
||||||
|
sim('PEMFuelCellSystem')
|
||||||
|
end
|
||||||
|
|
||||||
|
% Reuse figure if it exists, else create new figure
|
||||||
|
if ~exist('h1_PEMFuelCellSystem', 'var') || ...
|
||||||
|
~isgraphics(h1_PEMFuelCellSystem, 'figure')
|
||||||
|
h1_PEMFuelCellSystem = figure('Name', 'PEMFuelCellSystem');
|
||||||
|
end
|
||||||
|
figure(h1_PEMFuelCellSystem)
|
||||||
|
clf(h1_PEMFuelCellSystem)
|
||||||
|
|
||||||
|
plotIV(simlog_PEMFuelCellSystem)
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
% Plot fuel cell i-v curve
|
||||||
|
function plotIV(simlog)
|
||||||
|
|
||||||
|
% Get simulation results
|
||||||
|
i_cell = simlog.Membrane_Electrode_Assembly.i_cell.series.values('A/cm^2');
|
||||||
|
v_cell = simlog.Membrane_Electrode_Assembly.v_cell.series.values('V');
|
||||||
|
|
||||||
|
% Plot results
|
||||||
|
yyaxis left
|
||||||
|
plot(i_cell, v_cell, 'LineWidth', 1)
|
||||||
|
grid on
|
||||||
|
title('Fuel Cell I-V Curve')
|
||||||
|
ylabel('Cell Voltage (V)')
|
||||||
|
yyaxis right
|
||||||
|
plot(i_cell, i_cell.*v_cell, 'LineWidth', 1)
|
||||||
|
ylabel('Power Density (W/cm^2)')
|
||||||
|
xlabel('Current Density (A/cm^2)')
|
||||||
|
set(gca, 'LineWidth', 1)
|
||||||
|
|
||||||
|
end
|
||||||
|
|
@ -0,0 +1,66 @@
|
||||||
|
% Code to plot simulation results from PEMFuelCellSystem
|
||||||
|
%% Plot Description:
|
||||||
|
%
|
||||||
|
% This plot shows the electrical power produced by the fuel cell stack as
|
||||||
|
% well as the power consumed by the cathode air compressor and the coolant
|
||||||
|
% pump to maintain stable and efficient system operation. As a result, the
|
||||||
|
% net power produced by the system is a few percent less than the power
|
||||||
|
% produced by the stack. Note that this model assumes an isentropic
|
||||||
|
% compressor. Accounting for compressor efficiency would decrease net power
|
||||||
|
% gain by another couple of percent.
|
||||||
|
%
|
||||||
|
% This plot also shows the excess heat generated by the fuel cell stack,
|
||||||
|
% which must be removed by the cooling system. The maximum power produced
|
||||||
|
% by the fuel cell stack is 110 kW.
|
||||||
|
|
||||||
|
% Copyright 2020 The MathWorks, Inc.
|
||||||
|
|
||||||
|
% Generate simulation results if they don't exist
|
||||||
|
if ~exist('simlog_PEMFuelCellSystem', 'var')
|
||||||
|
sim('PEMFuelCellSystem')
|
||||||
|
end
|
||||||
|
|
||||||
|
% Reuse figure if it exists, else create new figure
|
||||||
|
if ~exist('h2_PEMFuelCellSystem', 'var') || ...
|
||||||
|
~isgraphics(h2_PEMFuelCellSystem, 'figure')
|
||||||
|
h2_PEMFuelCellSystem = figure('Name', 'PEMFuelCellSystem');
|
||||||
|
end
|
||||||
|
figure(h2_PEMFuelCellSystem)
|
||||||
|
clf(h2_PEMFuelCellSystem)
|
||||||
|
|
||||||
|
plotPower(simlog_PEMFuelCellSystem)
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
% Plot power produced and consumed
|
||||||
|
function plotPower(simlog)
|
||||||
|
|
||||||
|
% Get simulation results
|
||||||
|
t = simlog.Membrane_Electrode_Assembly.power_elec.series.time;
|
||||||
|
power_elec = simlog.Membrane_Electrode_Assembly.power_elec.series.values('kW');
|
||||||
|
power_dissipated = simlog.Membrane_Electrode_Assembly.power_dissipated.series.values('kW');
|
||||||
|
power_compressor = simlog.Oxygen_Source.Compressor.power.series.values('kW');
|
||||||
|
power_pump = simlog.Cooling_System.Pump.power.series.values('kW');
|
||||||
|
|
||||||
|
% Compute net power gain
|
||||||
|
power_net = power_elec - power_compressor - power_pump;
|
||||||
|
|
||||||
|
% Plot results
|
||||||
|
handles(1) = subplot(2, 1, 1);
|
||||||
|
plot(t, power_net, 'LineWidth', 1)
|
||||||
|
hold on
|
||||||
|
plot(t, power_elec, 'LineWidth', 1)
|
||||||
|
plot(t, power_compressor, 'LineWidth', 1)
|
||||||
|
plot(t, power_pump, 'LineWidth', 1)
|
||||||
|
grid on
|
||||||
|
legend('Net Gain', 'Electrical', 'Air Compressor', 'Coolant Pump', 'Location', 'best')
|
||||||
|
title('Power Produced and Consumed (kW)')
|
||||||
|
handles(2) = subplot(2, 1, 2);
|
||||||
|
plot(t, power_dissipated, 'LineWidth', 1)
|
||||||
|
grid on
|
||||||
|
title('Heat Dissipated (kW)')
|
||||||
|
xlabel('Time (s)')
|
||||||
|
|
||||||
|
linkaxes(handles, 'x')
|
||||||
|
|
||||||
|
end
|
||||||
|
|
@ -0,0 +1,78 @@
|
||||||
|
% Code to plot simulation results from PEMFuelCellSystem
|
||||||
|
%% Plot Description:
|
||||||
|
%
|
||||||
|
% This plot show the thermal efficiency of the fuel cell and its reactant
|
||||||
|
% utilization fraction. The thermal efficiency indicates the fraction of
|
||||||
|
% the hydrogen fuel's energy that the fuel cell has converted to useful
|
||||||
|
% electrical work. The theoretical maximum efficiency for a PEM fuel cell
|
||||||
|
% is 83%. However, actual efficiency is around 60% due to internal losses.
|
||||||
|
% Near maximum current, the efficiency drops to around 45%.
|
||||||
|
%
|
||||||
|
% The reactant utilization is the fraction of the reactants, H2 and O2,
|
||||||
|
% flowing into the fuel cell stack that has been consumed by the fuel cell.
|
||||||
|
% While higher utilization makes better use of the gases flowing through
|
||||||
|
% the fuel cell, it decreases the concentration of the reactants and thus
|
||||||
|
% reduces the voltage produced. Unconsumed O2 is vented to the environment,
|
||||||
|
% but unconsumed H2 is recirculated back to the anode to avoid waste.
|
||||||
|
% However, in practice, the H2 is periodically purged to remove
|
||||||
|
% contaminants.
|
||||||
|
|
||||||
|
% Copyright 2020 The MathWorks, Inc.
|
||||||
|
|
||||||
|
% Generate simulation results if they don't exist
|
||||||
|
if ~exist('simlog_PEMFuelCellSystem', 'var')
|
||||||
|
sim('PEMFuelCellSystem')
|
||||||
|
end
|
||||||
|
|
||||||
|
% Reuse figure if it exists, else create new figure
|
||||||
|
if ~exist('h3_PEMFuelCellSystem', 'var') || ...
|
||||||
|
~isgraphics(h3_PEMFuelCellSystem, 'figure')
|
||||||
|
h3_PEMFuelCellSystem = figure('Name', 'PEMFuelCellSystem');
|
||||||
|
end
|
||||||
|
figure(h3_PEMFuelCellSystem)
|
||||||
|
clf(h3_PEMFuelCellSystem)
|
||||||
|
|
||||||
|
plotEfficiencyUtilization(simlog_PEMFuelCellSystem)
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
% Plot efficiency and utilization
|
||||||
|
function plotEfficiencyUtilization(simlog)
|
||||||
|
|
||||||
|
% Get simulation results
|
||||||
|
t = simlog.Membrane_Electrode_Assembly.efficiency_HHV.series.time;
|
||||||
|
eta_HHV = simlog.Membrane_Electrode_Assembly.efficiency_HHV.series.values('1');
|
||||||
|
eta_LHV = simlog.Membrane_Electrode_Assembly.efficiency_LHV.series.values('1');
|
||||||
|
H2_consumed = simlog.Membrane_Electrode_Assembly.H2_consumed.series.values('g/s');
|
||||||
|
O2_consumed = simlog.Membrane_Electrode_Assembly.O2_consumed.series.values('g/s');
|
||||||
|
H2_in = simlog.Anode_Gas_Channels.Pipe_MA.mdot_g_A.series.values('g/s');
|
||||||
|
O2_in = simlog.Cathode_Gas_Channels.Pipe_MA.mdot_g_A.series.values('g/s');
|
||||||
|
|
||||||
|
% Compute utilization fraction
|
||||||
|
util_H2 = H2_consumed./H2_in;
|
||||||
|
util_O2 = O2_consumed./O2_in;
|
||||||
|
|
||||||
|
% Omit initial transients
|
||||||
|
idx = t >= 1;
|
||||||
|
|
||||||
|
% Plot results
|
||||||
|
handles(1) = subplot(2, 1, 1);
|
||||||
|
plot(t, eta_HHV*100, 'LineWidth', 1)
|
||||||
|
hold on
|
||||||
|
plot(t, eta_LHV*100, 'LineWidth', 1)
|
||||||
|
legend('Based on HHV', 'Based on LHV', 'Location', 'best')
|
||||||
|
grid on
|
||||||
|
title('Thermal Efficiency (%)')
|
||||||
|
handles(2) = subplot(2, 1, 2);
|
||||||
|
plot(t(idx), util_H2(idx)*100, 'LineWidth', 1)
|
||||||
|
hold on
|
||||||
|
plot(t(idx), util_O2(idx)*100, 'LineWidth', 1)
|
||||||
|
grid on
|
||||||
|
ylim([0 100])
|
||||||
|
legend('H_2', 'O_2', 'Location', 'best')
|
||||||
|
title('Reactant Utilization (%)')
|
||||||
|
xlabel('Time (s)')
|
||||||
|
|
||||||
|
linkaxes(handles, 'x')
|
||||||
|
|
||||||
|
end
|
||||||
|
|
@ -0,0 +1,66 @@
|
||||||
|
% Code to plot simulation results from PEMFuelCellSystem
|
||||||
|
%% Plot Description:
|
||||||
|
%
|
||||||
|
% This plot shows the temperatures at various locations in the system. The
|
||||||
|
% fuel cell stack temperature is maintained at a maximum of 80 degC by the
|
||||||
|
% cooling system. Fuel flowing to the anode is warmed by the recirculated
|
||||||
|
% flow. Air flowing to the cathode is warmed by the compressor.
|
||||||
|
%
|
||||||
|
% Maintaining an optimal temperature is critical to the operation of the
|
||||||
|
% fuel cell because higher temperatures lower the relative humidity which
|
||||||
|
% increases the membrane resistance. In this model, the cooling system is
|
||||||
|
% operated by a simple control of the coolant pump flow rate. The plot
|
||||||
|
% shows the temperature of the coolant after it has absorbed heat from the
|
||||||
|
% fuel cell stack and after it has rejected heat in the radiator.
|
||||||
|
|
||||||
|
% Copyright 2020 The MathWorks, Inc.
|
||||||
|
|
||||||
|
% Generate simulation results if they don't exist
|
||||||
|
if ~exist('simlog_PEMFuelCellSystem', 'var')
|
||||||
|
sim('PEMFuelCellSystem')
|
||||||
|
end
|
||||||
|
|
||||||
|
% Reuse figure if it exists, else create new figure
|
||||||
|
if ~exist('h4_PEMFuelCellSystem', 'var') || ...
|
||||||
|
~isgraphics(h4_PEMFuelCellSystem, 'figure')
|
||||||
|
h4_PEMFuelCellSystem = figure('Name', 'PEMFuelCellSystem');
|
||||||
|
end
|
||||||
|
figure(h4_PEMFuelCellSystem)
|
||||||
|
clf(h4_PEMFuelCellSystem)
|
||||||
|
|
||||||
|
plotTemperature(simlog_PEMFuelCellSystem)
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
% Plot temperature in fuel cell and coolant system
|
||||||
|
function plotTemperature(simlog)
|
||||||
|
|
||||||
|
% Get simulation results
|
||||||
|
t = simlog.MEA_Thermal_Mass.T.series.time;
|
||||||
|
T_stack = simlog.MEA_Thermal_Mass.T.series.values('degC');
|
||||||
|
T_anode = simlog.Anode_Gas_Channels.Pipe_MA.A.T.series.values('degC');
|
||||||
|
T_cathode = simlog.Cathode_Gas_Channels.Pipe_MA.A.T.series.values('degC');
|
||||||
|
T_coolant = simlog.Cooling_System.Fuel_Cell_Coolant_Channels.T_I.series.values('degC');
|
||||||
|
T_radiator = simlog.Cooling_System.Radiator.T_I.series.values('degC');
|
||||||
|
mdot_pump = simlog.Cooling_System.Pump.mdot_A.series.values('kg/s');
|
||||||
|
|
||||||
|
% Plot results
|
||||||
|
handles(1) = subplot(2, 1, 1);
|
||||||
|
plot(t, T_stack, 'LineWidth', 1)
|
||||||
|
hold on
|
||||||
|
plot(t, T_anode, 'LineWidth', 1)
|
||||||
|
plot(t, T_cathode, 'LineWidth', 1)
|
||||||
|
plot(t, T_coolant, 'LineWidth', 1)
|
||||||
|
plot(t, T_radiator, 'LineWidth', 1)
|
||||||
|
grid on
|
||||||
|
legend('Stack', 'Anode Inflow', 'Cathode Inflow', 'Coolant at Stack', 'Coolant at Radiator', 'Location', 'best')
|
||||||
|
title('Temperature (degC)')
|
||||||
|
handles(2) = subplot(2, 1, 2);
|
||||||
|
plot(t, mdot_pump, 'LineWidth', 1)
|
||||||
|
grid on
|
||||||
|
title('Coolant Pump Mass Flow Rate (kg/s)')
|
||||||
|
xlabel('Time (s)')
|
||||||
|
|
||||||
|
linkaxes(handles, 'x')
|
||||||
|
|
||||||
|
end
|
||||||
|
|
@ -0,0 +1,65 @@
|
||||||
|
% Code to plot simulation results from PEMFuelCellSystem
|
||||||
|
%% Plot Description:
|
||||||
|
%
|
||||||
|
% This plot shows the mass of hydrogen used during operation and the
|
||||||
|
% corresponding decrease in the hydrogen tank pressure. The energy of the
|
||||||
|
% consumed hydrogen fuel is converted to electrical energy.
|
||||||
|
|
||||||
|
% Copyright 2020 The MathWorks, Inc.
|
||||||
|
|
||||||
|
% Generate simulation results if they don't exist
|
||||||
|
if ~exist('simlog_PEMFuelCellSystem', 'var')
|
||||||
|
sim('PEMFuelCellSystem')
|
||||||
|
end
|
||||||
|
|
||||||
|
% Reuse figure if it exists, else create new figure
|
||||||
|
if ~exist('h5_PEMFuelCellSystem', 'var') || ...
|
||||||
|
~isgraphics(h5_PEMFuelCellSystem, 'figure')
|
||||||
|
h5_PEMFuelCellSystem = figure('Name', 'PEMFuelCellSystem');
|
||||||
|
end
|
||||||
|
figure(h5_PEMFuelCellSystem)
|
||||||
|
clf(h5_PEMFuelCellSystem)
|
||||||
|
|
||||||
|
plotEnergy(simlog_PEMFuelCellSystem, ...
|
||||||
|
getVariable(get_param('PEMFuelCellSystem', 'ModelWorkspace'), 'tank_V'))
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
% Plot energy produced and consumed
|
||||||
|
function plotEnergy(simlog, tank_V)
|
||||||
|
|
||||||
|
% Convert to m^3
|
||||||
|
tank_V = tank_V * 1e-3;
|
||||||
|
|
||||||
|
% Get simulation results
|
||||||
|
t = simlog.Hydrogen_Source.Fuel_Tank.p_I.series.time;
|
||||||
|
p_fuel = simlog.Hydrogen_Source.Fuel_Tank.p_I.series.values('MPa');
|
||||||
|
rho_fuel = simlog.Hydrogen_Source.Fuel_Tank.rho_I.series.values('kg/m^3');
|
||||||
|
power_elec = simlog.Membrane_Electrode_Assembly.power_elec.series.values('kW');
|
||||||
|
power_dissipated = simlog.Membrane_Electrode_Assembly.power_dissipated.series.values('kW');
|
||||||
|
|
||||||
|
M_consumed = (rho_fuel(1) - rho_fuel)*tank_V;
|
||||||
|
E_fuel = cumtrapz(t, power_elec + power_dissipated) / 3600;
|
||||||
|
E_produced = cumtrapz(t, power_elec) / 3600;
|
||||||
|
|
||||||
|
% Plot results
|
||||||
|
handles(1) = subplot(3, 1, 1);
|
||||||
|
plot(t, p_fuel, 'LineWidth', 1)
|
||||||
|
grid on
|
||||||
|
title('Fuel Tank Pressure (MPa)')
|
||||||
|
handles(2) = subplot(3, 1, 2);
|
||||||
|
plot(t, M_consumed, 'LineWidth', 1)
|
||||||
|
grid on
|
||||||
|
title('Hydrogen Consumed (kg)')
|
||||||
|
handles(2) = subplot(3, 1, 3);
|
||||||
|
plot(t, E_fuel, 'LineWidth', 1)
|
||||||
|
hold on
|
||||||
|
plot(t, E_produced, 'LineWidth', 1)
|
||||||
|
grid on
|
||||||
|
legend('Hydrogen Consumed', 'Energy Produced', 'Location', 'best')
|
||||||
|
title('Energy (kWh)')
|
||||||
|
xlabel('Time (s)')
|
||||||
|
|
||||||
|
linkaxes(handles, 'x')
|
||||||
|
|
||||||
|
end
|
||||||
Binary file not shown.
|
|
@ -0,0 +1,6 @@
|
||||||
|
<?xml version="1.0" encoding="UTF-8"?>
|
||||||
|
<MF0 version="1.1" packageUris="http://schema.mathworks.com/mf0/SlCache/19700101">
|
||||||
|
<slcache.FileAttributes type="slcache.FileAttributes" uuid="b42c146e-fe85-42df-9b99-23542a6c1f40">
|
||||||
|
<checksum>mbmopmGIoDRt8MVKbX70Ws8wIS5ByX3TRWfEs9xR1RggZQ3FJVfeveSmi2Xzru+1y/KC09UiQWui2s1zeVwxog==</checksum>
|
||||||
|
</slcache.FileAttributes>
|
||||||
|
</MF0>
|
||||||
Binary file not shown.
Binary file not shown.
|
|
@ -0,0 +1,6 @@
|
||||||
|
<?xml version="1.0" encoding="UTF-8"?>
|
||||||
|
<MF0 version="1.1" packageUris="http://schema.mathworks.com/mf0/SlCache/19700101">
|
||||||
|
<slcache.FileAttributes type="slcache.FileAttributes" uuid="d92eafd1-ffa3-4323-8aa7-03b8e13147c3">
|
||||||
|
<checksum>ekcxxpQ/6HRJC8J/tqoVG16wr9UZSgRqVNb+q5DgwRBcUIvB8a3lJ/vtwwx6/F4OqfWyoMon1TG+xKDQ9uMoYQ==</checksum>
|
||||||
|
</slcache.FileAttributes>
|
||||||
|
</MF0>
|
||||||
Binary file not shown.
|
|
@ -0,0 +1,23 @@
|
||||||
|
%Energy production
|
||||||
|
E_PV = 460000; %kWh
|
||||||
|
E_Grid = 140000; %kWh
|
||||||
|
E_ExternalGrid = 400000; %kWh
|
||||||
|
|
||||||
|
E_HospitalLoad = -1000000; %kWh
|
||||||
|
|
||||||
|
%add a units functionality?
|
||||||
|
|
||||||
|
%openExample('simscape/PEMElectrolysisSystemExample')
|
||||||
|
%Ready made electrolyser -
|
||||||
|
%Good news: it has a solar profile.
|
||||||
|
%Bad news: random constants - we want it at 128m3 hydrogen production a day
|
||||||
|
%I have attempted to find volume via the volumetric flow rate, however this
|
||||||
|
%uses "standard conditions". THIS BRINGS IT TO 12.8m3, JUST ADJUST (power?)
|
||||||
|
%SO ITS 10 TIMES BIGGER AND WE ARE GOOD!
|
||||||
|
|
||||||
|
%openExample('simscape/PEMFuelCellSystemExample')
|
||||||
|
%This one is also PEM so it should work?
|
||||||
|
|
||||||
|
%To do: Check if PEM is ok, change parameters for electrolyzer so that
|
||||||
|
%production is exactely 10times bigger, adjust fuel cell, create grid,
|
||||||
|
%ensure its power output is appropriate
|
||||||
Binary file not shown.
|
|
@ -0,0 +1,28 @@
|
||||||
|
%Energy production
|
||||||
|
E_PV = 460000; %kWh
|
||||||
|
E_Grid = 140000; %kWh
|
||||||
|
E_ExternalGrid = 400000; %kWh
|
||||||
|
|
||||||
|
E_HospitalLoad = -1000000; %kWh
|
||||||
|
|
||||||
|
%add a units functionality?
|
||||||
|
|
||||||
|
%openExample('simscape/PEMElectrolysisSystemExample')
|
||||||
|
%Ready made electrolyser -
|
||||||
|
%Good news: it has a solar profile.
|
||||||
|
%Bad news: random constants - we want it at 128m3 hydrogen production a day
|
||||||
|
%I have attempted to find volume via the volumetric flow rate, however this
|
||||||
|
%uses "standard conditions". THIS BRINGS IT TO 12.8m3, JUST ADJUST (power?)
|
||||||
|
%SO ITS 10 TIMES BIGGER AND WE ARE GOOD!
|
||||||
|
|
||||||
|
%openExample('simscape/PEMFuelCellSystemExample')
|
||||||
|
%This one is also PEM so it should work?
|
||||||
|
|
||||||
|
%To do: Check if PEM is ok, change parameters for electrolyzer so that
|
||||||
|
%production is exactely 10times bigger, adjust fuel cell, create grid,
|
||||||
|
%ensure its power output is appropriate
|
||||||
|
|
||||||
|
%%-------------------------------------------------------------------------
|
||||||
|
H2_leak = 0.001; %Fraction (%)
|
||||||
|
H2_cap = 128; %m3
|
||||||
|
E_H2_vol_h = 3000; %Wh/m3
|
||||||
Binary file not shown.
Binary file not shown.
|
|
@ -0,0 +1,28 @@
|
||||||
|
* text=auto
|
||||||
|
|
||||||
|
*.fig binary
|
||||||
|
*.mat binary
|
||||||
|
*.mdl binary diff merge=mlAutoMerge
|
||||||
|
*.mdlp binary
|
||||||
|
*.mex* binary
|
||||||
|
*.mlapp binary
|
||||||
|
*.mldatx binary
|
||||||
|
*.mlproj binary
|
||||||
|
*.mlx binary
|
||||||
|
*.p binary
|
||||||
|
*.sfx binary
|
||||||
|
*.sldd binary
|
||||||
|
*.slreqx binary merge=mlAutoMerge
|
||||||
|
*.slmx binary merge=mlAutoMerge
|
||||||
|
*.sltx binary
|
||||||
|
*.slxc binary
|
||||||
|
*.slx binary merge=mlAutoMerge
|
||||||
|
*.slxp binary
|
||||||
|
|
||||||
|
## Other common binary file types
|
||||||
|
*.docx binary
|
||||||
|
*.exe binary
|
||||||
|
*.jpg binary
|
||||||
|
*.pdf binary
|
||||||
|
*.png binary
|
||||||
|
*.xlsx binary
|
||||||
33
Shift Matlab Drive/Shift Matlab Drive/Simulations/Exisiting EST Model/est-model-main/.gitignore
vendored
Normal file
33
Shift Matlab Drive/Shift Matlab Drive/Simulations/Exisiting EST Model/est-model-main/.gitignore
vendored
Normal file
|
|
@ -0,0 +1,33 @@
|
||||||
|
# Autosave files
|
||||||
|
*.asv
|
||||||
|
*.m~
|
||||||
|
*.autosave
|
||||||
|
*.slx.r*
|
||||||
|
*.mdl.r*
|
||||||
|
|
||||||
|
# Derived content-obscured files
|
||||||
|
*.p
|
||||||
|
|
||||||
|
# Compiled MEX files
|
||||||
|
*.mex*
|
||||||
|
|
||||||
|
# Packaged app and toolbox files
|
||||||
|
*.mlappinstall
|
||||||
|
*.mltbx
|
||||||
|
|
||||||
|
# Deployable archives
|
||||||
|
*.ctf
|
||||||
|
|
||||||
|
# Generated helpsearch folders
|
||||||
|
helpsearch*/
|
||||||
|
|
||||||
|
# Code generation folders
|
||||||
|
slprj/
|
||||||
|
sccprj/
|
||||||
|
codegen/
|
||||||
|
|
||||||
|
# Cache files
|
||||||
|
*.slxc
|
||||||
|
|
||||||
|
# Cloud based storage dotfile
|
||||||
|
.MATLABDriveTag
|
||||||
|
|
@ -0,0 +1,8 @@
|
||||||
|
mu = 0.1; % friction coefficient (typical value 0.05-0.2)
|
||||||
|
eff_gen = 0.9; % generator efficiency (90%)
|
||||||
|
|
||||||
|
[Beta, nUse, PsellExtraction, PbuyExtraction] = dissipationBeta(Pdemand, massBlock, v, rho, Cd, areaBlock, n, mu, eta_gen);
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
@ -0,0 +1,7 @@
|
||||||
|
if PtoInjection <= n * P_required
|
||||||
|
PSell = nUse * P_required - PtoInjection; % Excess power to sell
|
||||||
|
PBuy = 0;
|
||||||
|
else
|
||||||
|
PBuy = PtoInjection - n * P_required; % Power deficit to buy
|
||||||
|
PSell = 0;
|
||||||
|
end
|
||||||
|
|
@ -0,0 +1,6 @@
|
||||||
|
*.asv
|
||||||
|
*.mat
|
||||||
|
.DS_Store
|
||||||
|
*.slxc
|
||||||
|
*.l
|
||||||
|
simulink_cache.xml
|
||||||
Binary file not shown.
|
|
@ -0,0 +1,21 @@
|
||||||
|
MIT License
|
||||||
|
|
||||||
|
Copyright (c) 2023 Energy Storage and Transport
|
||||||
|
|
||||||
|
Permission is hereby granted, free of charge, to any person obtaining a copy
|
||||||
|
of this software and associated documentation files (the "Software"), to deal
|
||||||
|
in the Software without restriction, including without limitation the rights
|
||||||
|
to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
|
||||||
|
copies of the Software, and to permit persons to whom the Software is
|
||||||
|
furnished to do so, subject to the following conditions:
|
||||||
|
|
||||||
|
The above copyright notice and this permission notice shall be included in all
|
||||||
|
copies or substantial portions of the Software.
|
||||||
|
|
||||||
|
THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
|
||||||
|
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
|
||||||
|
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
|
||||||
|
AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
|
||||||
|
LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
|
||||||
|
OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
|
||||||
|
SOFTWARE.
|
||||||
|
|
@ -0,0 +1,238 @@
|
||||||
|
# EST-model
|
||||||
|
This project contains the *Simulink model* for the **Energy Storage and Transport (EST)** project. This Simulink model contains a simplified version of a real-life energy storage and transport system, which describes the flow of energy in such a system. Supporting MATLAB files are provided which can be used to predefine parameters and to post-process data into figures.
|
||||||
|
|
||||||
|
## Getting started
|
||||||
|
To install the simulink model, just clone or download the entire repository (under `<> Code`) and open [EST.slx](EST.slx). If you download the zip file for the repository, make sure to properly unzip it, as otherwise errors occur.
|
||||||
|
|
||||||
|
### Version requirements
|
||||||
|
The simulink model [EST.slx](EST.slx) requires **Matlab R2022b** or newer, with the **Simulink** toolbox installed. Older, untested, versions are available in the [versions](versions/) directory. In order to run an older version, copy the Simulink file corresponding to your version into the main directory and make sure to restart both MATLAB and Simulink before running it.
|
||||||
|
|
||||||
|
### Overview of files and directory structure
|
||||||
|
* [EST.slx](EST.slx): Main file, containing the runnable Simulink model of the EST system.
|
||||||
|
* [preprocessing.m](preprocessing.m): Matlab script to define the model parameters and to read the supply and demand data files. Automatically executed by the Simulink model before running.
|
||||||
|
* [postprocessing.m](postprocessing.m): Matlab script to plot the model results. Automatically executed by the Simulink model after running.
|
||||||
|
* [data directory](data/): Directory from which the supply and demand data is to be read. By default, the directory contains example files for a storage system in a household with solar panels (see [Running the model](#running-the-model)), and the files for a charge/discharge cycle (see [Charge/discharge cycle](#chargedischarge-cycle)).
|
||||||
|
|
||||||
|
## Running the model
|
||||||
|
To run the model, open [EST.slx](EST.slx) and click the run button:
|
||||||
|
|
||||||
|
After running is complete, the following output is displayed:
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
> [!CAUTION]
|
||||||
|
> When the code gives an error, carefully check the following:
|
||||||
|
> - Did you install the complete repository? Only downloading [EST.slx](EST.slx) does not work.
|
||||||
|
> - If you downloaded the zip file, did you unpack it?
|
||||||
|
> - Is the Matlab current working directory set to the right folder, that is, the folder in which [EST.slx](EST.slx) is placed.
|
||||||
|
|
||||||
|
The results displayed above correspond to the default model settings, as configured in the [preprocessing.m](preprocessing.m) Matlab script:
|
||||||
|
|
||||||
|
```matlab
|
||||||
|
% Pre-processing script for the EST Simulink model. This script is invoked
|
||||||
|
% before the Simulink model starts running (initFcn callback function).
|
||||||
|
|
||||||
|
%% Load the supply and demand data
|
||||||
|
|
||||||
|
timeUnit = 's';
|
||||||
|
|
||||||
|
supplyFile = "SolarExample_supply.csv";
|
||||||
|
supplyUnit = "kW";
|
||||||
|
|
||||||
|
% load the supply data
|
||||||
|
Supply = loadSupplyData(supplyFile, timeUnit, supplyUnit);
|
||||||
|
|
||||||
|
demandFile = "SolarExample_demand.csv";
|
||||||
|
demandUnit = "kW";
|
||||||
|
|
||||||
|
% load the demand data
|
||||||
|
Demand = loadDemandData(demandFile, timeUnit, demandUnit);
|
||||||
|
|
||||||
|
%% Simulation settings
|
||||||
|
|
||||||
|
deltat = 5*unit("min");
|
||||||
|
stopt = min([Supply.Timeinfo.End, Demand.Timeinfo.End]);
|
||||||
|
|
||||||
|
%% System parameters
|
||||||
|
|
||||||
|
% transport from supply
|
||||||
|
aSupplyTransport = 0.01; % Dissipation coefficient
|
||||||
|
|
||||||
|
% injection system
|
||||||
|
aInjection = 0.1; % Dissipation coefficient
|
||||||
|
|
||||||
|
% storage system
|
||||||
|
EStorageMax = 10.*unit("kWh"); % Maximum energy
|
||||||
|
EStorageMin = 0.0*unit("J"); % Minimum energy
|
||||||
|
EStorageInitial = 2.0*unit("kWh"); % Initial energy
|
||||||
|
bStorage = 1e-6/unit("s"); % Storage dissipation coefficient
|
||||||
|
|
||||||
|
% extraction system
|
||||||
|
aExtraction = 0.1; % Dissipation coefficient
|
||||||
|
|
||||||
|
% transport to demand
|
||||||
|
aDemandTransport = 0.01; % Dissipation coefficient
|
||||||
|
```
|
||||||
|
|
||||||
|
## Theory and implementation
|
||||||
|
The EST system transports energy from the `Supply` to the `Demand`, both represented by a `block` in the Simulink model, possibly storing the energy in between. The EST model consists of five components (`blocks`), in the order of the energy flow:
|
||||||
|
1. `Transport from supply`: transports the energy from the supply site to the storage site.
|
||||||
|
2. `Injection`: inserts energy into the storage container.
|
||||||
|
3. `Storage`: container in which the energy is stored.
|
||||||
|
4. `Extraction`: extracts energy from the storage container.
|
||||||
|
5. `Transport to demand`: transports the energy from the storage site to the demand site.
|
||||||
|
|
||||||
|
The flow of energy between these components is managed by a `controller`, which ensures that the electrical power balance is satisfied through load balancing (buying or selling of energy).
|
||||||
|
|
||||||
|
Each subsystem in the model can be described by
|
||||||
|
|
||||||
|
$$\dot{E}=P_{\rm in}-P_{\rm out}-D$$
|
||||||
|
|
||||||
|
where $\dot{E}$ is the change of total energy in the subsystem, $P_{\rm in}$ the incoming power, $P_{\rm out}$ the outgoing power and $D$ the rate of dissipation in the subsystem.
|
||||||
|
|
||||||
|
For the `Transport from supply`, `Injection`, `Extraction`and `Transport to demand` components, the EST model assumes $\dot{E}=0$ and
|
||||||
|
|
||||||
|
$$D= a P_{\rm in}$$
|
||||||
|
|
||||||
|
where $a$ [-] is the subsystem dissipation coefficients. It then follows that
|
||||||
|
|
||||||
|
$$P_{\rm out} = P_{\rm in} - D = (1-a) P_{\rm in}$$
|
||||||
|
|
||||||
|
These relations are implemented in Simulink through `Matlab function blocks`, for example for the `Injection` component:
|
||||||
|
|
||||||
|
```matlab
|
||||||
|
function [PfromInjection, DInjection] = injection(PtoInjection, aInjection)
|
||||||
|
DInjection = aInjection * PtoInjection;
|
||||||
|
PfromInjection = PtoInjection - DInjection;
|
||||||
|
```
|
||||||
|
|
||||||
|
Note that $P_{\rm in}$ and $P_{\rm out}$ are here represented by `PtoInjection` and `PfromInjection`, respectively. The `Transport from supply` component follows the same implementation:
|
||||||
|
|
||||||
|
```matlab
|
||||||
|
function [PfromSupplyTransport, DSupplyTransport] = supplyTransport(PtoSupplyTransport, aSupplyTransport)
|
||||||
|
DSupplyTransport = aSupplyTransport * PtoSupplyTransport;
|
||||||
|
PfromSupplyTransport = PtoSupplyTransport - DSupplyTransport;
|
||||||
|
```
|
||||||
|
|
||||||
|
When $P_{\rm out}$ serves as the input of the system, for example for the `Extraction` component, the dissipation and power functions must be rewritten to
|
||||||
|
|
||||||
|
$$D = \frac{a}{1-a} P_{\rm out}$$
|
||||||
|
|
||||||
|
and
|
||||||
|
|
||||||
|
$$P_{\rm in} = P_{\rm out} + D = \frac{1}{1-a} P_{\rm out}$$
|
||||||
|
|
||||||
|
The implementation then follows as:
|
||||||
|
|
||||||
|
```matlab
|
||||||
|
function [PtoExtraction, DExtraction] = extraction(PfromExtraction, aExtraction)
|
||||||
|
DExtraction = aExtraction / (1-aExtraction) * PfromExtraction;
|
||||||
|
PtoExtraction = PfromExtraction + DExtraction;
|
||||||
|
```
|
||||||
|
|
||||||
|
Similarly, for the `Transport to demand` component, the implementation reads:
|
||||||
|
|
||||||
|
```matlab
|
||||||
|
function [PtoDemandTransport, DDemandTransport] = demandTransport(PfromDemandTransport, aDemandTransport)
|
||||||
|
DDemandTransport = aDemandTransport / (1-aDemandTransport) * PfromDemandTransport;
|
||||||
|
PtoDemandTransport = PfromDemandTransport + DDemandTransport;
|
||||||
|
```
|
||||||
|
|
||||||
|
For the `Storage` component, the dissipation model
|
||||||
|
|
||||||
|
$$D= b (E - E_{\rm min})$$
|
||||||
|
|
||||||
|
is assumed, where $E_{\rm min}$ is the minimum energy capacity of the system (by default set to 0) and $b$ [1/s] is the storage dissipation coefficient. This model essential states that the dissipation is proportional to the amount of energy stored.
|
||||||
|
|
||||||
|
Substitution of this dissipation model in the power balance results in the differential equation (DE)
|
||||||
|
|
||||||
|
$$\dot{E} + b E =P_{\rm in}-P_{\rm out}+b E_{\rm min}$$
|
||||||
|
|
||||||
|
In the Simulink model, this differential equation is integrated explicitly, meaning that $\dot{E}$ is computed based on the energy $E$ in the previous time step:
|
||||||
|
|
||||||
|
```matlab
|
||||||
|
function [DStorage, EdotStorage]= Storage(PtoStorage, PfromStorage, bStorage, aStorage , EStorageMin, EStorage)
|
||||||
|
DStorage = bStorage * (EStorage-EStorageMin);
|
||||||
|
EdotStorage = PtoStorage - PfromStorage - DStorage;
|
||||||
|
```
|
||||||
|
|
||||||
|
## Charge/discharge cycle
|
||||||
|
To illustrate the theory behind the model, we consider a single charge/discharge cycle. In this cycle, the energy system is charged for $T_{\rm charge}=3$ [h] with a power of $P_{\rm charge}=15$ [kW], after which the energy is stored for $T_{\rm store}=8$ [h], until the system is discharged for $T_{\rm discharge}=1$ [h] with a power of $P_{\rm discharge}=5$ [kW]. The corresponding supply and demand signals are stored in the [CycleExample_supply.csv](CycleExample_supply.csv) and [CycleExample_demand.csv](CycleExample_supply.csv) files. To read these files, the [preprocessing.m](preprocessing.m) file is configured as:
|
||||||
|
|
||||||
|
```matlab
|
||||||
|
% Pre-processing script for the EST Simulink model. This script is invoked
|
||||||
|
% before the Simulink model starts running (initFcn callback function).
|
||||||
|
|
||||||
|
%% Load the supply and demand data
|
||||||
|
|
||||||
|
timeUnit = 's';
|
||||||
|
|
||||||
|
supplyFile = "CycleExample_supply.csv";
|
||||||
|
supplyUnit = "kW";
|
||||||
|
|
||||||
|
% load the supply data
|
||||||
|
Supply = loadSupplyData(supplyFile, timeUnit, supplyUnit);
|
||||||
|
|
||||||
|
demandFile = "CycleExample_demand.csv";
|
||||||
|
demandUnit = "kW";
|
||||||
|
|
||||||
|
% load the demand data
|
||||||
|
Demand = loadDemandData(demandFile, timeUnit, demandUnit);
|
||||||
|
```
|
||||||
|
The time step size, `deltat`, used in the simulation is also specified in [preprocessing.m](preprocessing.m) and the final simulation time, `stopt`, is determined from the loaded time series:
|
||||||
|
```matlab
|
||||||
|
%% Simulation settings
|
||||||
|
|
||||||
|
deltat = 5*unit("min");
|
||||||
|
stopt = min([Supply.Timeinfo.End, Demand.Timeinfo.End]);
|
||||||
|
```
|
||||||
|
Finally, the dissipation coefficients are given by:
|
||||||
|
```matlab
|
||||||
|
%% System parameters
|
||||||
|
|
||||||
|
% transport from supply
|
||||||
|
aSupplyTransport = 0.01; % Dissipation coefficient
|
||||||
|
|
||||||
|
% injection system
|
||||||
|
aInjection = 0.1; % Dissipation coefficient
|
||||||
|
|
||||||
|
% storage system
|
||||||
|
EStorageMax = 40*unit("kWh"); % Maximum energy
|
||||||
|
EStorageMin = 0.0*unit("J"); % Minimum energy
|
||||||
|
EStorageInitial = 0.0*unit("J"); % Initial energy
|
||||||
|
bStorage = 5e-5/unit('s'); % Storage dissipation coefficient
|
||||||
|
|
||||||
|
% extraction system
|
||||||
|
aExtraction = 0.1; % Dissipation coefficient
|
||||||
|
|
||||||
|
% transport to demand
|
||||||
|
aDemandTransport = 0.01; % Dissipation coefficient
|
||||||
|
```
|
||||||
|
|
||||||
|
For this particular scenario, an exact solution to the model exists, which can be used to **verify the Simulink implementation**. During charging, the power balance differential equation for the storage container reads
|
||||||
|
|
||||||
|
$$\dot{E} + b E = c P_{\rm supply}$$
|
||||||
|
|
||||||
|
where $c = (1-a_{\rm supplyTransport}) (1-a_{\rm Injection})$. With the initial condition $E(0)=0$, the solution is given by
|
||||||
|
|
||||||
|
$$E = (1 - e^{-bt}) \frac{c}{b} P_{\rm supply} \qquad 0 \leq t < T_{\rm charge}$$
|
||||||
|
|
||||||
|
During storage, the power balance reads
|
||||||
|
|
||||||
|
$$\dot{E} + b E = 0$$
|
||||||
|
|
||||||
|
with the initial condition $E(T_{\rm charge}) = (1 - e^{-bT_{\rm charge}}) \frac{c}{b} P_{\rm supply} := E_{\rm charge}$. The solution during storage is given by
|
||||||
|
|
||||||
|
$$E = E_{\rm charge} e^{-b(t - T_{\rm charge})} \qquad T_{\rm charge} \leq t \leq \tau$$
|
||||||
|
|
||||||
|
where $\tau = T_{\rm charge} + T_{\rm store}$. Finally, during discharging, the differential equation reads
|
||||||
|
|
||||||
|
$$\dot{E} + b E = -d P_{\rm demand}$$
|
||||||
|
|
||||||
|
where $d = (1-a_{\rm Extraction})^{-1}(1-a_{\rm demandTransport})^{-1}$. With the initial condition $E(\tau)=E_{\rm charge} e^{-b T_{\rm store}} := E_{\rm store}$, the solution is given by
|
||||||
|
|
||||||
|
$$E = (e^{-b(t-\tau)} -1) \frac{d}{b} P_{\rm demand} + E_{\rm store} e^{-b(t-\tau)} \qquad \tau \leq t \leq \tau + T_{\rm discharge}$$
|
||||||
|
|
||||||
|
Comparison of this exact solution with the Simulink model conveys that the model solves the model equations as intended:
|
||||||
|
|
||||||
|
|
||||||
File diff suppressed because it is too large
Load diff
File diff suppressed because it is too large
Load diff
File diff suppressed because it is too large
Load diff
File diff suppressed because it is too large
Load diff
Binary file not shown.
|
After Width: | Height: | Size: 34 KiB |
|
|
@ -0,0 +1,162 @@
|
||||||
|
<?xml version="1.0" encoding="UTF-8" standalone="no"?>
|
||||||
|
<!-- Created with Inkscape (http://www.inkscape.org/) -->
|
||||||
|
|
||||||
|
<svg
|
||||||
|
width="157.17723mm"
|
||||||
|
height="143.70491mm"
|
||||||
|
viewBox="0 0 157.17723 143.70491"
|
||||||
|
version="1.1"
|
||||||
|
id="svg5"
|
||||||
|
inkscape:export-filename="background.png"
|
||||||
|
inkscape:export-xdpi="192"
|
||||||
|
inkscape:export-ydpi="192"
|
||||||
|
inkscape:version="1.2.1 (9c6d41e, 2022-07-14)"
|
||||||
|
sodipodi:docname="background.svg"
|
||||||
|
xml:space="preserve"
|
||||||
|
xmlns:inkscape="http://www.inkscape.org/namespaces/inkscape"
|
||||||
|
xmlns:sodipodi="http://sodipodi.sourceforge.net/DTD/sodipodi-0.dtd"
|
||||||
|
xmlns="http://www.w3.org/2000/svg"
|
||||||
|
xmlns:svg="http://www.w3.org/2000/svg"><sodipodi:namedview
|
||||||
|
id="namedview7"
|
||||||
|
pagecolor="#ffffff"
|
||||||
|
bordercolor="#000000"
|
||||||
|
borderopacity="0.25"
|
||||||
|
inkscape:showpageshadow="2"
|
||||||
|
inkscape:pageopacity="0.0"
|
||||||
|
inkscape:pagecheckerboard="0"
|
||||||
|
inkscape:deskcolor="#d1d1d1"
|
||||||
|
inkscape:document-units="mm"
|
||||||
|
showgrid="false"
|
||||||
|
inkscape:zoom="0.84096521"
|
||||||
|
inkscape:cx="264.57694"
|
||||||
|
inkscape:cy="250.30762"
|
||||||
|
inkscape:window-width="1440"
|
||||||
|
inkscape:window-height="777"
|
||||||
|
inkscape:window-x="0"
|
||||||
|
inkscape:window-y="25"
|
||||||
|
inkscape:window-maximized="1"
|
||||||
|
inkscape:current-layer="layer1"
|
||||||
|
showguides="false"><inkscape:grid
|
||||||
|
type="xygrid"
|
||||||
|
id="grid234"
|
||||||
|
originx="-7.9489565e-07"
|
||||||
|
originy="10.583334" /></sodipodi:namedview><defs
|
||||||
|
id="defs2"><rect
|
||||||
|
x="69.788693"
|
||||||
|
y="216.93353"
|
||||||
|
width="217.77436"
|
||||||
|
height="76.515314"
|
||||||
|
id="rect4694" /><marker
|
||||||
|
style="overflow:visible"
|
||||||
|
id="RoundedArrow"
|
||||||
|
refX="0"
|
||||||
|
refY="0"
|
||||||
|
orient="auto-start-reverse"
|
||||||
|
inkscape:stockid="RoundedArrow"
|
||||||
|
markerWidth="4"
|
||||||
|
markerHeight="3.2"
|
||||||
|
viewBox="0 0 6.1347524 5.9304951"
|
||||||
|
inkscape:isstock="true"
|
||||||
|
inkscape:collect="always"
|
||||||
|
preserveAspectRatio="xMidYMid"><path
|
||||||
|
transform="scale(0.7)"
|
||||||
|
d="m -0.21114562,-4.1055728 6.42229122,3.21114561 a 1,1 90 0 1 0,1.78885438 L -0.21114562,4.1055728 A 1.236068,1.236068 31.717474 0 1 -2,3 v -6 a 1.236068,1.236068 148.28253 0 1 1.78885438,-1.1055728 z"
|
||||||
|
style="fill:context-stroke;fill-rule:evenodd;stroke:none"
|
||||||
|
id="path1367" /></marker><marker
|
||||||
|
style="overflow:visible"
|
||||||
|
id="TriangleStart"
|
||||||
|
refX="0"
|
||||||
|
refY="0"
|
||||||
|
orient="auto-start-reverse"
|
||||||
|
inkscape:stockid="TriangleStart"
|
||||||
|
markerWidth="2"
|
||||||
|
markerHeight="2"
|
||||||
|
viewBox="0 0 5.3244081 6.1553851"
|
||||||
|
inkscape:isstock="true"
|
||||||
|
inkscape:collect="always"
|
||||||
|
preserveAspectRatio="xMidYMid"><path
|
||||||
|
transform="scale(0.5)"
|
||||||
|
style="fill:context-stroke;fill-rule:evenodd;stroke:context-stroke;stroke-width:1pt"
|
||||||
|
d="M 5.77,0 -2.88,5 V -5 Z"
|
||||||
|
id="path135" /></marker><rect
|
||||||
|
x="69.788693"
|
||||||
|
y="216.93353"
|
||||||
|
width="217.77436"
|
||||||
|
height="76.515314"
|
||||||
|
id="rect4913" /></defs><g
|
||||||
|
inkscape:label="Layer 1"
|
||||||
|
inkscape:groupmode="layer"
|
||||||
|
id="layer1"
|
||||||
|
transform="translate(-2.7948956e-6,10.583334)"><path
|
||||||
|
style="fill:#c81919;fill-opacity:1;stroke:#666666;stroke-width:0.1"
|
||||||
|
d="m 59.250838,119.02435 v 5.29167 h -11.90625 v 1.32292 l -3.96875,-3.96875 3.96875,-3.96875 v 1.32291 z"
|
||||||
|
id="path1406" /><rect
|
||||||
|
style="fill:#c81919;fill-opacity:0.3;stroke:none;stroke-width:4.99999;stroke-dasharray:none;stroke-opacity:1"
|
||||||
|
id="rect2686"
|
||||||
|
width="157.17723"
|
||||||
|
height="143.70491"
|
||||||
|
x="2.7948956e-06"
|
||||||
|
y="-10.583334" /><rect
|
||||||
|
style="fill:#e6e6e6;fill-opacity:1;stroke:#000000;stroke-width:0.41541;stroke-dasharray:none;stroke-opacity:1"
|
||||||
|
id="rect3703"
|
||||||
|
width="31.77051"
|
||||||
|
height="14.105489"
|
||||||
|
x="2.9574034"
|
||||||
|
y="55.991592" /><path
|
||||||
|
style="fill:#000000;fill-opacity:0.3;stroke:#c81919;stroke-width:3;stroke-dasharray:none;stroke-opacity:1;marker-end:url(#TriangleStart)"
|
||||||
|
d="m 4.1156944,60.250221 h 5.291664"
|
||||||
|
id="path3415"
|
||||||
|
sodipodi:nodetypes="cc" /><path
|
||||||
|
style="fill:#000000;fill-opacity:0.3;stroke:#000000;stroke-width:0.5;stroke-dasharray:none;stroke-opacity:1;marker-end:url(#RoundedArrow)"
|
||||||
|
d="m 4.1156944,66.696056 h 5.291664"
|
||||||
|
id="path3421"
|
||||||
|
sodipodi:nodetypes="cc" /><text
|
||||||
|
xml:space="preserve"
|
||||||
|
transform="matrix(0.26458333,0,0,0.26458333,-5.5922081,0.54175469)"
|
||||||
|
id="text4692"
|
||||||
|
style="font-size:13.3333px;white-space:pre;shape-inside:url(#rect4694);display:inline;fill:#000000;fill-opacity:1;stroke:none;stroke-width:1.88976;stroke-dasharray:none;stroke-opacity:1"><tspan
|
||||||
|
x="69.789062"
|
||||||
|
y="229.09797"
|
||||||
|
id="tspan1615">Energy flow</tspan></text><text
|
||||||
|
xml:space="preserve"
|
||||||
|
transform="matrix(0.26458333,0,0,0.26458333,-5.5493763,7.5130706)"
|
||||||
|
id="text4911"
|
||||||
|
style="font-size:13.3333px;white-space:pre;shape-inside:url(#rect4913);display:inline;fill:#000000;fill-opacity:1;stroke:none;stroke-width:1.88976;stroke-dasharray:none;stroke-opacity:1"><tspan
|
||||||
|
x="69.789062"
|
||||||
|
y="229.09797"
|
||||||
|
id="tspan1617">Data flow</tspan></text><path
|
||||||
|
style="fill:#c81919;fill-opacity:1;stroke:#666666;stroke-width:0.1"
|
||||||
|
d="m 43.433784,2.45037 v 5.2916667 h 11.90625 v 1.3229164 l 3.96875,-3.9687497 -3.96875,-3.9687499 V 2.45037 Z"
|
||||||
|
id="path512" /><path
|
||||||
|
style="fill:#c81919;fill-opacity:1;stroke:#666666;stroke-width:0.1"
|
||||||
|
d="m 98.252572,4.315956 v 3.1750001 h 11.906268 v 1.3229164 l 3.96875,-2.9104164 -3.96875,-2.9104166 V 4.315956 Z"
|
||||||
|
id="path1296"
|
||||||
|
sodipodi:nodetypes="cccccccc" /><path
|
||||||
|
style="fill:#c81919;fill-opacity:1;stroke:#666666;stroke-width:0.1"
|
||||||
|
d="m 141.90088,30.595528 h -2.95253 v 10.847933 h -1.32292 l 2.79919,3.96875 2.79918,-3.96875 h -1.32292 z"
|
||||||
|
id="path1394"
|
||||||
|
sodipodi:nodetypes="cccccccc" /><path
|
||||||
|
style="fill:#c81919;fill-opacity:1;stroke:#666666;stroke-width:0.1"
|
||||||
|
d="m 141.66492,81.395533 h -2.48061 v 10.847935 h -1.32292 l 2.56323,3.96875 2.56322,-3.96875 h -1.32292 z"
|
||||||
|
id="path1396"
|
||||||
|
sodipodi:nodetypes="cccccccc" /><path
|
||||||
|
style="fill:#c81919;fill-opacity:1;stroke:#666666;stroke-width:0.1"
|
||||||
|
d="m 114.67506,122.71267 v -2.24465 h -12.43544 v -1.32292 l -3.968776,2.44525 3.968776,2.44524 v -1.32292 z"
|
||||||
|
id="path1398"
|
||||||
|
sodipodi:nodetypes="cccccccc" /><path
|
||||||
|
style="fill:#c81919;fill-opacity:1;stroke:#666666;stroke-width:0.1"
|
||||||
|
d="m 86.867521,30.595528 h -2.95253 v 10.847933 h -1.32292 l 2.79919,3.96875 2.799177,-3.96875 h -1.322917 z"
|
||||||
|
id="path1402"
|
||||||
|
sodipodi:nodetypes="cccccccc" /><path
|
||||||
|
style="fill:#c81919;fill-opacity:1;stroke:#666666;stroke-width:0.1"
|
||||||
|
d="m 86.867521,81.924697 h -2.95253 V 92.77263 h -1.32292 l 2.79919,3.96875 2.799177,-3.96875 h -1.322917 z"
|
||||||
|
id="path1404"
|
||||||
|
sodipodi:nodetypes="cccccccc" /><path
|
||||||
|
style="fill:#c81919;fill-opacity:1;stroke:#666666;stroke-width:0.1"
|
||||||
|
d="m 52.104535,65.992483 v 2.95253 h 4.497933 v 1.32292 l 3.96875,-2.79919 -3.96875,-2.799177 v 1.322917 z"
|
||||||
|
id="path1408"
|
||||||
|
sodipodi:nodetypes="cccccccc" /><path
|
||||||
|
style="fill:#c81919;fill-opacity:1;stroke:#666666;stroke-width:0.1"
|
||||||
|
d="m 60.607968,57.525811 v 2.95253 h -4.497933 v 1.32292 l -3.96875,-2.79919 3.96875,-2.799177 v 1.322917 z"
|
||||||
|
id="path1410"
|
||||||
|
sodipodi:nodetypes="cccccccc" /></g></svg>
|
||||||
|
After Width: | Height: | Size: 7.3 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 246 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 614 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 587 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 435 KiB |
|
|
@ -0,0 +1,194 @@
|
||||||
|
% Post-processing script for the EST Simulink model. This script is invoked
|
||||||
|
% after the Simulink model is finished running (stopFcn callback function).
|
||||||
|
|
||||||
|
close all;
|
||||||
|
figure;
|
||||||
|
|
||||||
|
%% Supply and demand
|
||||||
|
subplot(2,2,1);
|
||||||
|
plot(tout/unit("day"), PSupply/unit("MW"));
|
||||||
|
hold on;
|
||||||
|
plot(tout/unit("day"), PDemand/unit("MW"));
|
||||||
|
xlim([0 tout(end)/unit("day")]);
|
||||||
|
grid on;
|
||||||
|
title('Supply and demand');
|
||||||
|
xlabel('Time [day]');
|
||||||
|
ylabel('Power [MW]');
|
||||||
|
legend("Supply","Demand");
|
||||||
|
|
||||||
|
%% Stored energy
|
||||||
|
subplot(2,2,2);
|
||||||
|
plot(tout/unit("day"), EStorage/unit("GWh"));
|
||||||
|
xlim([0 tout(end)/unit("day")]);
|
||||||
|
grid on;
|
||||||
|
title('Storage');
|
||||||
|
xlabel('Time [day]');
|
||||||
|
ylabel('Energy [GWh]');
|
||||||
|
|
||||||
|
%% Energy losses
|
||||||
|
subplot(2,2,3);
|
||||||
|
plot(tout/unit("day"), D/unit("MW"));
|
||||||
|
xlim([0 tout(end)/unit("day")]);
|
||||||
|
grid on;
|
||||||
|
title('Losses');
|
||||||
|
xlabel('Time [day]');
|
||||||
|
ylabel('Dissipation rate [MW]');
|
||||||
|
|
||||||
|
%% Load balancing
|
||||||
|
subplot(2,2,4);
|
||||||
|
plot(tout/unit("day"), PSell/unit("MW"));
|
||||||
|
hold on;
|
||||||
|
plot(tout/unit("day"), PBuy/unit("MW"));
|
||||||
|
xlim([0 tout(end)/unit("day")]);
|
||||||
|
grid on;
|
||||||
|
title('Load balancing');
|
||||||
|
xlabel('Time [day]');
|
||||||
|
ylabel('Power [MW]');
|
||||||
|
legend("Sell","Buy");
|
||||||
|
|
||||||
|
%% Pie charts
|
||||||
|
|
||||||
|
% integrate the power signals in time
|
||||||
|
EfromSupplyTransport = trapz(tout, PfromSupplyTransport);
|
||||||
|
EtoDemandTransport = trapz(tout, PtoDemandTransport);
|
||||||
|
ESell = trapz(tout, PSell);
|
||||||
|
EBuy = trapz(tout, PBuy);
|
||||||
|
EtoInjection = trapz(tout, PtoInjection);
|
||||||
|
EfromExtraction = trapz(tout, PfromExtraction);
|
||||||
|
EStorageDissipation = trapz(tout, DStorage);
|
||||||
|
EDirect = EfromSupplyTransport - ESell - EtoInjection;
|
||||||
|
ESurplus = EtoInjection-EfromExtraction-EStorageDissipation;
|
||||||
|
|
||||||
|
figure;
|
||||||
|
tiles = tiledlayout(1,2);
|
||||||
|
|
||||||
|
ax = nexttile;
|
||||||
|
pie(ax, [EDirect, EtoInjection, ESell]/EfromSupplyTransport);
|
||||||
|
lgd = legend({"Direct to demand", "To storage", "Sold"});
|
||||||
|
lgd.Layout.Tile = "south";
|
||||||
|
title(sprintf("Received energy %3.2e [GWh]", EfromSupplyTransport/unit('GWh')));
|
||||||
|
|
||||||
|
ax = nexttile;
|
||||||
|
pie(ax, [EDirect, EfromExtraction, EBuy]/EtoDemandTransport);
|
||||||
|
lgd = legend({"Direct from supply", "From storage", "Bought"});
|
||||||
|
lgd.Layout.Tile = "south";
|
||||||
|
title(sprintf("Delivered energy %3.2e [GWh]", EtoDemandTransport/unit('GWh')));
|
||||||
|
|
||||||
|
%% Number of charged towers and height of last partially charged tower (Robert's SSA10)
|
||||||
|
n_CH_dec = EStorage / Eg;
|
||||||
|
n_CH_wh = floor(n_CH_dec);
|
||||||
|
n_CH_frac = n_CH_dec - n_CH_wh;
|
||||||
|
h_last = n_CH_frac * h;
|
||||||
|
|
||||||
|
figure;
|
||||||
|
%Number of towers charged
|
||||||
|
subplot(2,1,1);
|
||||||
|
plot(tout/unit("day"), n_CH_dec);
|
||||||
|
xlim([0 tout(end)/unit("day")]);
|
||||||
|
grid on;
|
||||||
|
title('Towers Charged');
|
||||||
|
xlabel('Time [day]');
|
||||||
|
ylabel('Number of Towers Charged [-]');
|
||||||
|
|
||||||
|
%Height of last partially charged tower
|
||||||
|
subplot(2,1,2);
|
||||||
|
plot(tout/unit("day"), h_last);
|
||||||
|
xlim([0 tout(end)/unit("day")]);
|
||||||
|
ylim([(-0.2*h) (1.2*h)])
|
||||||
|
grid on;
|
||||||
|
title('Height of block on only non-fully charged/discharged tower');
|
||||||
|
xlabel('Time [day]');
|
||||||
|
ylabel('Partially charged height [m]');
|
||||||
|
|
||||||
|
%% Number of charging/discharging towers
|
||||||
|
n_CH_DCH = (PtoStorage - PfromStorage)/Ptower;
|
||||||
|
|
||||||
|
figure;
|
||||||
|
subplot(1,1,1);
|
||||||
|
plot(tout/unit("day"), n_CH_DCH);
|
||||||
|
xlim([0 tout(end)/unit("day")]);
|
||||||
|
grid on;
|
||||||
|
title('Number of towers charging (+ve) or discharging (-ve)');
|
||||||
|
xlabel('Time [day]');
|
||||||
|
ylabel('Number of Towers (dis)charging [-]');
|
||||||
|
%% Longest / shortest discharge duration (ChatGPT + Human commentary)
|
||||||
|
|
||||||
|
%PSupplyAvg = mean(PSupply)/unit("MW")
|
||||||
|
%PDemandAvg = mean(PDemand)/unit("MW")
|
||||||
|
|
||||||
|
% Logical vector where demand is greater than supply
|
||||||
|
condition = PDemand > PSupply; %Values (eg.: the avg Supply value) can be added on the PDemand side to ensure small jumps in PSupply are not treated as complete breaks in discharge)
|
||||||
|
condition = condition(:)'; % force condition to be a row vector
|
||||||
|
|
||||||
|
% Find start and end indices of stretches where condition is true
|
||||||
|
diffs = diff([0, condition, 0]);
|
||||||
|
start_idx = find(diffs == 1);
|
||||||
|
end_idx = find(diffs == -1) - 1;
|
||||||
|
|
||||||
|
% Lengths of those stretches
|
||||||
|
lengths = end_idx - start_idx + 1;
|
||||||
|
|
||||||
|
|
||||||
|
% Find the longest stretch
|
||||||
|
[~, longest_idx] = max(lengths);
|
||||||
|
longest_start = start_idx(longest_idx);
|
||||||
|
longest_end = end_idx(longest_idx);
|
||||||
|
longest_duration_steps = lengths(longest_idx);
|
||||||
|
longest_duration_seconds = longest_duration_steps * 900; % 900 seconds per interval
|
||||||
|
longest_duration_hours = longest_duration_seconds / 3600
|
||||||
|
|
||||||
|
|
||||||
|
% Find the shortest stretch (Human written)
|
||||||
|
[~, shortest_idx] = min(lengths);
|
||||||
|
shortest_start = start_idx(shortest_idx);
|
||||||
|
shortest_end = end_idx(shortest_idx);
|
||||||
|
shortest_duration_steps = lengths(shortest_idx);
|
||||||
|
shortest_duration_seconds = shortest_duration_steps * 900 % 900 seconds per interval
|
||||||
|
|
||||||
|
min_idces = [];
|
||||||
|
for i = 1:1:max(size(start_idx)) %Recording the indices of the shortest discharges
|
||||||
|
if lengths(i) == shortest_duration_steps
|
||||||
|
min_idces(end+1) = start_idx(i);
|
||||||
|
end
|
||||||
|
end
|
||||||
|
|
||||||
|
%% Average energy stored
|
||||||
|
EStorageAvg = mean(EStorage/unit("GWh")) %Average energy stored
|
||||||
|
|
||||||
|
%% Max, Min and Avg extraction power capacity
|
||||||
|
ExtCap = PDemand/unit("MW") - PSupply/unit("MW");
|
||||||
|
ExtCapFiltered = ExtCap(find(ExtCap>0)); %Filters negative values (charging rather than delivering)
|
||||||
|
%ExtCap = ExtCap(find(ExtCap>7.25)); %Filters negative values and small changes
|
||||||
|
ExtCapMax = max(ExtCapFiltered);
|
||||||
|
|
||||||
|
largest_idx = find(ExtCap == ExtCapMax); %Code finding the index of the largest discharge.
|
||||||
|
idx_before_largest = find(start_idx>=largest_idx); %All indeces of starting indeces of a dischage before the greatest discharge.
|
||||||
|
largest_length_seconds = lengths(idx_before_largest(end)) * 900;
|
||||||
|
largest_length_hours = largest_length_seconds / (60*60);
|
||||||
|
largest_lenght_all = [];
|
||||||
|
% for a = start_idx(idx_before_largest(end)):1:end_idx(idx_before_largest(end)) %All discharge values in said interval.
|
||||||
|
% largest_lenght_all(end+1) = ExtCap(a);
|
||||||
|
% end
|
||||||
|
% largest_lenght_avg = mean(largest_lenght_all)
|
||||||
|
%
|
||||||
|
% ExtCapMin = min(ExtCapFiltered)
|
||||||
|
% ExtCapAvg = mean(ExtCapFiltered)
|
||||||
|
|
||||||
|
% Max, Min and Avg extraction power capacity for the longest duration
|
||||||
|
for l = longest_start:1:longest_end
|
||||||
|
ExtCapsLongest(l) = PDemand(l)/unit("MW") - PSupply(l)/unit("MW");
|
||||||
|
end
|
||||||
|
ExtCapsLongestMax = max(ExtCapsLongest);
|
||||||
|
ExtCapsLongestMin = min(ExtCapsLongest);
|
||||||
|
ExtCapsLongestAvg = mean(ExtCapsLongest);
|
||||||
|
|
||||||
|
% Max, Min and Avg extraction power capacity for the shortest duration
|
||||||
|
ExtCapsShortest = [];
|
||||||
|
for s = 1:1:max(size(min_idces))
|
||||||
|
ExtCapsShortest(s) = ExtCap(min_idces(s));
|
||||||
|
end
|
||||||
|
ExtCapsShortestMax = max(ExtCapsShortest);
|
||||||
|
ExtCapsShortestMin = min(ExtCapsShortest);
|
||||||
|
ExtCapsShortestAvg = mean(ExtCapsShortest);
|
||||||
|
|
||||||
|
%disp(beta)
|
||||||
|
|
@ -0,0 +1,61 @@
|
||||||
|
% Pre-processing script for the EST Simulink model. This script is invoked
|
||||||
|
% before the Simulink model starts running (initFcn callback function).
|
||||||
|
|
||||||
|
%% Load the supply and demand data
|
||||||
|
|
||||||
|
timeUnit = 's';
|
||||||
|
|
||||||
|
supplyFile = "Team03_supply.csv";
|
||||||
|
supplyUnit = "MW";
|
||||||
|
|
||||||
|
% load the supply data
|
||||||
|
Supply = loadSupplyData(supplyFile, timeUnit, supplyUnit);
|
||||||
|
|
||||||
|
demandFile = "Team03_demand.csv";
|
||||||
|
demandUnit = "MW";
|
||||||
|
|
||||||
|
% load the demand data
|
||||||
|
Demand = loadDemandData(demandFile, timeUnit, demandUnit);
|
||||||
|
|
||||||
|
%% Simulation settings
|
||||||
|
|
||||||
|
deltat = 900*unit("s");
|
||||||
|
stopt = min([Supply.Timeinfo.End, Demand.Timeinfo.End]);
|
||||||
|
|
||||||
|
%% System parameters
|
||||||
|
|
||||||
|
% transmission line transport from supply (Robert SSA8)
|
||||||
|
etaTransformer = 0.99; % Transformer efficiency [-] (presentation)
|
||||||
|
Ls = 130e3; % Transmission line length [m] ("Egmond aan Zee"/"Luchterduinen" to Winterswwijk)
|
||||||
|
Rprime = 3e-5; % Resistance per unit length [Ohm/m] (presentation)
|
||||||
|
V = 132e3; % Transmission line voltage [V] (Industry standard: Carbon Trust, 2022. Reynard, 2025)
|
||||||
|
|
||||||
|
% injection system
|
||||||
|
g = 9.81; % gravitational acceleration [m/s^2]
|
||||||
|
h = 800; % !!! height of each tower [m] (Robert SSA10, to be changed after Phileines SSA10)
|
||||||
|
|
||||||
|
mu = 0.002; % friction coefficient [-] (typical value 0.05-0.2)
|
||||||
|
effMotor = 0.9;
|
||||||
|
theta = 15; % Angle of wrap on pulley
|
||||||
|
massBlock = 300000; % 500-5000[kg]/tower
|
||||||
|
v = 5; % 1-5 [m/s]
|
||||||
|
|
||||||
|
% extraction system
|
||||||
|
effGen = 0.9; % generator efficiency (9(grav0%)
|
||||||
|
mu_rolling = 0.001;
|
||||||
|
mu_bearing = 0.05;
|
||||||
|
n = 7646; % [-] no. of towers (now chosen by hand) %floor(EStorageMax/(massBlock*g*h));
|
||||||
|
|
||||||
|
|
||||||
|
% storage system
|
||||||
|
Eg = massBlock * g * h*unit("J"); % Maximum Gravitational Energy stored per tower [J]
|
||||||
|
Ptower = (Eg / deltat)*unit("W");
|
||||||
|
EStorageMax = Eg*n; % Maximum energy (now based on number of towers and their parameters) [J]
|
||||||
|
fprintf('EStorageMax = %.3f GWh',(EStorageMax./unit("GWh")));
|
||||||
|
|
||||||
|
EStorageMin = 0.0*unit("GWh"); % Minimum energy
|
||||||
|
EStorageInitial = 0.0*unit("GWh"); % Initial energy
|
||||||
|
bStorage = 0e-6/unit("s"); % Storage dissipation coefficient
|
||||||
|
|
||||||
|
% transmission line transport to demand (Robert SSA8)
|
||||||
|
Ld = 130e3; % Transmission line length [m] (Winterswwijk to Eindhoven)
|
||||||
|
|
@ -0,0 +1,29 @@
|
||||||
|
% Script defining constants, specifically the "units" dictionary.
|
||||||
|
|
||||||
|
global unit;
|
||||||
|
unit = containers.Map;
|
||||||
|
|
||||||
|
% time
|
||||||
|
unit("s") = 1.;
|
||||||
|
unit("min") = 60*unit("s");
|
||||||
|
unit("h") = 60*unit("min");
|
||||||
|
unit("day") = 24*unit("h");
|
||||||
|
unit("year") = 365*unit("day");
|
||||||
|
|
||||||
|
% energy
|
||||||
|
unit("J") = 1.;
|
||||||
|
unit("kJ") = 1000*unit("J");
|
||||||
|
unit("MJ") = 1000*unit("kJ");
|
||||||
|
unit("GJ") = 1000*unit("MJ");
|
||||||
|
|
||||||
|
% power
|
||||||
|
unit("W") = unit("J")/unit("s");
|
||||||
|
unit("kW") = 1000*unit("W");
|
||||||
|
unit("MW") = 1000*unit("kW");
|
||||||
|
unit("GW") = 1000*unit("MW");
|
||||||
|
|
||||||
|
% energy (Wh)
|
||||||
|
unit("Wh") = unit("W") *unit("h");
|
||||||
|
unit("kWh") = unit("kW")*unit("h");
|
||||||
|
unit("MWh") = unit("MW")*unit("h");
|
||||||
|
unit("GWh") = unit("GW")*unit("h");
|
||||||
|
|
@ -0,0 +1,10 @@
|
||||||
|
function Demand = loadDemandData(demandFile, timeUnit, demandUnit)
|
||||||
|
% LOADDEMANDDATA Load the demand data from the specified file.
|
||||||
|
% Demand = LOADDEMANDDATA(demandFile, timeUnit, demandUnit) loads the
|
||||||
|
% data in demandFile with the time in timeUnit and the power in
|
||||||
|
% demandUnit, returning a time series.
|
||||||
|
|
||||||
|
global unit;
|
||||||
|
demand = importdata(demandFile, ',');
|
||||||
|
Demand = timeseries(unit(demandUnit)*demand.data(:,2),unit(timeUnit)*demand.data(:,1));
|
||||||
|
end
|
||||||
|
|
@ -0,0 +1,10 @@
|
||||||
|
function Supply = loadSupplyData(supplyFile, timeUnit, supplyUnit)
|
||||||
|
% LOADSUPPLYDATA Load the supply data from the specified file.
|
||||||
|
% Supply = LOADSUPPLYDATA(supplyFile, supplyUnit) loads the
|
||||||
|
% data in supplyFile with the time in timeUnit and the power in
|
||||||
|
% supplyUnit, returning a time series.
|
||||||
|
|
||||||
|
global unit;
|
||||||
|
supply = importdata(supplyFile, ',');
|
||||||
|
Supply = timeseries(unit(supplyUnit)*supply.data(:,2),unit(timeUnit)*supply.data(:,1));
|
||||||
|
end
|
||||||
Binary file not shown.
Binary file not shown.
Binary file not shown.
|
|
@ -0,0 +1,12 @@
|
||||||
|
% Constants
|
||||||
|
massBlock = 10000; % [kg]
|
||||||
|
v = 2; % [m/s] (slower than extraction to reduce losses)
|
||||||
|
rho = 1.293; % [kg/m³]
|
||||||
|
Cd = 1.05; % Drag coefficient
|
||||||
|
areaBlock = 100; % [m²]
|
||||||
|
n = 10; % Number of towers
|
||||||
|
mu = 0.05; % Friction coefficient
|
||||||
|
eta_motor = 0.85; % Motor efficiency
|
||||||
|
|
||||||
|
% Run injection model
|
||||||
|
[Beta_inj, nUse, Pbuy, Psell] = injectionBeta(1e6, massBlock, v, rho, Cd, areaBlock, n, mu, eta_motor);
|
||||||
|
|
@ -0,0 +1,2 @@
|
||||||
|
<?xml version="1.0" encoding="UTF-8"?>
|
||||||
|
<MATLABProject xmlns="http://www.mathworks.com/MATLABProjectFile" xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" version="1.0"/>
|
||||||
|
|
@ -0,0 +1,93 @@
|
||||||
|
# EST Model
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
## Getting started
|
||||||
|
|
||||||
|
To make it easy for you to get started with GitLab, here's a list of recommended next steps.
|
||||||
|
|
||||||
|
Already a pro? Just edit this README.md and make it your own. Want to make it easy? [Use the template at the bottom](#editing-this-readme)!
|
||||||
|
|
||||||
|
## Add your files
|
||||||
|
|
||||||
|
- [ ] [Create](https://docs.gitlab.com/ee/user/project/repository/web_editor.html#create-a-file) or [upload](https://docs.gitlab.com/ee/user/project/repository/web_editor.html#upload-a-file) files
|
||||||
|
- [ ] [Add files using the command line](https://docs.gitlab.com/topics/git/add_files/#add-files-to-a-git-repository) or push an existing Git repository with the following command:
|
||||||
|
|
||||||
|
```
|
||||||
|
cd existing_repo
|
||||||
|
git remote add origin https://gitlab.tue.nl/p.b.r.arnaud.de.calavon/est-model.git
|
||||||
|
git branch -M main
|
||||||
|
git push -uf origin main
|
||||||
|
```
|
||||||
|
|
||||||
|
## Integrate with your tools
|
||||||
|
|
||||||
|
- [ ] [Set up project integrations](https://gitlab.tue.nl/p.b.r.arnaud.de.calavon/est-model/-/settings/integrations)
|
||||||
|
|
||||||
|
## Collaborate with your team
|
||||||
|
|
||||||
|
- [ ] [Invite team members and collaborators](https://docs.gitlab.com/ee/user/project/members/)
|
||||||
|
- [ ] [Create a new merge request](https://docs.gitlab.com/ee/user/project/merge_requests/creating_merge_requests.html)
|
||||||
|
- [ ] [Automatically close issues from merge requests](https://docs.gitlab.com/ee/user/project/issues/managing_issues.html#closing-issues-automatically)
|
||||||
|
- [ ] [Enable merge request approvals](https://docs.gitlab.com/ee/user/project/merge_requests/approvals/)
|
||||||
|
- [ ] [Set auto-merge](https://docs.gitlab.com/user/project/merge_requests/auto_merge/)
|
||||||
|
|
||||||
|
## Test and Deploy
|
||||||
|
|
||||||
|
Use the built-in continuous integration in GitLab.
|
||||||
|
|
||||||
|
- [ ] [Get started with GitLab CI/CD](https://docs.gitlab.com/ee/ci/quick_start/)
|
||||||
|
- [ ] [Analyze your code for known vulnerabilities with Static Application Security Testing (SAST)](https://docs.gitlab.com/ee/user/application_security/sast/)
|
||||||
|
- [ ] [Deploy to Kubernetes, Amazon EC2, or Amazon ECS using Auto Deploy](https://docs.gitlab.com/ee/topics/autodevops/requirements.html)
|
||||||
|
- [ ] [Use pull-based deployments for improved Kubernetes management](https://docs.gitlab.com/ee/user/clusters/agent/)
|
||||||
|
- [ ] [Set up protected environments](https://docs.gitlab.com/ee/ci/environments/protected_environments.html)
|
||||||
|
|
||||||
|
***
|
||||||
|
|
||||||
|
# Editing this README
|
||||||
|
|
||||||
|
When you're ready to make this README your own, just edit this file and use the handy template below (or feel free to structure it however you want - this is just a starting point!). Thanks to [makeareadme.com](https://www.makeareadme.com/) for this template.
|
||||||
|
|
||||||
|
## Suggestions for a good README
|
||||||
|
|
||||||
|
Every project is different, so consider which of these sections apply to yours. The sections used in the template are suggestions for most open source projects. Also keep in mind that while a README can be too long and detailed, too long is better than too short. If you think your README is too long, consider utilizing another form of documentation rather than cutting out information.
|
||||||
|
|
||||||
|
## Name
|
||||||
|
Choose a self-explaining name for your project.
|
||||||
|
|
||||||
|
## Description
|
||||||
|
Let people know what your project can do specifically. Provide context and add a link to any reference visitors might be unfamiliar with. A list of Features or a Background subsection can also be added here. If there are alternatives to your project, this is a good place to list differentiating factors.
|
||||||
|
|
||||||
|
## Badges
|
||||||
|
On some READMEs, you may see small images that convey metadata, such as whether or not all the tests are passing for the project. You can use Shields to add some to your README. Many services also have instructions for adding a badge.
|
||||||
|
|
||||||
|
## Visuals
|
||||||
|
Depending on what you are making, it can be a good idea to include screenshots or even a video (you'll frequently see GIFs rather than actual videos). Tools like ttygif can help, but check out Asciinema for a more sophisticated method.
|
||||||
|
|
||||||
|
## Installation
|
||||||
|
Within a particular ecosystem, there may be a common way of installing things, such as using Yarn, NuGet, or Homebrew. However, consider the possibility that whoever is reading your README is a novice and would like more guidance. Listing specific steps helps remove ambiguity and gets people to using your project as quickly as possible. If it only runs in a specific context like a particular programming language version or operating system or has dependencies that have to be installed manually, also add a Requirements subsection.
|
||||||
|
|
||||||
|
## Usage
|
||||||
|
Use examples liberally, and show the expected output if you can. It's helpful to have inline the smallest example of usage that you can demonstrate, while providing links to more sophisticated examples if they are too long to reasonably include in the README.
|
||||||
|
|
||||||
|
## Support
|
||||||
|
Tell people where they can go to for help. It can be any combination of an issue tracker, a chat room, an email address, etc.
|
||||||
|
|
||||||
|
## Roadmap
|
||||||
|
If you have ideas for releases in the future, it is a good idea to list them in the README.
|
||||||
|
|
||||||
|
## Contributing
|
||||||
|
State if you are open to contributions and what your requirements are for accepting them.
|
||||||
|
|
||||||
|
For people who want to make changes to your project, it's helpful to have some documentation on how to get started. Perhaps there is a script that they should run or some environment variables that they need to set. Make these steps explicit. These instructions could also be useful to your future self.
|
||||||
|
|
||||||
|
You can also document commands to lint the code or run tests. These steps help to ensure high code quality and reduce the likelihood that the changes inadvertently break something. Having instructions for running tests is especially helpful if it requires external setup, such as starting a Selenium server for testing in a browser.
|
||||||
|
|
||||||
|
## Authors and acknowledgment
|
||||||
|
Show your appreciation to those who have contributed to the project.
|
||||||
|
|
||||||
|
## License
|
||||||
|
For open source projects, say how it is licensed.
|
||||||
|
|
||||||
|
## Project status
|
||||||
|
If you have run out of energy or time for your project, put a note at the top of the README saying that development has slowed down or stopped completely. Someone may choose to fork your project or volunteer to step in as a maintainer or owner, allowing your project to keep going. You can also make an explicit request for maintainers.
|
||||||
|
|
@ -0,0 +1,31 @@
|
||||||
|
function [PfromInjection, DInjection, nUse, aInjection] = injection(Psupply, massBlock, v, rho, Cd, areaBlock, n, mu, eta_motor)
|
||||||
|
% Constants
|
||||||
|
g = 9.81; % Gravity [m/s²]
|
||||||
|
|
||||||
|
% --- Power to lift the block ---
|
||||||
|
P_lift = massBlock * g * v; % Mechanical power needed [W]
|
||||||
|
|
||||||
|
% --- Losses ---
|
||||||
|
P_drag = 0.5 * rho * Cd * areaBlock * abs(v)^3; % Air drag [W] (always opposes motion)
|
||||||
|
P_fric = mu * massBlock * g * abs(v); % Friction [W] (always dissipative)
|
||||||
|
Ploss = P_drag + P_fric;
|
||||||
|
|
||||||
|
% --- Total power required (incl. motor inefficiency) ---
|
||||||
|
P_required = (P_lift + Ploss) / eta_motor; % Electrical input power [W]
|
||||||
|
|
||||||
|
% --- Dynamic loss coefficient (replaces fixed aInjection) ---
|
||||||
|
aInjection = Ploss / P_required; % Ratio of losses to total input
|
||||||
|
|
||||||
|
% --- Dissipation and output ---
|
||||||
|
DInjection = aInjection * P_required; % Power lost [W]
|
||||||
|
PfromInjection = P_required - DInjection; % Net power after losses [W]
|
||||||
|
|
||||||
|
% --- Tower allocation logic ---
|
||||||
|
if Psupply == 0
|
||||||
|
nUse = 0; % No towers active
|
||||||
|
elseif Psupply <= n * P_required
|
||||||
|
nUse = ceil(Psupply / P_required); % Towers needed (rounded up)
|
||||||
|
else
|
||||||
|
nUse = n; % All towers active
|
||||||
|
end
|
||||||
|
end
|
||||||
|
|
@ -0,0 +1,2 @@
|
||||||
|
<?xml version='1.0' encoding='UTF-8'?>
|
||||||
|
<Info/>
|
||||||
|
|
@ -0,0 +1,2 @@
|
||||||
|
<?xml version='1.0' encoding='UTF-8'?>
|
||||||
|
<Info location="1" type="DIR_SIGNIFIER"/>
|
||||||
|
|
@ -0,0 +1,2 @@
|
||||||
|
<?xml version='1.0' encoding='UTF-8'?>
|
||||||
|
<Info/>
|
||||||
|
|
@ -0,0 +1,2 @@
|
||||||
|
<?xml version='1.0' encoding='UTF-8'?>
|
||||||
|
<Info location="SolarExample_supply.csv" type="File"/>
|
||||||
|
|
@ -0,0 +1,2 @@
|
||||||
|
<?xml version='1.0' encoding='UTF-8'?>
|
||||||
|
<Info/>
|
||||||
|
|
@ -0,0 +1,2 @@
|
||||||
|
<?xml version='1.0' encoding='UTF-8'?>
|
||||||
|
<Info location="CycleExample_demand.csv" type="File"/>
|
||||||
|
|
@ -0,0 +1,2 @@
|
||||||
|
<?xml version='1.0' encoding='UTF-8'?>
|
||||||
|
<Info/>
|
||||||
|
|
@ -0,0 +1,2 @@
|
||||||
|
<?xml version='1.0' encoding='UTF-8'?>
|
||||||
|
<Info location="CycleExample_supply.csv" type="File"/>
|
||||||
Some files were not shown because too many files have changed in this diff Show more
Loading…
Add table
Reference in a new issue