열화학적 CO₂ 메탄화 등온반응기의 수순환 냉각시스템 설계

열화학적 메탄화는 산업발전에서 발생되어 포집한 이산화탄소와 수전해로 생산된 수소를 반응시켜 메탄과 물을 생성하는 발열 반응이다. Sabatier에 의한 반응식은 다음과 같다⁶⁾.

촉매반응에 의해 메탄화가 이루어지고, Ni, Ru 등의 전이금속계열 촉매들이 사용된다. 본 연구에서는 발열에 의한 냉각시스템 설계가 목적이므로 촉매반응으로 인한 발열량을 사용하기 위해 다음과 같이 계산하였다.

촉매반응에 의한 기화량을 예측하기 위해 촉매가 충진된 튜브의 입출구의 실험값을 사용하였다. 식 (1)을 구성하는 CO₂, H₂, CH₄, H₂O의 입구와 출구의 엔탈피와 유량으로 냉각시스템으로의 방열량을 계산하였다. 입출구의 온도와 정압비열로 계산되는 엔탈피(h_s)는 조각다항식(piecewise-polynomial)으로 정의하여 입출구 온도(350℃)의 네제곱 형태로 구하고, 표준생성 엔탈피(standard enthalpy of formation, h_f)와 합하여 입출구 전엔탈피를 비교, 발열량을 계산하였다. 계산을 위한 CO₂, H₂, CH₄, H₂O의 입출구 조건을 Table 1에 나타내었다.

Table 1.

Enthalpy comparison between inlet & outlet

상용 CFD 코드인 ANSYS 사(Canonsburg, PA, USA)의 Fluent ver. 2020 R2를 통해 메탄화 반응기의 수순환 냉각시스템을 설계하였다. 압축성 점성 유동을 모사하기 위해 연속방정식, 모멘텀 방정식 및 에너지 방정식을 전산해석하고, 난류모델로 점성을 정의하였다. 상변화가 일어나지만 표면에서의 거동을 관찰할 목적이 아니며, 복잡한 형상과 소산에 대한 정확도를 고려하여 realizable k-ε 모델을 사용하였다. Realizable k-ε 모델은 standard k-ε 모델에서 난류 점성(μ_t)을 구성하는 점성계수(Cμ)를 변수로 정의하고 난류 소산(ε)의 수송방정식을 수정한 모델이다^7,8). 최적화 과정 종료 후 결정된 형상에 대해 uRANS 해석을 하여 시간에 따른 기화량과 온도를 관찰하였다.

다상유동의 모사를 위해 volume of fluid (VOF) 모델을 적용하였다. VOF 모델은 주로 자유표면이 있는 다상유동에서 사용하고, 초기조건에 독립적인 해석의 경우 정상해석이 가능해 해석 수가 많은 최적화 과정에 적합하다. 다상유동의 mass transfer는 기화-응축을 계산할 수 있는 Lee 모델을 사용하였다. Lee 모델은 증기 수송방정식(vapor transport equation)에 의해 지배되며, 각각의 격자에서 증기 체적분율에 해당하는 연속방정식으로 1차상(liquid)과 2차상(vapor)의 mass transfer rate를 계산한다. 1차상에서 2차상으로의 상변화를 $$ {\dot{m}}_{lv}$$, 반대의 상변화를 $$ {\dot{m}}_{vl}$$로 정의하면, 식(2)와 같이 증기 수송방정식이 성립한다.

v는 기체 상, α_v는 기체 체적분율, ρ_v는 기체 밀도, $$ {\overrightarrow{V}}_v$$는 기체의 속도이다. 이때, 1차상(liquid)의 온도가 포화온도보다 높으면 기화가 발생하고, 2차상(vapor)의 온도가 포화온도보다 낮으면 응축이 발생한다.

본 연구에서 사용한 등온반응기의 유동영역을 Fig. 1에 나타내었다. 전체 반응기의 높이는 1.8 m이고, 포화수는 y=0.25 m에서 직경 52.5 mm 배관으로 유입되고, 포화증기는 y=1.6 m에서 직경 62.7 mm 배관으로 배출된다. 메탄화 반응이 일어나는 튜브 내부는 고려하지 않고 수순환 냉각시스템과 맞닿은 튜브의 표면에 2.1.2절에서 계산된 발열량을 적용하였다. 열유속은 포화수의 증발잠열로 전달되고, 포화수의 기화가 발생한다. 수순환 냉각시스템은 기화된 기체가 상부로 배출되고, 기화량과 동일한 양의 포화수가 하부로 유입되는 구조이다. 총 74개의 튜브가 있고, 촉매는 Fig. 1에 표시된 만큼 충진되어 있다. 튜브 내 모든 촉매가 동시에 반응하지 않고, 상부에 위치한 촉매부터 반응하는 것을 실험을 통해 확인하여 발열이 발생하는 부분을 이에 맞게 고려하였다. 해석 시간과 정확성을 고려하여 대칭면의 한쪽 부분만 해석하였다.

Fig. 1.

CFD domain of methanation reactor and cooling system

격자는 fluent에서 제공하는 프로그램으로 형성하였다. 튜브의 간격이 좁기 때문에 proximity 조건을 주어 간격 간 최소 격자수를 정의하였다. 튜브의 곡률과 튜브 간 좁은 간격으로 발생할 수 있는 격자의 낮은 질을 polyhedral type 격자로 해소하였다. 기 구축된 CFD 해석 기법을 토대로 k-ε 난류모델을 사용하였기 때문에 경계층에 프리즘 격자를 형성하지 않았다. 격자의존성 해석의 결과로 nodes 총 21.7백만 개, cells 총 3.7백만 개를 해석에 사용하였다. 해석에 사용한 격자를 Fig. 2에 나타내었다.

Fig. 2.

Polyhedral mesh of CFD domain

CO2 메탄화 등온반응기의 냉각시스템은 입구로 포화수를 공급하고, 기화된 포화증기가 출구로 배출된다. 입구와 출구는 동일한 증기발생기(steam generator)에 연결되어 있고, 냉각 시스템에서 기화된 포화증기와 동일한 질량유량의 포화수가 입구로 공급된다. 따라서 입구와 출구의 유량은 다음 식과 같이 이전 iteration에서의 기화량을 적용하였다.

촉매반응에 의한 발열량은 튜브의 표면으로부터 포화수로 적용하였다. 위에서 언급한대로 튜브의 일부에서만 반응이 일어나기 때문에 발열을 적용하는 면적을 좁게 정의하였다. 해석의 효율성을 고려하여 대칭조건을 사용하였고, 발열이 발생하는 벽면을 제외하고 모두 단열조건을 적용하였다.

해석 영역은 중력가속도를 고려하였고, 출구면 기준 40 bar(g)의 고압 상태로 해석하였다. 본 연구에서 사용한 경계조건을 Table 2에 나타내었다.

Table 2.

Boundary conditions

본 연구에서는 기 설계된 CO₂ 메탄화 반응기의 수순환 냉각시스템 성능 향상을 위해 최적화 과정을 적용하였다. 튜브의 간격과 쉘 직경을 결정하는 변수 k와 쉘 직경의 크기를 조정하는 변수 Δd를 형상변수(factors)로 정의하였다. 기 설계된 반응기 튜브간 거리에 변수 k를 곱하였고, 쉘의 직경도 동일하게 k의 변수로 설정하였다. k로 결정된 쉘의 직경을 형상변수 Δd로 조정하여 형상을 구체화 하였다. 최적화 과정에서 생성되는 실험점의 형상에서 튜브와 쉘의 간섭으로 인한 오류를 없애기 위해 쉘 직경에 대해 두 개의 형상변수를 적용하였다. 형상변수 k와 Δd를 Fig. 3에 나타내었다.

Fig. 3.

Design variables; k & Δd

실험계획법(design of experiments) 및 최적화 과정을 Fig. 4에 나타내었다. 최적화를 진행할 문제를 정의한 후 최적화를 위한 변수를 정의한다. 본 연구에서는 동일한 경계조건에서 형상변수를 통한 영향을 연구하였다. 형상을 정의하는 모든 변수를 형상변수로 사용할 수 있지만, 각 변수의 민감도를 파악하기 어렵고, 해를 찾는 데 시간이 오래 걸리며 초기 실험점 배치(sampling method)에 따라 국소 최적해(local optimum)로 수렴할 가능성이 높다. 반응기 냉각시스템은 기존 설계의 변경 가능성, 형상변수의 민감도 파악을 고려하여 2개의 형상변수로 반응기의 물리적 규모를 결정하고자 하였다. 형상변수에 따른 냉각시스템의 결과는 튜브의 표면 온도로 관찰하였다. 튜브의 냉각이 적절히 이루어져야 하고, 위치별 온도 편차가 적어야 하며, 국부적으로 높아지는 온도를 해소하기 위해 목적함수와 구속조건을 다음과 같이 정의하였다.

Fig. 4.

CFD response surface optimization process

문제 정의와 변수 선정 이후 실험계획법을 통해 실험점(sample points)을 생성한다. 형상변수의 개수가 적고 초기 설계공간(design space)에서 변수의 영향을 파악하기 위해 중심합성계획(centered composite design)을 사용하였다. 변수 개수에 따른 중심합성계획법의 초기 실험점 개수는 다음 식 (5)를 따른다⁹⁾.

형상변수 개수가 2개이기 때문에 초기 실험점은 9개이다. 형상변수의 개수가 증가함에 따라 실험점의 개수가 가파르게 증가하기 때문에 중심합성계획법은 주로 형상변수의 개수가 적을 때 사용한다. Fig. 5에 설계공간을 나타내었다.

Fig. 5.

Design space of CCD model (DVs=2)

중심합성계획으로 생성한 실험점을 모두 해석한 후 반응면(response surface method)을 생성하였다. 설계공간에서 형상변수의 좌표와 목적함수의 값을 토대로 스플라인 곡면으로 반응면을 형성한다. 반응면을 생성하는 메타 모델은 크리깅(Kriging)을 사용하였다. 크리깅은 정규분포를 가정하지 않는 비모수 보간법 모델이고, 전산실험계획(design and analysis of computer experiments)으로도 알려져 있듯 전산유체역학에서 높은 근사 예측도를 갖는다^10,11).

크리깅 모델로 생성된 반응면을 토대로 다목적 유전 알고리즘(multi-objective genetic algorithm)을 수행하여 후보점(candidate points)을 도출하였다. 후보점은 초기 100개의 샘플과 1 iteration 당 100개의 샘플을 최대 20 iteration 수행하여 3개의 후보점을 선정하였다. 유전 알고리즘은 다음과 같이 식 (6) 교차(cross-over) 연산과 식 (7) 돌연변이(mutation) 연산을 통해 다음 세대를 도출한다^12-14).

C는 child, P는 parent로 C는 P의 다음 세대이고, δ는 다항분포로부터 계산된 분산이다.

생성된 후보점의 예측 결과와 CFD 해석 후 실제 결과의 차이가 0.1% 이하면 최적화 과정을 종료하였다. 만약 차이가 0.1% 이상이면 새로운 후보점을 추가한 뒤 반응면을 생성, 최적해 도출의 위 과정을 반복하였다.

9개의 초기 실험점을 해석한 후 22번의 후보점 생성 및 해석을 반복하고 최적화 과정을 종료하였다. 최적화 과정이 종료된 반응면을 Fig. 6에 나타내었다. 반응면은 형상변수 k와 Δd에 따른 튜브의 표면 평균 온도와 온도 편차이며, 31개의 실험점으로 형성한 반응면이다. Fig. 1의 tube (heat flux)면의 평균 온도인 t_ave (Fig. 6[a])는 영역 1에서 가장 낮게 나타났고, 영역 2-4에서 목표 평균 온도(330℃)가 분포하였다. Tube (heat flux)면의 온도 편차인 t_dev (Fig. 6[b])는 영역 1과 3에서 낮은 값을 보였다. 목적함수를 모두 만족하는 영역 3에서 파레토 최적해 중 하나를 선정하였고, 도출된 최적해를 반응면 위에 표시하였다.

Fig. 6.

Response surfaces; t_ave (top) & t_dev (bottom)

최적화 과정을 수행하기 전 모델을 Base로, 최적화 과정 이후의 모델을 Opt.로 명명하고 각각의 형상을 Fig. 7에 나타내었다. Opt. 형상은 Base 형상에 비해 튜브의 간격이 약 30.3% 감소하였고, 쉘 직경이 약 3.6% 감소하였다. 튜브 발열 부분의 목표 평균 온도를 만족하고 최고·저 온도차의 최소화를 위해 도출된 최적해의 물리적 분석을 수행하였다.

Fig. 7.

Base model (left) & Opt. model (right)

튜브의 발열에 의해 기화되는 포화수의 양을 토대로 본 연구에서 사용한 CFD 기법을 구축하였다. 고압의 포화수 증발잠열과 튜브의 발열량으로 기화량을 계산하여 CFD 결과와 비교하였고, 형상 변경 외에 경계조건은 동일하게 해석하였기 때문에 모든 해석에서 기화량은 동일하였다. 비정상 해석 시 해석 초기에는 포화수의 증발잠열로 적은 기화량을 보였고, 시간이 지남에 따라 기화량은 이론값과 일치하였다. 시간에 다른 기화량을 Fig. 8에 나타내었다. 그래프에 첨부된 그림은 대칭면에서의 증기 체적분율을 의미한다.

Fig. 8.

Base model (left) & Opt. model (right)

Base, Opt. 형상의 튜브 표면 온도분포를 Fig. 9와 같이 비교하였다. 두 형상 모두 반응기 중심에 위치한 튜브에서 높은 온도분포를, 바깥쪽에 위치한 튜브에서 낮은 온도분포를 보였다. 각각의 튜브 온도도 마찬가지로 반응기 중심쪽 면에서 높은 온도를 보였고, Base 형상의 튜브 표면 온도는 입출구가 위치한 +z축 방향에서 비교적 더 높은 분포를 보였다. Opt. 형상은 Base 형상에 비해 튜브의 간격이 감소하였음에도 불구하고 반응기 중심부의 국부적으로 높은 온도가 감소하였다. 반응기 중심에서 두 형상의 온도 분포가 큰 차이를 보인 것에 반해, 튜브와 쉘 사이의 온도 분포는 두 형상에서 큰 차이가 없었다.

Fig. 9.

Tube surface temperature (Base & Opt. model)

반응기 내부의 유동을 비교하여 형상 간 튜브의 온도차 이유를 확인하고자 하였다. 가장 높은 온도를 보이는 반응기 중앙에 위치한 튜브와 그 사이의 유속을 Fig. 10에 나타내었다. 튜브 표면에서 기화가 일어나면서 두 형상 모두 유속이 증가하였다. 두 형상의 기화량이 같기 때문에 입구로 유입되는 포화수의 양도 동일하지만 Opt. 형상의 튜브 간격이 좁아 유속이 더 높게 나타난 것으로 사료된다. 튜브의 발열부분 아래(y<1.15 m) 부분에서는 두 형상의 속도가 0.05-0.10 m/s로 비슷하였다. y>1.15 m부터 속도 차이가 점점 커지고, y=1.25 m 부분에서 가장 차이가 컸다. 튜브의 표면 온도는 Fig. 9와 마찬가지로 Base 형상에서 국부적으로 높은 온도를 보였으며, 온도의 편차가 컸다. 반면 Opt. 형상은 표면의 온도가 비교적 고르게 분포하였다. 포화증기로의 상변화가 일어나면 튜브 표면에서 열전달률이 감소하게 되는데, opt. 형상은 튜브 사이에서 비교적 빠른 유속으로 열전달이 활발히 발생하였다.

Fig. 10.

Tube surface temperature & velocity profiles between tubes

y=1.3 m 위치에서 y 방향 속도성분과 증기 체적분율을 비교하였다(Fig. 11). 속도성분의 방향을 비교하기 위해 양의 속도성분만 표시하였다. 기화가 일어나는 튜브의 사이에서 양의 y 방향 속도성분을 확인할 수 있었고, 튜브와 반응기 쉘 사이에서는 -y 방향의 속도성분을 보였다. No-slip condition이기 때문에 튜브의 표면에서 증기와 포화수의 속도는 모두 0이지만, 기화가 일어나는 튜브의 사이에서 높은 속도를 보였다. 증기 체적분율은 대부분 튜브 표면에서 보였으며, 튜브와 쉘 사이에서는 주로 포화수가 분포하였다. 튜브 간격이 줄어든 Opt. 형상에서도 튜브 간 증기 체적분율이 간섭하지 않아 튜브의 표면 온도가 유지된 것으로 사료된다. y 방향 속도성분과 체적분율 컨투어를 통해 반응기의 중앙부에서는 증기와 포화수가 y 방향으로 이동하고, 반응기 바깥 부분에서는 포화수가 -y 방향으로 이동하는 것을 확인하였다.

Fig. 11.

Velocity v & vapor volume fraction contour; Base (top) and Opt. (bottom)

반응기 대칭면에서의 증기 체적분율과 벡터장을 Fig. 12에 나타내었다. 반응기 상단 부분에서 두 형상 모두 포화증기 상태였고, 수면이 동일하였다. 발열이 발생하는 부분부터 상변화가 일어나기 시작하여 유체의 속도가 점차 증가하였다. 앞의 Fig. 11과 동일하게 튜브와 튜브 사이에서는 상변화가 활발히 일어나며 포화수와 증기가 +y 방향으로 이동하였다. 반면 쉘 부분에서는 주로 포화수가 분포하였고, 정체되거나 -y 방향으로의 낮은 속도분포를 보였다. 반응기 상단 부분에서는 포화증기가 낮은 속도로 와류를 형성한 후 출구 부분에서 속도가 상승하며 배출되었다. 출구 부분에서의 포화수는 출구 방향(+z)으로 이동하는 경향이었으나 좁아진 단면적을 고려하였을 때 속도가 매우 작았다. 전체적인 거동은 두 형상에서 비슷하였지만 튜브 사이에서의 속도가 Opt. 형상이 Base 형상에 비해 비교적 높았고, 증기 체적분율도 높았다.

Fig. 12.

Vapor volume fraction contour & velocity vector plot; Base (top) and Opt. (bottom)

Table 3에 형상변수와 목적함수, 구속조건의 결과를 기존 형상과 비교하였다. 형상변수 k는 Base 형상에 비해 약 30.3% 감소하였다. 튜브 간격은 이와 동일하게 감소하였고, 쉘 직경은 Δd (=-16.12 mm)에 의해 약 3.6% 감소하였다. 동일한 경계조건에서 변경된 형상에 의해 표면의 튜브 온도분포가 목적함수와 구속조건에 따라 변경되었다. 튜브 표면의 평균 온도는 330℃에 가까운 값을 목적함수로 정의하여 329.92℃의 결과가 나왔는데, 기존 형상 역시 평균 온도를 만족하도록 설계되었기 때문에 큰 차이는 없었다(약 0.3%). 촉매 반응 속도를 위하여 목적함수로 튜브 표면의 온도편차를 최소화하였고, 최고 온도의 상한을 구속조건으로 적용하였다. 그 결과, 튜브 표면에서의 최고·최저 온도차는 약 19.48% 감소하였으며, 최고 온도는 약 8.31% 감소하였다. 본 연구를 통해 기존 형상의 튜브 간격과 쉘 직경을 최적설계하여 튜브 간 온도차를 줄여 촉매반응이 고르게 일어날 수 있도록 하였다.

Table 3.

Optimization results