수직 원관에서의 메탄의 포화 유동 비등에 관한 수치해석적 연구

Ⅰ. 서 론

액화천연가스(Liquefied Natural Gas, LNG)는 SOx를 배출하지 않으며, 기존의 화석연료와 비교하여 친환경 연료로 각광 받고 있다. 따라서 조선, 해운 분야에서도 LNG를 선박의 연료로서 채택하고 있다. 육상에서는 가정용과 산업 분야의 발전용으로 LNG의 이용이 점차 증가하고 있고, LNG 연료추진 선박의 발주도 증가하고 있다(Son and Kim, 2020). LNG는 1 atm에서 –161 ℃의 낮은 비등점을 가지며, 운송을 위해서 액화하여 부피를 감소시킨다. 이때 LNG의 냉각을 위해서는 많은 양의 에너지가 투입되며, 이 에너지는 잠재적인 냉열로서 저장된다고 볼 수 있다(Choi and Lim, 2021). 이러한 냉열은 활용하지 않으면, 소산되어 사라지는 에너지의 형태이며, LNG의 냉열을 활용하는 다양한 열시스템이 개발되고 있다. LNG의 냉열을 활용하는 가장 기본적인 방법 중 하나는 열교환기를 통해서 직접 냉열을 회수하는 방안이다. 이 과정 중에서 LNG는 기화될 수 있는데, 상변화를 동반할 때의 열전달 특성은 단상과는 다르게 복잡한 메커니즘을 통해 나타난다. 이러한 비등 특성은 열교환기를 설계하기 위해서는 반드시 사전에 파악되어야 한다.

유동 비등 현상은 열교환기 및 원자력 발전소 등 다양한 산업 분야에서 발생하는 매우 중요한 물리적 현상이다. 하지만 물리적으로 복잡한 비등 양상의 예측은 많은 실험 연구와 함께 1차원적인 수치해석 연구가 주를 이루었고, 실제 유동현상을 잘 예측하지 못하고 있다. 특히 LNG와 메탄은 낮은 비등점과 폭발의 위험성을 가지고 있으므로 까다로운 실험조건으로 인하여 비등 특성에 대해 수행된 연구는 극히 드물다. 또한 이러한 실험 데이터를 기반으로 수행되는 수치해석적 연구 역시 많지 않은 실정이다. 산업 전반에서 급격히 증가하는 LNG의 사용에 비해서 열시스템을 설계하거나, 열교환기를 설계하기 위한 기초자료는 여전히 부족하다.

LNG의 비등 특성에 대한 연구는 Chen and Shi(2013)에 의해서 진행된 바 있다. 그들은 8 mm와 14 mm 직경의 수직 원관에서 LNG의 비등 현상을 실험을 통해서 연구하였으며, 연구 결과에 따르면 낮은 질량유속과 낮은 건도에서는 열전달 계수는 열유속에 의존하여 열유속이 증가할수록 열전달 계수는 증가하는 것으로 나타났다.

이러한 유동 비등 현상을 모사하기 위한 수치해석적인 연구는 신뢰도를 향상시키기 위하여 현재까지도 계속하여 연구중에 있다. 기본적으로 RPI 벽 비등 모델 및 열유속 분할 모델을 적용한다. Koncar et al.(2004)는 과냉각 비등을 모사하기 위한 수치해석적인 연구를 수행하였으며, 벽 비등 모델을 사용하여 실험결과와 비교하였으며, 그 결과 수치해석 모델의 신뢰성을 확보하였다. 최근 연구인 Setoodeh et al.(2021)은 기존의 벽 비등 모델에서 기포의 슬라이딩의 영향을 추가하여 분석할 수 있는 모델을 제안하였다.

LNG 및 메탄의 비등에 관한 연구는 주로 풀비등(pool boiling)에 집중되어 있는데, 그 이유는 LNG의 운송과정 중 발생하는 boil-off gas(BOG)를 확인하기 위해서이다. 따라서 상대적으로 유동 비등에 관한 연구는 부족한 편이다.

본 연구에서는 상용코드 ANSYS-CFX 18.2를 사용하여, LNG의 주성분인 메탄의 포화 유동 비등에 관한 수치해석적인 연구를 수행하였다. 상변화 모델과 모멘텀 전달 모델들을 사용하여 메탄의 포화 유동 비등 현상을 모사하였으며, 출구에서의 증기 건도를 이론적인 값과 비교하여, 전체적인 계산의 타당성을 입증하였으며, Chen and Shi(2013)의 실험에서 이용할 수 있는 열전달 계수의 실험 데이터와 비교, 분석하였다. 본 연구에서는 메탄의 포화 유동 비등을 모사하기 위한 모델을 분석하여 적용하였으며, 본 연구를 통해서 LNG 및 메탄의 열시스템과 열교환기를 설계할 시, 시뮬레이션을 위한 기초자료로 유용하게 사용될 것으로 판단된다.

Ⅱ. 수치해석 모델

Ⅲ. 수치해석 조건

1. 유동 도메인

본 연구에서는 Chen and Shi(2013)의 시험부 형상과 동일하게 내경 8 mm, 길이 1 m의 수직 상향 원관을 사용하였으며, 이때 바깥 표면에서 열이 공급된다. 수치해석에서는 계산시간의 단축을 위해서 1/4인 90°에 해당하는 영역만 계산하여 symmetry 조건을 적용하였다. <Table 1>은 LNG 터미널별 2017년 평균 가스 조성을 나타낸다(Yoo, 2018). <Table 1>과 같이 LNG의 경우 산지별로 조성이 다르기 때문에, 물성치를 정확하게 적용하기에는 어려운 측면이 있다. 따라서 본 연구에서는 LNG를 주성분인 메탄으로 단순화하여 수치해석에 사용하였다. 0.7 MPa에서의 포화액과 포화증기의 물성치는 NIST REPROP 9.1 프로그램을 이용하였으며, <Table 2>에 정리하였다.

<Table 1>

The average compositions by the different production

<Table 2>

Properties of saturated liquid and vapor at 0.7MPa

난류 모델은 유동해석에 일반적으로 사용되는 k-ε 모델을 적용하였다. 상변화 특성을 모사하기 위해서는 매우 작은 간격의 시간 스텝이 요구된다. 기준 시간 스텝에서 증가시키게 되면 발산하며, 너무 작은 시간 스텝에서는 수렴에 소요되는 시간이 급격히 증가하므로 본 연구에서는 물리적 시간 척도를 반복되는 계산을 통해 최적으로 축소하여 0.001초를 적용하였다. 반복계산은 기본적으로 5,000회 이후 도메인 전체의 증기 체적분율의 변화가 없었으며, 총 20,000번의 반복계산을 수행하여 최종적으로 수렴시켰다.

3. 격자

격자는 ICEM-CFD에 의해서 생성된 Hexahedral 격자계로 구성되었다. O-grid 기법을 사용하여 [Fig. 1]과 같이 일정한 모양의 격자가 적용되었다. 격자 의존성 검토를 위한 노드수는 1,100~40,300개로 설정하였다. k-ε 난류 모델과 벽비등 모델의 특성상 많은 격자를 요구하지 않는다. [Fig. 2]는 관의 길이 방향에 따른 증기의 체적분율을 노드수에 따라서 표시하였다. 입구부분에서는 격자수에 따른 차이가 다소 발생하였지만, 유동 진행방향에 따라서 격자 의존성이 크게 작용하지 않았다. 본 연구에서는 약 40,300개의 노드를 기준으로 계산을 수행하였다.

[Fig. 1]

Flow domain and hexahedral mesh.

[Fig. 2]

Mesh independence with respect to the number of nodes.

Ⅳ. 결과 및 고찰

1. 증기 건도

유동 비등에서 출구 건도 x는 에너지 평형으로부터 다음 식과 같이 이론적으로 계산된다.

$\[ x=\frac{1}{H_{fg}}\left[\frac{q}{\dot{m}}-c_p\left(T_{sat}-T_{in}\right)\right] \]$

(15)

여기서, 각 기호의 의미와 단위는 하기와 같다.

H_fg는 기체상과 액체상의 엔탈피의 차이[J/kg] 이며, $$ \dot{m}$$은 질량유속[kg/m²·s], q는 열유속[W/m²], C_p는 정압비열[J/kg·K], T_sat은 포화온도[K], T_in은 입구온도[K]를 의미한다.

본 연구는 포화 유동 비등으로, 원관의 입구로 포화액 상태로 유입되므로, (T_sat - T_in) = 0 이 된다.

[Fig. 3]은 식 (15)와 수치해석 결과를 비교하여 나타낸 것이다. 같은 열유속이 공급되었을 경우 식 (15)에서 계산되는 것과 같이 낮은 질량유속에서는 더 높은 증기 건도가 나타난다.

[Fig. 3]

Comparisons between the CFD data and theoretical values.

질량유속이 52.36 kg/m²·s일 때, 증기 건도는 약 0.14~0.25까지 계산되었으며, 질량유속이 108.06 kg/m²·s와 153.01 kg/m²·s일 때 증기 건도는 각각 0.07~0.20 그리고 0.04~0.22까지 계산되었다. 대체적으로 수치해석 결과는 식 (15)의 값에 비해서 건도를 약 0.01~0.02정도 낮게 예측하고 있다. 특히, 낮은 열유속일 때 식 (15)와 수치해석 결과의 차이가 크게 발생하였다. 질량유속 153.01 kg/m²·s와 열유속 7,950 W/m² 조건에서 식 (15)의 값과 수치해석 결과는 각각 0.06과 0.04로서 가장 큰 오차를 가진다. 여기에는 여러 가지 요인이 있을 수 있겠지만, 벽 비등 모델의 하위 모델인 기포이탈 빈도, 핵사이트 밀도 그리고 기포 이탈 직경 등의 세부 조정이 요구될 수 있다(Choi et al., 2016).

[Fig. 4]는 질량유속 153.01 kg/m²·s 에서의 길이 방향에 따른 유동단면적에서의 평균 증기 건도를 나타낸 것이다. 길이 방향에 따른 증기 건도의 변화는 대체로 선형적으로 나타나며, 이론적인 경향과 일치한다.

[Fig. 4]

Variations of vapor quality with respect to axial length.

일반적으로 열교환기의 설계에서 대부분의 관심사는 증기 건도 및 기공률(void fraction)의 분포와 열전달 계수의 예측에 있으므로, 본 연구에서는 전체적인 증기 건도와 열전달 계수의 확인을 목표로 하였다. 국소적으로 발생하는 기포의 거동 및 역학은 해석 도메인의 전체적인 질량전달과 출구 증기 건도에는 큰 영향을 미치지 않는 것으로 판단하여 이에 대한 상세한 계산은 수행하지 않았다. 물리적으로, 기포는 핵사이트에서 생성되어, 기포의 이탈빈도, 기포의 이탈직경 등이 주어지면, 벽면에서 이탈하여 포화액과 같이 유동한다. 이때 기포는 주변의 기포와 합쳐지거나 깨지기도 하는데, 작은 기포들이 유동하는 구간을 기포류(bubbly flow)라고 하며, 작은 기포들이 합쳐져서 하나의 큰 기포를 만드는데 이를 슬러그류(slug flow)라고 한다(Choi and Lim, 2015). 정상해석으로 수행되어, 기포의 거동은 모사되지 않았지만, 출구의 증기 건도를 보면 본 해석 조건에서의 결과는 기포류와 슬럭류에 해당되는 것으로 판단된다(Mishima, 1984).

[Fig. 5]는 질량유속 153.01 kg/m²·s와 열유속 31,858.9 W/m² 조건에서 출구 단면에서의 반경방향의 증기의 체적분율 분포를 나타낸 것이다. 벽 윤활 모델을 적용함으로써 벽면에서는 낮은 증기의 체적분율 분포를 확인할 수 있으며, 그 이후 가장 높은 체적분율 분포가 나타나며, 중심으로 갈수록 점차 감소한다. 이러한 체적분율 분포는 기존의 비등 현상에 관한 수치해석적 연구들(Koncar et al., 2004; Setoodeh et al., 2021)의 결과와 유사하다. 비항력 모델은 축방향의 체적분율 분포에는 크게 영향을 미치지 않으며, 반경방향의 체적분율 분포에 영향을 주는 것으로 보고되었다(Krepper et al., 2007).

[Fig. 5]

Variations of void fraction with respect to radial length.

2. 열전달 계수

일반적으로 실험에서의 열전달 계수 h는 식 (16)과 같이 계산된다.

$\[ h=\frac{q}{T_{wall}-T_{sat}} \]$

(16)

벽면의 온도 T_wall은 유체와 접하는 내부 벽면의 온도를 뜻하며, T_sat은 포화온도, q는 열유속을 말한다. 실험에서는 벽면 외부 온도를 측정하여 1차원 열전도 해석을 통해서 계산하였으며, 수치해석에서는 액체의 과열도를 이용하여 열전달 계수를 계산하였다.

[Fig. 6]은 Chen and Shi(2013)의 실험에서의 열전달 계수와 수치해석에 의한 열전달 계수를 비교하여 나타낸 것이다. 실험값은 LNG를 대상으로 수행되었기 때문에 물성치 차이에 의한 오차가 발생할 수 있다.

[Fig. 6]

Comparisons between the CFD and experimental data.

[Fig. 7]은 열유속에 따른 열전달 계수를 비교한 것이다. 수치해석 결과에 의하면, 본 연구의 수치해석 조건에서는 열전달 계수는 질량유속에는 큰 영향을 받지 않으며, 열유속에 크게 의존하는 것으로 나타났다. 이러한 비등 특성은 비등이 시작되어 기포류 및 슬러그류가 지배적인 영역에서 핵비등에 의한 열전달이 주를 이룬다고 볼 수 있다(Choi and Lim, 2016). 이것은 앞서 추측한 유동 패턴에 적절한 열전달 계수의 변화 경향으로 볼 수 있다.

[Fig. 7]

Variations of heat transfer coefficient with respect to heat flux.

본 연구에서 얻어진 메탄의 포화 유동 비등 특성 및 열전달 계수는 차후에 열교환기 설계에 직접적으로 이용될 수 있으며, 각종 열시스템 설계 및 해석의 기초자료로서 활용 가능할 것으로 판단된다.

Ⅴ. 결 론

본 연구에서는 메탄의 포화 유동 비등 현상에 대한 수치해석적인 연구를 수행하였으며, 수치해석 결과는 이론적인 값과 실험 데이터와 비교하여 분석하였으며, 다음과 같은 결론을 얻었다.

1. 출구 증기 건도는 수치해석은 이론값에 비해서 약 0.01~0.02 정도 낮게 예측하였으며, 오차는 0.05%에서 최대 0.26%로 나타났다. 이때 오차는 낮은 열유속 조건에서 크게 나타났으나, 전체적으로는 이론값과 유사한 경향을 나타내었다.

2. 출구에서의 체적분율은 벽면에서 낮게 나타나다가 가장 높은 값을 가지며, 내부로 갈수록 점차 감소한다. 이러한 경향은 기존의 수치해석 연구 결과와 일치한다.

3. 본 연구의 결과로 도출된 출구 증기 건도 0.25 이하에서는 기포류 및 슬러그류 영역에 해당되며, 이때는 질량유속보다는 열유속에 크게 의존한다. 이것은 핵비등이 지배적인 열전달 메커니즘으로 작용하는 것을 의미한다.

4. 추후의 연구에서는 본 연구에서 적용한 수치해석 기법을 토대로 다양한 조건에서의 시뮬레이션을 수행하고 신뢰성을 향상시켜, LNG 열교환기의 최적 설계 및 성능해석 등의 연구를 진행하고자 한다.

Acknowledgments

이 논문은 정부(과학기술정보통신부)의 재원으로 한구연구재단의 지원을 받아 수행된 연구임(No.2020R1G1A101062011).

References

• Antal SP, Lahey RT and Flaherty JE(1991). Analysis of phase distribution in fully developed laminar bubbly tow-phase flow, Int. J. Multiph. Flow, 17(5), 635~652. [https://doi.org/10.1016/0301-9322(91)90029-3]
• Burns AD, Frank T, Hamill I and Shi JM(2004). The favre averaged drag for turbulent dispersion in Eulerian multi-phase flow, ICMF. 4, 1~17.
• Chen D and Shi Y(2013). Experimental study on flow boiling heat transfer of LNG in a vertical smooth tube, Cryogenics, 57, 18~25. [https://doi.org/10.1016/j.cryogenics.2013.04.003]
• Choi YS and Lim TW(2015). Heat transfer characteristics and flow pattern investigation in micro-channels during two-phase flow boiling, JKOSME, 39(7), 696~701. [https://doi.org/10.5916/jkosme.2015.39.7.696]
• Choi YS and Lim TW(2016). A new correlation on heat transfer coefficient in horizontal multi channels, JFMSE, 28(5), 1388~1394. [https://doi.org/10.13000/JFMSE.2016.28.5.1388]
• Choi YS and Lim TW(2021). Design of multi-stage thermal cycle for recovery of LNG cold energy and waste heat, JFMSE, 33(2), 483~491. [https://doi.org/10.13000/JFMSE.2021.4.33.2.483]
• Choi YS, Kim YT, Lim TW(2016). CFD validation for subcooled boiling under low pressure, JKOSME, 40(4), 275~281. [https://doi.org/10.5916/jkosme.2016.40.4.275]
• Ishii M and Zuber N(1979). Drag coefficient and relative velocity in bubbly, droplet or particulate flows, AIChE. 25, 843~855. [https://doi.org/10.1002/aic.690250513]
• Koncar B, Kljenak L and Mavko B(2004). Modelling of local two-phase flow parameters in upward subcooled flow boiling at low pressure, Int. J. Heat Mass Transf. 47(6-7), 1499~1513. [https://doi.org/10.1016/j.ijheatmasstransfer.2003.09.021]
• Krepper E, Koncar B and Egorov Y(2007). CFD modelling of subcooled boiling-concept, validation and application to fuel assembly design, Nucl. Eng. Des. 237(7), 716~731. [https://doi.org/10.1016/j.nucengdes.2006.10.023]
• Mishima K(1984). Flow regime transition criteria for upward two-phase flow in vertical tubes, Int. J. Heat Mass Transf. 27(5), 723~737. [https://doi.org/10.1016/0017-9310(84)90142-X]
• Setoodeh H, Ding W, Lucas D and Hampel U(2021). Modelling and simulation of flow boiling with and Eulerian-Eulerian approach and integrated models for bubble dynamics and temperature-dependent heat partitioning, Int. J. Therm. Sci. 161, 106709. [https://doi.org/10.1016/j.ijthermalsci.2020.106709]
• Son H and Kim JK(2020). Energy-efficient process design and optimization of dual-expantion system for BOG (Boil-off gas) Re-liquefaction process in LNG fueled ship. Energy, 203(15), 117823. [https://doi.org/10.1016/j.energy.2020.117823]
• Tomiyama A, Tamai H, Zun I and Hosokawa S(2002). Transverse migration of single bubbles in simple shear flows, Chem. Eng. Sci. 57(11), 1849~1858. [https://doi.org/10.1016/S0009-2509(02)00085-4]
• Yoo, HS(2018). A study on the process analysis and the fuel gas heating temperature change for the improvement of dual fuel engine fuel gas methane number of LNGC, Korea Maritime and Ocean University, Thesis. http://kmou.dcollection.net/public_resource/pdf/200000105193_20211114033227.pdf