Porous Media Module

The porous module represents flow resistance and heat transfer in porous materials. Select one model:

Porous models
Input	Use
porousScreens	Pressure loss from porous resistance. Fluid storage is unchanged.
porousLTE	Porous heat transfer when fluid and solid share one temperature. Includes fluid volume fraction and solid heat storage.
porousLTNE	Porous heat transfer when fluid and solid temperatures differ. Solves a separate solid energy equation on the flow mesh.

Momentum Resistance

Load the module in the fluid control file, define a material, and assign it to a region. This example applies isotropic resistance to cells with the volume tag porous_region in the VOG file:

loadModule: porous
{
  // Existing flow setup and equation options

  // Porous resistance
  porousScreens: <
    material=DarcyForchheimer(darcy=400,forchheimer=40),
    screens=[volumeTag(name=porous_region,material=material)]
  >
}

Material names are arbitrary. Each screens entry selects a material by name and a region by volume tag or geometry. The resistance models, region selectors, and source controls below apply to porousScreens. Velocity is bulk (superficial) velocity, measured over the total flow area.

Inviscid cases (transport_model: none) require a nonzero initial velocity.

Darcy–Forchheimer Resistance

DarcyForchheimer combines linear viscous drag and quadratic inertial drag. For scalar coefficients, the force per unit volume is

\[\vec S_p=-\left(\mu d+\tfrac12\rho f\lVert\vec u\rVert\right)\vec u\]

Here \(\rho\) is density, \(\mu\) is dynamic viscosity, and \(\vec u\) is velocity. The inputs darcy and forchheimer specify \(d\) and \(f\). Supply at least one resistance option; all supplied coefficients must be finite.

Resistance coefficients
Option	SI unit	Default	Meaning
`darcy`	1/m²	0	Inverse permeability \(d\geq0\).
`permeability`	m²	Unset	Permeability \(K>0\); sets \(d=1/K\). Mutually exclusive with `darcy`.
`forchheimer`	1/m	0	Inertial resistance \(f\geq0\).

Bare values use SI units; permeability=25 cm*cm equals darcy=400. For the convention \(\mu/K+\rho\beta_F\lVert\vec u\rVert\), use darcy=1/K and forchheimer=2*beta_F.

Region Selection

Use volumeTag(name=porous_region,material=...) to select a named VOG volume tag, as in the setup example. Alternatively, select cells by geometry with circle, cylinder, or box. Each entry in screens names its material with material=...; a list may contain multiple entries.

Circular Screen

circle defines a planar circular screen.

Option	SI unit	Meaning
`radius`	m	Radius of the circular screen.
`center`	m	Center point of the circular screen.
`normal`	—	Normal vector to the plane of the circular screen.

porousScreens: <
  material=DarcyForchheimer(darcy=400,forchheimer=40),
  screens=[circle(radius=1.0, center=[0,0,0], normal=[0,1,0],
    material=material)]
>

Cylindrical Region

cylinder defines a cylindrical porous region.

Option	SI unit	Meaning
`radius`	m	Radius of the cylindrical screen.
`length`	m	Length of the cylindrical screen.
`center`	m	Center point of the cylindrical screen.
`axis`	—	Axis vector along the length of the cylinder.

porousScreens: <
  material=DarcyForchheimer(darcy=400,forchheimer=40),
  screens=[cylinder(radius=1.0, length=1.0, center=[0,0,0],
    axis=[0,1,0], material=material)]
>

Note

For a cylinder aligned with a 2D mesh’s extrusion direction, extend its ends slightly beyond both extrusion planes. Coincident ends can prevent cells from being tagged. For a mesh spanning z = -0.5 m to 0.5 m, use center=[0,0,0], axis=[0,0,1], and length=1.1.

Box Region

box defines a porous region by its center, orientation, and widths.

Option	SI unit	Meaning
`w1`	m	Width along `axis1`.
`w2`	m	Width along `axis2`.
`w3`	m	Width along the cross product of `axis1` and `axis2`.
`center`	m	Center point of the box screen.
`axis1`	—	First axis vector defining one edge of the box.
`axis2`	—	Second axis vector, perpendicular to `axis1`.

porousScreens: <
  material=DarcyForchheimer(darcy=400,forchheimer=40),
  screens=[box(center=[0,0,0], axis1=[1,0,0], axis2=[0,1,0],
    w1=1.0, w2=1.0, w3=1.0, material=material)]
>

ACC Screen Resistance

The Armour–Cannon–Cady model [ArCa1968] [Cady1973] specifies screen resistance from empirical parameters. Its force per unit volume is

\[\vec S_p=-\frac{B}{B_{\mathrm{CFD}}}\frac{Q}{\epsilon^2} \left(\alpha\mu a^2+\frac{\beta\rho\lVert\vec u\rVert}{D_p}\right)\vec u\]

All parameters are required:

ACC parameters
Option	SI unit	Meaning
`a`	1/m	Screen surface area per unit volume.
`B`	m	Physical screen thickness.
`BCFD`	m	Thickness of the modeled porous region.
`Dp`	m	Hydraulic pore diameter.
`Q`	—	Screen tortuosity factor.
`alpha`	—	Viscous resistance constant.
`beta`	—	Inertial resistance constant.
`epsilon`	—	Screen void fraction.

This example represents an 89 µm screen by a 1 mm thick porous box:

porousScreens: <
  mesh1=ACC(a=510235, B=8.9e-5 m, BCFD=0.001 m,
    Dp=5e-6 m, Q=1.3, alpha=3.2, beta=0.19, epsilon=0.245),
  screens=[box(center=[0,0,0], axis1=[1,0,0], axis2=[0,1,0],
    w1=0.001 m, w2=0.1 m, w3=0.1 m, material=mesh1)]
>

NDR Screen Resistance

The numerically determined resistance model (NDR) specifies quadratic resistance:

\[\vec S_p=-\frac{\rho K_q\lVert\vec u\rVert}{2\epsilon^2}\vec u\]

This is equivalent to DarcyForchheimer with darcy=0 and forchheimer=Kq/epsilon^2.

Both parameters are required:

NDR parameters
Option	SI unit	Meaning
`Kq`	1/m	Quadratic resistance coefficient \(K_q\).
`epsilon`	—	Screen void fraction.

Enter Kq as a bare value in inverse meters:

porousScreens: <
  mesh2=NDR(Kq=1742.782,epsilon=0.4),
  screens=[volumeTag(name=porous_region,material=mesh2)]
>

Directional Resistance

For directional resistance, supply three values for each coefficient, one per material axis. A scalar applies equally along all three axes. The force is

\[\vec S_p=-\left(\mu\mathbf D+\tfrac12\rho\lVert\vec u\rVert\mathbf F\right)\vec u\]

The entries of darcy and forchheimer are the principal values of \(\mathbf D\) and \(\mathbf F\). Both tensors share the material axes. Supply axis1 and axis2 together as finite, nonzero, perpendicular vectors. They are normalized; their cross product defines the third axis. The defaults are Cartesian x and y. Material axes are independent of the screen geometry’s axes.

Unequal principal coefficients require linearSolver=porousSGS in momentumEquationOptions and implicit resistance treatment, for stationary Cartesian incompressible flow. This example sets the first material axis at 30 degrees to x and includes the required solver settings:

// Flow properties
flowCompressibility: incompressible

// Time-stepping and fluxes
timeIntegrator: BDF
inviscidFlux: SOU

// Momentum and pressure-correction equations
momentumEquationOptions: <linearSolver=porousSGS>
momentumInterpolationOptions: <form=old>
pressureBasedMethod: SIMPLE

// Porous resistance
porousScreens: <
  material=DarcyForchheimer(darcy=[500,50000,500],forchheimer=[2,10,2],
    axis1=[0.8660254,0.5,0],axis2=[-0.5,0.8660254,0]),
  screens=[volumeTag(name=porous_cylinder,material=material)]
>

FOU may replace SOU. The default resistance treatment is implicit and conservative; standard assembly is also supported (see Source Treatment). The conservativePressure interpolation modifier is optional. Adjust relaxationFactor and maxIterations in momentumEquationOptions as needed for convergence.

Translational periodic boundaries are supported. Rotating frames, rotational periodicity, moving meshes, and interface boundaries are unsupported. With twoDimensions, resistance must not couple in-plane and out-of-plane velocity components. These restrictions also apply when porousSGS is selected for an isotropic material.

Source Treatment

porositySourceOptions selects the resistance treatment and spatial assembly. The defaults are implicit and conservative.

Source options
Option	Effect
`implicit`	Add resistance to the momentum diagonal.
`explicit`	Evaluate resistance using the current iterate and place it on the momentum right-hand side.
`conservative`	Assemble face contributions into adjacent cells.
`standard`	Evaluate the force at each cell center.

Choose implicit or explicit and conservative or standard; omitted choices retain their defaults.

See Porous Source Treatment in the Theory Guide for the source-split diagram and derivation.

Fitting Pressure-Loss Measurements

fit_porous_resistance fits nonnegative darcy and forchheimer coefficients. It requires Python, NumPy, and a CSV file with columns velocity_m_s and pressure_drop_pa. Include at least two distinct nonzero speed magnitudes. Run the installed utility:

fit_porous_resistance measurements.csv \
  --length 0.2 --density 1.0 --viscosity 0.01 --output fit.json

--length is the porous-layer thickness \(L\); --density and --viscosity are fluid density and dynamic viscosity. Use SI units, superficial velocity, and \(\Delta p=p(0)-p(L)\), with positive velocity from 0 to \(L\). Exclude entrance, exit, and pipe losses. Fit each principal direction separately and specify material axes independently.

The command prints a material definition and writes fitted coefficients, measured and predicted pressure losses, and pressure-error estimates to fit.json.

Thermal Model Setup

LTE and LTNE use one isotropic material. Set volumeTag to a VOG volume tag, or omit it to fill the mesh. An empty selection is an error. Cells outside the material retain ordinary fluid properties.

Use the fluid mesh, boundary conditions, initial conditions, and transport properties of a single-species compressible case. The examples below show the porous material and equation settings to add to that control file; retain its grid, fluid conditions, transport, and time-step entries.

Shared LTE and LTNE requirements
Setting	Required value
Fluid model	Single-species `chemEOS`.
`flowRegime` / `flowCompressibility`	`laminar` / `compressible`.
`pressureBasedMethod`	`SIMPLE` or `SIMPLEC` with the default compressible mass form of pressure correction.
`timeIntegrator`	`BDF`. `SIMPLEC` requires `CNBDFBlend: 0`.
`inviscidFlux`	`FOU` or `SOU`.
`gridCoordinates`	`cartesian` with a stationary mesh.
`momentumInterpolationOptions`	`form=old`.
`energyEquationOptions`	`form=totalEnthalpy`.
`drhodt_predictor`	`zero`.

Species transport, gravity, rotating frames, periodic boundaries, conditions at mesh interfaces, farfield, totalPressureInlet, fixedMassOutlet, and double flux are unsupported. Do not set numCorrectorIterations. An internal material interface in a conforming mesh is supported; see Porous/Free-Fluid Interface.

Use these conventions for thermal-model inputs and results:

Density is the fluid EOS density. Velocity is intrinsic fluid velocity; superficial velocity is porosity times intrinsic velocity.
Resistance coefficients use superficial velocity. Prescribed mass and volume fluxes and wall heat fluxes use total area.
Constant effective conductivities and heat sources already include phase volume fractions. Do not scale them again.
Solid heat capacity per total volume is (1-porosity)*solidDensity*solidSpecificHeat.

Thermal models accept scalar darcy and forchheimer coefficients. The porousScreens material definitions and geometry selectors do not apply to LTE or LTNE.

Local Thermal Equilibrium

LTE assumes that the fluid and stationary solid share one temperature. This example assigns the material to the porous_region volume tag:

loadModule: porous
{
  // Flow properties
  flowRegime: laminar
  flowCompressibility: compressible

  // Time-stepping
  timeIntegrator: BDF
  CNBDFBlend: 0

  // Momentum and pressure-correction equations
  momentumInterpolationOptions: <form=old,conservativePressure>
  pressureBasedMethod: SIMPLEC

  // Energy equation
  energyEquationOptions: <form=totalEnthalpy,linearSolver=SGS,
    relaxationFactor=.7,maxIterations=20>

  // Porous material
  porousLTE: <volumeTag=porous_region,porosity=.5,
    solidDensity=2500 kg/m/m/m,solidSpecificHeat=800 J/kg/K,
    effectiveConductivity=5 W/m/K,darcy=1e8,forchheimer=1000>

  // Plot variables
  plot_output: t,porousLTEPorosity
}

Set the common initial temperature through initialCondition and thermal boundary conditions through boundary_conditions.

LTE material options
Option	SI unit	Default	Meaning
`porosity`	—	Required	Fluid volume fraction, greater than zero and at most one.
`solidDensity`	kg/m³	Required	Positive density of the solid material.
`solidSpecificHeat`	J/(kg K)	Required	Positive, constant solid specific heat.
`effectiveConductivity`	W/(m K)	Required	Nonnegative conductivity of the combined material, per total area.
`darcy`	1/m²	0	Nonnegative inverse permeability.
`forchheimer`	1/m	0	Nonnegative inertial resistance, using the one-half convention in Darcy–Forchheimer Resistance.
`heatSource`	W/m³	0	Heat supplied per total volume; negative values remove heat.
`volumeTag`	—	Whole mesh	VOG volume tag occupied by the material.

For conduction-dominated cases, use full energy relaxation and enough scalar-solver iterations to converge each time step.

Local Thermal Nonequilibrium

LTNE solves separate fluid and solid energy equations on the flow mesh. The stationary solid stores and conducts heat and exchanges heat with the fluid at \(H_v(T_s-T_f)\) per total volume, positive into the fluid. Use LTNE when the two phases cannot be assumed to share one temperature.

Set linearSolver=porousSGS in energyEquationOptions and use backward Euler (BDF with CNBDFBlend: 0). AMR is unsupported. This example fills the mesh with porous material and starts the solid at 350 K:

loadModule: porous
{
  // Flow properties
  flowRegime: laminar
  flowCompressibility: compressible

  // Time-stepping
  timeIntegrator: BDF
  CNBDFBlend: 0

  // Momentum and pressure-correction equations
  momentumInterpolationOptions: <form=old,conservativePressure>
  pressureBasedMethod: SIMPLEC

  // Energy equation
  energyEquationOptions: <form=totalEnthalpy,linearSolver=porousSGS,
    relaxationFactor=.7,maxIterations=200,linearSolverTolerance=1e-9>

  // Porous material
  porousLTNE: <porosity=.5,
    solidDensity=2500 kg/m/m/m,solidSpecificHeat=800 J/kg/K,
    solidEffectiveConductivity=5 W/m/K,
    volumetricHeatTransferCoefficient=1e5 W/m/m/m/K,
    initialSolidTemperature=350 K>

  // Plot variables
  plot_output: t,porousSolidTemperature,porousSolidMask
}

LTNE material options
Option	SI unit	Default	Meaning
`porosity`	—	Required	Fluid volume fraction, greater than zero and at most one.
`solidDensity`	kg/m³	Required	Positive density of the solid material.
`solidSpecificHeat`	J/(kg K)	Required	Positive, constant solid specific heat.
`solidEffectiveConductivity`	W/(m K)	Required	Nonnegative solid-phase conductivity per total area.
`volumetricHeatTransferCoefficient`	W/(m³ K)	Required	Nonnegative distributed interphase coefficient; zero removes volumetric exchange.
`fluidConductivityModel`	—	`constant` with a specified conductivity; otherwise `porosityScaled`	Fluid conductivity model; see Fluid Conductivity.
`fluidEffectiveConductivity`	W/(m K)	Required for `constant`	Finite, nonnegative fluid-phase conductivity per total area; incompatible with `porosityScaled`.
`fluidHeatSource`, `solidHeatSource`	W/m³	0	Heating of each phase per total volume; negative values remove heat.
`initialSolidTemperature`	K	Fluid initial temperature	Positive solid temperature on a fresh run. Restart uses the saved solid temperature.
`volumeTag`	—	Whole mesh	VOG volume tag occupied by the material.
`solidTermination`	—	Required with `volumeTag`	`insulated` or `coupled()` at internal porous/free-fluid faces; see Porous/Free-Fluid Interface for resistance options.
`darcy`, `forchheimer`	1/m², 1/m	0	Nonnegative resistance coefficients, as for LTE.

With porosity=1, set solid conductivity, interphase exchange, and solid heating to zero. Cells outside the material tag have no solid temperature.

Solid Boundary Conditions

porousSolidBoundaryConditions assigns solid conditions by the same names used in boundary_conditions. Omitted solid conditions are adiabatic:

// Fluid boundary conditions
boundary_conditions: <inlet=subsonicInlet(massFlux=.2 kg/m/m/s,T=300 K),
  outlet=fixedPressureOutlet(p=1e5 Pa),sides=symmetry>

// Solid boundary conditions
porousSolidBoundaryConditions: <inlet=temperature(T=350 K),
  outlet=heatFlux(q=-100 W/m/m),sides=adiabatic()>

Supported solid laws are adiabatic(), temperature(T=...), and heatFlux(q=...). Temperature is in kelvin; heat flux is in W/m² of total area and is positive outward. Negative heat flux supplies heat to the solid. A fluid temperature or heat-flux boundary acts only on the fluid. Assign or split heater power explicitly between the phases.

These conditions apply to external mesh boundaries. Use solidTermination for the solid surface inside a mesh containing both porous and clear-fluid cells.

Fluid Conductivity

Choose fluidConductivityModel according to the available material data:

porosityScaled computes porosity times the fluid conductivity from transport_model at the current fluid state. It includes no tortuosity, thermal dispersion, or turbulent contribution.
constant uses a measured or calibrated bulk coefficient. For example, set fluidConductivityModel=constant,fluidEffectiveConductivity=.1 W/m/K.

porosityScaled supports const_viscosity, sutherland, powerLaw, transportDB, database, and chemkin. Module-defined and tabular transport providers are unsupported with this selection.

Fluid temperature and heat-flux walls use the selected effective conductivity. Nonzero imposed heat flux requires positive conductivity. Zero conductivity is allowed for insulated paths.

Porous/Free-Fluid Interface

Use one conforming mesh containing the porous material and surrounding fluid. Set volumeTag to the porous region and choose solidTermination:

Solid termination
Option	Effect
`insulated`	No heat passes through the solid surface at the material interface.
`coupled()`	The porous solid exchanges heat with the pore fluid and adjacent clear fluid at the interface.

Both options allow fluid to cross the interface and conduct heat. This example enables heat transfer from the solid surface:

porousLTNE: <volumeTag=porous_region,solidTermination=coupled(),
  porosity=.5,solidDensity=2500 kg/m/m/m,solidSpecificHeat=800 J/kg/K,
  solidEffectiveConductivity=5 W/m/K,
  volumetricHeatTransferCoefficient=1e5 W/m/m/m/K,darcy=1e8>

To add thermal resistance at the solid surface, use, for example, coupled(Rwall=.002 m*m*K/W). Rwall is finite and nonnegative, in m² K/W of total interface area, and defaults to zero. Zero resistance does not force neighboring cell temperatures to be equal; conduction through the adjacent material still limits heat transfer.

Zero solid conductivity removes heat transfer from the solid surface. Zero porous-fluid conductivity removes conductive exchange between the pore fluid and clear fluid. Volumetric exchange within porous cells remains controlled by volumetricHeatTransferCoefficient.

coupled() requires porosity < 1 and at least one internal face between the tagged material and clear fluid. It applies to every such face. External solid boundaries retain porousSolidBoundaryConditions. Separate meshes, nonmatching interfaces, radiation, and empirical slip or pressure-jump laws are not provided by this option.

Specify coolant supply at an external fluid inlet using the ordinary boundary_conditions. The solver transports mass and enthalpy through shared internal faces; no interface mass source is needed. Normal mass flux uses total area and remains continuous. Intrinsic velocity can change with porosity. Resolve the interface region and check mesh convergence of local velocities and temperatures. LTE uses its single temperature and combined conductivity at the same material boundary; it has no solidTermination option.

Convergence

In energyEquationOptions, maxIterations limits the coupled linear iterations and linearSolverTolerance sets their relative tolerance. The global linearSolverAbsTolerance sets the absolute tolerance in watts. For conduction-dominated cases, increase maxIterations; full energy relaxation can also improve convergence.

Set the solid nonlinear tolerances in porousLTNE:

Solid convergence options
Option	SI unit	Default	Meaning
`solidAbsoluteTolerance`	W	1e-8	Positive absolute tolerance for a cell’s integrated solid-energy residual.
`solidRelativeTolerance`	—	1e-8	Nonnegative relative tolerance.

Nonlinear acceptance requires the fluid convergence criteria and the solid limit solidAbsoluteTolerance + solidRelativeTolerance * initialResidual. Here initialResidual is the largest absolute cell residual in the solid at the start of the timestep. Reaching a linear or nonlinear iteration limit without convergence is an error. LTNE: log entries report the solid residual and its acceptance ratio; a ratio at most one passes.

Output and Restart

Select fields with plot_output. The ordinary plot_freq and restart_freq settings control plot and checkpoint output. Fluid temperature uses t; the module adds the following fields:

Porous plot variables
Variable	Model	SI unit	Meaning
`screenID`	Screens	—	Screen region identifier.
`porousLTEPorosity`	LTE	—	Fluid volume fraction.
`porousLTECapacity`	LTE	J/(m³ K)	Solid heat capacity per total volume.
`porousLTEConductivity`	LTE	W/(m K)	Combined material conductivity.
`porousSolidTemperature`	LTNE	K	Solid temperature; zero outside the solid.
`porousTemperatureDifference`	LTNE	K	\(T_s-T_f\); zero outside the solid.
`porousInterphaseHeatRate`	LTNE	W	Interphase heat rate per cell, positive into the fluid.
`porousSolidMask`	LTNE	—	One where a solid exists, zero elsewhere. Use to identify valid solid temperatures.
`porousFluidConductivity`	LTNE	W/(m K)	Effective fluid conductivity.

Standard mass and energy volume integrals account for fluid porosity and solid energy in both thermal models. Use the budget files below to compare stored mass and energy with imposed sources and outward boundary flows. Boundary fluid energy flows include enthalpy transport, heat conduction, and viscous work.

Restart

Follow the ordinary Stream restart procedure. Keep LTE and LTNE material properties and region selections unchanged.

LTNE writes solid temperature, the solid-cell mask, and a material record with the Stream checkpoint. Retain these files together; missing solid data or changed material or domain settings are errors. Keep solidTermination and Rwall unchanged. With porosityScaled, transport parameters and property data must also remain unchanged. The MPI rank count may change. An LTE checkpoint cannot initialize an LTNE restart.

LTE Budget

output/porous_lte_budget.dat is written at print_freq:

LTE budget columns
Column	Name	SI unit	Meaning
1	`cycle`	—	Cycle number.
2	`time_s`	s	Simulation time.
3	`fluid_mass_kg`	kg	Fluid mass.
4	`total_energy_J`	J	Combined fluid/solid total energy.
5	`heat_input_W`	W	Imposed bulk heating, positive into the domain.
6	`outward_mass_kg_s`	kg/s	Outward mass flow.
7	`outward_energy_W`	W	Outward total-energy flow.

LTNE Budget

output/porous_ltne_budget.dat records the initial state and every accepted timestep, including the final step, independently of output frequencies:

LTNE budget columns
Column	Name	SI unit	Meaning
1	`cycle`	—	Cycle number.
2	`time_s`	s	Simulation time.
3	`mass_kg`	kg	Fluid mass.
4	`fluid_energy_J`	J	Fluid total energy.
5	`solid_energy_J`	J	Solid energy.
6	`fluid_input_W`	W	Imposed fluid heating, positive into the fluid.
7	`solid_input_W`	W	Imposed solid heating, positive into the solid.
8	`exchange_to_fluid_W`	W	Interphase heat rate, positive into the fluid and out of the solid.
9	`outward_mass_kg_s`	kg/s	Outward mass flow.
10	`outward_fluid_energy_W`	W	Outward fluid energy flow.
11	`outward_solid_energy_W`	W	Outward solid heat flow.
12	`interface_mass_kg_s`	kg/s	Mass flow through coupled interfaces, positive from porous to clear fluid.
13	`interface_fluid_heat_W`	W	Conductive heat transfer between pore fluid and clear fluid, positive into clear fluid.
14	`interface_solid_heat_W`	W	Surface heat transferred from the porous solid to the pore fluid and clear fluid.

Columns 12–14 are zero without coupled interfaces. Add columns 8 and 14 to the fluid energy balance and subtract their sum from the solid balance. Column 13 redistributes energy within the fluid and contributes no net domain heating. These interface quantities are internal transfers; they are not included in the external boundary totals.

References

[ArCa1968]

1. Armour, J. N. Cannon, Fluid flow through woven screens 1968.

[Cady1973]

1. Cady, Study of Thermodynamic Vent and Screen Baffle Integration for Orbital Storage and Transfer of Liquid Hydrogen, 1973.