Liquid Hydrogen Storage and Transportation

This example shows how to model a cryogenic tank by using Simscape™ Fluids™ blocks. Aviation and aerospace applications commonly use liquid hydrogen storage instead of compressed gas storage. Engineering challenges include minimizing liquid boil-off, managing tank internal pressure, and designing a robust storage tank.

Liquid hydrogen is stored at approximately 20 K. Tank pressure can increase due to liquid boil-off caused by heat ingress, hydrogen isomer reaction, and slosh energy during tank transport. To mitigate boil-off, the tank structure typically consists of three layers: a low thermal conductivity inner layer, a vacuum layer, and a low thermal conductivity outer layer. Hydrogen embrittlement at low temperatures influences the material choice and thickness, which contribute to the tank weight and energy density. This example demonstrates verifying the tank mechanical design against the hydrogen boil-off and tank pressure increase.

Define System

Liquid hydrogen is typically stored and transported in tanks. This example simulates a tank fill process and transport conditions. The image shows the cryogenic tank with an inlet port, a liquid outlet port, and a gas vent. A pump fills the hydrogen from a reservoir. During filling, the liquid outlet and vent are closed. As pressure rises, the vent opens to release gas and maintain a specified pressure within the tank. The Two-Phase Fluid Properties (2P) block is parameterized with hydrogen properties by using the twoPhaseFluidTables function.

Initialize the system model:

open_system("LiquidHydrogenStorageAndTransportation")
run("LiquidHydrogenStorageAndTransportationParameters")

Specify Requirements

The tank is filled and transported in the span of a few hours. The tank is near empty when filling starts and the fill liquid ortho-isomer ratio is 0.1. It takes 20 minutes to fill the tank with liquid hydrogen and 1.5 hours to transport by road. The tank then spends the next 22 hours at rest. The key design requirements are:

• The pressure inside the tank must remain below 1 MPa.

• Maximum boil-off rate must be less than 1% mass per day.

SysParams.PressureMax    = 1;     % Maximum tank pressure, MPa
SysParams.MaxBoilOffRate = 1/100; % Maximum boil-off rate

Specify Tank Geometry

In this example, a cylindrical tank is used for storing liquid hydrogen. Specify the tank geometry:

InputParam.TankLen = 4.15; % Tank length, m
InputParam.TankDia = 1.6;  % Tank diameter, m
InputParam.PipeDia = 0.11; % Diameter of pipe connections, m

Update model parameters based on user-defined inputs:

Tank.CrossSectionalArea = 0.25*pi*InputParam.TankDia^2;
Tank.Volume = Tank.CrossSectionalArea*InputParam.TankLen;

The tank is modeled using the Receiver Accumulator (2P) block in Simscape Fluids. The liquid and vapor occupy separate volumes inside the block. As the liquid warms, some of it evaporates and transfers to the vapor volume. Because vapor is less dense, this increases the pressure in the tank, requiring periodic venting to remain below maximum pressure. This process is called boil-off losses.

Model Boil-off Losses

The thermal network connected to port H of the Receiver Accumulator (2P) block models the hydrogen physics related to boil-off. The Boil-Off Model subsystem consists of three mechanisms: slosh model, ambient heat transfer model, and hydrogen isomer conversion model.

Model Heat Exchange with Environment

The impact of environment is modeled using thermal networks for heat conduction, convection, and radiation. The three layers of the tank wall are modeled as the inner wall, the outer wall, and the insulation layer. The tank outer wall gains energy from the environment through convective and radiative heat transfer.

Specify the tank outer vessel properties:

Tank.OuterWall.Thickness = 3e-2; % Thickness, meters
Tank.OuterWall.ThermalK = 5e-2;  % Thermal conductivity, W/m*K
Tank.OuterWall.Density = 2400;   % Material density, kg/cu.m
Tank.OuterWall.SpHeat = 800;     % Heat Capacity, J/kg.K

Specify the tank insulation layer properties:

Tank.Insulation.Thickness = 5e-3; % Thickness, meters
Tank.Insulation.ThermalK = 1e-4;  % Thermal conductivity, W/m*K
Tank.Insulation.Density = 1.25;   % Material density, kg/cu.m
Tank.Insulation.SpHeat = 1;       % Heat Capacity, J/kg.K

Specify the tank inner vessel properties:

Tank.InnerWall.Thickness = 2e-2; % Thickness, meters
Tank.InnerWall.ThermalK = 5e-2;  % Thermal conductivity, W/m*K
Tank.InnerWall.Density = 800;    % Material density, kg/cu.m
Tank.InnerWall.SpHeat = 400;     % Heat Capacity, J/kg.K

There are two pipe connections to the tank: the liquid pipe and the vapor vent pipe. The liquid filling line and liquid withdrawal line in the tank both connect to the liquid pipe. The safety vent line connects to the vapor vent pipe. The pipes are typically not as well shielded and present a heat transfer path into the tank. The figure below shows the heat transfer model schematic:

Specify the vapor vent pipe parameters:

Tank.VentPipe.Area = 0.25*pi*InputParam.PipeDia^2; % Tank pressure vent cross-sectional area, m^2
Tank.VentPipe.Thickness = 0.006;                   % Vent pipe thickness, m
Tank.VentPipe.Length = 0.1;                        % Vent pipe length, m
Tank.VentPipe.ThermalK = 5;                        % Vent pipe thermal conductivity, W/m*K

Specify the liquid tapping pipe parameters:

Tank.LiqPipe.Area = 0.25*pi*InputParam.PipeDia^2; % Tank liquid tapping point cross-sectional area, m^2
Tank.LiqPipe.Thickness = 0.006;                   % Liquid tapping pipe thickness, m
Tank.LiqPipe.Length = 0.1;                        % Liquid tapping pipe length, m
Tank.LiqPipe.ThermalK = 5;                        % Liquid tapping pipe thermal conductivity, W/m*K

The Simscape model of the tank walls consists of three layers: the inner vessel in contact with hydrogen, the insulating or vacuum layer, and the outer vessel in contact with the environment. The pipes conduct heat to all three tank wall layers. The subsystem SegmentedPipeFittings models the three pipes, each of which is segmented based on the tank wall thickness. The pipe metal conducts heat into the tank and to the tank vessel. This is represented in the Simscape diagram below:

Model Isomer Conversion Heat Generation

Hydrogen has two isomers: ortho and para. At room temperature, hydrogen composition is 75% ortho and 25% para isomer. The two isomers have slightly different thermal properties, but almost the same chemical properties. At temperatures at low as 20 K, liquid hydrogen exists predominantly as para isomer. Energy is released when ortho-hydrogen converts to para-hydrogen and energy is absorbed during the reverse process. The uncatalyzed conversion of ortho-hydrogen to para-hydrogen is a slow process. During liquefaction, if the conversion is not complete, there could be heat release during storage from the ortho-hydrogen converting to para equilibrium at low temperatures. Different components may have different isomer ratio and the flow of fluid through tanks, storage, pipes might lead to isomer heat generation. A Simscape custom component models the isomer conversion heat generation.

Specify the ortho-isomer fraction for the liquid in the tank and for the liquid flowing into the tank separately. The following plot shows the specific enthalpy of isomer conversion versus temperature:

figure
plot(Hydrogen.TempVec, Hydrogen.EnthalpyVec, LineWidth = 1)
grid on
xlabel("Temperature (K)")
ylabel("Isomer Conversion Heat (kJ/kg)")
title("Isomer Conversion Heat")

The rate of isomer conversion is modeled using first-order dynamics based on the reaction rate constant $\mathit{k}$. Doing a mass balance of ortho-hydrogen in the tank gives differential equation for the ortho-hydrogen fraction in the tank:

$$\frac{\mathit{d}{\mathit{f}}_{\mathrm{tank}}}{\mathit{dt}}{\mathit{M}}_{\mathrm{liq}}=\mathit{k}\left({\mathit{f}}_{\mathrm{equil}}-{\mathit{f}}_{\mathrm{tank}}\right){\mathit{M}}_{\mathrm{liq}}+{\dot{\mathit{m}}}_{\mathrm{inflow}}\left({\mathit{f}}_{\mathrm{inflow}}-{\mathit{f}}_{\mathrm{tank}}\right)$$

where ${\mathit{f}}_{\mathrm{tank}}$, ${\mathit{f}}_{\mathrm{equil}}$, and ${\mathit{f}}_{\mathrm{inflow}}$ are the mass fractions of ortho-hydrogen in the tank, at equilibrium, and for inflow during the filling process, respectively, ${\mathit{M}}_{\mathrm{liq}}$ is the liquid mass in the tank, and ${\dot{\mathit{m}}}_{\mathrm{inflow}}$ is the liquid mass flow rate into the tank during the filling process.

The rate of heat generated by the conversion is

$$\mathit{Q}=-\mathit{h}\left({\mathit{T}}_{\mathrm{liq}}\right)\mathit{k}\left({\mathit{f}}_{\mathrm{equil}}-{\mathit{f}}_{\mathrm{tank}}\right){\mathit{M}}_{\mathrm{liq}}$$

where $\mathit{h}\left({\mathit{T}}_{\mathrm{liq}}\right)$ is the specific enthalpy of the isomer conversion as a function of the liquid temperature. The Simscape language file LiquidHydrogenIsomerConversion.ssc implements the equations for the custom component.

Model Impact of Slosh

The liquid slosh during transportation increases the kinetic energy of fluid and fluid mixing. This increase leads to a more uniform temperature in the fluid, changes in heat transfer, and heat generation due to the fluid impact on walls. The vehicle transporting the tank moves between time t = 30 minutes and t = 2 hours. After that time, the tank is at rest. This example assumes that the rate of heat addition due to slosh can be derived from experimental measurements:

load LiquidHydrogenSloshImpact.mat
figure
plot(SloshImpact.Time/3600, squeeze(SloshImpact.Data), LineWidth = 1)
grid on
xlabel("Time (hr)")
ylabel("Power (W)")
title("Impact of Slosh")

Evaluate Design

Set the average temperature during the day:

SysParams.AmbientT = 300;   % Ambient temperature, K
SysParams.AmbientTvar = 10; % Ambient temperature variation amplitude during the day, K

Set the tank fill time:

SysParams.TankFillTime = 1200; % Time for which tank is filled, s
SysParams.TankFillRate = 0.4;  % Tank fill rate, kg/s

Set the initial ortho-isomer fraction in the tank:

Hydrogen.OrthoInit = 0.1; % Initial ortho-isomer mass fraction in tank (post liquefaction)

Set the ortho-isomer fraction in the incoming liquid flow:

Hydrogen.OrthoInflow = 0.13; % Ortho-isomer mass fraction in incoming flow (post liquefaction)

Enable slosh modeling:

SysParams.SloshModelOpt = 1; % Set 1/0 to model/not-model slosh impact

When fluid flow through a two-phase fluid node changes direction, the internal energy advected through the node changes between the former and new upstream values. This change in internal energy is smoothed for numerical robustness, but the smoothing manifests as a slight heat transfer through the node analogous to conduction. Because the heat transfer rate for boil-off loss modeling is already very low, the additional heat transfer due to the numerical smoothing effect may not be negligible. Therefore, the smoothing parameter located in the Two-Phase Fluid Properties (2P) block is reduced from the default value of 0.01 to 1e-4:

SysParams.SmoothRevThres = 1e-4; % Dynamic pressure threshold for smoothed flow reversal, Pa

Run simulation for tank filling, subsequent transport (1.5 hrs), and storage (~22 hrs) for a full day.:

simResult = sim("LiquidHydrogenStorageAndTransportation", "StopTime", "3600*24");
data = simResult.logsout.extractTimetable;
time = hours(data.Time);
time_s = seconds(data.Time);

Plot Results

Plot the heat ingress due to different sources:

result.totalTime = time_s(end) - time_s(1);
result.heatGenMech = ["Isomer Reaction"; "Heat Ingress - Walls"; "Heat Ingress - Pipes"; "Slosh Effect"];
result.avgPow = [
    trapz(time_s, data.("Heat.isomer"));
    trapz(time_s, data.("Heat.vesselWalls"));
    trapz(time_s, data.("Heat.pipeFittings"));
    trapz(time_s, data.("Heat.slosh"))] / result.totalTime;
result.avgPowPercent = round(result.avgPow*100/sum(result.avgPow), 1);
result.avgPowPlot = table(result.heatGenMech, result.avgPow, result.avgPowPercent);
result.avgPowPlot.Properties.VariableNames = ["Mechanism", "Power Added (W)", "% Contribution"];
figure
piechart(result.avgPow, result.heatGenMech)
title("Average Power Dissipated for Entire Duration")

Display the data for heat ingress due to different sources:

disp(result.avgPowPlot)

          Mechanism           Power Added (W)    % Contribution
    ______________________    _______________    ______________

    "Isomer Reaction"             256.36              34.4     
    "Heat Ingress - Walls"        293.12              39.3     
    "Heat Ingress - Pipes"        194.99              26.2     
    "Slosh Effect"                 1.154               0.2     

Plot the tank pressure changes:

figure
plot(time, data.("pressure"), LineWidth = 1)
hold on
plot([time(1) time(end)], [1 1]*SysParams.PressureMax, "k--")
hold off
grid on
xlabel("Time (hr)")
ylabel("Pressure (MPa)")
title("Tank Pressure")
legend("Pressure in Tank", "Pressure Limit", Location = "best")

Calculate the boil-off rate by integrating the gas outflow rate to obtain the mass lost and integrating the liquid inflow rate to obtain the mass supplied:

result.boilOffRate = trapz(time_s, data.("Outflow.g")) ...
    / (data.("mass")(1) + trapz(time_s, data.("inflow")));
if result.boilOffRate > SysParams.MaxBoilOffRate
    result.dispMessage = "Tank boil-off rate (" + result.boilOffRate*100 + " percent per day) does not meet the requirements (" + SysParams.MaxBoilOffRate*100 + " percent per day)";
else
    result.dispMessage = "Tank boil-off rate (" + result.boilOffRate*100 + " percent per day) meets the requirements (" + SysParams.MaxBoilOffRate*100 + " percent per day)";
end
disp(result.dispMessage)

Tank boil-off rate (0.99194 percent per day) meets the requirements (1 percent per day)

Plot the amount of vented hydrogen gas:

figure
plot(time, cumtrapz(time_s, data.("Outflow.g")), LineWidth = 1)
grid on
xlabel("Time (hr)")
ylabel("Mass (kg)")
title("Cumulative Mass of Vented Gas")

Plot the tank surface temperature:

figure
plot(time, data.("TankTemp.innerVessel"), LineWidth = 1)
hold on
plot(time, data.("TankTemp.insulation"), LineWidth = 1)
plot(time, data.("TankTemp.outerVessel"), LineWidth = 1)
hold off
grid on
xlabel("Time (hr)")
ylabel("Temperature (K)")
title("Tank Temperature")
legend("Inner Vessel", "Insulation", "Outer Vessel", Location = "best")

Plot the fluid temperature inside the tank:

figure
plot(time, data.("H2Temp.liq"), LineWidth = 1)
hold on
plot(time, data.("H2Temp.vap"), LineWidth = 1)
hold off
grid on
xlabel("Time (hr)")
ylabel("Temperature (K)")
title("Fluid Temperature")
legend("Liquid", "Vapor", Location = "best")

Related Topics

• Hydrogen Refueling Station