극저온 LNG 저장탱크에서 벽면 거칠기가 열성층화 및 자가 가압 현상에 미치는 영향에 대한 수치해석

2.1 해석 모델 및 물성치

본 연구에서는 선행 연구 Ferrín and Pérez (2020)에서 제시된 LNG 저장탱크 해석 모델을 기반으로, 동일한 물성, 초기 조건 및 경계 조건을 적용하여 수치해석 기법의 정확성을 검증하였다. 해석 모델은 축 대칭 조건을 가정한 총 부피 0.5m³의 원통형 Type C 저장탱크이며, 내조의 높이는 1.5m, 반지름은 0.3315m이다. 탱크 내부는 초기 시점에 충전율 80%인 경우를 가정하였으며, 해당 모델의 형상은 Fig. 1에 주요 구조와 치수가 제시되어 있다. 초기 온도는 포화 조건 하의 113K로 설정하였고, 이에 상응하는 초기 압력은 1.13bar로 설정하였다. 탱크 외부 열유입에 의해 자연대류 및 증발 현상이 유발된다는 전제 하에, 초기 유동은 정지 상태로 가정되었다. 고체 벽 구조는 진공 단열재인 펄라이트 층과 두 개의 304 스테인리스강 내조와 외조로 구성되며, 각 강재층의 두께는 3mm이고, 단열재 두께는 11mm이다. 펄라이트는 낮은 열전도도와 불연성으로 우수한 단열 성능을 가지며, 304 스테인리스강은 극저온에서도 구조적 안정성을 유지할 수 있는 특징을 가진다. 이와 같은 다층 단열 구조의 상세도는 Fig. 2에 제시되어 있다.

Fig. 1

Computational domain and mesh

Fig. 2

Cross section of the solid wall

LNG 저장탱크 내 유체는 다수의 선행 연구 (Ferrín and Pérez, 2020; Roh et al., 2013; Saleem et al., 2018)를 참고하여 메탄 단일 조성으로 단순화하였다. 이는 LNG의 주성분이 메탄(90% 이상)임을 고려하였으며, 이러한 가정은 내부 유동 및 열성층화 거동의 거시적 구조에는 큰 영향을 미치지 않으나, 혼합조성에 따라 발생할 수 있는 롤오버 현상이나 증발 특성의 세부 차이 등은 모사되지 않는다는 한계가 존재한다 (Kulitsa and Wood, 2018). 액체상은 온도 차가 크지 않아 물성 변화가 미미하므로 밀도, 열전도도, 비열, 점성 계수를 Table 1 같이 일정한 값으로 가정하였다.

Table 1

Properties for liquid and vapor phase

액체의 비압축성을 유지하면서도 자연대류를 묘사하기 위해 밀도 변화를 부력 항에서만 반영하는 부시네스크 근사를 적용하였다. 액체의 체적력은 식 (1)에 제시된 관계식을 따른다.

여기서 열팽창 계수 β = 3.54×10^-3 K^-1, 기준온도 T_ref는 초기 온도인 113K로 설정하였다.

기체는 큰 온도 구배가 발생하므로, 식 (2), (3)에서 점성 계수 μ_v와 열전도도 k_v는 Lemmon et al. (2025)의 데이터를 기반으로 온도의 선형 함수로 표현하였다.

여기서 T는 기체의 절대 온도(K)이다.

기체의 밀도는 상변화 과정에서 발생하는 자가 가압 현상을 정밀하게 반영하기 위해 압력과 온도의 함수로 정의하였으며, 이를 위해 Peng–Robinson 상태방정식을 적용하였다. 해당 방정식은 이상 기체 상태 방정식에 비해 극저온 조건에서 나타나는 비이상적 열역학 거동을 보다 정확하게 모사할 수 있어, LNG 저장탱크 내부의 압력–온도 관계 예측에 유리하다. 해당 방정식은 식 (4)-(8)에 제시하였다.

V_m은 몰부피 (m³/mol), R은 기체 상수이며, 계수 a,b,α,Z는 임계온도(T_c), 임계압력(p_c) 및 이심 인자 (ω_ac = 0.0114)에 기반하여 정의된다.

증발 잠열과 포화 온도는 상평형 조건에서 열역학적 특성이 압력 변화에 따라 변동하므로 각각 압력의 함수로 표현된다. 식 (9), (10)은 Lemmon et al. (2025)의 데이터를 기반으로 한다.

펄라이트는 온도에 따라 열전도 특성이 크게 변화하므로, 저장탱크 고체벽 내부의 열전달을 정밀하게 모사하기 위해 온도의 함수로 고려된다. 해당 열전도도는 Hofmann (2006)의 실험 결과를 기반으로 회귀된 경험식으로 표현된다.

외조는 외부 대기 온도(293K), 내조는 저장 유체와 열적으로 접촉된 상태(113K)이다. 이에 따라 각 강재에는 해당 온도 조건에서의 정적 열 물성치가 부여되었으며, 물성값은 NIST Cryogenics Technology Group (1995)에서 제시한 데이터에 기반하여 Table 2에 정리하였다.

Table 2

Properties for solid materials

2.2 지배 방정식

극저온 LNG 저장탱크 내부의 유동장은 Unsteady Reynolds-Averaged Navier-Stokes(URANS) 방정식을 기반으로 수치적으로 해석되었으며, 다상 유동을 해석을 위해 Volume of Fluid (VOF) 기법을 적용하였다. VOF 기법은 비혼합성 유체 사이의 상경계를 고정된 오일러 격자 상에서 추적하는 기법으로 각 유체의 체적 분율을 계산 셀 내에 정의함으로써 상경계를 추적한다. 각 제어 체적 내에서는 모든 상의 체적 분율의 합이 항상 1이 되도록 한다. 이를 통해 계면을 정밀하게 추적할 수 있다.

질량 보존 방정식:

운동량 보존 방정식 :

에너지 보존 방정식 :

여기서 k_eff는 유효 열전도도이며, S는 상변화로 인한 소스 항이다.

상 간 경계의 추적은 각 상의 체적 분율 변화에 대한 연속 방정식을 통해 수행된다. 본 연구에서는 액상의 체적 분율에 대해 다음과 같은 이송방정식이 적용되었다.

여기서 하첨자 l은 액상을 나타내며 $$ {\dot{m}}_{lv}$$는 액상과 기상 사이의 상변화에 의한 질량 전달 항으로 정의된다. 기상의 체적 분율 α_v는 α_l + α_v = 1의 조건을 만족하므로 별도의 수송 방정식은 고려하지 않는다.

상변화에 의한 증발 질량 전달을 모델링하기 위해 Lee 모델(Ansys Inc, 2021b)을 적용하였으며, 해당 모델은 사용자 정의 함수(user-defined functions, UDF)를 통해 수치적으로 구현하였다.

여기서 γ는 완화 계수(s^-1)이다. 해석 모델의 경우 탱크 내부로의 열유입으로 인한 액체의 온도 상승이 미미하여 비등을 유도하기에 부족하므로 액체-기체 경계면에서의 증발 과정만을 고려하였다. 선행 연구에서 γ = 0.1, 1, 100의 세 가지 값에 대해 민감도 분석이 수행되었으며, 그 결과 장기적으로 BOR 차이는 16% 이내로 수렴하고 기체 압력 상승은 유사하게 나타났다. 이에 본 연구에서는 γ = 0.1을 적용하였다.

액체는 증발하는 과정에서 증발 잠열을 흡수하므로 이에 대한 소스 항을 에너지 보존 방정식에서 고려하였다.

여기서 H_vap는 증발 잠열이다.

고체 영역에서는 다음과 같은 열전달 방정식을 통해 온도 분포를 계산하였다.

여기서 밀도, 비열, 열전도도 값들은 고체 재료의 물성치에 해당한다.

2.3 경계 조건 및 초기 조건

해석 영역은 외벽-외기 경계, 내벽-유체 경계, 중심축 대칭 경계로 구성되며, 각 경계의 물리적 특성에 따라 적절한 조건을 적용하였다. 외벽은 대류 열전달 경계 조건을 따르며, 외기온도는 293K, 대류 열전달 계수는 7W/m²·K로 설정하였다. 고체 벽면을 통한 열전달을 모사하기 위해 다음과 같은 열전달 관계식을 적용하였다.

내벽과 유체가 접하는 경계에는 no-slip 조건을 적용하였으며, 탱크 중심축에 해당하는 경계에는 축 대칭 조건을 부여하였다.

초기 유체 영역에서는 액상과 기상이 모두 정지한 상태를 가정하였으며, 초기 온도는 113K, 초기 압력은 1.13bar로 설정되었다. 초기 내부 액체의 충전율은 80%이며, 난류 변수 K = 0 m²/s², ω = 1s^-1로 부여하였다. 단열 시스템에 해당하는 고체 영역의 초기 온도 분포는 외기 열 경계 조건과 유체 초기 온도인 113K를 내벽 온도로 설정하여 정상 상태 해석을 수행한 뒤, 그 결과를 초기 조건으로 적용하였다.

2.4 수치해석 기법

본 연구의 수치해석은 Ansys Fluent 2021 R1을 기반으로 수행되었다. 압력-속도 결합 알고리즘으로 SIMPLEC(Semi-Implicit Method for Pressure Linked Equation Consistent)를 적용하였다. 운동량 방정식에서의 압력 보간을 위해 PRESTO! 스킴을 적용하였다. 대류항과 확산항의 이산화에는 각각 QUICK 스킴과 Least Squares Cell-Based 스킴을 적용하여 수치 확산을 최소화하고 고차 정확도를 확보하였다. 시간 이산화는 first-order implicit 스킴을 적용하였으며, 액체-기체 경계면의 포착은 Geo-Reconstruct 방법을 적용하여 경계면에서의 수치 확산을 억제하였다. 본 해석 조건에서 산정된 Grashof 수는 2.67×10¹⁰, Rayleigh 수는 5.79×10¹⁰이다. Enayati et al. (2016)은 Rayleigh 수 ≈ 5×10⁹을 고난류의 시작 기준으로 제시한 바 있으며, 본 연구 조건은 이를 상회한다. 따라서 난류 해석에는 k - ω SST RANS 모델을 적용하였다. 이 모델은 자유 유동 영역에서는 k - ϵ 모델을 사용하고 경계층 근처에서는 k - ω 모델을 조합하여 사용하므로 k - ω 모델이 벽면 근처에서 우수한 성능을 발휘하는 장점을 유지할 수 있다. 초기 단계에서는 시간 간격 0.001s로 20단계를 수행하여 수치적 안정성을 확보하였으며, 이후 전 시간 구간에서는 시간 간격 0.01s로 해석을 진행하여 계산 효율성과 해석의 수렴성을 동시에 확보하도록 하였다.

2.5 수치해석 검증

액체-기체 경계면에서는 온도 및 속도의 급격한 구배가 형성되므로, Fig. 1에 제시된 바와 같이 해당 영역에 대해 고밀도 격자 배치하여 상변화 및 유동 특성을 정밀하게 해석할 수 있도록 하였다. 또한, 벽 근처 유동을 보다 정확히 모사하기 위해 벽 인접 첫 번째 셀 중심이 점성저층 내에 위치하도록 구성하였으며, y⁺는 1 미만으로 유지하였다.

본 연구에 적용된 격자 수는 총 42,110개로, 이는 선행 연구(Ferrín and Pérez, 2020)에서 사용된 격자 수와 유사하며, 격자 민감도 분석을 통해 수치적 신뢰성이 확보된 수준이다.

수치해석 검증은 액체 80% 충전 조건에서 증발로 인한 내부 압력 상승을 선행 연구 결과와 비교하였다. Fig. 3에 도시된 결과에서 확인되듯, 본 연구의 해석 결과는 4000s 동안 선행 연구에서 보고된 압력 상승 거동과 매우 유사하며, 최대 오차는 2% 이내로 나타났다. 이는 본 연구에서 사용된 지배방정식, 격자 구성 및 물성치가 선행 해석 결과를 안정적으로 재현함을 입증하며, 수치해석 기법의 타당성과 재현성이 충분히 확보되었음을 나타낸다.

Fig. 3

Comparison of pressure rise with previous study for validation

2.6 표면 거칠기 모델

CFD 기반 해석에서 표면 거칠기의 영향을 반영하기 위해, 난류 벽 함수 영역에서 벽면 근처 속도 분포에 대한 보정을 수행하는 방법이 일반적으로 사용된다. 본 연구에서는 Ansys Fluent 2021 R1에서 제공하는 Standard Law-of-the-Wall Modified for Roughness 모델 (Ansys Inc, 2021a)을 적용하였다. 해당 모델은 난류 유동의 로그 속도 분포에서 같은 기울기를 가지는 절편 감소를 통해 거칠기 효과를 반영하며, 벽 법칙의 형태는 다음과 같다.

여기서 $$ u^{\mathrm{*}}=C_{\mu }^{1/4}{k}^{1/2}$$는 난류 마찰 속도, κ는 폰 카르만 상수, E는 실험적 상수이며, ΔB는 다음과 같이 정의된다.

여기서 f_r은 거칠기 효과로 인해 절편이 이동하는 정도를 정량적으로 나타내는 거칠기 함수이다. ΔB는 Cebeci and Bradshaw (1977)가 제안한 경험식을 바탕으로 무차원 거칠기 높이 $$ K^{+}=\rho K_s{u}^{\mathrm{*}}/\mu $$와 거칠기 상수 C_s에 따라 결정되며, 이는 다음과 같은 세 가지 유동 영역으로 구분된다.

Hydrodynamically smooth regime

Transitional regime

Fully rough regime

여기서 C_s는 거칠기 형상에 따라 상이하게 결정되며 일반적으로 0.5에서 1 사이의 값을 가진다. 본 연구에서는 거칠기 높이에 따른 차이를 비교하기 위해 모든 조건에서 일정한 값인 1을 적용하였다. 거칠기 높이 K_s는 0.01m 및 0.05m의 두 조건으로 설정하였다. 해당 값은 물리적으로 실제 탱크 내벽에 부착이 가능하며, 구조적으로 실현 가능한 범위를 고려하여 선정되었다. 각 조건은 서로 다른 거칠기 수준에서의 유동 저항 변화를 비교 분석할 수 있을 만큼 충분한 차이를 가지며, 해석 결과에 유의미한 영향을 미칠 수 있도록 설계되었다. 아울러, 거칠기 조건 모두 y⁺ > 30을 만족하여 난류 벽 함수가 적용 가능한 영역에 해당하며, 해당 조건 하에 Standard Law-of-the-Wall Modified for Roughness 모델을 통해 표면 거칠기 효과가 안정적으로 반영되도록 하였다.

3.1 자연대류 유동 및 열성층화

저장탱크 벽면을 통한 열유입은 내부 유체의 온도 구배로 인한 밀도 차이를 유발하여 부력에 의해 자연대류 순환 유동을 형성한다. 이러한 유동은 벽면을 따라 상승하고 중심부를 따라 하강하는 전형적인 순환 구조를 보인다. Fig. 4에 제시된 매끄러운 벽 조건의 유동 함수 분포를 보면, 초기 100s 시점에서는 벽면에서 가열된 유체가 부력에 의해 상승하여 자유수면 근처에서 탱크 축 방향으로 흐르고, 이 과정에서 증발로 열을 손실하여 온도가 낮아지고 밀도가 증가한 유체가 탱크 축을 따라 하강하는 단일 순환이 관찰된다. 시간이 경과함에 따라 단일 순환 구조는 유지되면서 자유수면 인근에서 이차 순환이 형성되는데, 이는 벽면 열 경계층과 벌크 유체간의 상호 작용 및 역류 흡입층의 결합에 기인한다. 이와 같은 유동 구조는 Fu et al. (2014)에서도 관찰된 바 있으며 시간이 지남에 따라 그 구조가 더욱 뚜렷해지고 전체 유동은 점차 안정화된다. 이러한 자연대류 순환은 밀도가 높은 저온 유체는 하부로 가라앉고 상대적으로 밀도가 낮은 고온 유체는 상부에 분포하는 열성층화를 유도한다. Fig. 5에서 확인되듯, 초기 단계에서는 벽면 인근과 자유수면 부근에서 국부적인 온도 구배가 형성된다. 자연대류가 지속되는 과정에서 상·하부 간 온도 구배가 확대되고, 1000s 시점에서는 상부 고온층과 하부 저온층이 명확히 분리된 열성층화 구조가 확립된다. 이와 같이 발달한 열성층화는 자유수면 인근에서의 증발을 촉진하고, 그 결과 극저온 액체화물 탱크에서의 BOG 발생과 자가 가압을 증가시키는 요인으로 작용한다.

Fig. 4

Stream function contours of the liquid phase for a smooth wall

Fig. 5

Temperature contours of the liquid phase for a smooth wall

3.2 거칠기에 의한 속도 분포

Fig. 6은 저장탱크 높이 y = 0.7m에서 반지름 방향 r에 따른 수직 속도 분포를 나타내며, 벽면 거칠기가 자연대류 유동 특성에 미치는 영향을 분석한 결과이다. 매끄러운 벽과 거칠기 높이 0.01m, 0.05m의 세 조건에 대해, 1000s, 4000s 시점에서의 유속 분포를 비교하였다. 모든 조건에서 벽면 근처(r/D ≈ 1)에서 최대 유속이 관찰되었으며, 이는 탱크 외벽을 통한 열유입이 해당 위치에서 가장 큰 온도 구배를 형성하고, 이에 따른 부력 효과가 지배적으로 작용했기 때문이다. 반면, 중심축 방향으로 멀어질수록 온도 구배와 부력이 감소하여 유속이 저하되는 양상을 보인다. 1000s 시점에서 벽면 인근 최대 유속은 매끄러운 벽 조건에서 약 1.47×10^-2m/s, K_s = 0.01m 조건에서 1.11×10^-2m/s, K_s = 0.05m 조건에서 0.84×10^-2m/s로 나타났으며, 거칠기 높이가 증가함에 따라 최대 유속은 매끄러운 벽 대비 각각 24.3%, 42.2% 감소하였다. 이러한 감소는 거칠기 요소가 벽면 전단 응력을 증가시켜 벽면 인근 유동 저항이 커지고, 자연대류에 의해 유도되는 상향 흐름을 효과적으로 억제한 결과이다. 4000s 시점에서는 탱크 내부 온도 구배가 완화되면서 부력 기반 자연대류의 구동력이 약화되어 전체 유속이 모든 조건에서 감소하였다. 거칠기에 의한 벽면 근처 고온 유체의 상승 유속 저하는 유지되었으며, 이는 열성층화의 형성을 지연시키고 자유수면으로의 열 전달을 억제하여 열성층화의 심화를 효과적으로 완화하는 메커니즘으로 작용한다.

Fig. 6

Vertical velocity distribution along the r-direction in the liquid phase (y = 0.7m)

3.3 거칠기에 따른 열성층화 억제 효과

Fig. 7은 1000s 시점에서 매끄러운 벽과 거칠기 조건에 따른 저장탱크 내부 액체의 온도 분포를 나타낸다. 모든 조건에서 상부 고온층과 하부 저온층으로 구분되는 열성층화가 형성되었으나, 거칠기 조건에서는 벽면 근처 상승 유속 저하의 영향으로 상부 온도가 낮고 고온층 발달이 억제되는 양상이 나타났다. 특히 거칠기 높이가 증가할수록 이러한 억제 효과가 강화되어 열성층화 강도가 완화되었다.

Fig. 7

Temperature contours at t = 1000 s for smooth and rough wall conditions

Fig. 8은 액체상의 r = 0.15m 위치에서의 수직 온도 분포를 통해 이러한 거동을 정량적으로 뒷받침한다. 1000s 시점에서 자유수면 온도는 매끄러운 벽 조건에서 113.42K로 가장 높았으며, K_s = 0.01m 및 K_s = 0.05m 조건은 각각 113.28K, 113.25K로 나타나, 거칠기 높이 증가에 따라 액체 상부 온도가 점진적으로 감소함을 보여준다. 4000s 시점에서는 상·하부 간 온도 구배 증가와 액체 상부 온도 상승으로 열성층화가 심화되었다. 그러나 거칠기 조건에서는 상부 온도가 매끄러운 벽 대비 여전히 낮게 유지되어, 거칠기가 열성층화를 억제하는 메커니즘이 장기적으로 지속됨을 확인할 수 있다.

Fig. 8

Vertical temperature distribution in the liquid phase at r = 0.15 m for smooth and rough wall conditions (H : filling ratio 80%)

3.4 액체상 탱크 측벽 입열

Fig. 9는 시간에 따른 액체상의 벽면 평균 열유속 변화를 도시한 것이다. 모든 조건에서 평균 열유속은 초기 과도 구간 이후 약 23.5 W/m²로 수렴하였으며, 조건 간 편차는 ±1% 이내로 매우 작았다. 이는 벽면 거칠기 적용 여부와 관계없이 벽면을 통한 액체상으로의 입열은 동일하게 유지됨을 의미한다.

Fig. 9

Average heat flux at the liquid phase side wall for smooth and rough wall conditions

벽면 거칠기는 입열에는 직접적인 영향을 미치지 않으나, 탱크내부의 대류 구조를 변화시켜 자유수면 온도를 낮추는 효과를 유발한다. 이와 같은 온도 저하는 액체 증발 거동에 차이를 유발하는 요인으로 작용하며, 동일한 입열 조건에서도 증발 거동이 달라질 수 있음을 시사한다.

3.5 탱크 내부 압력 상승 비교

극저온 저장탱크에서는 벽면을 통한 열 유입이 제한적이므로, 내부 압력 상승은 비등보다는 기체-액체 경계면에서의 증발이 주된 원인으로 작용한다. 앞선 결과에서 확인된 바와 같이, 거칠기 조건에서는 열성층화의 완화를 통해 자유수면 온도가 매끄러운 벽 조건에 비해 지속적으로 낮게 유지되었으며, 이는 자유수면에서의 증발량 감소로 이어져 탱크 내부 압력 상승을 억제하였다.

Fig. 10의 결과에 따르면, 500s 이전 구간에서는 세 조건 간 압력 상승의 차이가 두드러지지 않았다. 이는 초기 구간에서 열성층화가 충분히 형성되지 않아, 거칠기의 영향이 자유수면 온도 및 증발량 변화에 유의미하게 작용하지 않았기 때문이다. 그러나 1000s 이후 열성층화가 발달하면서 조건별 압력 상승률의 차이가 점차 확대되었다. 매끄러운 벽 조건에서는 상부 고온층의 온도가 상대적으로 높게 유지되고, 이에 따른 증발량 증가로 높은 압력 상승을 보였다. 반면, 거칠기 조건에서는 열성층화가 완화되면서 압력 상승 속도가 상대적으로 낮게 유지되었다.

Fig. 10

Pressure rise in the tank for smooth and rough wall

4000s에서 매끄러운 벽 조건은 1.232bar로 가장 높은 압력을 기록한 반면, K_s = 0.01m와 K_s = 0.05m 조건에서는 각각 1.197 bar 및 1.188bar로 더 낮은 압력을 나타냈다. 이 값은 매끄러운 벽 대비 압력 상승량이 34.3%, 43.1% 낮은 수준으로, 거칠기 높이가 클수록 내부 압력 상승 억제 효과가 뚜렷해짐을 보여준다.

3.6 Boil-off gas 증가량 비교

Fig. 11은 시간에 따른 BOG 증가량 변화를 나타내며, 벽면 거칠기가 탱크 내 증발 거동에 미치는 영향을 비교한 결과이다. 모든 조건에서 초기 200s 구간에는 탱크 내벽과 액체 사이의 큰 온도 구배로 인해 BOG가 급격히 증가하는 양상을 보였다. 시간이 경과함에 따라 매끄러운 벽 조건에서는 액체의 증발이 뚜렷하게 발생하였으나, 거칠기 조건에서는 열성층화 억제에 따라 BOG 발생량이 매끄러운 벽 대비 현저히 낮게 나타났다. 1000s 시점에서 매끄러운 벽 조건의 BOG 질량은 약 0.0016kg, K_s = 0.01m 조건은 약 0.0005kg, K_s = 0.05m 조건은 약 0.00037kg으로 나타났다. 이는 각각 매끄러운 벽 대비 69%, 77% 감소한 수준으로, 거칠기 증가에 따라 증발 억제 효과가 뚜렷하게 향상됨을 확인할 수 있다. 이러한 차이는 앞선 결과에서 확인된 열성층화 완화 및 액체 상부 온도 저하에 기인하며, 벽면 거칠기가 LNG 저장탱크 내 증발 및 자가 가압을 효과적으로 저감할 수 있음을 보여준다.

Fig. 11

Boil-off gas generation in the tank for smooth and rough wall conditions