1. 서 론
2. 수치 해석 방법
2.1 해석 모델
2.2 CO2 흡착 등온선 모델 및 속도론 모델
2.3 등량흡착열 모델
2.4 질량 보존 방정식
2.5 운동량 보존 방정식
2.6 에너지 보존 방정식
2.7 벽면 경계조건
2.8 재료 물성치 모델
2.9 모델 파라미터 및 해석 조건
3. 모델 검증 데이터
4. 결과 및 분석
4.1 격자 민감도 분석
4.2 실험 파과곡선 및 온도 응답 검증
4.3 3차원 CO2 농도장 및 온도장 분석
4.4 흡착률 분포를 통한 추가 흡착 영역 분석
5. 결 론
1. 서 론
최근 산업 활동의 고도화와 에너지 소비 증가로 인해 대기 중 이산화탄소 농도가 지속적으로 상승하면서, CO2 배출 저감은 탄소중립 실현을 위한 핵심 과제로 부상하고 있다[1]. 특히 화석연료 기반 발전, 제철, 시멘트 및 석유화학 공정과 같은 대규모 고정 배출원은 고농도 CO2 배출의 주요 원인으로 지적되며, 이들 배출원으로부터 CO2를 효과적으로 분리·회수하기 위한 탄소 포집 기술의 중요성이 지속적으로 확대되고 있다[2]. 현재 아민 흡수 기반 공정은 높은 기술 성숙도와 상용화 가능성을 확보한 대표적인 CO2 포집 기술이지만, 흡수제 재생에 요구되는 높은 에너지 소비, 용매 열화, 장치 부식, 휘발성 손실 및 운전 비용 증가와 같은 한계를 수반한다[3,4]. 이러한 한계를 보완하기 위해 고체 흡착제를 이용한 흡착 기반 CO2 포집 공정이 유망한 대안으로 주목받고 있으며, 이는 낮은 재생 에너지 요구량, 단순한 공정 구성, 우수한 흡착제 재사용성 및 다양한 운전 조건에 대한 적용 가능성 측면에서 기존 액상 흡수 공정을 보완할 수 있는 잠재력을 가진다[5].
고체 흡착제 중 금속-유기 골격체(Metal-Organic Framework, MOF)는 높은 비표면적, 조절 가능한 기공 구조, 다양한 금속 노드와 유기 리간드 조합을 바탕으로 CO2 선택 흡착에 유리한 소재로 평가된다[6,7,8]. 특히 유연 금속-유기 골격체(flexible MOF)는 일반적인 강체 다공성 흡착제와 달리 압력, 온도 또는 흡착질 분압 변화에 따라 골격 구조가 가역적으로 변화하며, 특정 압력 구간에서 흡착량이 급격히 증가하는 gate-opening 또는 sigmoidal adsorption 거동을 나타낸다. 이러한 구조 응답성은 좁은 압력 범위에서도 큰 흡착량 차이를 유도할 수 있어, 압력변동흡착(PSA) 및 진공압력변동흡착(VPSA) 기반 CO2 분리 공정에서 높은 working capacity와 향상된 재생 효율을 기대하게 한다[9].
본 연구에서 대상으로 하는 Elastic Layer-structured Metal-organic framework-11(ELM-11)은 CO2 분압 변화에 따라 gate-opening 거동을 보이는 대표적인 flexible MOF이다. ELM-11은 낮은 압력 영역에서는 제한적인 흡착 거동을 나타내지만, 특정 분압 이상에서 구조적 전이가 발생하면서 흡착량이 급격히 증가하는 비선형 sigmoidal 등온선을 보인다. 이러한 특성은 CO2 포집 공정에서 선택적 흡착과 재생 효율 향상에 유리하게 작용할 수 있으나, 동시에 hysteresis와 복잡한 gate-adsorption 거동을 동반하기 때문에 흡착 평형, 흡착속도 및 고정층 파과 거동의 예측을 어렵게 한다[10]. 특히 gate-opening 압력 이하에서는 흡착이 충분히 진행되지 않은 상태로 기체가 흡착층을 통과하는 slipping-off 현상이 발생할 수 있으며, 이는 일반적인 S자형 파과곡선이 아닌 계단형 파과곡선과 복잡한 흡착 전선 이동을 유발한다. 따라서 ELM-11 기반 CO2 포집 공정을 정확히 해석하기 위해서는 gate-opening에 따른 흡착 평형과 속도론적 거동을 동시에 고려할 수 있는 모델링 접근이 필요하다.
기존 고정층 흡착 해석은 주로 1차원 축방향 모델 또는 lumped parameter 기반 공정 모델을 활용하여 출구 농도 변화와 전체 흡착량을 예측하는 데 초점을 두었다. 이러한 접근은 계산 효율성과 공정 최적화 측면에서 장점을 가지지만, 충전층 내부에서 발생하는 국부적인 CO2 농도 구배, 분압 변화, 흡착량 분포, 흡착열에 의한 온도 상승 및 반경방향 비균일성을 직접적으로 해석하는 데에는 한계가 있다[11]. 실제 흡착 거동은 입구 조건만으로 결정되지 않으며, 충전층 내부 각 위치에서의 CO2 분압, 온도, 흡착 포화도 및 물질전달 저항에 의해 지배된다. 특히 ELM-11과 같이 gate-opening 특성을 갖는 흡착제의 경우, 국부 CO2 분압과 온도 변화가 구조 전이와 흡착속도를 좌우하므로, 단순한 출구 파과곡선만으로는 내부 흡착 전선의 형성 및 이동, 국부 포화 영역, 열적 비균일성을 충분히 설명하기 어렵다. 따라서 ELM-11 기반 흡착 공정의 정밀한 해석과 설계를 위해서는 농도장, 분압장, 온도장 및 흡착량장을 공간적으로 동시에 고려할 수 있는 3차원 비정상 수치해석이 요구된다.
이에 본 연구에서는 pelletized ELM-11이 충전된 고정층 흡착탑을 대상으로 3차원 비정상 흡착 수치모델을 구축하고, CO2 흡착 및 파과 거동을 해석하였다. ELM-11의 gate-opening에 따른 급격한 흡착량 변화와 비선형 흡착 특성을 정량적으로 반영하기 위해 계단형 흡착 전이를 고려한 흡착평형식과 속도론적 흡착 모델을 적용하였으며, 비정상 해석을 통해 출구 CO2 파과곡선뿐만 아니라 충전층 내부의 CO2 농도장, 온도장 및 흡착량 분포를 함께 분석하였다. 이를 통해 기존 1차원 파과곡선 중심 해석에서 확인하기 어려운 흡착 전선의 3차원 이동, 국부 포화 영역의 형성, slipping-off에 따른 계단형 파과 특성 및 흡착열에 의한 온도 변화 메커니즘을 규명하고자 하였다.
2. 수치 해석 방법
2.1 해석 모델
본 연구의 해석 대상은 참조 논문[12]에서 실험적으로 사용된 원통형 고정층 흡착탑이다. 해당 흡착탑은 Elastic Layer-structured Metal-organic framework-11(ELM-11) 흡착제가 충전된 구조로, CO2/N2 혼합기체가 흡착제 베드를 통과하는 동안 발생하는 CO2 흡착 및 파과 거동을 모사하기 위해 사용되었다.
해석 형상은 Fig. 1과 같이 ELM-11이 충전된 원통형 흡착제 베드로 구성하였다. 참조 논문[12]의 실험 조건과 직접 비교하기 위해, 실험에서 사용된 흡착제 베드의 내경 7.75 mm, 길이 170 mm, 벽 두께 0.889 mm를 동일하게 반영하여 3차원 형상을 구축하였다. 실제 펠릿 입자의 개별 형상과 배열은 직접 구현하지 않고, 흡착제 충전층을 평균화된 다공성 매질(porous media)로 가정하였다. 이를 통해 복잡한 입자 스케일 구조를 단순화하면서도, 충전층 내부의 CO2 농도장, 온도장 및 흡착량 분포를 공간적으로 해석할 수 있도록 하였다.
수치 모델은 질량, 운동량 및 에너지 보존 방정식을 기반으로 구성하였다. CO2 흡착 거동을 반영하기 위해 Ansys Fluent[13]와 사용자 정의 함수(User Defined Function, UDF)를 연계하였으며, 흡착평형식으로부터 계산된 평형흡착량과 선형구동력(Linear Driving Force, LDF) 모델을 이용해 시간에 따른 흡착속도를 산정하였다. 계산된 흡착속도는 각 시간 단계에서 질량 및 에너지 보존 방정식의 소스항으로 적용되며, 이를 통해 유동장, 농도장, 온도장 및 흡착량장의 상호작용을 비정상적으로 해석하였다. 본 모델은 원통형 고정층뿐만 아니라 다양한 형상과 다차원 흡착 해석에도 확장 적용할 수 있다.
본 연구의 수치해석에는 다음과 같은 가정을 적용하였다. 기체상은 이상기체 거동을 따르며, 흡착탑 내부 유동은 낮은 Reynolds 수 조건(Re < 10)에 따라 층류로 가정하였다. 흡착제 베드는 균질한 다공성 매질로 취급하였으며, 공극률, 입자 특성 및 열물성은 해석 영역 내에서 일정하다고 가정하였다. 또한 흡착 과정에서의 물질전달 저항은 LDF 모델로 표현하였고, 흡착에 따른 CO2 제거량과 흡착열 효과는 각각 질량 및 에너지 방정식의 소스항으로 반영하였다.
2.2 CO2 흡착 등온선 모델 및 속도론 모델
ELM-11의 CO2 흡착 평형 특성은 참조 논문[12]에서 제시된 실험 흡착 등온선 데이터를 기반으로 정의하였다. ELM-11은 일반적인 강체 흡착제와 달리 CO2 분압 변화에 따라 gate-opening 거동을 나타내며, 이로 인해 특정 압력 구간에서 흡착량이 급격히 증가하는 비선형 흡착 특성을 가진다. 따라서 본 연구에서는 이러한 흡착 거동을 반영하기 위해 일부 수정된 weighted dual-site Langmuir(wDSL) 모델을 적용하였다[14].
wDSL 모델은 ELM-11의 CO2 흡착 거동을 gate-closing 상태와 gate-opening 상태의 조합으로 표현한다. 평형 흡착량 는 gate-closing 상태에 해당하는 흡착량 과 gate-opening 상태에 해당하는 흡착량 를 가중함수 로 조합하여 다음과 같이 나타낼 수 있다.
여기서 , , , 는 등온식 매개변수이다. 는 CO2의 분압을 의미한다. 과 는 각각 이상기체상수와 온도이다. 가중함수 는 step pressure, 즉 gate-opening 또는 gate-closing 압력을 나타내는 와 step pressure 부근에서 등온선의 기울기를 나타내는 γ를 사용하여 다음 식으로 표현된다.
여기서 역시 등온선의 형상, 즉 등온선이 시작되는 압력과 기울기 모두에 영향을 주는 매개변수이다. Step pressure 는 다음 식으로 표현된다.
여기서 는 기준 온도 = 273.15 K에서의 등온선 step pressure이고, 는 상전이 엔탈피이다. 흡착 및 탈착 가지 각각에 대한 , , , , , γ, , 매개변수는 참조 논문[12]의 실험 등온선을 바탕으로 최소자승법을 사용하여 추정하였다.
흡착 속도론은 선형 구동력(Linear Driving Force, LDF) 모델을 사용하여 표현하였다. LDF 모델은 실제 흡착량 가 평형 흡착량 에 접근하는 속도를 1차 지연 형태로 나타내며, 다음과 같이 정의된다.
여기서 는 시간에 따라 변화하는 실제 흡착량, 는 해당 위치의 온도와 CO2 분압에서 계산되는 평형 흡착량을 의미한다. 는 흡착 속도 상수로, 기체가 흡착제 입자 외부 경막, 거대기공 및 미세기공을 통과하면서 발생하는 전체 물질전달 저항을 포함하는 유효 계수이다.
본 연구에서는 를 wDSL 모델로부터 계산하고, LDF 모델을 이용하여 각 셀에서의 시간에 따른 흡착량 변화를 계산하였다. 또한 는 참조 실험의 CO2 파과곡선 형태를 재현할 수 있도록 curve fitting을 통해 역으로 추정하여 적용하였다. 이를 통해 ELM-11의 gate-opening에 따른 비선형 평형흡착 특성과 고정층 내 동적 흡착속도 효과를 동시에 반영하였다.
추정된 wDSL 모델의 적합성을 확인하기 위해, 각 온도 조건에서의 실험 흡착 등온선과 wDSL 모델 계산 결과를 비교하였다. Fig. 2는 273.15–323.15 K 범위에서 측정된 ELM-11의 CO2 흡착 등온선과 wDSL 모델의 예측 결과를 나타낸다. 그림에서 파란색 사각형은 실험값을, 빨간색 실선은 wDSL 모델 계산값을 의미한다. 전체 압력 범위에서 wDSL 모델은 실험 흡착량 변화를 효과적으로 재현하였으며, 특히 ELM-11의 특징적인 gate-opening 구간에서 나타나는 급격한 흡착량 증가를 적절히 표현하였다. 이러한 적합성은 정량적 오차 평가에서도 확인되었으며, 전체 흡착 등온선에 대한 RMSE는 0.134 mol kg-1, R2는 0.991로 나타났다. 이는 추정된 wDSL 모델이 실험 흡착 등온선과 높은 일치도를 가지며, ELM-11의 비선형 흡착 거동을 적절히 설명할 수 있음을 의미한다.
또한 온도가 증가함에 따라 gate-opening이 발생하는 CO2 분압이 높은 압력 방향으로 이동하는 경향이 확인되었다. 이는 ELM-11의 구조전이가 온도와 CO2 분압에 강하게 의존함을 의미하며, 식 (6)의 항을 통해 step pressure의 온도 의존성을 반영한 wDSL 모델이 ELM-11의 비선형 흡착 특성을 설명하는 데 적합함을 보여준다. 따라서 본 연구에서는 Fig. 2에서 검증된 wDSL 모델을 이용하여 각 계산 셀의 국부 온도와 CO2 분압에 따른 평형흡착량 를 계산하였다.
2.3 등량흡착열 모델
ELM-11의 CO2 흡착은 gate-opening에 의해 특정 흡착량 구간에서 급격한 흡착 거동 변화를 보인다. 따라서 본 연구에서는 흡착열을 단일 상수로 가정하지 않고, 흡착량에 따라 변화하는 등량흡착열(isosteric heat of adsorption) 로 표현하였다. 이를 위해 먼저 wDSL 흡착 등온선 모델로부터 각 온도에서의 평형흡착량을 계산하고, 동일 흡착량 조건에서의 CO2 평형분압을 역산하였다. 이후 Clausius–Clapeyron 관계식을 이용하여 흡착량에 따른 를 계산하였다.
먼저, wDSL 모델로부터 주어진 온도 T와 CO2 분압 에서의 평형흡착량은 다음과 같이 표현된다.
여기서 는 wDSL 모델로 계산되는 CO2 평형흡착량이다. 등량흡착열을 계산하기 위해서는 동일한 흡착량 에서 온도 변화에 따른 평형분압 의 변화를 알아야 하므로, 각 온도 에 대해 다음의 비선형 방정식을 만족하는 압력을 수치적으로 계산하였다.
이를 통해 특정 흡착량 q에서의 평형분압 를 얻을 수 있다. 이후 Clausius–Clapeyron 식을 적용하면, 등량흡착열은 다음과 같이 계산된다.
여기서 R은 이상기체상수이며, 는 동일 흡착량 에서의 CO2 평형분압이다. 실제 계산에서는 여러 온도 조건에서 얻은 를 이용하여 와 사이의 선형 회귀를 수행하였다.
계산된 는 Fig. 3과 같이 흡착량에 따라 비선형적으로 변화하였다. 낮은 흡착량 영역에서는 가 급격히 증가한 뒤 약 21 kJ mol−1 수준의 완만한 plateau를 형성하였다. 이후 높은 흡착량 영역에서는 gate-opening에 따른 흡착 상태 변화의 영향으로 가 다시 증가하여 약 33 kJ mol−1 수준의 피크를 나타내었다. 이러한 결과는 ELM-11의 CO2 흡착이 단순한 Langmuir형 흡착이 아니라, 낮은 흡착량 영역, gate-opening 전후 영역, 고흡착량 포화 영역에서 서로 다른 열적 특성을 가짐을 의미한다.
수치해석에 적용하기 위해서는 계산된 데이터를 연속적인 함수 형태로 표현할 필요가 있다. 본 연구에서는 Clausius–Clapeyron 식으로 계산된 등량흡착열을 세 개의 sigmoid 함수로 fitting하여 흡착량 의존 흡착엔탈피 모델을 구성하였다.
여기서 는 fitting된 등량흡착열, 는 기준 흡착열, , , 는 각 sigmoid 항의 크기를 의미한다. , , 는 각 전이 구간의 기울기를 조절하는 매개변수이며, , , 는 각 변화가 발생하는 흡착량 위치를 나타낸다. 첫 번째 sigmoid 항은 저흡착량 영역에서의 급격한 흡착열 증가를, 두 번째 sigmoid 항은 gate-opening 이후 고흡착량 영역에서의 흡착열 증가를, 세 번째 항은 포화 영역에서 나타나는 완만한 흡착열 감소를 표현하기 위해 도입하였다.
Fig. 3은 Clausius–Clapeyron 관계식으로 계산한 ELM-11의 CO2 등량흡착열과 three-sigmoid 함수로 fitting한 결과를 나타낸다. 파란색 기호는 각 흡착량 조건에서 계산된 Clausius–Clapeyron 기반 값을 의미하며, 빨간색 실선은 이를 식 (12)의 three-sigmoid 함수로 근사한 결과이다. 그림에서 확인할 수 있듯이, fitting 함수는 저흡착량 영역에서의 급격한 증가, 중간 흡착량 영역의 plateau, 그리고 고흡착량 영역에서의 추가적인 증가 및 포화 부근의 완만한 감소 경향을 전반적으로 잘 재현하였다. 또한 fitting 결과의 정량적 오차를 평가한 결과, RMSE는 0.106 kJ mol-1, R2는 0.999로 나타나 Clausius–Clapeyron 기반 계산값과 three-sigmoid fitting 결과가 높은 일치도를 보임을 확인하였다.
2.4 질량 보존 방정식
질량 보존 방정식은 다음과 같이 모든 화학종의 대류-확산 방정식을 풀어 각 화학종의 국소 질량 분율을 계산한다.
여기서 𝜌는 가스 밀도, 는 화학종 (CO2)의 질량 분율, 는 속도 벡터, 𝜖은 베드 공극률(bed porosity), 는 화학종 의 질량 분산 계수(mass dispersion coefficient), 는 화학종 의 분자량, 는 성분 의 흡착량, 는 시간이다. 우변(식 (13))의 두 번째 항은 흡착 소스(source)를 의미한다. 전체 질량 수지(overall mass balance)는 다음과 같다.
2.5 운동량 보존 방정식
가스 유동의 운동량 방정식은 다음과 같이 표현될 수 있다.
여기서 는 사용 압력(operating pressure), 는 응력 텐서(stress tensor), 는 중력 가속도 벡터, 는 다공성 매질 내의 운동량 소스 항(momentum source term)이다. 이는 관성 및 점성 저항을 포함하는 Ergun 방정식으로부터 다음과 같이 계산한다.
여기서 𝜇는 가스의 동점성계수(dynamic viscosity), 는 다공성 매질의 점성 저항(viscous resistance), 는 관성 저항(inertial resistance), 는 속도 벡터의 크기, 는 방향의 속도 성분이다.
2.6 에너지 보존 방정식
CO2 분리를 위한 에너지 방정식은 흡착제 베드 내부에 저장된 에너지와 대류 유동, 압력 일(pressure work), 열 확산 및 이류(advection)로 인한 에너지 변화, 그리고 흡착 과정에 의해 방출되거나 소비되는 에너지 간의 평형을 나타낸다. 이는 다음과 같이 작성할 수 있다.
여기서 는 총 가스 에너지, 는 총 흡착제 에너지이다. 는 흡착열로써 물질에 따라 다르며 흡착된 가스량의 변화에 따라 달라지기도 한다. 는 평형 온도, 는 현열 엔탈피(sensible enthalpy), 그리고 는 가스 성분 의 확산 플럭스(diffusion flux)이다.
매개변수 는 흡착제 베드의 유효 열전도도(effective conductivity)이며 다음과 같이 표현될 수 있다.
여기서 와 는 각각 가스 혼합물과 흡착제의 열전도도이다.
금속 벽면의 경우, 벽을 통과하는 에너지 방정식은 다음과 같이 벽에 저장된 열과 벽을 통해 확산되는 열의 평형에 따라 결정된다.
여기서 는 벽면 밀도, 는 벽면 열용량, 는 벽면의 국소 온도, 그리고 는 벽면 재질의 열전도도이다.
2.7 벽면 경계조건
흡착 과정에서 발생하는 열은 흡착제 베드 내부의 온도 분포와 평형흡착량에 영향을 미치므로, 본 연구에서는 에너지 방정식을 이용하여 비등온 흡착 해석을 수행하였다. 흡착탑 벽면과 외부 대기 사이의 열전달은 대류 및 복사 열전달을 함께 고려하여 다음과 같이 표현할 수 있다.
여기서 는 벽면의 열전도도, 는 벽면 온도, 는 외기 온도, 은 벽면에 수직한 방향을 의미한다. 는 외부 대류 열전달계수, σ는 Stefan-Boltzmann 상수이다. 식 (20)의 우변 첫 번째 항은 벽면과 외부 공기 사이의 자연대류 열전달을 나타내며, 두 번째 항은 벽면과 주변 환경 사이의 복사 열전달을 의미한다.
외부 대류 열전달계수 는 Nusselt 수 상관식을 이용하여 계산할 수 있다. 원통형 벽면 주위의 자연대류 열전달에 대해 다음과 같은 상관식이 사용된다.
여기서 는 흡착탑 외부 직경, 는 공기의 열전도도, 는 공기의 Prandtl 수, 는 Rayleigh 수이다. 공기의 열물성은 벽면 온도와 외기 온도의 평균온도에서 평가된다. 식 (21)를 통해 계산된 Nusselt 수로부터 외부 대류 열전달계수 를 산출할 수 있다.
흡착제 베드와 벽면 사이의 계면에서는 고체 벽면과 다공성 흡착층 사이의 열전달이 발생한다. 이때 계면에서의 열전달은 벽면 측 열전도와 흡착제 베드 측 유효 열전도가 서로 연속된다는 조건을 이용하여 다음과 같이 나타낼 수 있다.
여기서 는 벽면과 흡착제 베드 사이의 내부 대류 열전달계수를 의미한다. 식 (22)은 벽면과 흡착제 베드 계면에서 열유속이 연속적으로 보존됨을 나타내며, 본 연구에서는 수치해석 프로그램에서 계산되는 wall-interface 열전달 조건을 이용하여 이를 반영하였다.
2.8 재료 물성치 모델
본 연구에서 사용된 물질은 흡착제인 ELM-11, 기체 성분인 CO2와 N2, 그리고 흡착탑 벽면 재료로 구성된다. 수치해석에서는 CO2/N2 혼합기체가 ELM-11이 충전된 다공성 베드를 통과하면서 발생하는 유동, 종수송, 흡착 및 열전달 현상을 해석하였다. 이를 위해 기체 혼합물의 밀도, 비열, 점도 및 열전도도와 같은 물성은 국부 온도와 조성 변화에 따라 계산되도록 설정하였다.
혼합기체의 밀도는 이상기체 거동을 따른다고 가정하였으며, CO2와 N2의 조성 변화에 따른 평균 분자량을 반영하였다. 혼합기체의 비열은 각 성분의 몰분율을 고려한 혼합법칙을 적용하였고, 점도와 열전도도는 CO2와 N2 각각의 물성을 바탕으로 혼합기체의 유효 물성으로 계산하였다. 흡착제 베드의 열적 특성은 기상, 고체상 ELM-11, 그리고 흡착된 CO2가 가지는 열용량 효과를 함께 고려하여 정의하였다. 다만 N2의 흡착량은 CO2에 비해 매우 작다고 가정하였으므로, 흡착상 열용량 및 흡착열 효과는 N2를 제외한 CO2 흡착에 의해 지배되는 것으로 보았다.
베드 내부에서 CO2와 N2의 질량전달은 대류뿐 아니라 분자확산 및 기계적 분산에 의해 영향을 받는다. 특히 고정층 흡착탑에서는 기체가 pellet 사이의 공극을 통과하면서 축방향 및 국부적인 농도 분산이 발생하므로, 종수송 방정식에는 유효 질량 분산계수 를 적용하였다. 낮은 Reynolds 수의 층류 조건에서는 다음의 상관식을 이용하여 각 성분의 질량 분산계수를 계산할 수 있다[15].
여기서 는 성분 의 질량 분산계수, 는 기체의 겉보기 속도, 는 ELM-11 pellet의 평균 입자 직경, 는 pellet 직경을 기준으로 정의한 Reynolds 수이다.
한편, ≥0.3 조건을 만족하는 경우에는 분자확산계수와 베드 공극률의 영향을 함께 고려한 다음의 상관식을 적용할 수 있다[16].
여기서 는 성분 의 분자확산계수, 𝜖은 베드 공극률, 는 Schmidt 수이다. Reynolds 수와 Schmidt 수는 각각 혼합기체의 밀도, 점도, 겉보기 속도, pellet 직경 및 분자확산계수를 이용하여 계산하였다.
2.9 모델 파라미터 및 해석 조건
본 연구에서 사용한 ELM–11 pellet의 물성값은 Table 1에 정리하였다. 입자 반경, 공극률, 충전 밀도 및 비열은 고정층 내부의 유동 저항, 흡착 소스항 및 에너지 방정식의 열용량 항을 계산하는 데 사용하였다. 특히 흡착에 의한 CO2 질량 소스항과 에너지 소스항은 ELM–11의 충전 밀도에 비례하도록 구성하였으며, pellet의 비열과 충전 밀도는 흡착열에 따른 고정층 온도 응답을 계산하는 데 적용하였다.
Table 1.
Physical and thermal properties of ELM-11 pellets
흡착 평형은 weighted dual-site Langmuir(wDSL) 모델을 이용하여 계산하였으며, 이에 사용된 주요 등온선 파라미터는 Table 2에 나타내었다. 본 연구에서는 문헌의 정적 등온선 기반 파라미터를 기본으로 사용하되, gate-open 상태의 포화 흡착량을 나타내는 는 breakthrough curve 기반 보정값을 적용하였다. 이에 따라 는 정적 등온선 기반 값인 3.635 mol/kg에서 3.110 mol/kg으로 수정하였다. 이러한 보정은 정적 등온선 측정 조건과 반복적인 fixed-bed breakthrough 조건에서 ELM–11의 유효 흡착 평형 특성이 완전히 동일하지 않을 수 있기 때문이다. 선행 연구에서는 반복 breakthrough test 과정에서 slipping-off plateau가 변화하며, 이를 ELM–11의 gate-opening pressure가 낮아지는 현상으로 해석하였다[12]. 또한 는 breakthrough curve fitting 과정에서 함께 재추정되어 반복 실험 중 관찰되는 흡착 평형 특성 변화를 반영하는 값으로 사용되었다[12]. 따라서 본 연구에서도 고정층 내 실제 CO2 breakthrough 및 온도 응답을 보다 현실적으로 재현하기 위해 보정된 값을 적용하였다.
3. 모델 검증 데이터
본 연구에서는 구축한 3차원 비정상 CO2 흡착 수치모델의 예측성을 평가하기 위해 참조 문헌[12]에 보고된 ELM-11 pellet 고정층 breakthrough 실험 데이터를 검증 데이터로 사용하였다. 해당 실험은 20% CO2/80% N2 혼합가스를 원통형 고정층 흡착탑에 공급하고, 시간에 따른 출구 CO2 농도와 흡착층 내부 85 mm 지점의 온도 변화를 측정한 것이다.
모델 검증에는 Table 3에 정리된 바와 같이 압력, 온도 및 선속도가 서로 다른 Run 1–8 조건을 사용하였다. 이 조건들은 수치모델의 운전 조건별 적용성을 평가하기 위한 validation dataset으로, CO2 분압, 기체 체류시간, 흡착 전선 이동 속도 및 흡착열 발생 조건의 변화를 포함한다. 따라서 본 검증은 단일 운전 조건에 대한 단순 재현이 아니라, 다양한 열·유동 조건에서 모델이 ELM-11의 gate-opening 기반 비선형 흡착 거동을 일관되게 예측할 수 있는지를 확인하는 데 목적이 있다.
각 조건에서 계산된 출구 CO2 breakthrough curve와 85 mm 위치의 온도 응답을 실험 결과와 비교함으로써, 흡착평형식, LDF 기반 흡착속도 모델 및 흡착열 소스항이 결합된 3차원 비정상 흡착 모델의 타당성을 검증하였다.
Table 3.
Operating conditions for fixed-bed CO2 adsorption breakthrough simulations
본 연구에서는 각 검증 case에 대해 시간에 따른 입구 대비 출구 CO2 농도비 C/C0와 흡착층 내부 온도비 T/Tamb를 비교 대상으로 설정하였다. CO2 breakthrough curve는 흡착 평형 및 물질전달 모델의 적절성을 평가하는 지표로 사용하였고, 온도 응답은 흡착열 및 에너지 소스항 구현의 타당성을 평가하는 지표로 사용하였다.
특히 ELM–11은 gate-opening 특성에 의해 S-shaped 흡착등온선과 slipping-off breakthrough 거동을 나타내므로, 단순한 breakthrough time뿐만 아니라 초기 plateau 구간, 급격한 농도 상승 구간, 그리고 흡착열에 의한 온도 피크의 재현성이 중요하다. 따라서 본 연구에서는 각 case에 대해 수치모델이 C/C0 곡선의 전체적인 형태와 T/Tamb의 피크 높이 및 감소 경향을 동시에 재현할 수 있는지를 중점적으로 검토하였다.
4. 결과 및 분석
4.1 격자 민감도 분석
수치해석 결과의 격자 의존성을 확인하기 위해 433,152, 1,489,600, 2,552,000개의 셀로 구성된 정렬 격자계를 사용하여 격자 민감도 분석을 수행하였다. 격자 민감도 평가는 출구 CO2 농도비(C/C0)와 흡착층 내부 온도비(T/Tamb)를 기준으로 수행하였다. 두 변수는 각각 흡착 평형 및 물질전달 모델에 의해 결정되는 파과 거동과, 흡착열 및 열전달 모델에 의해 결정되는 열적 응답을 대표하므로, 본 연구에서 구축한 비정상 흡착 모델의 격자 의존성을 평가하는 주요 지표로 사용하였다.
Fig. 4는 Run 2 조건에서 세 가지 격자계에 대해 계산된(T/Tamb)와 (C/C0)의 시간 변화를 비교한 결과이다. 전체적으로 모든 격자계에서 온도 응답과 CO2 파과곡선은 동일한 경향을 보였다. 격자 민감도가 크지 않은 주된 이유는 본 해석 조건이 낮은 Reynolds 수 영역에 해당하기 때문이다. 본 연구의 고정층 흡착탑 내부 유동은 (Re<10)의 층류 조건에서 진행되며, 유동장은 난류 와류나 강한 관성 효과보다는 다공성 매질 내 점성 저항과 축방향 대류에 의해 지배된다. 따라서 유동 구조가 복잡한 난류 문제에 비해 격자 해상도 변화가 속도장 및 압력장에 미치는 영향이 상대적으로 작다. 또한 본 모델에서는 실제 pellet 간 미세 유로를 직접 해상하지 않고, 흡착층을 평균화된 다공성 매질로 처리하였다. 이에 따라 입자 스케일의 유동보다는 셀 단위에서 계산되는 CO2 분압, 온도, 평형흡착량 및 LDF 기반 흡착속도가 전체 파과 거동을 결정하게 되며, 충분한 축방향 및 반경방향 해상도가 확보되면 격자 추가 세분화에 따른 변화는 제한적으로 나타난다. 따라서 본 연구에서는 계산 비용과 해석 정확도를 함께 고려하여 433,152 cells의 정렬 격자계를 최종 해석 격자로 선정하였다.
4.2 실험 파과곡선 및 온도 응답 검증
Fig. 5는 Run1–Run8 조건에 대해 문헌 실험 데이터와 본 연구의 3차원 비정상 수치해석 결과를 비교한 것이다. 각 조건에서 상단 그래프는 흡착층 내부 온도 응답을 나타내는 T/Tamb, 하단 그래프는 출구 CO2 농도비를 나타내는 C/C0이다. 전체적으로 수치해석 결과는 실험에서 관찰되는 CO2 파과곡선의 계단형 증가 거동과 흡착열에 따른 온도 상승 및 감소 경향을 재현하였다. 이는 본 연구에서 구축한 wDSL 기반 흡착 평형 모델과 LDF 흡착 속도 모델이 ELM–11 고정층 흡착의 주요 동적 거동을 합리적으로 모사하고 있음을 보여준다.
CO2 파과곡선에서 가장 두드러지는 특징은 일반적인 Langmuir형 흡착제에서 나타나는 단일 S자형 파과와 달리, 초기 농도 상승 이후 일정한 plateau를 거친 뒤 다시 급격히 증가하는 계단형 곡선이 나타난다는 점이다. 본 연구에서는 이러한 거동을 1차 파과, plateau 구간, 2차 파과로 구분하여 해석하였다. 초기 단계에서는 CO2가 흡착층 입구로 유입되지만, 낮은 CO2 분압 영역에서는 ELM–11의 gate-opening이 충분히 발생하지 않는다. 따라서 일부 CO2가 흡착되지 못하고 출구로 통과하면서 C/C0가 0에서 일정 수준까지 먼저 상승하는 1차 파과가 형성된다.
이후 흡착층 내부의 CO2 분압이 증가하면 ELM–11의 flexible 구조 확장이 진행되고, gate-open 상태로 전환되면서 평형 흡착량이 급격히 증가한다. 이 구간에서는 유입된 CO2가 추가적으로 흡착되기 때문에 출구 농도의 상승이 일시적으로 억제되며, 파과곡선에서 plateau가 형성된다. 즉, plateau 구간은 ELM–11의 구조 확장에 의해 새로운 흡착 용량이 활성화되는 구간으로 해석할 수 있다. 마지막으로 gate-open 상태의 흡착 용량이 점차 포화에 가까워지면, 흡착층 내부의 CO2 농도 전선이 출구 방향으로 이동하고 출구 C/C0가 다시 급격히 증가한다. 이 구간이 2차 파과에 해당한다.
온도 응답의 경우, 실험과 수치해석 모두 CO2 흡착 초기에 T/Tamb가 상승한 뒤 시간이 지남에 따라 서서히 감소하는 경향을 보였다. 이는 CO2 흡착에 의해 발생한 흡착열이 고정층 내부에 축적되고, 이후 유입 기체와 벽면 열손실에 의해 점차 완화되기 때문이다. 특히 온도 상승은 CO2 흡착이 활발하게 발생하는 흡착 전선 부근에서 나타나므로, 온도 피크의 발생 시점은 파과곡선의 plateau 및 2차 파과 구간과 밀접하게 연결된다.
Run별 비교 결과, 수치해석은 각 조건에서 계단형 파과곡선의 전체적인 형태와 온도 응답의 상승 및 감소 경향을 대체로 재현하였다. 유속이 큰 조건에서는 기체 체류 시간이 짧아져 파과가 상대적으로 빠르게 나타났으며, 유속이 작은 조건에서는 흡착 전선의 이동이 지연되어 2차 파과 시점이 늦어지는 경향을 보였다. 또한 압력 및 온도 조건 변화에 따라 CO2 분압과 gate-opening 조건이 달라지면서 plateau 구간과 온도 피크의 크기에도 차이가 나타났다. 이러한 경향은 ELM–11의 흡착 거동이 단순히 유속에 의해서만 결정되는 것이 아니라, 국부 CO2 분압, 온도, 평형 흡착량, 그리고 구조 확장 조건이 복합적으로 작용한 결과임을 의미한다.
추가적으로 Fig. 5에 나타난 실험값에 대한 수치해석 결과의 신뢰성을 정량적으로 평가하기 위해, 각 Run 조건에서 출구 CO2 농도비(C/C0)와 온도 응답(T/Tamb)에 대한 RMSE를 계산하였다. 오차 계산 시 수치해석 결과는 실험 데이터의 시간점으로 선형 보간하여 동일한 시간 기준에서 비교하였다. Table 4에 나타낸 바와 같이, C/C0에 대한 RMSE는 2.095×10-2 – 5.151×10-2 범위로 나타났으며, T/Tamb에 대한 RMSE는 3.781×10-2 – 9.666×10-2 범위로 계산되었다. 이러한 결과는 본 모델이 다양한 운전 조건에서 ELM–11 고정층의 계단형 CO₂ 파과 거동과 흡착열에 따른 온도 응답을 정량적으로도 합리적인 수준에서 재현하고 있음을 보여준다.
Table 4.
RMSE-based quantitative validation of CO₂ breakthrough curves and temperature responses under different operating conditions
그러나 대부분의 조건에서 온도 피크의 크기는 실험과 수치해석 사이에 차이를 보였다. 이는 등량흡착열 산정 과정의 불확실성에서 기인한 것으로 판단된다. 본 연구에서 적용한 등량흡착열은 온도별 흡착 등온선으로부터 Clausius–Clapeyron 관계를 이용해 계산되었기 때문에, 동일 흡착량 조건에서의 압력 보간 오차, 그리고 gate-opening 구간에서의 급격한 흡착량 변화에 민감하게 영향을 받을 수 있다. 특히 ELM-11은 hysteresis를 동반하는 flexible MOF이므로, 평형 등온선 기반으로 도출된 등량흡착열이 실제 비정상 breakthrough 과정에서 발생하는 순간적인 구조 전이와 이에 따른 열 발생 특성을 완전히 반영하지 못할 가능성이 있다. 그럼에도 불구하고 본 모델은 CO2 출구 농도와 온도 응답을 동시에 비교했을 때, ELM–11 고정층 흡착에서 나타나는 계단형 파과 특성과 흡착열에 따른 열적 응답을 전반적으로 재현하였다. 따라서 이후 절에서는 출구 응답만으로는 확인하기 어려운 고정층 내부의 CO2 농도 전선, 온도장, 흡착률 및 평형 흡착량 분포를 3차원 시각화를 통해 추가적으로 분석하였다.
4.3 3차원 CO2 농도장 및 온도장 분석
Fig. 6는 계단형 파과곡선과 plateau 구간이 가장 뚜렷하게 나타난 Run 4 조건에서 시간에 따른 3차원 CO2 농도비(C/C0)와 온도장(T)을 나타낸 것이다. 3차원 내부장 분석은 출구 breakthrough curve만으로는 확인하기 어려운 흡착 전선의 공간적 이동, gate-opening에 따른 국부 흡착 영역의 변화, 그리고 흡착열에 의한 열적 응답을 직접적으로 보여준다.
초기 25 s에서는 유입부 근처에서 높은 CO2 농도 영역이 형성되며, 축방향으로 갈수록 CO2 농도가 감소하는 분포가 나타난다. 이는 ELM-11의 gate-opening이 충분히 진행되기 전, 일부 CO2가 흡착되지 못한 채 충전층을 통과하는 slipping-off 현상과 관련된다. 이때 출구에 도달한 CO2는 breakthrough curve에서 관찰되는 1차 파과를 형성한다. 이후 150 s에서는 고농도 영역이 충전층 내부로 확장되고, 중간 영역에서 뚜렷한 농도 구배가 형성된다. 이 구간은 출구 파과곡선에서 plateau가 나타나는 시점과 대응된다. 즉, 초기 slipping-off에 의해 출구 CO2 농도는 먼저 상승하지만, 동시에 충전층 내부에서는 gate-opening에 따른 추가 흡착이 진행되면서 출구 농도의 급격한 증가가 억제된다. 따라서 plateau 구간은 초기 CO2 통과와 구조 전이에 따른 추가 흡착이 동시에 작용한 결과로 해석된다.
300 s 이후에는 CO2 농도 전선이 점차 출구 방향으로 이동하며, 고농도 영역이 충전층 후단부까지 확장된다. 450 s와 600 s에서는 대부분의 흡착층이 높은 CO2 농도 상태에 도달하고, 출구 C/C0 역시 1에 근접한다. 이는 breakthrough curve에서 2차 파과가 진행되는 구간과 일치한다. 따라서 3차원 농도장 결과는 Run 4에서 나타나는 계단형 파과곡선이 흡착층 내부의 농도 전선 이동과 ELM-11의 gate-opening 기반 구조 확장 거동에 의해 형성됨을 보여준다.
온도장에서도 농도 전선과 유사한 공간적 이동 특성이 확인된다. 초기에는 유입부 근처에서 CO2 흡착이 집중적으로 발생하면서 국부적인 온도 상승이 나타난다. 시간이 경과함에 따라 유입부 영역은 점차 흡착 포화 상태에 가까워지고, 흡착이 활발하게 일어나는 영역이 충전층 내부 및 출구 방향으로 이동한다. 이에 따라 고온 영역 역시 농도 전선과 함께 축방향으로 전파된다. 이러한 열적 응답은 CO2 흡착이 발열 과정이며, 흡착 전선의 이동이 열 발생 위치를 결정한다는 점을 나타낸다. 결과적으로 Fig. 6의 온도장은 실험 및 수치해석에서 관찰된 T/Tamb 피크와 이후 완만한 감소 경향이 흡착 전선의 이동, 국부 포화 진행, 그리고 충전층 내 열전달 과정의 복합적 결과임을 설명한다.
4.4 흡착률 분포를 통한 추가 흡착 영역 분석
Fig. 7은 Run 4 조건에서 시간에 따른 흡착률 와 평형흡착량 의 3차원 분포를 나타낸 것이다. Fig. 7의 좌측에 제시된 분포는 실제 CO2 흡착이 고정층 전체에서 균일하게 발생하지 않고, 농도 전선이 형성되는 제한된 영역에 집중됨을 보여준다. 이는 이미 CO2 흡착이 상당히 진행된 포화 영역이나, 아직 CO2가 충분히 도달하지 않은 미도달 영역에서는 흡착률이 낮고, 국부 CO2 분압이 gate-opening 조건에 도달하면서 추가 흡착이 활발하게 일어나는 전이 영역에서 흡착률이 가장 크게 나타나기 때문이다.
초기 25 s에서는 높은 흡착률이 유입부 근처에 국한된다. 이 시점에서는 유입부 인근의 CO2 분압이 먼저 증가하면서 일부 영역에서 gate-opening이 유도되고 강한 흡착이 발생하지만, 고정층 내부 대부분의 영역에서는 아직 CO2 분압이 충분히 상승하지 않아 흡착률이 낮게 유지된다. 이는 초기 CO2의 일부가 충분히 흡착되지 못하고 충전층을 통과하는 slipping-off 현상과 연관되며, 출구 breakthrough curve에서 관찰되는 1차 파과와 대응된다.
150 s와 300 s에서는 높은 영역이 유입부에서 고정층 중앙부로 이동한다. 이 구간은 출구 파과곡선에서 plateau가 형성되거나 2차 파과로 전환되는 시점에 해당한다. 출구 CO2 농도는 비교적 완만하게 증가하거나 일시적으로 제한된 상승을 보이지만, 3차원 흡착률 분포는 충전층 내부에서 흡착 전선이 계속 이동하며 추가 흡착이 활발히 진행되고 있음을 보여준다. 즉, plateau 구간은 내부 흡착이 정지된 상태가 아니라, gate-opening에 의해 새롭게 활성화된 영역에서 CO2가 지속적으로 흡착되는 동적 과정으로 해석된다.
450 s에서는 높은 흡착률 영역이 고정층 후단부로 이동하며, 이는 흡착 전선이 출구에 근접했음을 의미한다. 이후 600 s에서는 대부분의 영역에서 가 감소하는데, 이는 고정층 전반이 높은 흡착 포화도에 도달하면서 추가 흡착 여력이 감소했기 때문이다. 따라서 분포는 계단형 파과곡선의 plateau와 2차 파과가 고정층 내부에서 국부적으로 이동하는 흡착 전선과 ELM-11의 gate-opening 기반 추가 흡착 과정에 의해 형성됨을 보여준다.
Fig. 7의 우측은 시간에 따른 평형흡착량 분포를 나타낸다. 는 국부 CO2 분압과 온도 조건에서 ELM-11이 도달할 수 있는 평형 흡착 용량을 의미하며, 실제 순간 흡착량 와는 구별된다. 본 연구에서는 흡착등온선 특성에 근거하여 가 2.5 mol kg−1 이상인 영역을 ELM-11의 gate-opening 흡착 상태가 유도될 수 있는 구조 확장 가능 영역으로 해석하였다. 따라서 분포는 실제 흡착이 완료된 위치를 나타내기보다는, 국부 열역학 조건이 gate-open 상태를 형성할 수 있는지를 판단하는 지표로 사용된다.
25 s에서는 높은 영역이 유입부 근처에 제한적으로 형성된다. 이는 CO2가 먼저 유입부에 도달하면서 해당 영역의 국부 분압이 gate-opening 조건에 가까워졌기 때문이다. 150 s와 300 s에서는 > 2.5 mol kg−1인 영역이 고정층 내부로 점차 확장되며, 이는 CO2 농도 전선의 이동에 따라 gate-opening이 가능한 열역학적 영역도 함께 전파됨을 의미한다. 450 s에서는 높은 영역이 고정층의 상당 부분을 차지하고, 600 s에서는 거의 전체 흡착층으로 확장된다. 이는 고정층 대부분이 gate-open 흡착 상태를 유도할 수 있는 조건에 도달했으며, 추가 흡착 여력이 감소함에 따라 출구 C/C0가 1에 근접하는 2차 파과 이후의 거동과 연결된다.
다만 는 평형 상태에서의 흡착 가능량이므로, > 2.5 mol kg−1인 영역을 실제 구조 확장이 완료된 영역으로 직접 해석하는 것은 적절하지 않다. 실제 흡착 진행 정도는 순간 흡착량 와 흡착률 를 함께 고려해야 한다. 이에 본 연구에서는 분포를 통해 gate-opening이 유도될 수 있는 열역학적 영역을 확인하고, 분포를 통해 해당 영역 중 실제로 추가 흡착이 활발하게 발생하는 위치를 식별하였다. 이러한 결합 해석을 통해 ELM-11 고정층에서 계단형 파과곡선이 형성되는 내부 메커니즘을 공간적으로 규명할 수 있다.
5. 결 론
본 연구에서는 pelletized ELM-11이 충전된 원통형 고정층 흡착탑을 대상으로 3차원 비정상 CO2 흡착 수치모델을 구축하고, CO2/N2 혼합기체 조건에서의 파과 거동과 열적 응답을 해석하였다. ELM-11의 gate-opening 및 비선형 흡착 특성을 반영하기 위해 흡착평형식, LDF 기반 흡착속도 모델, 흡착열 소스항을 질량 및 에너지 보존 방정식에 결합하였다. 구축된 수치모델은 압력, 온도 및 속도 변화에 따른 출구 CO2 농도 변화와 흡착층 내부 온도 응답을 전반적으로 재현하였으며, 다양한 운전 조건에서 ELM-11의 비선형 흡착 및 계단형 파과 특성을 예측할 수 있음을 확인하였다.
3차원 농도장 분석을 통해 ELM-11 고정층의 계단형 파과곡선은 단순한 유동 지연이 아니라, 초기 slipping-off, gate-opening에 따른 추가 흡착, 흡착 전선 이동이 결합되어 형성됨을 확인하였다. 특히 plateau 구간에서는 출구 농도 변화가 제한적으로 나타나지만, 고정층 내부에서는 흡착 전선이 지속적으로 이동하며 추가 흡착이 진행되고 있었다. 온도장 분석에서는 CO2 흡착에 따른 발열 영역이 농도 전선과 함께 유입부에서 출구 방향으로 이동하는 것을 확인하였다. 이는 실험 및 해석에서 관찰되는 온도 피크가 흡착 전선 이동, 국부 포화, 열전달 과정의 상호작용에 의해 형성됨을 의미한다.
또한 와 분포를 통해 실제 흡착이 활발히 발생하는 영역과 gate-opening이 가능한 열역학적 영역을 구분할 수 있었다. 는 실제 추가 흡착이 집중되는 흡착 전선 위치를 나타내며, 는 국부 조건에서 구조 확장이 유도될 수 있는 영역을 보여준다. 결과적으로 본 연구는 기존 1차원 파과곡선 중심 해석으로는 확인하기 어려운 ELM-11 고정층 내부의 흡착 전선 이동, 계단형 파과 형성 메커니즘 및 흡착열에 따른 열적 응답을 3차원적으로 규명하였다. 제안된 수치모델은 향후 flexible MOF 기반 CO2 포집 공정의 설계 및 운전 조건 최적화에 활용될 수 있다.









