NUMERICAL ANALYSIS FOR A SELF-PRESSURIZED LIQUEFIED HYDROGEN STORAGE TANK USING OPENFOAM

S. Ahn; S. Park

doi:10.6112/kscfe.2026.31.3.174

Preview

Original Article

Journal of Computational Fluids Engineering. 30 September 2026. 174-189
https://doi.org/10.6112/kscfe.2026.31.3.174

NUMERICAL ANALYSIS FOR A SELF-PRESSURIZED LIQUEFIED HYDROGEN STORAGE TANK USING OPENFOAM

OpenFOAM을 활용한 액화수소 저장탱크 자가가압 수치해석

S. Ahn¹

S. Park¹²^†

안 서연¹

박 선호¹²^†

¹Department of Convergence Study on the Ocean Science and Technology, Korea Maritime & Ocean University

²Department of Ocean Engineering, Korea Maritime & Ocean University

¹국립한국해양대학교 해양과학기술전문대학원 해양과학기술융합학과

²국립한국해양대학교 해양공학과

^{*Corresponding Author}

License (open-access, http://creativecommons.org/licenses/by-nc/4.0):

This is an Open-Access article distributed under the terms of the Creative Commons Attribution Non-Commercial License (http://creativecommons.org/licenses/by-nc/4.0) which permits unrestricted noncommercial use, distribution, and reproduction in anymedium, provided the original work is properly cited.

ABSTRACT

As liquefied hydrogen (LH₂) emerges as a key future energy carrier, reliable prediction of its phase-change behavior during self-pressurization is essential for the design and operation of LH₂ storage and receiving terminals. In this study, the phase change within a self-pressurized LH₂ tank was numerically investigated using OpenFOAM. The simulation modeled the self-pressurizing tank tested at the NASA Multipurpose Hydrogen Test Bed (MHTB). Two solvers were compared: an Euler-Euler two-fluid solver, which solves governing equations for each phase separately, and a Volume of Fluid (VOF) solver, which treats both phases with a single set of equations while capturing the interface. To examine turbulence effects, results from the k-ω SST, buoyant k-ε, and LES-WALE models were evaluated against those of laminar flow. For phase change, the Lee and Schrage models were assessed. The Euler-Euler solver slightly overpredicted the tank pressure, whereas the VOF solver demonstrated closer agreement with the experimental data. Among the turbulence treatments, both the laminar assumption and LES captured the experimental trends well, while the Schrage model exhibited superior performance among the phase- change models. Overall, the internal pressure of the MHTB was predicted with a mean absolute error of approximately 1.8 kPa, and key phenomena, including boil-off gas (BOG), boil-off rate (BOR), thermal stratification, and internal flow fields, were accurately reproduced.

Keywords

오픈폼(OpenFOAM)

전산유체역학 (Computational Fluid Dynamics, CFD)

저장탱크(Storage tank)

액화수소(Liquefied hydrogen)

극저온 유체(Cryogenic fluid)

상변화(Phase change)

MAIN

1. 서 론
2. 수치해석 기법
2.1 지배방정식
2.2 상변화 모델
3. 해석 모델 및 방법론
3.1 해석 모델 및 경계조건
3.2 격자 의존도 평가
4. 결과 및 고찰
4.1 2상유동 솔버의 영향
4.2 난류모델의 영향
4.3 상변화 모델의 영향
4.4 결과 분석
5. 결 론

1. 서 론

국제에너지기구(International Energy Agency, IEA)의 Net Zero 2050 시나리오에서는 상업용 수소 및 수소 기반 연료 수요의 20% 이상이 2030년까지 국가 간에 교역될 것으로 전망하였다[1].이러한 대규모 인프라 확장을 뒷받침하기 위해서는 액화수소(Liquefied Hydrogen, LH₂)의 운송과 저장이 이루어지는 인수기지의 안정적 운영 기술을 확보하는 것이 중요하다. 그러나 액화수소의 장기 저장과 이송에 관한 실험 및 수치해석 연구는 아직 충분하지 않으며 점차 확대되는 단계에 있다. 특히 육상 대용량 액화수소 저장탱크 관련 기술은 향후 국제 수소 교역의 핵심을 차지할 것으로 예상되므로 이에 대한 연구가 시급한 실정이다. 액화수소 저장탱크의 자가 가압(self-pressurization) 및 열성층화 거동에 관한 실험 연구는 주로 항공우주 분야를 중심으로 지속적으로 수행되어 왔다. 초기 연구에서는 소형 구형 탱크를 이용하여 열유속과 충전율(fill level)에 따른 자가 가압 특성을 실험적으로 규명하였고[2], 이후 NASA는 K-Site 설비에서 4.89 m³규모의 비행체급 액화수소 탱크를 대상으로 2.0~3.5 W/m² 수준의 저열유속 조건에서 자가 가압 및 열성층화 실험을 수행하였다[3]. 또한, 대형화된 다목적 수소 시험대(Multipurpose Hydrogen Test Bed, MHTB)를 이용하여 가변밀도 다층 단열재(variable density multilayer insulation)와 분사 바(spray bar) 기반 열역학적 벤트 시스템(Thermodynamic Vent System, TVS) 조건에서의 저장 성능을 평가하였다[4].

이러한 실험 데이터를 기반으로 자가 가압 거동을 재현하기 위한 전산유체역학(Computational Fluid Dynamics, CFD) 후속 연구가 진행되었다. Kassemi and Kartuzova[5]는 체적 분율(Volume of Fluid, VOF) 기법을 기반으로 계면 열 및 물질 전달 모델을 결합한 상변화 모델을 개발하였다. 이들은 해당 모델의 결과를 NASA K-Site[3] 및 MHTB 액화수소 실험[4]과 비교하였으며, 상변화 계수와 계면 및 기상 영역의 난류가 압력 및 온도 예측에 미치는 영향을 분석하였다[5]. 최근에는 다양한 난류 모델, 충전율(29, 49, 83%) 및 벽면 열유속 조건에서 K-Site 탱크의 자가 가압 거동을 검증하는 연구로 확장되었다[6]. 그러나 이들 연구는 모두 항공우주선용 액화수소 탱크, 즉 저열유속 및 구조물 질량 대비 저용량 조건을 대상으로 하였으며 후속 수치해석 연구 또한 항공우주 분야에 편중되어 있다.

항공우주 산업 이외에 모빌리티 분야에서도 연료용 액화수소 저장탱크에 관한 설계 및 수치해석 연구가 수행되어 왔다. 대표적으로 대형 트럭용 액화수소 연료탱크의 구조 및 열해석 기반 설계 절차가 제안되었으며[7], 육상 차량용 액화수소 연료탱크의 단열 성능과 열적 자립 시간(thermal autonomy)을 규정하는 국제 표준(ISO 13985)이 마련되어 있다[8]. 또한 민수용 대형 저장을 겨냥하여 능동 냉각을 적용한 액화수소 저장탱크에서의 압력 거동을 수치적으로 분석한 연구도 보고되었다[9]. 그러나 이러한 다양한 액화수소 관련 연구 중에서도 해상 운송 및 육상 인수기지와 같은 대용량 저장 시스템에 대해서는 시나리오 기반 설계 방법론만 제시되었을 뿐, 실험적·수치해석적 연구는 여전히 부족한 실정이다[10]. 특히 육상 액화수소 인수기지를 대상으로 한 실험적 연구는 막대한 비용과 시간 그리고 안전성 문제를 동반하므로 수치해석적 접근이 반드시 선행되어야 한다.

본 연구에서는 신뢰성이 확보된 NASA MHTB 자가 가압 실험 결과[4]를 기준으로 OpenFOAM 기반 수치해석의 예측 성능을 검증하였다. 먼저 다상 유동 해석 기법으로 2상을 각각 계산하는 오일러-오일러(Euler-Euler) 2상 유동 모델과 2상을 단일 유체로 풀고 계면을 추적하는 체적 분율(VOF) 기반 2상 유동 모델을 각각 적용하여 두 기법의 압력 예측 성능을 비교하였다.이어서 선정된 2상 유동 모델에 대해 난류 모델, 상변화 모델, 충전율(fill level) 조건을 변화시키며 실험값과의 정합성을 체계적으로 비교 및 분석하였다. 정확도가 높은 수치해석 기법들을 조합하여 최적 조건을 도출하고, 해당 해석 결과에 KGS 코드에서 규정하는 산정 방식을 적용하여 증발가스 발생량(Boil-Off Gas, BOG) 및 증발률(Boil-Off Rate, BOR)을 도출·검증하였다. 최종적으로 확립된 수치해석 기법 및 결과 도출 체계는 향후 실제 육상 액화수소 저장탱크 및 이송 시스템의 열유동 해석을 위한 기초 자료로 활용될 수 있을 것으로 기대된다.

2. 수치해석 기법

본 연구에서는 저장탱크 내부의 자가 가압과 열성층화를 해석하기 위해 서로 다른 두 가지 다상 유동 해석 솔버(solver)를 검토하였다. 첫 번째 솔버는 각 상마다 별도의 보존식을 푸는 오일러-오일러 2상유체(two-fluid) 방식의 솔버(이하 ‘오일러-오일러 솔버’)인 chtMultiRegionTwoPhaseEulerFoam이고, 다른 솔버는 하나의 방정식을 풀어 2상을 정의하고 체적분율 방정식을 풀어 계면을 정의하는 솔버(이하 ‘VOF 솔버’)인 multiRegionPhaseChangeFoam이다. 오일러-오일러 솔버는 ESI OpenFOAM v2312의 표준 솔버이며, VOF 솔버는 동일 버전 기반으로 TwoPhaseFlow 라이브러리를 컴파일하여 사용하였다. 두 솔버 모두 유체–고체 영역 간 복합열전달(conjugate heat transfer)과 상변화를 고려할 수 있다. 극저온 유체 해석에서 예측 정확도를 좌우하는 또 다른 핵심 요소는 상변화 모델로, 본 연구에서는 다수의 선행 연구에서 사용된 상변화 모델들을 적용 및 비교하여 실험 결과를 가장 잘 재현하는 모델을 채택하였다. NASA MHTB 검증을 위해 사용한 두 종류의 솔버와 두 종류의 상변화 모델에 대한 세부 사항은 아래와 같다.

2.1 지배방정식

2.1.1 오일러-오일러 솔버

오일러-오일러 솔버는 액상(liquid, l)과 기상(gas, g)을 상호 침투하는 연속체로 취급하며, 각 상에 대해 체적 평균된 보존 방정식을 개별적으로 푼 뒤 상간 교환항으로 결합한다[11,12]. 각 상이 고유한 속도장과 보존식을 가진다는 점이 오일러-오일러 솔버의 본질적인 특징이다.

오일러-오일러 솔버의 질량보존 방정식은 식 (1)과 같다. 각 상 $k$ 의 질량보존은 체적분율( $α_{k}$ ) 의 수송방정식으로 표현되며, 상변화에 따른 상간 질량 전달을 포함한다.

(1)

\frac{\partial (α_{k} ρ_{k})}{\partial t} + \nabla \cdot (α_{k} ρ_{k} \vec{u_{k}}) = {\dot{m}}_{k}

여기서 $ρ_{k}$ 는 상 $k$ 의 밀도, $U_{k}$ 는 상 $k$ 의 속도, ${\dot{m}}_{k}$ 는 상변화에 의해 상 $k$ 로 전달되는 단위 체적당 질량 전달율이다. 이때 각 상은 $α_{l} + α_{g} = 1, {\dot{m}}_{l} = - {\dot{m}}_{g}$ 을 만족한다.

이어서 각 상 $k$ 의 운동량 보존 방정식은 식 (2)와 같다.

(2)

\frac{\partial (α_{k} ρ_{k} \vec{u_{k}})}{\partial t} + \nabla \cdot (α_{k} ρ_{k} \vec{u_{k} u_{k}}) = - α_{k} \nabla p + \nabla \cdot (α_{k} τ_{k}) + α_{k} ρ_{k} \vec{g} + \vec{M_{k}}

여기서 $p$ 는 두 상이 공유하는 압력, $τ_{k}$ 는 점성 및 난류 응력을 포함하는 응력텐서, $g$ 는 중력가속도, $M_{k}$ 는 상간 모멘텀 전달항으로서 항력, 양력, 벽면 윤활(wall lubrication), 가상 질량, 난류 소산 등을 포함한다. 본 연구에서 항력은 Ishii–Zuber 상관식[13], 양력은 Tomiyama 모델[14]을 적용하여 모델링하였다.

각 상 $k$ 의 에너지보존은 엔탈피( $h_{k}$ )에 대한 수송방정식으로 표현되며, 이는 식 (3)과 같다.

(3)

\frac{\partial (α_{k} ρ_{k} h_{k})}{\partial t} + \nabla \cdot (α_{k} ρ_{k} \vec{u_{k}} h_{k}) = α_{k} \frac{D p}{D t} + \nabla \cdot (α_{k} k_{e f f, k} \nabla T_{k}) + Q_{k}

여기서 $k_{e f f, k}$ 는 각 상의 유효 열전도도, $T_{k}$ 는 온도, $Q_{k}$ 는 계면 열전달과 상변화 잠열을 포함하는 상간 에너지 교환항이다. 우변의 $D p / D t$ 항은 압력일(pressure work)항으로, 밀폐 상태인 저장탱크의 가압에 따른 증기 압축 가열을 나타낸다.

유체 영역과 달리 외조와 단열재로 구성된 고체 영역에서는 열전도 방정식만을 풀며, 해당 지배 방정식은 식 (4)와 같다.

(4)

ρ_{s} c_{p, s} \frac{\partial T_{s}}{\partial t} = \nabla \cdot (k_{s} \nabla T_{s})

유체와 고체가 맞닿는 벽면은 연성 경계 조건으로 처리되며, 온도 및 열유속 연속 조건을 통해 두 영역의 에너지 방정식이 결합되어 해석된다.

2.1.2 VOF 솔버

체적분율 방정식을 계산하여 계면을 포착하는 multiRegionPhaseChangeFlow은 두상을 단일 속도장( $u$ ), 압력장( $p$ )으로 다루는 단일 유체방식이며, 각 물성은 계면 셀에서 체적분율에 따라 가중 평균된다[15]. 두 상의 계면은 체적분율(𝛼)의 변화를 통해 추적되며, 기하학적 VOF 기법인 isoAdvector[16]와 plicRDF 계면 재구성 기법[17]을 적용하여 정밀하게 포착된다. 각 상의 질량보존으로부터 유도된 체적분율 수송방정식(VOF transport equation)은 식 (5)와 같다[18,19].

(5)

\frac{\partial α}{\partial t} + \nabla \cdot (\vec{u} α) - {\dot{α}}_{p c} = α \nabla \cdot \vec{u} + α (1 - α) (\frac{ψ_{g}}{ρ_{g}} - \frac{ψ_{l}}{ρ_{l}}) \frac{D p}{D t}

여기서 ${\dot{α}}_{p c}$ 는 상변화 소스항이고 𝜓는 압축 계수이며, 우변은 가열 및 압축에 따른 체적 변화를 나타낸다. 이어서 운동량 보존 방정식은 식 (6)과 같이 단일 방정식 형태로 계산되며, 표면장력은 연속표면력(continuum surface force, CSF) 형태의 계면항으로 포함된다[20].

(6)

\frac{\partial (ρ \vec{u})}{\partial t} + \nabla \cdot (ρ \vec{u u}) - \nabla \cdot [μ_{e f f} (\nabla \vec{u} + (\nabla \vec{u})^{T})] = - \nabla p_{r g h} + (\vec{g} \cdot \vec{x}) (ρ_{l} - ρ_{g}) \hat{n_{s}} δ_{s} + \vec{f}

여기서 $p_{r g h} = p - (\vec{g} \cdot \vec{x}) ρ$ 이고, ${\hat{n}}_{s}$ 은 계면 단위 법선이고, $δ_{s}$ 는 계면에서 법선 방향에 대한 적분값이 1인 계면 델타함수이다. 이는 $\hat{n_{s}} δ_{s} \approx \nabla α$ 로 정확히 근사되며, 밀도가 급변하는 계면에 작용하는 부력 항만을 고려하기 위함이다. $f$ 는 표면장력항이다.

에너지 보존 방정식은 식 (7)과 같이 액상과 기상 각각의 온도장에 대해 계산된다.

(7)

\frac{\partial (α_{i} ρ_{i} c_{p, i} T_{i})}{\partial t} + \nabla \cdot (α_{i} ρ_{i} c_{p, i} T_{i} \vec{u}) = \nabla \cdot (λ_{i} \nabla T_{i}) + \dot{q_{p h}} + \dot{q}

여기서 $\dot{q_{p h}}$ 는 상변화 에너지 소스항이며, $\dot{q}$ 는 압력일 효과를 반영하는 항이다.

두 솔버의 지배방정식에서 알 수 있듯이 오일러-오일러 솔버와 VOF 솔버는 다음과 같은 점에서 차이가 있다. 오일러-오일러 솔버는 각 상이 고유한 속도장을 가지며 식 (2)와 같이 각 상에 대한 운동량 보존 방정식을 개별적으로 계산한다. 반면 VOF 솔버는 두 상이 하나의 속도장을 공유하여 식 (6)과 같이 단일 운동량 보존 방정식을 풀며, 밀도 및 점성 등의 물성은 계면 셀에서 체적 분율로 가중 평균한 혼합값을 사용한다. 계면을 표현하는 방식 또한 다른데, 오일러-오일러 솔버는 두 상을 상호 침투하는 연속체로 취급하여 뚜렷한 계면 없이 분산상으로 다루는 반면, VOF 솔버는 체적 분율의 기하학적 재구성을 통해 계면을 명확히 정의한다. 이러한 수치적 특성은 상간 결합 방식의 차이로 이어진다.오일러-오일러 솔버에서는 항력 및 양력 등의 상간 계면력( $M_{k}$ )과 누셀트 수(Nusselt number) 상관식 기반의 계면 열전달을 별도의 모델로 정의하지만, VOF 솔버는 상간 항력 항이 존재하지 않고 단일 유동장을 통해 운동량이 자동으로 연성되며 표면장력은 계면에서만 작용한다. 이러한 뚜렷한 차이점이 존재하므로, 본 연구에서는 두 솔버의 해석 결과를 실험 결과와 비교하여 적합한 솔버를 선정한 후 예측 정확도 검증을 수행하였다.

2.1.3 상태방정식

각 상의 밀도를 압력 및 온도와 연결하기 위해 상태방정식이 필요하며, 각 상에 적용된 상태 방정식은 해석 결과에 직접적인 영향을 미친다. 본 수치해석에서는 기상에는 이상기체 상태방정식을, 액상에 Boussinesq 근사를 적용하였다.

기상은 다음과 같이 이상기체 상태방정식으로 통해 압축성을 반영한다.

(8)

ρ_{g} = \frac{p}{R_{g} T}

여기서 $R_{g} = R_{u} / M_{g}$ 은 기체상수( $R_{u}$ :보편기체상수, $M_{g}$ :몰질량)이고 이상기체 가정은 증기의 압축성을 반영하여 밀폐된 저장탱크의 자가가압을 모사할 수 있도록 한다.

액상은 다음과 같이 Boussinesq 근사를 적용하여 밀도를 온도만의 함수로 취급한다.

(9)

ρ_{l} = ρ_{0} [1 - β (T - T_{0})]

여기서 $ρ_{0}$ 은 기준온도 $T_{0}$ 에서의 기준 밀도이고, 𝛽는 열팽창계수이다. 액상을 비압축으로 간주하되 온도에 따른 미소한 밀도 변화를 부력항에 반영함으로써 ,벽면 가열에 의한 자연대류와 열성층화를 물리적으로 정확히 재현할 수 있다[21,22].

2.2 상변화 모델

극저온 저장탱크의 자가가압은 벽면을 통한 열 침입이 액상 과열을 유발하고, 액상과열이 증발을 일으키며, 이로 인한 증발이 증기압력 상승으로 이어지는 경로로 진행되므로 증발 및 응축에 따른 질량 전달률 산정 방법이 예측 정확도를 결정한다. 본 연구에서는 대표적인 두 가지 상변화 모델인 Lee 모델[23]과 Schrage 모델[24]을 비교 검토하여 실험 결과를 잘 모사하는 모델을 최종적으로 선정하였다.

2.2.1 Lee 모델

Lee 모델은 국소 온도와 포화 온도의 차이에 비례하여 상변화에 따른 질량 전달을 산정하는 모델이다[23]. 본 연구에서는 OpenFOAM v2512에 구현된 Lee 모델의 구조를 참고하여, TwoPhaseFlow 라이브러리 솔버에서 구동될 수 있도록 코드를 직접 추가하여 구현하였다. Lee 모델에서 단위체적당 상변화 질량전달률은 다음과 같이 정의된다.

(10)

{\dot{m}}_{e v a p} = C ρ α \frac{(T - T_{s a t})}{T_{s a t}}, {\dot{m}}_{c o n d} = - C ρ α \frac{(T_{s a t} - T)}{T_{s a t}}

증발(evaporation)은 $C > 0, T > T_{s a t}$ 일 때 계산되며, 응축(condensation)은 $C < 0, T < T_{s a t}$ 일 때 계산된다. 여기서 $T_{s a t}$ 는 포화 온도이고, $C$ 는 질량전달의 강도를 결정하는 경험적 계수이다. $C$ 가 클수록 계면 과열이 빠르게 해소되어 상변화가 촉진된다. Lee 모델의 계수 값의 변화는 압력 상승률 및 증발가스 발생량 발생 예측에 직접적 영향을 미친다.

2.2.2 Schrage 모델

Schrage 모델은 TwoPhaseFlow 라이브러리에서 계면 열저항(interface heat resistance) 모델로 구현되어 있다[24]. 이 모델은 계면이 포화온도( $T_{s a t}$ )에 도달해 있다고 가정하며, 각 셀의 온도와 포화온도의 차로부터 계면 열유속을 산정한다.

(11)

q_{p c} = \frac{T_{g} - T_{sat}}{R_{int}} + \frac{T_{1} - T_{sat}}{R_{int}}

계면 열저항( $R_{int}$ )은 기체 분자 운동론에 기반한 관계식으로 주어지며, 수용계수(accommodation coefficient)인 $C_{a c c}$ 를 포함한다.

(12)

R_{int} = \frac{2 - C_{acc}}{2 C_{acc}} \frac{T_{sat}^{3 / 2} \sqrt{2 π R_{gas}}}{ρ_{g} L^{2}}

여기서 $R_{g a s}$ 는 기체상수, $L$ 은 잠열이다. $C_{a c c} (0 < C_{a c c} \leq 1)$ 가 클수록 $R_{int}$ 이 작아져 상변화가 빠르게 일어난다. 상변화 질량전달률은 $\dot{m} = q_{p c} / L$ 로 계산되며, 계면에 직접 부과 할 때 발생하는 압력 진동을 방지하기 위해 Hardt & Wondra [25]방식을 적용하여 소스항을 계면 인접 영역으로 평활화한다.

3. 해석 모델 및 방법론

3.1 해석 모델 및 경계조건

MHTB는 Fig. 1과 같이 높이 3.05 m, 직경 3.05 m의 원통형 몸체에 상·하부 2:1 타원형 경판(dome)이 결합된 형태의 액화수소 저장탱크로, 내부 체적은 약 18.09 m³이다[4,5]. 단열재는 탱크 표면에 도포된 스프레이 발포 단열재(spray-on foam insulation, SOFI) 위에 45겹의 다층 단열(Multi-layer insulation, MLI) 블랭킷을 적용한 구조이며, MLI는 고온측 경계에서 저온측(SOFI) 경계로 갈수록 층 밀도가 변하는 가변 밀도(variable-density) 방식으로 구성되어 있다[1,6]. 실험에서는 열 유입 조건을 정밀하게 제어하여 측정 오차를 최소화하기 위해 단열재로 감싼 저장탱크를 진공 챔버 내에 설치하고, 온도가 제어되는 크라이오슈라우드(cryoshroud)를 통해 일정한 열 유입을 인가하였다.

본 연구에서는 Table 1에 제시된 각 충전율(filling level) 조건에서의 열 유입량(heat leak)과 초기 압력 및 초기 온도를 그대로 반영하여 해석 조건을 설정하였다.

각 충전율 조건에서의 물성은 NIST Chemistry WebBook에서 포화 상태의 값을 추출하여 사용하였으며, 추출된 물성치는 Table 2와 같다. 이 물성치들은 각 상의 상태 방정식 및 열전달·수송 물성 모델의 입력값으로 사용된다. 경계 조건은 다음과 같이 구성하였다. 압축성 부력 유동을 해석하기 위해 본 연구에서는 정수압 성분을 분리한 수정 압력 변수인 $p_{r g h} = p - ρ (\vec{g} \cdot \vec{x})$ 를 풀며, $p_{r g h}$ 로부터 재구성된 전체 압력 $p$ 가 상태 방정식에 반영되도록 하였다. 이를 위해 압력장( $p$ )에는 zeroGradient 조건을, $p_{r g h}$ 에는 fixedFluxPressure 조건을 적용하여 벽면에서의 압력 구배가 중력 등의 체적력과 균형을 이루고 지정된 속도 플럭스와 일치하도록 설정하였다. 속도장( $u$ )에 대해서는 모든 고체 벽면에 no-slip 조건을 부여하였다. 온도 경계 조건으로는 externalWallHeatFluxTemperature를 적용하여, 총 열 침입량(heat leak)을 탱크 표면적으로 나눈 균일 열유속을 벽면에 부과하였다. 또한, 2차원 축대칭 형상을 고려하여 전면과 후면 경계에는 wedge 경계 조건을 적용하였다.

https://cdn.apub.kr/journalsite/sites/kscfe/2026-031-03/N0500310311/images/jkscfe_2026_313_174_F1.jpg

Fig. 1.

Geometric configuration and dimensions of the NASA MHTB tank

Table 1.

Initial conditions for different filling levels

Filling levels [%]	Heat leak [W]	T [ K]	P [Pa]
90	54.1	20.599	111500
50	51.0	20.71	111600

Table 2.

Thermophysical properties of LH₂ and H₂

	Density [kg/m³]	Specific heat [J/kg·K]	Conductivity [W/(m·K)]
90 Filling H₂	1.4596	6472.9	0.017083
90 Filling LH₂	70.450	9928.2	0.10090
50 Filling H₂	1.4508	6462.2	0.017168
50 Filling LH₂	70.466	9961.3	0.10091

3.2 격자 의존도 평가

격자는 OpenFOAM의 기본 유틸리티인 blockMesh를 사용하여 정렬 육면체(hexahedral) 격자로 생성하였다. 벽면 근처의 유동 및 열전달 특성이 해석 결과에 큰 영향을 미치므로, 벽면 인접 셀이 경계면에 정렬되도록 구성하여 벽면 수직 방향의 급격한 온도 및 속도 구배를 정확히 포착하도록 하였다. 또한 계면에서의 열전달과 상변화는 저장탱크의 자가 가압 거동을 지배하는데, 이때 계면을 통한 열 및 물질 전달과 계면 기하학적 재구성의 정확도는 국소 격자 크기에 매우 민감하다. 따라서 계면 영역에 격자를 조밀하게 집중 배치함으로써 계면에서의 급격한 온도 구배와 상변화율을 정밀하게 모사하였으며, 이러한 계면 집중 격자계는 액화수소 저장탱크의 다상 유동 해석에서 널리 채택되는 방식이다[6,26].

격자 독립성(mesh independence)을 검토하기 위해 Very coarse, Coarse, Medium, Fine의 4단계 격자를 구성하여 비교를 수행하였다. 각 격자의 총 셀 수는 각각 2,848개, 5,502개, 10,386개, 20,211개로, 격자 세분화 단계마다 셀 수가 약 1.9배씩 증가하도록 설정하였다. 각 격자 시스템의 형상 및 분포는 Fig. 2에 나타내었다.

https://cdn.apub.kr/journalsite/sites/kscfe/2026-031-03/N0500310311/images/jkscfe_2026_313_174_F2.jpg

Fig. 2.

Computational meshes

Fig. 3은 4가지 격자 시스템에서 시간에 따라 예측된 (a) 탱크 내부 압력과 (b) 증발가스량을 비교하여 나타낸 것이다. 증발가스량는 액체 질량 감소량으로부터 산정한 것으로 계면에서의 상변화를 직접 반영하는 계면 민감 지표이다. 압력과 함께 증발가스량의 격자 수렴성을 검토함으로써 계면 인접 영역의 해상도에 대한 격자 의존도까지 평가하였다. Fig. 3(a)에서 4가지 격자에 대한 압력예측 결과는 전반적으로 잘 일치하였으며 특히 Medium과 Fine 격자 간 3,000 s 시점의 최종 압력 차이는 약 86 Pa(상대 오차 약 0.07%)에 불과하였다. 또한 Fig. 3(b)에 나타난 증발가스량는 전 구간에 걸쳐 네 격자 간 차이가 크지 않았다. 이를 상세히 확인하기 위해 600–1,200 s 구간을 확대하여 함께 제시하였으며, 확대 구간에서 네 격자의 증발가스량 곡선이 서로 근접함을 확인할 수 있다. 특히 Medium과 Fine 격자의 증발가스량 차이가 작아, 증발가스량 예측 또한 Medium 격자 수준에서 격자 의존도가 충분히 낮음을 확인하였다. 따라서 계산 효율성과 예측 정확도를 함께 고려하여, 약 1만 개의 셀로 구성된 Medium 격자를 본 연구의 기본 격자로 최종 선정하였다. 또한 Medium 메쉬에 대해 벽면 경계층의 해상도는 난류 모델의 예측 결과에 중요한 영향을 미치므로 벽 인접 셀의 무차원 벽거리( $y^{+}$ )를 확인하였다. $y^{+}$ 는 약 56 정도로 30~100 범위에 분포하여 벽함수(wall function)를 사용하는 난류 모델(k–ω SST, buoyant k–ε)의 적용에 적합한 수준임을 확인하였다.

https://cdn.apub.kr/journalsite/sites/kscfe/2026-031-03/N0500310311/images/jkscfe_2026_313_174_F3.jpg

Fig. 3.

Pressure and boil-off gas(BOG) for mesh convergence assessment

4. 결과 및 고찰

NASA MHTB 조건에 대해 오일러-오일러 솔버와 VOF 솔버를 각각 적용하여 수치해석을 수행하고, 탱크 내부의 시간에 따른 압력 변화를 비교·평가하였다. 압력 예측 거동 비교를 통해 적합한 이상유동 솔버를 먼저 선정한 후, 난류 모델에 따른 예측 차이를 비교하고 그 메커니즘을 분석하였다. 저장탱크의 자가 가압은 벽면 가열에 의해 유도되는 부력, 자연 대류 및 열성층화 거동에 의해 지배되며, 이러한 열유동 특성은 적용된 난류 모델에 따라 민감하게 반응하므로 모델별 영향을 체계적으로 검토하였다. 고려된 난류 모델은 k–ω SST, buoyant k–ε, LES WALE(Wall-Adapting Local Eddy-viscosity) 모델[27]이며, 이를 층류 유동 해석 결과와 상호 비교하였다.마지막으로 상변화 모델에 따른 압력 및 온도장 예측 거동을 비교·분석하였다.

4.1 2상유동 솔버의 영향

Fig. 4는 두 솔버를 이용하여 NASA MHTB 내부 압력을 예측한 결과이며, Table 3은 실험값과의 정량적 오차를 비교하여 나타낸 것이다. VOF 솔버의 평균 절대 오차(MAE)는 1.82 kPa, 평균 제곱근 오차(RMSE)는 2.16 kPa로 오일러-오일러 솔버에 비해 낮게 나타났다. 평균 부호 오차(Mean Signed Error, MSE)의 경우, VOF 솔버는 +0.29 kPa로 특정 방향으로의 편향 없이 전 해석 구간에서 실험값을 우수하게 추종하였다. 반면 오일러-오일러 솔버는 평균 부호 오차가 +2.25 kPa로 계산되어 전반적으로 압력을 과대 예측하는 경향을 보였다. 이러한 과대 예측 경향은 시간이 지날수록 누적되어, 20,000 s 시점에서 VOF 솔버는 실험 대비 오차가 -0.49 kPa로 높은 일치도를 보인 반면 오일러-오일러 솔버는 +4.40 kPa의 유의미한 편차를 나타내었다.

https://cdn.apub.kr/journalsite/sites/kscfe/2026-031-03/N0500310311/images/jkscfe_2026_313_174_F4.jpg

Fig. 4.

Pressure for different solvers

Table 3.

Quantitative error of each solver against the NASA MHTB experiment

Solver	MAE [kPa]	RMSE [kPa]	Mean bias [kPa]	Max. error [kPa]	Final error [kPa]
VOF solver	1.82	2.16	0.29	3.42	-0.49
Euler–Euler solver	2.25	2.68	2.25	4.65	+4.40

Table 3에 나타난 예측 정확도의 차이는 두 솔버가 지닌 수치 모델링 특성의 차이로 설명할 수 있다. 자가 가압 저장탱크 내부에서는 액상과 기상이 명확하고 비교적 안정된 자유 수면을 형성하며, 증기 압력 상승은 이 계면을 통한 열 및 물질 전달과 액상 영역의 열성층화 거동에 의해 지배된다. VOF 솔버는 체적 분율의 기하학적 재구성을 통해 계면을 얇고 선명하게 포착하므로, 뚜렷한 계면에서 발생하는 상변화와 열전달 현상을 물리적으로 타당하게 분리하여 해석한다. 반면 오일러-오일러 솔버는 두 상을 상호 침투하는 연속체, 즉 분산상으로 취급하며, 상간 열 및 물질 전달을 항력 계수 및 누셀트 수 상관식과 같은 경험적 모델에 의존한다. 그러나 이러한 모델들은 주로 기포나 액적이 분산된 유동 조건에 맞춰 정식화되었기 때문에, 정지된 평면 계면 조건에는 물리적 적합성이 떨어진다. 그 결과 계면에서의 열 및 물질 전달량이 과대평가되어 내부 압력 상승률을 지속적으로 과대 예측한 것으로 판단된다. 따라서 명확한 자유 계면이 존재하는 극저온 저장탱크의 자가 가압 해석에는 VOF 솔버가 물리적으로 더욱 타당하다고 판단되어, 이후의 난류 모델 및 상변화 모델 민감도 해석에는 VOF 솔버를 채택하였다.

4.2 난류모델의 영향

난류 모델로는 k–ω SST, buoyant k–ε 그리고 LES WALE 모델[27]을 적용하였으며, 그 결과를 층류 유동 해석 결과와 상호 비교하였다. 자가 가압 저장탱크 내부 유동은 낮은 벽면 열 유입에 따른 미약한 부력으로 구동되어 유속이 매우 느리고, 밀도 차에 의한 열성층화가 난류를 억제하므로 실제로는 층류 또는 천이 영역에 머무를 수 있다. 그럼에도 불구하고 난류 확산에 기인한 열전달 및 성층화 촉진 효과를 평가하기 위해 해당 난류 모델들을 선정하였다. 선행 연구에서는 k–ω SST 모델이 난류 확산을 과대평가하여 계면 열전달과 성층화를 과도하게 예측하는 반면, 오히려 층류 가정이 자가 가압률을 실험과 더 우수하게 일치시키는 것으로 보고된 바 있다[5]. 본 연구에서는 선행 연구[5]에서 다루지 않은 buoyant k–ε 모델을 추가로 도입하였다.자가 가압 거동은 부력 구동 자연 대류와 열성층화에 의해 지배되므로, 부력 생성항이 난류 운동에너지에 미치는 영향을 직접 반영하는 buoyant k–ε 모델이 물리적으로 타당한 대안이 될 수 있기 때문이다. 아울러 LES 기법을 도입하여 큰 와류는 직접 해석하고 격자 크기 이하의 와류는 부격자 스케일(Sub-Grid Scale, SGS) 모델로 처리하였으며, 부격자스케일 모델로는 WALE 모델[27]을 적용하였다. WALE 모델은 벽면 근처에서 와점성이 적절히 감쇠하고 순수 전단 및 층류 영역에서 와점성이 0으로 수렴하는 물리적 특성을 지녀, 층류·천이·난류 영역이 공존하는 자연 대류 유동 해석에 적합하다[27]. 이들 4가지 유동 모델에 따른 압력 예측 결과를 NASA MHTB 실험 데이터와 비교하여 Fig. 5에 나타내었다.

https://cdn.apub.kr/journalsite/sites/kscfe/2026-031-03/N0500310311/images/jkscfe_2026_313_174_F5.jpg

Fig. 5.

Pressure for different turbulence models

k–ω SST 모델은 전 해석 구간에 걸쳐 압력을 크게 과소 예측하였으며, buoyant k–ε모델은 해석 초기에는 실험값을 잘 추종하였으나 후반부로 갈수록 압력 상승률이 둔화되며 과소 예측하는 경향을 보였다. 반면 층류 유동과 LES WALE 모델은 실험 압력 상승 곡선을 매우 우수하게 재현하였다.

밀폐 저장탱크의 압력 상승은 계면 인근 및 기상 영역의 열성층화 거동에 의해 결정된다. 이때 난류 혼합이 강할수록 계면과 벽면에서 유입된 열이 액상 벌크 영역으로 과도하게 확산·수송되어 열성층화가 약화되고, 계면 온도 상승이 둔화되면서 압력 상승률이 낮아지게 된다. RANS 기반의 k–ω SST 및 buoyant k–ε 모델은 실제보다 높은 난류 점성을 생성하여 열을 과도하게 혼합시킴으로써 성층화를 과소평가하고, 결과적으로 압력을 낮게 예측한다. 반면 실제 자가 가압 유동은 미약한 부력, 극저속, 밀도 성층화가 복합적으로 작용하여 난류 생성이 강하게 억제된다.따라서 인위적인 난류 확산을 배제하는 층류 가정과, 층류 및 저전단 영역에서 와점성이 0으로 자연스럽게 수렴하는 WALE 모델[27] 기반의 LES 기법이 열성층화 및 압력 상승 거동을 물리적으로 정확하게 포착한다. 이는 층류 유동 가정이 자가 가압률 예측에서 실험과 더 우수한 정합성을 보인다는 선행 연구 결과와도 잘 부합한다[5]. 층류 유동과 LES WALE 모델은 실험값과의 일치도가 유사하였으나, LES 기법은 와류 구조를 직접 해상해야 하므로 계산 비용이 대폭 증가한다. 따라서 계산 효율성과 예측 정확도를 동시에 확보할 수 있는 층류 가정을 이후 상변화 모델 민감도 해석 등 후속 연구의 기본 유동 모델로 채택하였다.

4.3 상변화 모델의 영향

두 상변화 모델인 Lee 모델[23]과 Schrage 모델[24]을 적용하여 탱크 내부의 자가가압 거동을 해석하였다. Fig. 6은 적용된 상변화 모델에 따른 탱크 내부 압력 예측 결과의 비교를 나타낸다.

https://cdn.apub.kr/journalsite/sites/kscfe/2026-031-03/N0500310311/images/jkscfe_2026_313_174_F6.jpg

Fig. 6.

Pressure for different phase change models

Lee 모델[23]의 증발 계수(C)에 대해서는 10⁴과 10¹⁰의 두 가지 값을 고려하였다. 두 조건 모두 시간에 따른 압력 상승률을 크게 과대 예측하였으며, 압력이 시간에 따라 과도하게 선형적으로 상승하는 거동을 보였다. 본 해석에서는 계수 값의 차이가 매우 큼에도 불구하고 두 경우 모두 유사한 수준으로 압력을 과대 예측하였는데, 이는 계수 조정만으로는 실험 결과와 일치하는 해를 도출하는 데 한계가 있음을 보여준다. 이러한 결과는 Lee 모델[23]이 상변화를 계면 영역에 국한하지 않고 과열된 액상 영역 전반에서 국소 과열도에 비례하여 증발을 발생시키기 때문에, 기상 영역으로 유입되는 증기 발생량을 과도하게 산정하여 압력을 선형적으로 과대 예측하기 때문인 것으로 판단된다. 반면 Schrage 모델[24]은 실험 압력 상승 곡선을 매우 우수하게 재현하였다. Schrage 모델[24]은 계면 열저항 이론에 기반하여 계면에서의 열전달을 통해 상변화 질량 유속을 산정하므로, 증발량이 계면으로 유입되는 열유속에 의해 물리적으로 제한된다. 따라서 불필요한 체적 증발을 배제하고 증기 생성을 물리적으로 타당하게 모사함으로써, 실험에서 나타난 자가 가압 거동을 높은 정확도로 포착하였다. 이상의 결과를 바탕으로 본 연구에서는 극저온 액화수소 저장탱크의 자가 가압 수치해석을 위한 최종 상변화 모델로 Schrage 모델[24]을 선정하여 후속 검증 해석을 수행하였다.

4.4 결과 분석

앞선 민감도 해석 결과, 다상 유동 솔버로는 VOF 솔버, 유동 가정은 층류유동, 그리고 상변화 모델로는 Schrage 모델[24]이 각각 실험 결과를 가장 우수하게 재현하는 것으로 확인되었다. 이에 따라 이들 세 가지 기법을 결합한 최적 수치 조건을 수립하여 최종 해석을 수행하였다. 예측된 압력 상승 거동을 NASA MHTB 실험 데이터와 정량적으로 비교하고, 이를 바탕으로 증발가스 발생량과 증발률을 산정하였다. 아울러 탱크 내부의 열성층화 형성 과정과 유동장을 시각화하여 자가 가압 거동의 물리적 메커니즘을 규명하였다.

Fig. 7은 최종 수치 조건을 적용하여 해석한 시간에 따른 탱크 내부 압력 및 증발가스 발생량을 나타낸 것이다. Fig. 7에서 확인할 수 있듯이, 해석을 통해 도출된 시간에 따른 압력 변화는 NASA MHTB 실험 계측값과 매우 높은 일치도를 보였다. 또한 증발량은 시간에 따라 점진적으로 증가하다가, 약 15,000 s 이후에는 상승한 증기 압력에 의해 계면에서 유도된 응축 현상까지 물리적으로 타당하게 재현하였다.

Table 4는 해석 시간에 따라 산정된 기화량(Standard Liters Per Minute, SLPM)과 질량 증발률을 나타낸 것이다. 해석 초기에는 SLPM과 증발률이 모두 양(+)의 값으로 나타나 순 증발(boil-off)이 우세하였으나, 15,000 s 이후에는 음(-)의 값으로 전환되어 순 응축이 지배적으로 발생하는 것으로 나타났다. 이는 밀폐 저장탱크 자가 가압 과정의 열역학적 메커니즘으로 설명할 수 있다. 해석 초기에는 벽면 열 침입에 의해 액상이 가열·증발하면서 기상 영역으로 증기가 유입되므로 양의 순 증발률을 보인다. 그러나 가압이 진행됨에 따라 기상 압력(증기압)이 상승하면, 상승한 압력에 대응하는 포화 온도가 계면 액체 온도보다 높아지게 된다. 그 결과 계면에서 증기 응축이 유도되어 순 증기 발생률이 음의 값으로 역전된다. 이러한 경향은 Fig. 7에 나타난 누적 증발량 곡선이 약 15,000~16,000 s 부근에서 극댓값에 도달한 후 점차 감소하는 거동과 일치한다.

https://cdn.apub.kr/journalsite/sites/kscfe/2026-031-03/N0500310311/images/jkscfe_2026_313_174_F7.jpg

Fig. 7.

Pressure and BOG

Table 4.

Calculated SLPM and BOR over time

Time [sec]	SLPM [L/min]	BOR [%/day]
5000	+1.07	+0.17
10000	+0.22	+0.04
15000	-2.03	-0.33
20000	-3.43	-0.55

아울러 전 해석 구간에서 SLPM 및 증발률의 절대량이 매우 작거나 음수로 유지된다는 점은, 밀폐 자가 가압 조건에서는 외부로 배출(venting)되는 증발가스가 발생하지 않으며 유입된 열량이 주로 기상 영역의 가열 및 압력 상승에 소비됨을 의미한다. 즉, MHTB와 같이 단열 성능이 매우 우수한 극저온 탱크의 밀폐 저장 조건에서는 순 증발량이 극히 제한적이며, 본 수치해석 기법은 이러한 자가 가압 및 계면 응축 거동을 정량적으로 정확히 재현하였다.

Fig. 8은 탱크 내부의 온도 분포와 열성층화 거동을 나타낸 것이다. 액상 영역의 온도는 약 20.5~21 K 범위로 비교적 균일하게 유지되는 반면, 액면 상부의 기상 영역에서는 상단으로 갈수록 온도가 급격히 상승하여 상부 돔 부근에서 약 24~25 K에 이르는 뚜렷한 열성층화가 관찰된다. 즉, 열성층화는 주로 기상 영역에 집중적으로 발생하며, 액상은 얇은 계면 인접 층을 제외하면 거의 등온 상태를 유지한다.

Fig. 9는 액상 및 기상 영역에서의 수직 방향 속도 분포를 나타낸 것이다. 이러한 온도 성층화 구조는 유동장의 수직 속도 분포를 통해 명확히 설명할 수 있다. 액상 영역에서는 벽면 가열에 따른 부력 구동 자연 대류가 활발하게 발달한다. 이러한 자연 대류 순환이 액상 벌크 영역을 지속적으로 혼합하여 균일한 온도를 형성하게 하며, 가열된 벽면을 따라 상승한 온난 유체가 계면 인근으로 열을 수송한다. 반면 기상 영역의 유속은 액상에 비해 전반적으로 매우 느리며, 국소 유동 역시 계면 및 상부 벽면 부근에 국한된다. 이로 인해 기상 영역에서는 대류에 의한 열 혼합이 제한되어 상부로 유입된 열이 전체로 확산되지 못하고 상단부에 축적되며, Fig. 8에서 확인된 강한 기상 열성층화를 유발한다.

https://cdn.apub.kr/journalsite/sites/kscfe/2026-031-03/N0500310311/images/jkscfe_2026_313_174_F8.jpg

Fig. 8.

Temperature contours and vertical thermal stratification in the LH₂ storage tank

https://cdn.apub.kr/journalsite/sites/kscfe/2026-031-03/N0500310311/images/jkscfe_2026_313_174_F9.jpg

Fig. 9.

Vertical velocity contours in the LH₂ storage tank

결과적으로 기상 상부에 축적된 고온 증기층과 계면 인근의 국소 열역학적 상태가 탱크 내부의 증기 압력 상승을 지배하며, 본 수치해석 기법은 이러한 액상 자연 대류와 기상 열성층화의 복합적 상호작용을 물리적으로 타당하게 재현하였다. 이는 앞서 층류 유동이 자가 가압 거동을 가장 정확하게 예측했던 결과와도 일치하는데, 미약한 부력 및 극저속 조건에서 기상 영역의 열성층화가 인위적인 난류 확산 없이 온전히 유지되어야 실제 압력 상승률을 정확히 포착할 수 있기 때문이다.

5. 결 론

본 연구는 육상 액화수소 인수기지 저장탱크의 열유동 및 상변화 해석을 위한 기초 수치해석 플랫폼을 구축하고자 수행되었다. 이를 위해 오픈소스 전산유체역학 라이브러리인 OpenFOAM을 활용하여 NASA MHTB 자가 가압 실험을 대상으로 수치해석 체계를 수립하고 정밀한 검증을 진행하였다. 먼저 오일러-오일러 2상 유동 솔버와 VOF 솔버를 적용하여 밀폐 액화수소 탱크 내부의 압력 변화 거동을 예측하고 실험 결과와 비교·평가하였다. 또한 난류 유동의 영향을 체계적으로 검토하기 위해 k–ω SST, buoyant k–ε, LES WALE 모델[27]을 고려하였으며, 이를 층류 유동의 결과와 상호 비교하였다. 상변화 모델로는 Lee 모델[23]과 Schrage 모델[24]을 적용하여 압력 예측 성능을 평가하고 최적의 수치 모델 조합을 도출하였다. 다만 본 연구는 솔버, 난류 모델, 상변화 모델을 순차적으로 고정·선정하는 실용적 절충안을 택하였으며, 이로 인해 인자 간 상호작용 효과(예: 오일러-오일러 솔버에 대해서는 다른 상변화 모델이 최적일 가능성)를 완전히 탐색하지는 못하였다. 향후 결합 최적화 방식을 통해 이러한 상호작용을 규명할 필요가 있다. 선정된 최적 수치 조건을 통해 NASA MHTB의 시간에 따른 압력 상승 거동을 매우 우수하게 재현하였으며, 압력 이력뿐만 아니라 증발가스 발생량, 증발률, 기상 및 액상의 열성층화, 그리고 탱크 내부 자연 대류 유동장까지 물리적 타당성을 확보하며 정확히 예측함을 확인하였다. 아울러 본 연구는 상변화 해석 기법에 초점을 두어 수행되었으므로, 향후 실질적인 설계 자료 도출을 위해서는 자가가압 및 상변화에 영향을 미칠 수 있는 지지대(support)를 해석에 포함할 필요가 있다. 본 연구를 통해 확립된 수치해석 기법과 검증 결과는 향후 대용량 육상 액화수소 저장탱크의 실제 설계 제원 및 단열 환경을 반영한 신뢰성 높은 열유동 시뮬레이션을 수행하는 데 핵심적인 기초 자료로 활용될 수 있을 것으로 기대된다.

Acknowledgements

본 연구는 기후에너지환경부(MCEE)와 한국에너지기술평가원(KETEP)의 지원을 받아 수행한 연구 과제입니다. (RS-2025-02311765)

References

International Energy Agency (IEA), “Global Hydrogen Review 2023,” IEA, Paris, 2023.

1967, Aydelott, J.C., “Normal Gravity Self-Pressurization of 9-Inch (23 cm) Diameter Spherical Liquid Hydrogen Tankage,” NASA TN D-4171,

1991, Hasan, M.M., Lin, C.S. and Van Dresar, N.T., “Self-Pressurization of a Flightweight Liquid Hydrogen Storage Tank Subjected to Low Heat Flux,” NASA TM-103804.

2003, Hastings, L.J., Flachbart, R.H., Martin, J.J., Hedayat, A., Fazah, M., Lak, T., Nguyen, H. and Bailey, J.W., “Spray Bar Zero-Gravity Vent System for On-Orbit Liquid Hydrogen Storage,” NASA/TM-2003-212926.

2016, Kassemi, M. and Kartuzova, O., “Effect of interfacial turbulence and accommodation coefficient on CFD predictions of pressurization and pressure control in cryogenic storage tank,” Cryogenics, Vol.74, pp.138-153.

10.1016/j.cryogenics.2015.10.018

2025, Kartuzova, O., Kassemi, M. and Hauser, D., “CFD validation of K-site tank self-pressurization under varying fill levels and heat fluxes with different turbulence models,” Cryogenics, Vol.152, 104210,

10.1016/j.cryogenics.2025.104210

2022, Choi, Y., Kim, J., Park, S., Park, H. and Chang, D., “Design and analysis of liquid hydrogen fuel tank for heavy duty truck,” Int. J. Hydrog. Energy, Vol.47, No.32, pp.14687-14702.

10.1016/j.ijhydene.2022.02.210

2006, International Organization for Standardization, “ISO 13985:2006 — Liquid hydrogen — Land vehicle fuel tanks.”

2023, Wan, C., Zhu, S., Shi, C., Bao, S., Zhi, X., Qiu, L. and Wang, K., “Numerical simulation on pressure evolution process of liquid hydrogen storage tank with active cryogenic cooling,” Int. J. Refrig., Vol.150, pp.47-58.

10.1016/j.ijrefrig.2023.01.012

2021, Kim, J., Park, S. et al., “Operation scenario-based design methodology for large-scale storage systems of liquid hydrogen import terminal,” Int. J. Hydrog. Energy, Vol.46, No.79, pp.40262-40277.

10.1016/j.ijhydene.2021.09.218

2002, Weller, H.G., “Derivation, modelling and solution of the conditionally averaged two-phase flow equations,” Technical Report TR/HGW/02, OpenCFD Ltd.

2002, Rusche, H., “Computational Fluid Dynamics of Dispersed Two-Phase Flows at High Phase Fractions,” Ph.D. dissertation, Imperial College London.

1979, Ishii, M. and Zuber, N., “Drag coefficient and relative velocity in bubbly, droplet or particulate flows,” AIChE J., Vol.25, No.5, pp.843-855.

10.1002/aic.690250513

2002, Tomiyama, A., Tamai, H., Žun, I. and Hosokawa, S., “Transverse migration of single bubbles in simple shear flows,” Chem. Eng. Sci., Vol.57, No.11, pp.1849-1858.

10.1016/S0009-2509(02)00085-4

2023, Scheufler, H. and Roenby, J., “TwoPhaseFlow: A framework for developing two phase flow solvers in OpenFOAM,” OpenFOAM J., Vol.3, pp.200-224.

10.51560/ofj.v3.80

2016, Roenby, J., Bredmose, H. and Jasak, H., “A computational method for sharp interface advection,” R. Soc. Open Sci., Vol.3, No.11, 160405.

10.1098/rsos.16040528018619PMC5180117

2019, Scheufler, H. and Roenby, J., “Accurate and efficient surface reconstruction from volume fraction data on general meshes,” J. Comput. Phys., Vol.383, pp.1-23.

10.1016/j.jcp.2019.01.009

2013, Miller, S.T., Jasak, H., Boger, D.A., Paterson, E.G. and Nedungadi, A., “A pressure-based, compressible, two-phase flow finite volume method for underwater explosions,” Comput. Fluids, Vol.87, pp.132-143.

10.1016/j.compfluid.2013.04.002

2014, Jadidi, M., Tembely, M., Moghtadernejad, S. and Dolatabadi, A., “A coupled level set and volume of fluid method with application to compressible two-phase flow,” In Proc. 22nd Annu. Conf. CFD Society of Canada, Toronto, ON, Canada, pp.1-4.

1992, Brackbill, J.U., Kothe, D.B. and Zemach, C., “A continuum method for modeling surface tension,” J. Comput. Phys., Vol.100, No.2, pp.335-354.

10.1016/0021-9991(92)90240-Y

1903, Boussinesq, J., Théorie Analytique de la Chaleur, Vol.2. Paris: Gauthier-Villars.

1960, Spiegel, E.A. and Veronis, G., “On the Boussinesq approximation for a compressible fluid,” Astrophys. J., Vol.131, pp.442-447.

10.1086/146849

1980, Lee, W.H., “A pressure iteration scheme for two-phase flow modeling,” in Multiphase Transport: Fundamentals, Reactor Safety, Applications, T. N. Veziroglu, Ed. Washington, DC: Hemisphere, pp.407-431.

1953, Schrage, R.W., A Theoretical Study of Interphase Mass Transfer. New York: Columbia University Press.

2008, Hardt, S. and Wondra, F., “Evaporation model for interfacial flows based on a continuum-field representation of the source terms,” J. Comput. Phys., Vol.227, No.11, pp.5871-5895.

10.1016/j.jcp.2008.02.020

2026, Fernandes, S., Korsukova, E.V., Ellis, C.D., Ambrose, S. and Eastwick, C.N., “A CFD comparison of interfacial phase change models for boil-off, self-pressurisation and thermal stratification in liquid hydrogen storage tanks,” Int. J. Heat Mass Transf., Vol.256, 128067.

10.1016/j.ijheatmasstransfer.2025.128067

1999, Nicoud, F. and Ducros, F., “Subgrid-scale stress modelling based on the square of the velocity gradient tensor,” Flow Turbul. Combust., Vol.62, No.3, pp.183-200.

10.1023/A:1009995426001

Journal of Computational Fluids Engineering ISSN:1598-6071(Print) 3022-6252(Online) 한국전산유체공학회지

Preview

NUMERICAL ANALYSIS FOR A SELF-PRESSURIZED LIQUEFIED HYDROGEN STORAGE TANK USING OPENFOAM

ABSTRACT

MAIN

(1)

(2)

(3)

(4)

(5)

(6)

(7)

(8)

(9)

(10)

(11)

(12)

Fig. 1.

Geometric configuration and dimensions of the NASA MHTB tank

Table 1.

Initial conditions for different filling levels

Table 2.

Thermophysical properties of LH2 and H2

Fig. 2.

Computational meshes

Fig. 3.

Pressure and boil-off gas(BOG) for mesh convergence assessment

Fig. 4.

Pressure for different solvers

Table 3.

Quantitative error of each solver against the NASA MHTB experiment

Fig. 5.

Pressure for different turbulence models

Fig. 6.

Pressure for different phase change models

Fig. 7.

Pressure and BOG

Table 4.

Calculated SLPM and BOR over time

Fig. 8.

Temperature contours and vertical thermal stratification in the LH2 storage tank

Fig. 9.

Vertical velocity contours in the LH2 storage tank

Acknowledgements

References

Thermophysical properties of LH₂ and H₂

Temperature contours and vertical thermal stratification in the LH₂ storage tank

Vertical velocity contours in the LH₂ storage tank