Main Content

Liquid Hydrogen Storage and Transportation

In this example, you will learn how to use Simscape™ Fluids™ libraries to model a cryogenic tank. In the aviation and aerospace industries, hydrogen storage in a liquid form is preferred over compressed gas storage. The key engineering challenges involve minimizing the liquid boil-off, managing tank internal pressure, and ensuring a safe storage tank design.

Liquid hydrogen is stored at around 20 K. The pressure in tank may increase due to heat ingress and subsequent liquid boil-off. Hydrogen isomer reaction and slosh energy during tank transport can also cause liquid boil-off. To mitigate boil-off, the tank structure typically consists of three layers: a low thermal conductivity inner layer, a vacuum, and a low thermal conductivity outer layer. Hydrogen embrittlement at low temperatures impact the material choice and thickness, which contribute to the overall storage tank weight and energy density. In this example, you will learn how to verify your tank mechanical design against the boil-off requirements and safety due to pressure build-up inside the tank.

Define System

Liquid hydrogen is typically stored and transported in tanks. In this example, you will learn how to simulate a tank fill process and conditions during transport. As shown in the picture above, a cryogenic tank has an inlet port, an outlet port for liquid, and a vent for releasing gaseous hydrogen. 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.

LiquidHydrogenStorageTransport_IMG05.png

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:

LiquidHydrogenStorageTransport_IMG03A.pngLiquidHydrogenStorageTransport_IMG03B.png

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")

Figure contains an axes object. The axes object with title Isomer Conversion Heat, xlabel Temperature (K), ylabel Isomer Conversion Heat (kJ/kg) contains an object of type line.

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

dftankdtMliq=k(fequil-ftank)Mliq+m˙inflow(finflow-ftank)

where ftank, fequil, and finflow are the mass fractions of ortho-hydrogen in the tank, at equilibrium, and for inflow during the filling process, respectively, Mliq is the liquid mass in the tank, and m˙inflow is the liquid mass flow rate into the tank during the filling process.

The rate of heat generated by the conversion is

Q=-h(Tliq)k(fequil-ftank)Mliq

where h(Tliq) 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")

Figure contains an axes object. The axes object with title Impact of Slosh, xlabel Time (hr), ylabel Power (W) contains an object of type line.

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")

Figure contains an object of type piechart. The chart of type piechart has 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.35              34.4     
    "Heat Ingress - Walls"        293.04              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")

Figure contains an axes object. The axes object with title Tank Pressure, xlabel Time (hr), ylabel Pressure (MPa) contains 2 objects of type line. These objects represent Pressure in Tank, Pressure Limit.

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.99514 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")

Figure contains an axes object. The axes object with title Cumulative Mass of Vented Gas, xlabel Time (hr), ylabel Mass (kg) contains an object of type line.

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")

Figure contains an axes object. The axes object with title Tank Temperature, xlabel Time (hr), ylabel Temperature (K) contains 3 objects of type line. These objects represent Inner Vessel, Insulation, Outer Vessel.

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")

Figure contains an axes object. The axes object with title Fluid Temperature, xlabel Time (hr), ylabel Temperature (K) contains 2 objects of type line. These objects represent Liquid, Vapor.

Related Topics