NUMERICAL STUDY ON MASS TRANSFER AND HYDROGEN BUBBLE GENERATION IN WATER ELECTROLYSIS

J. Park; G. Son

doi:10.6112/kscfe.2024.29.3.173

Preview

Research Article

Journal of Computational Fluids Engineering. 30 September 2024. 173-185
https://doi.org/10.6112/kscfe.2024.29.3.173

NUMERICAL STUDY ON MASS TRANSFER AND HYDROGEN BUBBLE GENERATION IN WATER ELECTROLYSIS

수전해에서 물질 전달과 수소 기포 생성에 대한 수치적 연구

J. Park¹

G. Son¹^†

박 정주¹

손 기헌¹^†

¹Department of Mechanical Engineering, Sogang 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

In alkaline water electrolysis(AWE), hydrogen bubbles are generated through electrochemical reactions involving the interaction between continuous and dispersed phases. In this work, a multiphase numerical model for AWE is developed to simulate the evolution of hydrogen bubble associated with the transport of dissolved hydrogen in liquid water. The numerical results show that bubble nucleation under supersaturation conditions of dissolved hydrogen gas is a critical parameter in determining gas evolution in AWE. The effects of ambient temperature and current density on bubble evolution are investigated. As current density increases and bubble resistance rises, the ohmic overpotential is observed to dominate over the activation overpotential.

Keywords

Alkaline water electrolysis

Hydrogen production

Concentration

Two-Fluid model

키워드

알칼라인 수전해

수소 생산

농도

MAIN

1. 서 론
2. 수치해석 기법
2.1 지배방정식
2.2 경계조건
3. 결과 및 고찰
3.1 수전해의 기포 유동에 대한 의 효과
3.2 수전해의 기포 유동에 대한 온도 효과
3.3 수전해의 기포 유동에 대한 전류 밀도 효과
4. 결 론

1. 서 론

알칼라인 수전해(Alkaline Water Electrolysis, AWE)는 물을 전기 분해하여 수소를 생산하는 친환경적인 공정으로, 높은 효율성과 비용 절감 효과로 인해 대규모 생산 공정에서 널리 활용되고 있다. 수전해를 통하여 수소를 생산하는 과정에서, 전극 부근의 화학반응으로 인해 수소는 기포로 형성되며, 이는 전해질과 전극의 접촉을 방해하고 전해질 내 이온의 이동을 저해하는 저항으로 작용한다[1]. 이러한 저항은 전류 밀도와 온도가 높아짐에 따라 전압이 평형 상태 이상으로 크게 증가하는 과전압(overpotential)으로 나타나며, 수소 생산에 필요한 에너지를 증가시키는 문제를 초래한다[2]. 따라서 수전해 과정에서 발생하는 기포를 포함한 다상유동에 대한 해석은 수소 생산 효율을 향상시키기 위한 방안을 도출하는 데 매우 중요하다.

수전해에 대한 기존의 해석적인 연구에서는 전해질을 연속 상(continuous phase)으로 기포를 분산 상(dispersed phase)으로 구분하여 해석하는 two-fluid 방법을 제시하고 있다. Mat 등[3]은 수직 및 수평 수전해 셀에서의 유동 특성과 수소 발생에 의한 기포 분포를 해석하기 위해 2상 유동 모델을 적용하였다. 그러나 이 연구에서는 수전해 실험에서 관찰되는 난류 효과에 의한 기포 상승에 수직한 방향으로 분산되는 힘을 고려하지 않았다. El-Askary 등[4]은 수전해 과정에서의 수소 기포 유동을 예측하기 위해 난류 분산력이 추가된 다상 모델을 제시하였으며, 전류 밀도가 증가하고 양극과 음극 사이의 간격이 감소함에 따라 수소 생산량이 증가함을 확인하였다. Zarghami 등[5], Li 등[6], Cao 등[7]은 수소 기포 유동 해석에서 수소뿐만 아니라 액상의 물로부터 증발되는 수증기를 포함시킴으로써 기포 분포를 보다 정확하게 예측하였다. Park and Son[8]은 two-fluid 방법을 기반으로 전기화학 반응으로 생성된 수소량을 패러데이 법칙에 따라 모델링하고, 전류 밀도와 기포 발생이 과전압에 미치는 영향을 분석하였다. 대부분의 기존 연구에서는 전기화학 반응으로 생성된 수소가 전부 기포로 변환된다고 가정하였으며, 이는 용존 수소가 포화 농도에 도달하지 않았음에도 기포가 형성된다는 가정을 바탕으로 하고 있어 물리적으로 적절하지 않다. 최근에는, 수전해 과정에서 용존 수소 농도와 기포 발생 조건을 함께 고려하는 해석기법이 제시되고 있다. Lee 등[9]은 mixture 모델을 기반으로 농도장 방정식에 페러데이 법칙을 적용하여 용존 수소 생성항을 추가하고, 기포 생성이 시작되는 임계 농도에 대한 모델링을 추가하였다. 그러나 그들의 연구에서는 난류 효과와 기포 내부의 수증기 효과가 포함되지 않았다. Moshin 등[10]은 two-fluid 방법과 VOF(volume-of-fluid) 방법을 결합하고 용존 수소 가스의 농도장 방정식을 추가하는 포괄적인 해석기법을 제안하였다. 이와 같은 수전해에 대한 많은 연구에도 불구하고, 수소 생산 효율을 저해하는 기포 발생에 대한 연구는 전기화학과 다상역학 즉 화학과 기계 분야의 다학제적 요소로 인해 체계적인 연구는 상당히 미진한 실정이다. 특히 수소 기포의 발생과 성장을 포함한 다상유동에 대한 적절한 해석기법의 개발은 여전히 부족하다.

본 연구에서는 수전해 과정의 수소 기포의 발생과 성장에 집중하여, two-fluid 모델을 기반으로 전해질 영역에서 용존 수소 농도의 해석과 과포화 조건에서의 기포 발생을 포함한 다상유동에 대한 보다 체계적인 해석 연구를 수행하고자 한다.

2. 수치해석 기법

본 연구에서는 OpenFOAM(Open-source Field Operator and Manipulation) 중 two-fluid 모델 기반인 reactingMultiphaseEulerFoam[11]에 용존 수소의 농도장 방정식을 추가한다. 해석 영역은 2차원으로, 수소 발생 영역만을 고려하여, Fig. 1과 같이 음극과 유체 영역으로 구성된다. 여기서 $H_{f}$ 는 채널의 높이, $H_{c}$ 는 전극의 길이, $W_{f}$ 는 채널의 폭을 의미한다. 본 연구의 전해질은 1M KOH 수용액이며, 전극에 인접한 유체 영역에서 전기화학 반응을 통해 수소가 발생하며, 일부는 전해질 내부로 용해되고, 용존 수소의 농도가 포화 농도보다 큰 경우 기포로 발생된다고 가정한다. 유체 영역 전체는 등온 조건이며 기상과 액상 모두 비압축성으로 가정한다. 기포는 수소뿐만 아니라 액상의 물에서 증발하는 수증기의 혼합물로 구성되며, 열역학적 평형 상태에 있다고 가정한다. 기포 내 수증기의 분압은 포화압력으로 나타낼 수 있으며, 수소의 분압은 대기압에 수증기 분압을 차감하여 구할 수 있다. 유체 영역의 온도가 일정하다고 가정하면, 수증기와 수소의 분압비는 일정하게 유지되며, 이에 따라 기포의 조성비 또한 일정함으로 밀도가 변하지 않는다고 가정할 수 있다.

https://cdn.apub.kr/journalsite/sites/kscfe/2024-029-03/N0500290313/images/jkscfe_2024_293_173_F1.jpg

Fig. 1.

Schematics of an alkaline water electrolysis cell

2.1 지배방정식

유체 영역에서의 액상과 기상에 대한 질량 및 운동량 보존 방정식은 다음과 같이 표현된다.

(1)

\frac{\partial α_{g} ρ_{g}}{\partial t} + \nabla ∙ α_{g} ρ_{g} u_{g} = S_{g}

(2)

\frac{\partial α_{l} ρ_{l}}{\partial t} + \nabla ∙ α_{l} ρ_{l} u_{l} = - S_{g}

(3)

\frac{\partial α_{i} ρ_{i} u_{i}}{\partial t} + \nabla \cdot α_{i} ρ_{i} u_{i} u_{i} = - α_{i} \nabla p + \nabla ∙ α_{i} τ_{i} + α_{i} ρ_{i} g + f_{i, m}

여기서, $S_{g}$ 는 기포 생성항이며 전극에 인접한 셀에서만 적용한다. 운동량 보존 방정식에서의 하첨자 i는 전해질(l)과 기상(g)을 의미하며 $f_{i, m}$ 은 각 상 사이에 작용하는 계면 운동량 전달항이고 각 상에 대한 점성응력( $τ_{i}$ )은 식 (4)과 같다.

(4)

τ_{i} = μ_{i, c f f} [(\nabla u_{i} + \nabla u_{i}^{T}) - \frac{2}{3} (\nabla \cdot u_{i}) I]

기존 연구에서는 전극에서 발생한 수소가 모두 기포로 변환된다고 가정한 반면, 본 연구에서는 전극에서 발생한 수소 중 일부는 전해질에 용존되고 포화 상태를 초과하는 조건에서만 수소 기포로 변환된다는 모델링을 적용하였다. 전해질 내에 존재하는 화학종은 수소( $H_{2}$ )와 물( $H_{2} O$ )로 구성되며, 액상에서 용해된 수소의 질량 분율을 $Y_{l}$ 로 정의하면, 액상에서 용존 가스의 농도장 방정식은 다음과 같다.

(5)

\frac{\partial α_{l} ρ_{l} Y_{l}}{\partial t} + \nabla ∙ α_{l} ρ_{l} Y_{l} u_{1} = \nabla ∙ α_{1} \frac{μ_{1, eff}}{{Sc}_{t}} \nabla Y_{1} + S_{1}

2.1.1 기포 생성과 용존 수소 가스의 생성 항

식 (1)의 기포 생성항( $S_{g}$ )과 식 (5)의 용존 수소 가스 생성항( $S_{l}$ )은 전극에서 전기화학 반응으로 발생한 수소량에 의해 결정되며, 이는 전극에 인접한 계산 셀에서만 적용된다. 발생한 수소 중 일부는 기포로 변환되고, 나머지는 전해질에 용존되는 것을 모델링하기 위해, 본 연구에서는 포화 수소 농도( $Y_{l, s a t}$ ) 대비 기포가 발생하는 수소 농도( $Y_{l, n u c}$ )의 비를 𝜉로 정의한다. $Y_{l, s a t}$ 는 KOH 수용액의 포화압력( $p_{w, K O H}$ )과 Henry 법칙에 의해 용존 수소의 포화 몰 분율( $x_{l, s a t}$ )을 구하고 질량비로 아래와 같이 변환하여 구할 수 있다.

(6)

χ_{l, s a t} = \frac{p_{a t m} - p_{w, K O H}}{H}

(7)

Y_{l, s a t} = \frac{M_{H_{2}} χ_{l, s a t}}{M_{H_{2} O} (1 - χ_{l, s a t}) + M_{H_{2}} χ_{l, s a t}}

여기서, 기포는 대기압( $p_{a t m}$ )으로 가정하고, $H$ 는 Henry 상수, $M_{H_{2}}$ 와 $M_{H_{2} O}$ 는 수소와 수용액의 분자량이다.

기포 생성항( $S_{g}$ )과 용존 수소 가스 생성항( $S_{l}$ )는, 전극 부근에서 생성된 수소의 농도가 포화 농도 대비 𝜉 이하이면 기포 발생없이( $S_{g}$ =0) 전해질에 용해되고, 𝜉이상이면 전해질에 용해되지 않고( $S_{l}$ =0) 기포로 변환된다는 가정하에, 다음과 같이 모델링하였다. 용존 수소의 농도가 과포화 상태에 도달하여 수소 기포가 생성될 경우 기포 주변의 물이 함께 증발하여 수증기가 기포에 혼합되며, 이때 수소 대비 수증기 발생량은 기포 내 수소 대비 수증기의 질량비 𝛽와 같다고 가정하였다. 이에 따라 혼합 기포 생성항( $S_{g}$ )은 𝛽를 포함하여 식 (8)과 같이 나타낼 수 있다.

(8)

S_{g} = \{\begin{cases} \begin{array}{cl} 0 & if Y_{l} \leq ξ Y_{l, s a t} \end{array} \\ \begin{array}{cl} (1 + β) \frac{| i |}{n_{ϵ} F w} M_{H_{2}} & if Y_{l} > ξ Y_{l, s a t} \end{array} \end{cases}

(9)

S_{l} = \{\begin{cases} \frac{| i |}{n_{c} F w} M_{H_{2}} & if Y_{l} \leq ξ Y_{l, s a t} \\ 0 & if Y_{l} > ξ Y_{l, s a t} \end{cases}

여기서, 전류 밀도 i는 외부 변수로 취급하며, $n_{e}$ 는 화학반응에서 수소 대비 전자의 몰비, F는 Faraday 상수, $w$ 는 계산 셀의 폭, 기포 내의 수소 대비 수증기의 질량비 𝛽는 기포 성장 과정에서 기상과 액상 사이에 열역학적 평형이 지속적으로 유지된다고 가정하면 다음과 같이 표현할 수 있다.

(10)

β = \frac{p_{w, K O H} M_{H_{2} O}}{(p_{a t m} - p_{w, K O H}) M_{H_{2}}}

2.1.2 계면 운동량 전달 항

운동량 보존 방정식에서 각 상에 작용하는 생성 항은 다음과 같다.

(11)

f_{g, m} = - f_{1, m} = f_{drag} + f_{td}

(12)

f_{d r a g} = - \frac{3}{4} \frac{ρ_{l} α_{g} α_{l}}{d_{b}} C_{D} | u_{g} - u_{l} | (u_{g} - u_{l})

(13)

f_{t d} = - \frac{3}{4} \frac{C_{D}}{d_{b}} α_{g} | u_{g} - u_{l} | \frac{μ_{l, t}}{S c_{t}} (\frac{1}{α_{g}} + \frac{1}{α_{l}}) \nabla α_{g}

계면에서 작용하는 힘( $f_{g, m}, f_{l, m}$ )은 일반적으로 항력(drag)과 난류 분산력(turbulent dispersion force), 벽 윤활력(wall lubrication), 양력(lift), 가상 질량력(virtual mass)으로 구성된다. Haug 등[12]의 연구에 따르면, 전류 밀도가 1000 $A / m^{2}$ 이상인 경우 기포의 평균 직경이 100 $μ m$ 로 수렴하는 것으로 보고되었다. 본 연구의 알칼라인 수전해에서 인가된 전류 밀도는 최소 1500 $A / m^{2}$ 이므로, 기포의 직경을 100 $μ m$ 로 일정하게 가정하였다. 이 때 기포에 관한 레이놀즈 수( $R e_{b} = d_{b} | u_{g} - u_{l} | / ν_{l}$ )는 대략 11.4로 계산되며 생성된 기포와 유속의 크기가 매우 작기 때문에 시간에 따른 유동 변화는 미미하다고 볼 수 있다[4,5,8]. 따라서 식 (11)과 같이 벽 윤활력, 양력, 가상 질량력 효과를 무시하고 유체 영역에서 움직이는 기포에 의하여 유발되는 항력과 난류 분산력만 고려한다. 여기서, $C_{D}$ 는 Schiller & Naumann 모델을, 난류 분산력의 경우에는 Burns 모델을 적용하였다.

2.1.3 난류 효과

유체 영역에서의 난류 효과는 standard k-ε 모델을 사용하였다.

(14)

\frac{\partial α_{i} ρ_{i} k_{i}}{\partial t} + \nabla \cdot α_{i} ρ_{i} u_{i} k_{i} = α_{i} ρ_{i} (P_{i, k} - ε_{i}) + \nabla ∙ α_{i} (μ_{i} + \frac{μ_{i, t}}{σ_{k}}) \nabla k_{i}

(15)

\frac{\partial α_{i} ρ_{i} ε_{i}}{\partial t} + \nabla ∙ α_{i} ρ_{i} u_{i} ε_{i} = α_{i} ρ_{i} (c_{ε, 1} P_{i, k} - c_{ε, 2} ε_{i}) + \nabla ∙ α_{i} (μ_{i} + \frac{μ_{i, t}}{σ_{ε}}) \nabla ε_{i}

(16)

μ_{i, t} = ρ_{i} C_{μ} \frac{k_{i}^{2}}{ε_{i}}

식에서의 하첨자 i는 기상(g)과 액상(l)을 나타내고 $P_{i, k}$ 는 난류 생성 항이며 모델 계수는 기본값을 사용하였다.

2.2 경계조건

채널의 하단부에는 일정 유입 속도 조건을, 상단부는 압력 출구 조건을 설정하였다. 전극의 표면, 채널의 벽과 분리막을 포함한 모든 벽에 표준 벽함수 난류 조건과 점착 조건을 적용하였다. 유체 영역에 식 (5), (1)과 (2), (3)의 순서로 지배 방정식을 풀게 되고, 운동량 보존식에는 PISO 알고리즘을 적용하였다. 시간 항은 implicit Euler 기법을 사용하였다. 식 (5)에 대해서는 Gauss limitedLinear 기법을, 식 (1)과 (2)의 대류 항은 Gauss vanLeer 기법을, 나머지 식에 대해서는 QUICK 기법을 사용하였고 확산 항은 중앙 차분 기법을 이용하였다.

2.2.1 격자 의존도

본 연구의 해석에 대한 격자의 의존성 검증을 위해 수직 방향에 대한 단위 격자의 수를 600개로 고정하고 수평 방향에 대한 단위 격자의 개수를 80,160,320 세 가지로 수행하였다. 전극의 상단부(y=0.5m)를 기준으로 수평 방향에 따른 $α_{g}$ 를 측정하였다. 결과는 Fig. 2와 같으며 $∆_{x} = W_{f} / 160$ 160과 $∆_{x} = W_{f} / 320$ 320일 때의 결과가 거의 일치함을 볼 수 있다. 따라서 단위 격자의 크기는 $0.05 m m \times 1 m m$ 로 고정한 후 수전해 셀 내부 2상 유동을 해석하였다.

https://cdn.apub.kr/journalsite/sites/kscfe/2024-029-03/N0500290313/images/jkscfe_2024_293_173_F2.jpg

Fig. 2.

Dependence of grid size on the prediction of $α_{g}$ profile at $y = 0.5 m$ with $i = 1500 A / m^{2}$

3. 결과 및 고찰

계산 영역은 2차원이며 Riegel[13]의 실험 조건을 모사하여 $i = 1500 A / m^{2}$ , 채널의 폭은 $8 m m$ , 높이는 $0.6 m$ 로 구성된다. 전극의 길이는 $0.4 m$ 이며 중심점은 채널 높이의 중심점( $y = 0.3 m$ )과 일치한다. 수전해 셀 내부 온도는 일정하게 유지되며, 작동 압력은 대기압으로, 전해질은 1M KOH를 사용하였으며, 계산에 사용된 물성은 Table 1과 같다[14]. 입구에서는 일정한 속도( $u_{y, i n} = 0.69 m / s$ )를 부여하였으며, 수증기와 수소가 결합된 기포가 음극(cathode)에서 생성되어 출구로 배출되는 과정을 계산하였다. 이 때 기포의 반경은 100 $μ m$ 로 고정하였으며, 식 (8)과 (9)에서 셀의 폭(w)는 기포의 반경과 동일하게 설정되었다.

Table 1.

Physical properties used in AWE calculation

Properties	Value
$ρ_{l} (k g / m^{3})$	1026
$ρ_{g} (k g / m^{3})$	0.25
$μ_{l} (P a ∙ s)$	0.45×10^-3
$μ_{g} (P a ∙ s)$	1.27×10^-5
$Y_{l, s a t}$	9.91×10^-7
𝛽	4.24
$M_{H_{2}} (k g / m o l)$	0.002
$i (A / m^{2})$	1500

3.1 수전해의 기포 유동에 대한 𝜉의 효과

용존 상태의 수소의 농도가 포화 상태에 도달하여 기포로 생성되는 경우(𝜉=1)를 고려하기에 앞서, 전극 표면에서 발생된 수소가 기포로 전혀 변환되지 않는 경우(𝜉>10)와 모두 기포로 변환되는 경우(𝜉=0)에 대한 계산을 수행하였으며, 정상상태에서의 결과는 Fig. 3과 같다. 𝜉>10의 경우 단상 유동이며, 이로 인해 속도, 난류 에너지, 난류 동점성 계수 필드는 거의 대칭을 이루는 것을 확인할 수 있다. 용존 수소의 질량비의 경우 포화 상태에 비해 4.6배까지 상승하는 것을 볼 수 있다. 반면 𝜉=0의 경우 기포가 발생하면서 속도, 난류 에너지, 난류 역학 점도 필드는 비대칭을 이루게 된다. 전극 표면에서 형성된 기포에 의하여 벽쪽의 속도 구배를 감소시키기 때문에 난류 에너지가 𝜉>10의 경우 대비 기포가 존재하는 영역에서 작아지며 전해질의 속도는 증가함을 볼 수 있다.

https://cdn.apub.kr/journalsite/sites/kscfe/2024-029-03/N0500290313/images/jkscfe_2024_293_173_F3.jpg

Fig. 3.

Effect of 𝜉 on the flow fields in AWE

Fig. 4는 𝜉에 따른 전극과 전해질 경계면에서의 용존 수소 질량비( $Y_{l}$ )와 기상 체적비( $α_{g}$ )를 나타낸다. Fig. 4(a)에 따르면, 전극이 존재하는 구간( $0.1 m \leq y \leq 0.5 m$ )에서 𝜉>10일 경우 용존 수소의 질량비가 포화 상태에 비해 4.6배까지 상승하는 반면, 𝜉=1일 경우 용존 수소의 질량비가 거의 포화 상태로 유지된다. Fig. 4(b)에서는, 𝜉=1일 때 기포가 차지하는 영역이 𝜉=0에 비해 상당히 작음을 볼 수 있다. 이는 전극에서 발생한 수소 중 상당 부분이 기포로 변환되지 않고 전해질에 용해되어 있음을 의미한다.

https://cdn.apub.kr/journalsite/sites/kscfe/2024-029-03/N0500290313/images/jkscfe_2024_293_173_F4.jpg

Fig. 4.

Effect of 𝜉 on $Y_{l}$ and $α_{g}$ at the electrode-fluid interface in AWE

3.2 수전해의 기포 유동에 대한 온도 효과

Figs. 5와 6은 온도에 따른 $α_{g}$ 분포와 전극 상단부(y=0.5m)에서 채널의 폭방향에 따른 $α_{g}$ 의 분포를 나타내며, Riegel 등[13]의 실험 데이터는 경향성을 비교하기 위해 포함하였다. 온도가 70°C에서 80°C로 증가하면, $p_{w, K O H}$ 가 32kPa에서 50kPa로 증가한다. 전류 밀도가 일정한 상태에서 온도가 중가하면, 식 (8)의 수소 생성항은 일정하지만 기포 내의 수소 대비 수증기의 질량비인 𝛽가 4.24에서 9로 증가하여 수소와 수증기가 혼합된 수소 기포 발생량이 증가하고, 결과적으로 $α_{g}$ 가 증가한다. 온도가 증가함에 따라 $Y_{l, s a t}$ 는 9.91×10^-7에서 7.36×10^-7로 감소하고, 이로 인해 𝜉=0과 𝜉=1의 $α_{g}$ 의 차이는 온도가 증가할수록 감소하는 경향을 보인다. 본 연구의 해석 결과는 80°C에서 Riegel 등[13]의 실험 데이터와 잘 일치함을 볼 수 있다. Riegel 등의 실험 데이터는 50°C에서 얻어진 것으로 보고되었으나, 이후 다수의 해석 연구[5,6]에서는 80°C의 경우에 해당 데이터가 더 잘 일치한다고 보고하고 있다. 기포 내 수소의 양을 온도에 따라 비교하면, Fig. 7에 나타난 바와 같이 𝜉=0인 경우에는 거의 변화가 없으며, 𝜉=1인 경우에는 온도가 증가함에 따라 약간 증가함을 볼 수 있다. 이는 $Y_{l, s a t}$ 의 감소로 기포로 변환되지 않고 전해질에 용해되어 있는 용존 수소의 양이 상대적으로 감소하기 때문이다.

https://cdn.apub.kr/journalsite/sites/kscfe/2024-029-03/N0500290313/images/jkscfe_2024_293_173_F5.jpg

Fig. 5.

Effect of temperature on $α_{g}$ fields in AWE. The white solid lines represent $α_{g}$ =0.01 contours

https://cdn.apub.kr/journalsite/sites/kscfe/2024-029-03/N0500290313/images/jkscfe_2024_293_173_F6.jpg

Fig. 6.

Effect of temperature on $α_{g}$ profiles at y=0.5m in AWE

https://cdn.apub.kr/journalsite/sites/kscfe/2024-029-03/N0500290313/images/jkscfe_2024_293_173_F7.jpg

Fig. 7.

Effect of temperature on $α_{H_{2}}$ profiles at y=0.5m in AWE

3.3 수전해의 기포 유동에 대한 전류 밀도 효과

Figs. 8과 9는 $T$ =80℃에서 인가된 전류 밀도가 각각 3250, 6250 $A / m^{2}$ 일 때 기포 유동을 나타낸다. 전극 길이( $H_{c}$ )는 Riegel 등[13]의 실험조건을 고려하여 $0.4 m$ 에서 $0.12 m$ 로 변경하였다. 전류 밀도가 3250 $A / m^{2}$ 에서6250 $A / m^{2}$ 로 증가함에 따라 $α_{g}$ 경계층 두께의 평균값이 𝜉=0일 때 17%, 𝜉=1일 때 21% 증가하는 것을 볼 수 있다. 𝜉=0과 𝜉=1인 경우 평균 $\bar{α_{g}}$ 의 차이를 비교하면 3250 $A / m^{2}$ 인 경우 10%, 6250 $A / m^{2}$ 의 경우 4%로 작아지며, 이는 포화 상태의 용존 수소 농도는 7.36×10^-7로 동일하지만, 기포로 변환되는 양이 전류 밀도의 증가와 함께 상대적으로 증가되기 때문이다. 전류 밀도에 따른 본 해석 결과는 Riegel 등[13]의 실험 데이터와 잘 일치함을 볼 수 있다.

https://cdn.apub.kr/journalsite/sites/kscfe/2024-029-03/N0500290313/images/jkscfe_2024_293_173_F8.jpg

Fig. 8.

Effect of current density on $α_{g}$ fields. The white solid lines represent $α_{g}$ =0.01 contours

https://cdn.apub.kr/journalsite/sites/kscfe/2024-029-03/N0500290313/images/jkscfe_2024_293_173_F9.jpg

Fig. 9.

Effect of current density on $α_{g}$ profiles at y=0.5m in AWE

Fig. 10은 𝜉=1 조건에서 전류 밀도를 15000 $A / m^{2}$ 까지 증가시켰을 때의 결과를 나타낸다. 전류 밀도가 증가함에 따라 전극 상단부에서의 $α_{g}$ 는 0.85까지 증가하며, 기포가 많아짐에 따라 부력에 의한 전해질의 유동방향 속도가 전극 부근에서 두드러지게 증가하는 것을 볼 수 있다.

Fig. 11은 𝜉=1 조건에서 전류 밀도에 따른 활성화 과전압 $η_{c, a c t}$ 과 저항 과전압 $η_{c, o h m}$ 의 크기의 변화를 보여준다. 여기서, 점선은 기포 없이 전해질만 존재하는 경우이며, x방향에 대하여 1차원 해를 구하였다[8]. 기포에 따른 과전압의 변화는 다음과 같이 구할 수 있다[15,16].

(17)

η_{c, a c t} = \frac{1}{b_{c}} \ln (\frac{i}{α_{l} i_{c, 0}})

(18)

η_{c, ohm} = - W_{f} \frac{i}{κ_{f}}

(19)

κ_{f} = κ_{l} {(1 - α_{g})}^{3 / 2}

여기서, $i_{c, 0}$ 는 음극의 교환 전류 밀도로 3.192 $A / m^{2}$ , $b_{c}$ 는 음극 Tafel 기울기로 –15.15 $V^{- 1}$ , $W_{f}$ 는 채널의 폭, 기포를 포함한 전해질의 전기 전도도 κ_f는 순수 전해질의 전기 전도도 (κ_l38.95 $S / m$ [17])와 기포의 효과를 반영하기 위해 Bruggeman 방정식[18]을 활용하였다. 전류 밀도가 증가하면, 기포 생성과 함께 전극 표면에서의 활성화(activation) 과전압은 미세하게 증가하여 $A / m^{2}$ 일 때 전해질만 존재하는 경우보다 0.09 $V$ 정도 크며, 옴(ohmic) 과전압은 0.5 $V$ 까지 상대적으로 크게 증가하는 것을 볼 수 있다.

https://cdn.apub.kr/journalsite/sites/kscfe/2024-029-03/N0500290313/images/jkscfe_2024_293_173_F10.jpg

Fig. 10.

Effect of current density on $α_{g}$ and streamwise velocity profile at $y = 0.5 m$

https://cdn.apub.kr/journalsite/sites/kscfe/2024-029-03/N0500290313/images/jkscfe_2024_293_173_F11.jpg

Fig. 11.

Effect of current density on overpotentials

4. 결 론

본 연구에서는 수전해 과정의 수소 기포를 포함한 다상유동 해석을 위해, two-fluid 모델을 기반으로 전해질 영역에서 용존 수소 농도의 해석과 과포화 조건에서의 기포 발생을 포함한 모델링을 추가하였다. 또한, 기포는 수소뿐만 아니라 액상의 물에서 증발하는 수증기의 혼합물로 가정하는 열역학적 모델도 적용하였다. 포화 수소 농도 대비 기포가 발생하는 수소 농도의 비인 𝜉로 도입하여, 전극 표면에서 발생된 수소가 모두 기포로 변환되는 경우(𝜉=0)와 용존 수소 농도가 포화 상태에 도달하여 기포로 생성되는 경우(𝜉=1)에 대한 결과를 분석하고 비교하였다. 𝜉=1일 때 기포가 차지하는 영역이 𝜉=0에 비해 상당히 작으며, 전극 상부에서 배출되는 기포의 양이 적음을 확인할 수 있었다. 기포 분포의 해석 결과는 문헌에 보고된 실험 데이터와 잘 일치함을 확인하였다. 온도와 전류밀도가 증가함에 따라 기포 발생량이 증가하지만, 𝜉에 따른 영향은 상대적으로 감소하였다. 전류 밀도가 일정한 경우, 기포 내 수소의 양을 온도에 따라 비교하면, 𝜉=0인 경우에는 거의 변화가 없으며, 𝜉=1인 경우에는 온도가 증가함에 따라 약간 증가함을 확인할 수 있었다. 전류 밀도가 증가할수록, 기포가 전극 표면의 활성화 과전압보다 전해질에서의 옴(저항) 과전압에 더 크게 영향을 미침을 확인하였다.

Acknowledgements

이 성과는 정부(과학기술정보통신부)의 재원으로 한국연구재단의 지원을 받아 수행된 연구임(RS-2024-00350110).

References

2014, Wang, M., Wang, Z., Gong, X. and Guo, Z., "The Intensification Technologies to Water Electrolysis for Hydrogen Production-A Review," Renew. Sustain. Energy Rev., Vol.29, pp.573-588.

10.1016/j.rser.2013.08.090

1984, Krenz, M., Untersuchung des elektrodennahen Raumes gasentwickelnder Elektroden. Diss.

2004, Mat, M.D., Aldas, K. and Ilegbusi, O.J., "A two-phase flow model for hydrogen evolution in an electrochemical cell," Int. J. Hydrogen Energy, Vol.29, pp.1015-1023.

10.1016/j.ijhydene.2003.11.007

2015, El-Askary, W.A., Sakr, I.M., Ibrahim, K.A. and Balabel, A., "Hydrodynamics characteristics of hydrogen evolution process through electrolysis: Numerical and experimental studies," Energy, Vol.90, pp.722-737.

10.1016/j.energy.2015.07.108

2023, Zarghami, A., Deen, N.G. and Vreman, A.W., "CFD modeling of multiphase flow in an alkaline water electrolyzer," Chem. Eng. Sci., Vol.227, 115926.

10.1016/j.ces.2020.115926

2024, Li, J., Liang, H., Rehan, M.A. and Li, G., "Numerical simulation and analysis of two-phase flow model considering bubble coverage for alkaline electrolytic water," Appl. Therm. Eng., Vol.245, 122890.

10.1016/j.applthermaleng.2024.122890

2024, Cao, X., Zhao, N., Zhang, S., Zhou, L., Hu, Y. and Yun, J., "Investigation of the hydrogen bubble effect on the overpotential in an alkaline water electrolyzer," Int. J. Hydrogen Energy, Vol.49, pp.47-57.

10.1016/j.ijhydene.2023.06.242

2023, Park, S. and Son, G., "Numerical study on hydrogen bubble production and overpotential in water electrolysis," KSCFE., Vol.28, pp.64-71.

10.6112/kscfe.2023.28.3.064

2022, Lee, J., Alam, A., Park, C., Yoon, S. and Ju, H., "Modeling of gas evolution processes in porous electrodes of zero-gap alkaline water electrolysis cells," Fuel, Vol.315, 123273.

10.1016/j.fuel.2022.123273

2024, Mohsin, H.M., Zhuo, Y. and Shen, Y., "Eulerian-Eulerian-VOF multifluid modelling of liquid-gas reacting flow for hydrogen generation in an alkaline water electrolyser," Fuel, Vol.373, 132164.

10.1016/j.fuel.2024.132164

OpenFOAM: User Guide v2112, Available online: https://www.openfoam.com/documentation/guides/latest/doc/ index.html.

2017, Haug, P., Kreitz, B., Koj, M. and Turek, T., "Process modelling of an alkaline water electrolyzer," Int. J. Hydrogen Energy, Vol.42, pp.15689-15707.

10.1016/j.ijhydene.2017.05.031

1998, Riegel, H., Mitrovic, J. and Stephan, K., "Role of mass transfer on hydrogen evolution in aqueous media," J. Appl. Electrochem., Vol.28, pp.10-17.

10.1023/A:1003285415420

2010, Guo Y.j., Xu, H.B., Fen, G.U.O., Zheng, S.L. and Zhang, Y., "Density and viscosity of aqueous solution of K2CrO4/KOH mixed electrolytes," TNMSC., Vol.20, pp.32-36.

10.1016/S1003-6326(10)60007-6

2017, MacA, E. "Modelling and simulation of an alkaline electrolyser cell," Energy, Vol.138, pp.316-331.

10.1016/j.energy.2017.07.053

2022, Colli, A.N. and Bisang, J.M., "Current and Potential Distribution in Two-Phase (Gas Evolving) Electrochemical Reactors by the Finite Volume Method," JES., Vol.169, 034524.

10.1149/1945-7111/ac5d90

2007, Gilliam, R.J., Graydon, J.W., Kirk, D.W. and Thorpe, S.J. "A review of specific conductivities of potassium hydroxide solutions for various concentrations and temperatures," Int. J. Hydrogen Energy, Vol.32, No.3, pp.359-364.

10.1016/j.ijhydene.2006.10.062

2024, Karimi-Sibaki, E., Varkhrushev, A., Wu, M., Bohacek, J. and Kharicha, A., "A numerical study on effects of current density distribution, turbulence, and magnetohydrodynamics (MHD) on electrolytic gas flow with application to alkaline water electrolysis (AWE)," Chem. Eng. Res. Des., Vol.208, pp.731-739.

10.1016/j.cherd.2024.07.042

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

Preview

NUMERICAL STUDY ON MASS TRANSFER AND HYDROGEN BUBBLE GENERATION IN WATER ELECTROLYSIS

ABSTRACT

MAIN

Fig. 1.

Schematics of an alkaline water electrolysis cell

(1)

(2)

(3)

(4)

(5)

(6)

(7)

(8)

(9)

(10)

(11)

(12)

(13)

(14)

(15)

(16)

Fig. 2.

Dependence of grid size on the prediction of αg profile at y=0.5m with i=1500A/m2

Table 1.

Physical properties used in AWE calculation

Fig. 3.

Effect of 𝜉 on the flow fields in AWE

Fig. 4.

Effect of 𝜉 on Yl and αg at the electrode-fluid interface in AWE

Fig. 5.

Effect of temperature on αg fields in AWE. The white solid lines represent αg=0.01 contours

Fig. 6.

Effect of temperature on αg profiles at y=0.5m in AWE

Fig. 7.

Effect of temperature on αH2 profiles at y=0.5m in AWE

Fig. 8.

Effect of current density on αg fields. The white solid lines represent αg=0.01 contours

Fig. 9.

Effect of current density on αg profiles at y=0.5m in AWE

(17)

(18)

(19)

Fig. 10.

Effect of current density on αg and streamwise velocity profile at y=0.5m

Fig. 11.

Effect of current density on overpotentials

Acknowledgements

References

Dependence of grid size on the prediction of $α_{g}$ profile at $y = 0.5 m$ with $i = 1500 A / m^{2}$

Effect of 𝜉 on $Y_{l}$ and $α_{g}$ at the electrode-fluid interface in AWE

Effect of temperature on $α_{g}$ fields in AWE. The white solid lines represent $α_{g}$ =0.01 contours

Effect of temperature on $α_{g}$ profiles at y=0.5m in AWE

Effect of temperature on $α_{H_{2}}$ profiles at y=0.5m in AWE

Effect of current density on $α_{g}$ fields. The white solid lines represent $α_{g}$ =0.01 contours

Effect of current density on $α_{g}$ profiles at y=0.5m in AWE

Effect of current density on $α_{g}$ and streamwise velocity profile at $y = 0.5 m$