모델링 방법에 따른 기체확산층 변형 비교 연구

1. 서 론

고분자 전해질 연료전지(Proton exchange membrane fuel cell, PEMFC)는 수소와 대기 중의 산소의 전기화학반응에 의하여 전기에너지를 발생시키는 친환경적 발전장치로써, 부산물로 물과 열을 배출한다.¹⁾ 고분자 전해질 연료전지의 성능은 주요 구성 요소들인 기체확산층(Gas diffusion layer, GDL), 분리판(Bipolar plate, BP), 막전극접합체(Membrane electrode assembly, MEA)의 상호 작용에 의하여 영향을 받는다. 기체확산층은 분리판을 통하여 공급된 수소와 산소를 전기화학 반응이 발생하는 막전극접합체로 전달하는 통로이며, 촉매층에서 수소와 산소의 전기화학반응으로 인해 발생하는 수분을 배출하는 역할을 한다. 이러한 주 기능 외에도 기체확산층은 전류와 열의 전도 역할을 수행하며, 기계적으로 약한 강도를 가지는 전해질막이 분리판으로부터 직접적인 손상을 입지 않도록 하는 완충 역할까지 수행한다.^2,3)

다공성 구조의 기체확산층에서는 볼트 체결시 분리판의 유동 채널과 리브의 구조로 인한 불균일한 응력 분포와 과도한 변형이 발생할 수 있고 그 결과로 인한 공극체적의 감소는 수분 및 물질의 이동을 방해 할 수 있다.^4,5) 따라서 연료전지 스택의 체결 시 체결압은 성능과 내구성에 영향을 미치는 중요한 변수로 작용하며 체결압에 따른 기체확산층의 변형거동 예측은 최적의 연료전지 성능을 구현하는데 필요한 요소이다.^6,7) 또한 기체확산층의 변형은 볼트 체결시 기체확산층 강성과 막전극접합체와의 접촉면에서의 상호 작용에 영향을 받는다.⁸⁾

이러한 기체확산층의 변형에 대한 해석적, 수치적 방법을 적용한 연구들은 많은 연구자들에 의하여 진행되어 왔다.^9-11) 기체확산층의 기계적/전기적 특성은 제조과정, 체결 조립 과정, 작동, 작동 후 손상 정도 등에 대하여 많은 영향을 미친다.⁷⁾ 이러한 특성 및 변형과정에 대한 문제는 주로 2차원 모델이나 실험적 방법에 의하여 수행되어 왔으며 유한 요소법과 같은 수치 방법이 적용되었다.^2,4) 해석모델링에 필요한 재료모델 및 재료 물성에 대한 연구도 진행되었다. Kleemann 등¹²⁾은 기체확산층의 변형에 따른 전기적 특성에 대한 연구에서 수치적 모델을 제시하고 실험결과와 비교를 하였다.

대다수의 연구에서는 연료전지 형상 특성을 고려한 평면 변형 모델과 균일 압력 하중을 작용 하중으로 사용하였다.^2,4,8) 그러나, 실제 기체확산층의 변형 특성은 불균일한 섬유 배치에 의하여 등방성 모델과는 차이가 있으며 체결에 의한 압축 하중을 받지만 형상 특성에 의한 불균일한 하중 상태를 나타내어 평면 변형 모델보다는 3차원 모델을 사용하는 경우가 실제 변형을 잘 구현할 수 있다.^7,13)

볼트 체결에 의한 압축하중하의 고분자 전해질 연료전지의 구성품 간의 기계적 거동, 전기적, 열적 특성, 질량 유동 현상에 대한 연구는 기존의 2차원 모델을 사용하여 수행하기에는 어려움이 따를 수 있다. 그러므로 연구자들은 모델의 일부가 아닌 3차원 연료전지모델을 활용하여 연구를 진행하였다.^7,13) 3차원 모델은 실제 변형을 대표하며 다양한 연구 목적에 적용할 수 있으나 효율적인 리소스 사용에 제약이 있다. 적절한 모델 크기의 제어가 가능한 3차원 서브모델의 사용은 이러한 제약을 제거하여 연료전지 구성품의 상세 해석을 가능하게 한다.

본 연구에서는 기체확산층의 변형거동을 모사하기 위하여 직교이방성(Orthotropic) 재료 거동모델을 가정하였고, 2차원 평면변형모델(2D plane strain model), 3차원 글로벌 모델(3D global model), 3차원 서브 모델(3D sub-model)을 사용하였으며 제시한 3 모델간의 변형 거동의 비교, 고찰을 통하여 최적의 연료전지 거동 모사에 적합한 모델을 제시하고자 한다.

직교이방성 재료모델은 상용유한요소해석 프로그램인 ABAQUS¹⁴⁾의 사용자재료서브루틴(UMAT)을 사용하여 구현하였고, 모델링 전처리 프로그램으로는 HYPERMESH, 유한요소해석은 ABAQUS를 사용하여 연구를 진행하였다.

2. PEMFC GDL 구조해석

2.1 유한요소모델

유한요소해석에 사용된 3차원 연료전지 모델은 분리판, 기체확산층, 막전극접합체로 구성된 단위 셀을 적층하여 양 끝판(End plate)에 볼트를 체결하여 구성하였다(Fig. 1).

Fig. 1

Assembled PEMFC global 3D model

효율적인 해석을 위하여 2차원 모델은 분리판, 기체확산층, 막전극접합체로 구성된 단위 셀을 평면 변형모델(2D plane strain model)로 이상화하였다(Fig. 2). 2차원 모델은 283,713개의 절점과 281,239개의 요소로 구성되어있으며, 모든 요소의 크기를 2μm로 모델링하였다. ABAQUS에서 제공하는 평면변형요소인 CPE4 요소를 사용하였고 Fig. 3과 같이 두께 방향을 y축으로 설정하여 모델링 하였으며 리브의 상단에 같은 크기의 압력 하중을 가하였다. 분리판, 기체확산층, 막전극접합체간의 접촉면에서는 접촉조건을 적용하였다. 반복적인 분리판 채널형상과 기체확산층의 적층구조 때문에 해석모델은 Fig. 3과 같은 대칭모델을 사용하였다.

Fig. 2

2D FEA model i llustrating width of channel(1 mm) with enlarged fillet area of bipolar plate

Fig. 3

2D FEA model with boundary conditions

연료전지 스택 모델은 두께가 얇은 구성품으로 이루어져 있기 때문에 해석 모델의 크기는 증가하게 된다.¹⁵⁾ 스택 전체를 모델링한 3차원 글로벌 모델은 2,648,322개의 절점과 2,840,316개의 요소로 구성되어있으며, 반복적인 분리판의 채널, 리브형상의 폭을 2차원 모델과 동일한 1 mm로 설정하였다(Fig. 4). 모델의 두께 방향은 z축으로 설정하였다(Fig. 5). 사용된 요소는 ABAQUS에서 제공하는 C3D6, C3D8이며 Fig. 5와 같이 양 끝판(End plate)에 경계조건과 실제로 적용될 볼트체결하중을 주어 해석을 수행하였다. 2차원 모델과 같이 분리판, 기체확산층, 막전극접합체간의 접촉면에서는 접촉조건을 적용하였다.

Fig. 4

Bipolar plate of 3D global model illustrating width of channel(1 mm)

Fig. 5

3D global FEA model with boundary conditions

3차원 글로벌 모델의 최 하단 단위 셀(분리판, 기체확산층, 막전극접합체로 구성)의 일부를 절단하여 구성한 3차원 서브 모델은 132,477개의 절점과 110,612개의 요소로 구성되어있다. ABAQUS에서 제공하는 C3D6, C3D8의 요소를 사용하였으며 Fig. 6과 같이 두께 방향을 z축으로 설정하고 서브 모델의 절단면에 글로벌 모델의 경계 조건을 부여하여 해석을 수행하였다.

Fig. 6

3D sub FEA model with boundary conditions

해석에 사용한 2차원 모델과 3차원 연료전지 스택 모델의 두께 방향 설정과 각 구성품의 두께 치수는 Table 1에 정리하였다. 2차원 평면 변형 모델에서의 두께 방향은 Y 방향으로 설정하였고 3차원 모델에서는 Z 방향이 두께 방향으로 설정되었다.

Table 1

Set-up of thickness direction and dimension in thickness direction of model

2.2 재료 모델

탄소섬유가 얽혀있는 다공성 구조의 기체확산층은 유연하여 체결압 작용 시 두께 방향으로 불균일한 공극 분포를 나타낼 수 있다. 이는 원활한 수분 배출을 어렵게 하여 연료전지 성능 저하의 원인이 될 수 있다.

기체확산층의 변형에 대한 기존 수치모사 연구들은 주로 2차원 모델을 사용하였으며 수치모사에 필요한 재료 모델을 해석적 방법이나 실험 결과값의 Curve-fitting 함수를 이용하였다.^16-18) 대부분의 연구에서 기체확산층 자체의 형상 및 내부 구조의 복잡성에 의하여 기체확산층의 재료 모델링은 간단한 등방성을 가정하여 연구가 진행되어 왔지만, 기체확산층은 탄소섬유의 배열로 인하여 거시적 관점에서 직교 이방성의 성질을 나타낸다고 할 수 있다.¹²⁾

일반적인 재료의 응력-변형률 관계는 아래와 같이 Hooke의 법칙(Hooke’s law)으로 나타낼 수 있다.

$\[ {\sigma }_{ij}=C_{ijkl}{ϵ}_{kl} \]$

(1)

여기서, 탄성텐서 행렬 C_ijkl는 36개의 재료상수가 필요한 반면에 이방성 재료는 21개의 재료상수의 결정이 필요하다. 직교이방성을 나타내는 기체확산층의 경우에는 9개의 재료상수의 결정이 필요하며, 단면 수직 방향으로 긴 형상을 지닌 연료전지 셀은 평면변형 모델로 이상화가 가능하다. 이상화된 2차원 직교이방성 평면변형 모델의 Hooke 법칙은 아래 식 (2)와 같으며 4개의 재료상수 (E_x,E_y,G_xy,ν_xy)가 필요하다.¹²⁾

$\[ \left\{\begin{array}{c}{\sigma }_x\\ {\sigma }_y\\ {\tau }_{xy}\end{array}\right\}=\left[\begin{array}{ccc}\frac{E_x}{1-{\nu }_{xy}{\nu }_{yx}}& \frac{{\nu }_{xy}{E}_y}{1-{\nu }_{xy}{\nu }_{yx}}& 0\\ \frac{{\nu }_{xy}{E}_y}{1-{\nu }_{xy}{\nu }_{yx}}& \frac{E_y}{1-{\nu }_{xy}{\nu }_{yx}}& 0\\ 0& 0& G_{xy}\end{array}\right]\left\{\begin{array}{c}{ϵ}_x\\ {ϵ}_y\\ {\gamma }_{xy}\end{array}\right\} \]$

(2)

기체확산층에 압축 하중이 작용하면 발생하는 부피의 감소는 기체확산층의 다공성 구조와 탄소섬유의 높은 강성으로 인하여 기공(Pore)에 의한 부피감소로 생각할 수 있어 포아송 비는 0으로 가정할 수 있다. 실험으로 측정한 포아송비 ν_xy도 0에 가까운 값을 나타내었다.¹²⁾ Fig. 3의 x, z, 방향으로의 재료 변형은 실제 다소간의 차이를 나타내며 단면상의 평면변형상태와 탄소섬유의 배열로 인한 높은 면내 평면 강성(In-plane stiffness)으로 인하여 x, z 방향의 변형은 리브(Rib)의 영향에 의한 국소 대변형을 제외하고 낮은 변형률을 나타낸다. 그러므로 x, z 방향에서의 변형은 서로 유사하며 선형탄성으로 가정할 수 있다.¹⁹⁾ 이러한 가정들을 식 (2)에 대입하면 기체확산층 모델에 적용할 수 있는 아래 식 (3)과 같은 2차원 직교이방성 모델 구성방정식을 구성할 수 있다.

$\[ \left\{\begin{array}{c}{\sigma }_x\\ {\sigma }_y\\ {\tau }_{xy}\end{array}\right\}=\left[\begin{array}{ccc}{E}_x& 0& 0\\ 0& E_y& 0\\ 0& 0& G_{xy}\end{array}\right]\left\{\begin{array}{c}{ϵ}_x\\ {ϵ}_y\\ {\gamma }_{xy}\end{array}\right\} \]$

(3)

연료전지를 체결하는 압축 과정에서 분리판 리브 부분과 접촉하는 기체확산층 부분은 처음 두께에 대하여 약 10-40 %의 변형을 나타낸다.¹⁹⁾ 이러한 두께 방향(2차원 모델에서 y 방향)의 대변형에 의한 비선형 거동은 두께 방향의 탄성계수 E_y를 변형률의 함수로 가정하여 나타낼 수 있다. 여기서는, 아래 식 (4)와 같이 두께 방향 변형률의 함수로 나타내었다.

$\[ \left\{\begin{array}{c}{\sigma }_x\\ {\sigma }_y\\ {\tau }_{xy}\end{array}\right\}=\left[\begin{array}{ccc}{E}_x& 0& 0\\ 0& E_y\left({ϵ}_y\right)& 0\\ 0& 0& G_{xy}\end{array}\right]\left\{\begin{array}{c}{ϵ}_x\\ {ϵ}_y\\ {\gamma }_{xy}\end{array}\right\} \]$

(4)

Table 2에는 수치모사에 사용한 기체확산층 두께방향의 비선형 탄성계수 함수 E(ϵ_y)를 정리하였다. 몇 가지 유형의 기체확산층에 대한 실험을 통하여 얻어진 실험값의 Curve-fitting 함수이고 구역은 3 부분으로 나누어져 있다. 소규모의 경화(Small strain hardening)를 나타내는 초기영역, 대규모의 경화(Large strain hardening)를 나타내는 마지막 영역, 그리고 2 영역을 이어주는 상수 영역(Constant modulus)으로 구성되어있다.¹⁹⁾

Table 2

Piecewise polynomial fitting model of the thickness direction non-linear behavior of GDL19)

일반적인 3차원 상태의 직교이방성 재료에 대한 Hooke의 법칙은 다음과 같다.

$\[ \left[\begin{array}{c}{\sigma }_{xx}\\ {\sigma }_{yy}\\ {\sigma }_{zz}\\ {\sigma }_{yz}\\ {\sigma }_{zx}\\ {\sigma }_{xy}\end{array}\right]=\left[\begin{array}{cccccc}\frac{1-v_{yz}{v}_{zy}}{E_y{E}_z\mathrm{\Delta }}& \frac{v_{yx}+v_{zx}{v}_{yz}}{E_y{E}_z\mathrm{\Delta }}& \frac{v_{zx}+v_{yx}{v}_{zy}}{E_y{E}_z\mathrm{\Delta }}& 0& 0& 0\\ \frac{v_{xy}+v_{xz}{v}_{zy}}{E_z{E}_x\mathrm{\Delta }}& \frac{1-v_{zx}{v}_{xz}}{E_z{E}_x\mathrm{\Delta }}& \frac{v_{zy}+v_{zx}{v}_{xy}}{E_z{E}_x\mathrm{\Delta }}& 0& 0& 0\\ \frac{v_{xz}+v_{xy}{v}_{yz}}{E_x{E}_y\mathrm{\Delta }}& \frac{v_{yz}+v_{xz}{v}_{yx}}{E_x{E}_y\mathrm{\Delta }}& \frac{1-v_{xy}{v}_{yx}}{E_x{E}_y\mathrm{\Delta }}& 0& 0& 0\\ 0& 0& 0& 2G_{yz}& 0& 0\\ 0& 0& 0& 0& 2G_{zx}& 0\\ 0& 0& 0& 0& 0& 2G_{xy}\end{array}\right]\left[\begin{array}{c}{\epsilon }_{xx}\\ {\epsilon }_{yy}\\ {\epsilon }_{zz}\\ {\epsilon }_{yz}\\ {\epsilon }_{zx}\\ {\epsilon }_{xy}\end{array}\right] \]$

(5)

여기서

$\[ \Delta =\frac{1-{\nu }_{xy}{\nu }_{yx}-{\nu }_{yz}{\nu }_{zy}-{\nu }_{zx}{\nu }_{xz}-2{\nu }_{xy}{\nu }_{yz}{\nu }_{zx}}{E_x{E}_y{E}_z} \]$

기체확산층을 3차원 직교이방성 재료 모델을 사용하여 해석하는 경우에는 재료 특유의 다공성 구조의 고려가 필요하고 이로 인하여 포아송비들을 0으로 가정할 수 있다. 그 결과로 식 (5)는 다음과 같이 식 (6)으로 간략화할 수 있다.⁷⁾

$\[ \left\{\begin{array}{c}{\sigma }_{xx}\\ {\sigma }_{yy}\\ {\sigma }_{zz}\\ {\sigma }_{yz}\\ {\sigma }_{zx}\\ {\sigma }_{xy}\end{array}\right\}=\left[\begin{array}{cccccc}{E}_x& 0& 0& 0& 0& 0\\ 0& E_y& 0& 0& 0& 0\\ 0& 0& E_z& 0& 0& 0\\ 0& 0& 0& G_{yz}& 0& 0\\ 0& 0& 0& 0& G_{zx}& 0\\ 0& 0& 0& 0& 0& G_{xy}\end{array}\right]\left\{\begin{array}{c}{ϵ}_{xx}\\ {ϵ}_{yy}\\ {ϵ}_{zz}\\ {\gamma }_{yz}\\ {\gamma }_{xz}\\ {\gamma }_{xy}\end{array}\right\} \]$

(6)

여기서, 2차원 모델과 마찬가지로 두께 방향(3차원 모델에서 z 방향)의 비선형 변형 거동은 두께 방향의 변형률의 함수(E_z(ϵ_z))로 나타내며 실험 결과값을 Curve-fitting 한 함수를 이용한다.

2차원 모델에서 유추한 것과 같이 x, y 방향에서의 재료 거동은 유사한 선형 탄성으로 가정하며 xy 면에서의 전단계수 G_xy는 잘 알려진 다음 관계식 (7)을 사용하여 계산할 수 있다.¹³⁾

$\[ {\left(G_{GDL}\right)}_{ xy}=\frac{{\left(E_{GDL}\right) }_x}{2\left(1+{\mathrm{v}}_{ GDL}\right)} \]$

(7)

2차원 평면 변형모델의 변형 거동 수치모사에는 식 (4)를 사용하였고 3차원 글로벌 모델과 서브 모델에는 간략화한 직교이방성 구성방정식(Simplified orthotropic constitutive equation)을 적용하여 체결압에 따른 각각의 거동을 비교하였다. 수치모사에 사용한 물성데이터는 Table 3에 정리하였다.

Table 3

Material properties used for simulation4,19)

3. 유한요소해석결과

상기한 재료 모델들을 상용유한요소해석프로그램인 ABAQUS에서 제공되는 사용자 재료서브루틴인 UMAT을 이용하여 구현하였다.

연료전지 단위 셀의 평면 변형만을 고려한 2차원모델의 경우에는 볼트 체결 하중을 압력 하중으로 이상화하였지만 3차원 모델의 경우에는 실제 볼트 체결 하중을 적용하였다.

2차원 모델과 3차원 모델의 비교 및 검토를 수행하기 위하여 2차원 모델의 적용하중인 압력하중을 동일한 하중 크기의 3차원 모델의 볼트 체결력으로 변환하는 과정이 필요하다. 압력하중을 같은 크기의 볼트 체결력으로 변환하여 적용하기 위하여 볼트 몸통부(Shank)의 단면적과 분리판과 기체확산층 사이의 접촉면적의 비가 필요하다. 계산된 비를 볼트 단면적을 기준으로 계산한 볼트 하중에 곱하여 3차원 모델의 볼트 체결력으로 사용하였으며 이 체결력을 등가 볼트하중이라 정의한다. 사용 모델에 대하여 계산된 접촉면적의 비는~5이고 이때 사용된 볼트의 직경은 10 mm이다.

2차원 압력하중에 대하여 동일한 압력 결과를 나타내는 변환된 등가볼트하중은 다음과 같은 식으로 나타낼 수 있다.

여기서 (F_bolt)_eq: 등가 볼트하중, P_pr: 2차원 압력하중, A_bolt: 볼트 단면적, C_conυert: 변환계수 (~5)

식 (8)을 적용하여 계산된 등가 볼트하중을 Table 4에 정리하였다.

Table 4

Applied 2D pressure vs. equivalent bolt load

수치모사에 사용된 유한요소 모델은 Fig. 7과 같다. 2D 평면변형모델은 가장 단순한 형태로 분리판, 기체확산층, MEA로 구성된 대칭 형상이다. 3차원 글로벌 모델은 실제 연료전지 스택 형상이며 서브 모델은 3차원 글로벌 모델 단위 셀에서 채널을 포함하여 일부분을 잘라낸 형상이다.

Fig. 7

The models used in the simulation: (a) 2D plane strain model with an inset showing mesh density, (b) 3D global, (c) 3D sub-model

기체확산층의 변형에 대한 연구는 2차원 모델을 사용한 경우가 대다수였고 실제 연료전지 스택 모델을 사용한 경우는 드물다. 실제 연료전지 스택의 체결은 볼트를 사용하므로 2차원 모델의 경우와는 다른 결과를 나타낼 수 있다. 특히 볼트를 스택 모서리에서 체결하므로 이상화한 2차원 모델에서의 경우처럼 균일한 접촉부 압력을 기대할 수 없다. 그러므로 접촉부에 균등한 체결압을 작용하는 2차원 모델과 동일한 하중 조건의 3차원 모델을 구현하면 모델 사이의 차이점과 어느 모델이 기체확산층 변형 연구에 유용한지의 판단이 가능할 것이다.

3차원 글로벌 모델에서의 기체확산층의 변형 형상이나 응력 분포에 대한 관찰은 조밀한 요소를 지닌 2차원 모델과 비교하여 요소망 조밀도에 의한 영향으로 인해 결과값의 큰 차이가 나타났다. 그러므로, 신뢰할 수 있는 결과를 도출하기 위하여 3차원 서브 모델을 구성하였다. 서브 모델은 스택 첫 번째 셀 중심부의 일부분을 잘라낸 형상을 사용하였고, 3차원 글로벌 모델의 변형 형상보다 2차원 평면 모델에 더 유사한 경향을 나타내었다.

Fig. 8에서는 2차원 모델에서의 압력하중과 동일한 크기의 볼트 하중을 적용한 3차원 서브 모델의 응력 분포를 나타내었다. 여기서, Cross direction 응력은 모델 단면 수평방향(x-dir)이고, Thickness direction은 모델 단면 수직 방향(z-dir)을 의미한다.

Fig. 8

Contour plot of stresses of the sub-model: (a) Von-mises stress, (b) Cross direction stress, (c) Thickness direction stress

Figs. 9~11에서는 2차원 평면 모델, 3차원 글로벌 모델, 3차원 서브 모델의 Von-mises 응력 분포를 나타내었다. 0.5 MPa에서 2.5 MPa까지 압축 하중이 증가함에 따라 기체확산층의 볼록한 정도는 심화되며 응력의 증가를 관찰할 수 있다. 2차원 평면 모델과 3차원 서브 모델의 채널 양 끝단에서 응력 집중이 관찰되었다(Fig. 9, Fig. 11). 반면에 3차원 글로벌 모델에서는 같은 위치에서 응력 집중 영역이 나타나지 않았다(Fig. 10).

Fig. 9

Contour plot of Von-mises stress for 2D plane strain model (1st: 0.5 (MPa), 2nd: 1 (MPa), 3rd: 1.5 (MPa), 4th: 2 (MPa), 5th: 2.5 (MPa))

Fig. 10

Contour plot of Von-mises stress for 3D global model (1st: 0.5 (MPa), 2nd: 1 (MPa), 3rd: 1.5 (MPa), 4th: 2 (MPa), 5th: 2.5 (MPa))

Fig. 11

Contour plot of Von-mises stress for 3D sub-model (1st: 0.5 (MPa), 2nd: 1 (MPa), 3rd: 1.5 (MPa), 4th: 2 (MPa), 5th: 2.5 (MPa))

기체확산층의 채널침입량은 기체확산층의 변형 정도를 나타내며 Fig. 12에 표시한 적색 원 사이의 수직 거리로 정의한다. Fig. 13에는 볼트 체결압의 증가에 따른 채널침입량의 변화량을 각 모델들과 문헌에서 인용한 실험결과 값에 대해서 도시하였다. 시뮬레이션에서 사용한 체결압력은 실험결과에서 사용한 것과 같은 크기를 사용하였다.

Fig. 12

The schematic of GDL intrusion illustrating vertical distance(distance between two red hollow circle)

Fig. 13

GDL intrusion vs compression pressure2,19)

실험 결과와 3 모델들이 모두 하중이 증가함에 따라서 채널 칩입량은 증가하였지만 요소망의 조밀도가 가장 낮은 3차원 글로벌 모델은 다른 모델들에 비하여 낮은 침입량을 나타내었다(Fig. 13). 조밀한 요소망을 사용한 모델은 실험 결과에 더 근접하는 결과를 나타내었으며 이는 서브 모델 구현시 현 모델보다 더 조밀한 요소망을 필요로 함을 나타낸다.

Fig. 14에는 동등한 하중 상태(2차원 모델: 2.5 (MPa) 압력하중, 3차원 모델: 981.75 (N) 볼트 하중)에서 각 모델들의 응력 분포 차이를 비교하고자 같은 최대/최소 응력 구간(Legend)를 설정하여 나타내었다. (a)에서는 Von-mises 응력의 분포, (b)에서는 단면 수평방향(Cross direction) 응력 분포, (c)에서는 단면 수직 방향(Thickness direction) 응력 분포를 나타내었다. 2차원 모델과 서브 모델의 경우 리브 Fillet의 유무의 차이에 의하여 다소 형상의 차이는 나타나지만 유사한 응력 분포와 분리판 접촉부 모서리에서의 응력 집중을 관찰할 수 있었다.

Fig. 14

Contour plot of stresses to observe differences among simulated 3 model cases with same level of stress: (a) Von-mises stress, (b) Cross direction stress, (c) Thickness direction stress

Fig. 15에서는 2차원, 3차원 글로벌, 서브 모델에 대해서 작용하중 증가에 따른 응력값의 변화를 도시하였다. 이때, 2차원 모델의 응력 결과는 경계조건에서의 영향을 최소화하기 위하여 리브의 모서리 부근과 대칭 경계조건 사이의 중간 위치에서의 값을 취하였다. 3차원 모델은 2차원 모델과 같은 최 하단 위치에서의 응력값을 취하였다. 3차원 글로벌 모델에 비하여 2차원 모델과 3차원 서브 모델은 유사한 경향을 나타내었지만 요소망의 조밀도 때문에 다소의 응력값 차이는 나타났다.

Fig. 15

Compression pressure vs. stress of GDL in the sub model: (a) Von-mises stress, (b) Cross direction stress, (c) Thickness direction stress

Fig. 16에서는 3차원 서브 모델의 채널길이 방향(y 방향)의 변형률 분포를 나타내었다. Fig. 15의 응력값을 취한 위치에서의 변형률 값은 거의 0에 근접하여 2차원 평면변형의 가정인 채널 길이 방향의 변형률이 0인 조건을 만족하였다. 그러므로 2차원 모델과 3차원 서브모델의 거시적인 변형 차이는 모델 요소망의 차이에 기인한 것으로 판단된다.

Fig. 16

Contour plot of elastic strain in through-plane direction

4. 결 론

고분자 전해질 연료전지의 주요 구성품인 기체확산층의 변형거동을 적절히 구현할 수 있는 재료모델에 대한 비교 연구를 수행하였다. 2차원, 3차원 모델을 구성하였고 재료의 거동은 직교이방성을 가정하였다. 특히, 서브 모델링기법을 사용하여 기체 확산층의 3차원 변형 양상에 대한 연구를 시도하였다. 2차원 모델에서 가정한 압력하중과 동등한 볼트하중을 3차원 모델에 적용하였으며 동일한 위치에서 응력값을 취하여 모델의 적합성을 확인하였다. 이러한 모델간의 비교를 통하여 본 연구에서는 다음과 같은 결과를 도출할 수 있었다.

• 1) 연료전지 압축 체결 과정에 대한 수치 모사 방법을 2차원 평면 변형모델, 3차원 글로벌 모델, 3차원 서브 모델에서 구현하여 모델의 적합성을 확인하였다.
• 2) 2차원 평면 변형모델의 적용 압력 하중과 3차원 글로벌 모델, 3차원 서브 모델에서 사용하는 볼트 하중의 선형관계모델을 제안하였다.
• 3) 조밀한 요소망을 적용한 2차원 평면 변형모델을 기준으로 결과를 비교하였으며 요소망의 차이 때문에 3차원 글로벌 모델보다는 3차원 서브 모델에서의 응력 결과값들이 유효한 경향을 나타내었다.
• 4) 거시적 관점에서는 3차원 모델의 채널 길이 방향(Through-plane direction)의 변형률은 대부분의 위치에서 0에 가까운 값을 나타내어 2차원 평면 변형 모델의 유효함을 나타내었으나 응력이 집중하는 리브 모서리 부분에서는 1-3.7 %의 변형률을 나타내어 기체확산층 변형에 대한 3차원 효과가 있을 것으로 판단된다.
• 5) 3차원 글로벌모델과 적절한 모델 크기를 지닌 서브 모델의 적용은 연료전지 스택의 기체확산층 변형을 잘 구현할 수 있을 것으로 판단되고 다양한 후속 연구에 응용될 수 있을 것이다.

Acknowledgments

이 논문은 안동대학교 기본연구지원사업에 의하여 연구되었음.

References

• H. Jeong, J. Kim, S. Lee, C. Lim, B. Ahn and C. Kim, “Analysis of Mass Transport in PEMFC GDL,” Transactions of the Korean Society of Mechanical Engineers-B, Vol.36, No.10, pp.979-988, 2012. [https://doi.org/10.3795/KSME-B.2012.36.10.979]
• T. H. Lee and C. H. Yang, “A Study on the Design Parameters Affecting GDL Deformation Behavior,” Transactions of the KSAE, Vol.26, No.6, pp.764-772, 2018. [https://doi.org/10.7467/KSAE.2018.26.6.764]
• C. J. Oh and Y. T. Lee, “Measurement of In-plane Gas Permeability of Gas Diffusion Layers in Proton Exchange Membrane Fuel Cells Under Compressive Strain,” Korean Journal of Air-Conditioning and Refrigeration Engineering, Vol.28, No.9, pp.367-372, 2016. [https://doi.org/10.6110/KJACR.2016.28.9.367]
• Y. Lai, P. A. Rapaport, C. Ji and V. Kumar, “Channel Intrusion of Gas Diffusion Media and the Effect on Fuel Cell Performance,” Journal of Power Sources, Vol.184, pp.120-128, 2008. [https://doi.org/10.1016/j.jpowsour.2007.12.065]
• P. Chippar, O. Kyeongmin, K. Kang and H. Ju, “A Numerical Investigation of the Effects GDL Compression and Intrusion in Electrolyte Fuel Cells(PEFCs),” International Journal of Hydrogen Energy, Vol.37, No.7, pp.6326-6338, 2012. [https://doi.org/10.1016/j.ijhydene.2011.04.154]
• J. S. Kim and J. B. Kim, “Effect of Gas Diffusion Layer Compression and Inlet Relative Humidity on PEMFC Performance,” Applied Chemistry for Engineering, Vol.32, No.1, pp.68-74, 2021.
• Y. Chen, C. Jiang and C. D. Cho, “An Investigation of a Three-Dimensional Constitutive Model of Gas Diffusion Layers in Polymer Electrolyte Membrane Fuel Cells,” World Academy of Science Engineering and Technology, International Journal of Energy and Power Engineering, Vol.12, No.12, pp.906-912, 2018.
• T. H. Lee and C. H. Yang, “A Parametric Study on the Deformation of Gas Diffusion Layer in PEM Fuel Cell,” Journal of Mechanical Science and Technology, Vol.34, pp.259-268, 2020. [https://doi.org/10.1007/s12206-019-1227-8]
• J. Yoon and J. Park, “Structural Analysis of Gasket and GDL for Enhanced Performance of PEMFC,” Journal of the Korean Society for Aeronautical and Space Science, Vol.36, No.7, pp.642-650, 2008. [https://doi.org/10.5139/JKSAS.2008.36.7.642]
• P. Chippar, K. O. K. Kang and H. Ju, “A Numerical Investigation of the Effects of GDL Compression and Intrusion in Polymer Electrolyte Fuel Cells (PEMFCs),” International Journal of Hydrogen Energy, Vol.37, No.7, pp.6326-6338, 2012. [https://doi.org/10.1016/j.ijhydene.2011.04.154]
• P. A. Gigos, Y. Faydi and Y. Mayer, “Mechanical Characterization and Analytical Modeling of Gas Diffusion Layer under Cyclic Compression,” International Journal of Hydrogen Energy, Vol.40, No.17, pp.5958-5965, 2015. [https://doi.org/10.1016/j.ijhydene.2015.02.136]
• J. Kleemann, F. Finsterwalder and W. Tillmetz, “Characterization of Mechanical Behaviour and Coupled Electrical Properties of Polymer Electrolyte Membrane Fuel Cell Gas Diffusion Layers,” Journal of Power Sources, Vol.190, pp.92-102, 2009. [https://doi.org/10.1016/j.jpowsour.2008.09.026]
• B. Liu, L. F. Liu, M. Y. Wei and C. W. Wu, “Vibration Mode Analysis of the Proton Exchange Membrane Fuel Cell Stack,” Journal of Power Sources, Vol.331, No.1, pp.229-307, 2016. [https://doi.org/10.1016/j.jpowsour.2016.09.056]
• ABAQUS User’s Manual, Ver.6.4, HKS, Providence, RI, 2003.
• D. Kim, J. Kim and H Kim, “Assembly Analysis for Evaluation of Sealing in PEMFC Stack,” Transactions of the KSAE, Vol.18, No.5, pp.68-75, 2010.
• H. S. Chu, C. Yeh and F. Chen, “Effects of Porosity Change of Gas Diffuser on Performance of Proton Exchange Membrane Fuel Cell,” Journal of Power Sources, Vol.123, No.1, pp.1-9, 2003. [https://doi.org/10.1016/S0378-7753(02)00605-5]
• Z. Y. Su, C. T. Liu, H. P. Chang, C. H. Li, K. J. Huang and P. C. Sui, “A Numerical Investigation of the Effects of Compression Force PEM Fuel Cell Performance,” Journal of Power Sources, Vol.183, No.1, pp.182-192, 2008. [https://doi.org/10.1016/j.jpowsour.2008.04.060]
• V. Norouzifard and M. Bahrami, “Deformation of PEM Fuel Cell Gas Diffusion Layers under Compressive Loading: An Analytical Approach,” Journal of Power Sources, Vol.264, pp.92-99, 2014. [https://doi.org/10.1016/j.jpowsour.2014.04.057]
• P. A. Garcia-Salaberri, M. Vera and R. Zaera, “Nonlinear Orthotropic Model of the Inhomogeneous Assembly Compression of PEM Fuel Cell Gas Diffusion Layer,” International Journal of Hydrogen Energy, Vol.36, No.18, pp.11856-11870, 2011. [https://doi.org/10.1016/j.ijhydene.2011.05.152]