Parametric Studies on Hydrogen Embrittlement in Liquified Hydrogen Tank using Molecular Dynamics Simulation

Song-Hyun Cha; Hyun-Seok Kim; Seonho Cho

doi:10.7734/COSEIK.2022.35.6.325

Preview

Research Paper

Journal of the Computational Structural Engineering Institute of Korea. 31 December 2022. 325-331
https://doi.org/10.7734/COSEIK.2022.35.6.325

Parametric Studies on Hydrogen Embrittlement in Liquified Hydrogen Tank using Molecular Dynamics Simulation

분자동역학을 이용한 액화수소 연료탱크의 수소취성화 파라메터 연구

Song-Hyun Cha¹

Hyun-Seok Kim²^†

Seonho Cho³

차 송현¹

김 현석²^†

조 선호³

¹Graduate Student, Department of Naval Architecture and Ocean Engineering, Seoul National University, Seoul, 08826, Korea

²Research Engineer, Alternative Fuels and Power System Research Division, KRISO, Daejeon, 34103, Korea

³Professor, Department of Naval Architecture and Ocean Engineering, Seoul National University, Seoul, 08826, Korea

¹서울대학교 조선해양공학과 박사과정

²선박해양플랜트연구소 친환경연료추진연구본부 선임연구원

³서울대학교 조선해양공학과 교수

^{†Corresponding Author}

ABSTRACT

Hydrogen embrittlement in metals has been a serious issue with regard to structural safety. In this study, molecular dynamics simulations revealed that the aggregation of hydrogen atoms at the crack tips suppresses the dislocation emission and thus results in cleavage fracture. A series of molecular dynamics simulations were performed considering factors such as the concentration of hydrogen atoms, loading rate, and diffusion coefficient. We investigated the conditions that minimize hydrogen embrittlement. The simulation results were consistent with the experimental results and used to quantify hydrogen embrittlement.

Keywords

molecular dynamics

hydrogen embrittlement

dislocation emission

cleavage fracture

금속의 취성화는 수소와 접촉하는 구조물을 안정적으로 설계하는데 있어서 큰 문제가 되어왔다. 본 논문에서는 분자동역학 해석을 통해 균열선단 주변에 모인 수소원자들이 전위 이동 현상을 억제하고, 이로 인해 벽개 파괴 현상이 발생하는 것을 확인하였다. 다양한 수소 농도, 하중 속도, 수소 확산 속도 등을 바꾸어가며 분자동역학 해석을 수행하였고, 이에 따른 수소 취성화를 최소화시킬 수 있는 조건들을 조사하였다. 분자동역학 해석 결과는 기존의 실험결과와 잘 일치하였으며 이를 바탕으로 수소 취성화 현상을 정량화하여 평가하였다.

키워드

분자동역학

수소 취성화

전위 이동

벽개 파괴

MAIN

1. 서 론
2. 본 론
2.1 분자동역학 시뮬레이션
2.2 NVT 앙상블
2.3 EAM 퍼텐셜
2.4 분자동역학 모델
2.5 시뮬레이션 조건
2.6 파라미터
2.7 시뮬레이션 결과
3. 결 론

1. 서 론

수소에 의한 금속의 취성화는 수소와 관련된 산업 현장에서 항상 중요한 문제가 되어왔다. 특히 최근 환경문제로 수소 에너지의 수요가 증가하여 액화수소 탱크에서 그 중요성이 점점 커지고 있다. 기존에 알려진 바에 따르면(Zhang, 2016), 수소가 금속 표면에서 영향을 주는 방식은 크게 두 가지가 있는데, ‘기포 발생’, ‘수소 취성화’가 있다. 기포 발생은 용기 내의 빈 공간에 수소 기체가 누적되어 기포가 커지는 현상을 의미한다. 수소 취성화는 Fig. 1에서와 같이 수소 입자가 철 입자 사이사이에 침투하여 취성을 갖는 물질을 형성하고, 여기에 하중이 가해지면서 균열이 발생하게 된다.

https://cdn.apub.kr/journalsite/sites/jcoseik/2022-035-06/N0040350601/images/Figure_jcoseik_35_06_01_F1.jpg

Fig. 1.

Hydrogen embrittlement(Zhang, 2016)

이미 존재하는 마이크로 균열 주변에 하중이 가해졌을 때, 수소가 존재하는 경우와 수소가 존재하지 않는 경우를 서로 비교하여 수소 취성화를 설명하는 메카니즘(Song and Curtin, 2011)이 개발되었다. 수소가 존재하지 않는 마이크로 균열은 Fig. 2(a)와 같이 격자구조가 어긋난 전위(dislocation)를 방출하거나 흡수하면서 균열의 끝부분이 뭉툭하게 변하게 된다, 이러한 마이크로 균열은 진전이 되지 않고 따라서 액화수소 탱크의 파괴에 큰 영향을 주지 않는다. 그러나 마이크로 균열 주변에 수소가 충분히 존재할 경우, Fig. 2(b)와 같이 균열 끝부분의 응력 분포로 인해서 수소 원자가 모이게 된다. 모인 수소 원자들은 전위의 흡수와 방출을 막고, 수화물과 같은 수소와 철이 결합한 물질을 형성하게 되어 여기에서 벽개(cleavage)가 발생한다.

https://cdn.apub.kr/journalsite/sites/jcoseik/2022-035-06/N0040350601/images/Figure_jcoseik_35_06_01_F2.jpg

Fig. 2.

Effect of hydrogen in crack propagation(Song and Curtin, 2011)

본 논문에서는 이러한 현상을 원자 단위로 분석하기 위해서 분자동역학 해석을 도입하였으며(Matsumoto et al., 2009; Taketomi et al., 2010), 여기에 사용된 수치 모델은 Fig. 3과 같다.

https://cdn.apub.kr/journalsite/sites/jcoseik/2022-035-06/N0040350601/images/Figure_jcoseik_35_06_01_F3.jpg

Fig. 3.

Atomistic configuration of simulation model

분자동역학 해석은 LAMMPS 코드를 활용하여 수행하였으며, 확산을 통해 균열 선단에 도달하는 수소 원자를 모사하기 위하여 수소 원자는 균열 선단 주위의 반경 10nm 원통형 영역 내 사면체의 틈새 자리에 피코초(ps)당 1개의 비율로 무작위로 배치하였다. 100개의 수소원자가 추가될 때마다 수소가 분리되어 안정된 위치를 갖도록 4나노초(ns) 동안 분자동역학 시스템을 완화시켜 주었다.

2. 본 론

이 장에서는 분자동역학을 활용한 수소 취성화 모델을 소개한 후, 분자동역학 모사 과정에서 얻은 결과와 실험 결과를 비교하고, 일치함을 확인하고자 한다.

2.1 분자동역학 시뮬레이션

역학 문제를 해결할 때, 연속체 스케일의 현상에 대해서는 유한요소법이 주로 사용되었지만, 나노 스케일의 현상에 대해서는 원자 수준의 물체를 해석하는 데에 있어서 유한요소법은 충분하지 않다. 나노 스케일에서 물체의 특성과 역학적 거동을 구하기 위해서 분자동역학 전산모사는 가장 흔히 사용되는 방법이다(Cha et al., 2014). 원자구조는 각 원자 사이의 상호작용에 의해 지배를 받으므로 원자 사이의 퍼텐셜이 정의될 필요가 있다. 원자들 사이의 모든 상호작용을 고려하게 되면, 퍼텐셜에 대한 계산비용이 커지게 되므로 이를 극복하기 위해 주기 경계 조건(periodic boundary condition)을 도입하였다.

2.2 NVT 앙상블

입자의 개수가 N, 부피가 V, 온도가 T로 고정되어 있는 시스템에서 분자동역학 시뮬레이션을 수행하기 위해 물리적 시스템이 열 저장조에 접촉해 있다고 가정한다. 이 확장된 시스템의 해밀토니안을 구하면 식 (1)과 같다.

(1)

H (\tilde{p}, q, {\tilde{p}}_{s}, s) = \sum_{i = 1}^{3 N} \frac{{\tilde{p}}_{i}^{2}}{2 m_{i} s^{2}} + \frac{{\tilde{p}}_{s}^{2}}{2 Q} + V_{p o t} (q) + z l n s

여기서, $q, \tilde{p}, V_{p o t}$ 은 각각 원자의 위치, 운동량, 퍼텐셜 에너지를 의미한다. $s, {\tilde{p}}_{s}, Q$ 는 각각 열원(heat bath)의 자유도를 나타내는 파라메터, 운동량, 질량을 의미한다. $z$ 는 NVT 앙상블의 샘플링 인수이다. $p$ 위에 표시된 틸다는 다른 척도의 시간 𝜏에 대한 변수를 의미한다. 위의 해밀토니안을 변환하여 운동방정식을 구하면 식 (2), (3), (4)와 같다.

(2)

{\ddot{q}}_{i} = \frac{F_{i}}{m_{i}} - {\dot{q}}_{i} ζ

(3)

\dot{ζ} = \frac{1}{Q} (\sum_{i = 1}^{3 N} m_{i} {\dot{q}}_{i}^{2} - z)

(4)

z = 3 N k_{B} T

여기서, $F_{i}$ 는 원자 $i$ 가 받는 힘, $m_{i}$ 는 원자 $i$ 의 질량, 𝜁는 $\dot{s} / s$ 를 나타내는 값으로, 마찰 계수로 부른다. $N$ 은 원자의 개수, $k_{B}$ 는 볼츠만 상수, $T$ 는 절대온도를 의미한다.

2.3 EAM 퍼텐셜

EAM(Embedded Atom Method) 퍼텐셜은 금속 시스템에서 결합을 나타내기 위해 국부적인 전자 밀도를 고려하고 있기 때문에 가장 적절한 다체(many body) 퍼텐셜로 알려져 있다. EAM 퍼텐셜을 식으로 나타내면 식 (5)와 같다.

(5)

E_{i} = F_{i} ({\bar{ρ}}_{i}) + \frac{1}{2} \sum_{j \neq i} ϕ_{i j} (r_{i j})

여기서, $F_{i}$ 는 전자밀도의 함수로 주어지는 삽입(embedding) 에너지를 의미하고, 전자 밀도는 식 (6)과 같이 주어진다.

(6)

{\bar{ρ}}_{i} = \sum_{j \neq i} ρ_{j} (r_{i j})

𝜙는 쌍(pair) 퍼텐셜에 의한 상호작용을 의미하고 식 (7)과 같다.

(7)

ϕ_{i j} (r_{i j}) = z_{i} (r_{i j}) z_{j} (r_{i j}) / r_{i j}

$z_{i}$ 는 원자 $i$ 의 유효 전하 함수이다.

2.4 분자동역학 모델

분자동역학 모델은 Fig. 3과 같으며 시뮬레이션 상자의 크기는 485Å×480Å×21Å이다. 원기둥의 반지름은 43.1Å이며 원기둥의 높이는 21Å로서 시뮬레이션 상자와 높이가 동일하여 주기적 경계조건으로 z축 방향의 격자구조가 이어져 있다. 균열의 길이는 반지름보다 약간 짧은 37Å이며 격자 구조의 방향성(orientation)은 $(111) [1 \bar{1} 0]$ 이다.

원자간 퍼텐셜(interatomic potential)은 기존에 개발된 EAM/ fs 퍼텐셜(Ramasubramaniam et al., 2009)을 사용하였다. 그러나 이를 그대로 사용하면 서로 이웃한 수소 입자 간의 상호작용이 너무 강해서 실제보다 너무 잘 뭉치게 된다. 이를 보정하기 위해 수소 입자와 수소 입자 사이의 척력에 해당하는 항을 더해 주었으며 식 (8)과 같다(Song and Curtin, 2013).

(8)

f (r_{i j}, λ, k, C_{0}, B_{0}, r_{0}, r_{1}) = \{\begin{matrix} C_{0} {[\frac{(r_{i j} - r_{0})}{λ}]}^{k - 1} \exp [- {(\frac{r_{i j} - r_{0}}{λ})}^{k} - B_{0} {(r_{i j} - r_{1})}^{2}] if r_{i j} \geq r_{1} \\ C_{0} {[\frac{(r_{i j} - r_{0})}{λ}]}^{k - 1} \exp [- {(\frac{r_{i j} - r_{0}}{λ})}^{k}] if r_{0} \leq r_{i j} \leq r_{1} \\ 0 if r_{i j} < r_{0} \end{matrix}

$r_{i j}$ 는 원자 $i$ 와 원자 $j$ 사이의 거리를 의미하고, 나머지 변수인 $λ, k, C_{0}, B_{0}, r_{0}, r_{1}$ 은 퍼텐셜을 나타내기 위한 상수이다.

2.5 시뮬레이션 조건

분자동역학 시뮬레이션은 NVT 앙상블로 300K의 일정한 온도에서 수행되었고, 시간 간격은 0.4펨토초(fs)로 설정하였다. 하중은 원기둥 모델의 바깥 경계에 있는 각 원자들에 식 (9)와 같이 변위를 가하였다. 균열의 끝부분을 중심으로 입자가 위치한 r과 θ가 정해지면, 각 입자별로 변위가 결정된다.

(9)

\{\begin{matrix} u \\ v \end{matrix}\} = \frac{K_{I}}{G} \sqrt{\frac{r}{2 π}} \{\begin{matrix} \cos \frac{θ}{2} \\ \sin \frac{θ}{2} \end{matrix}\} (β - \cos^{2} \frac{θ}{2}) = \frac{K_{I}}{G} \sqrt{\frac{r}{2 π}} \frac{1}{4} \{\begin{matrix} (4 β - 3) \cos \frac{θ}{2} - \cos \frac{3 θ}{2} \\ (4 β - 1) \sin \frac{θ}{2} - \sin \frac{3 θ}{2} \end{matrix}\}

$K_{I}$ 은 응력확대계수, $G$ 는 전단 탄성계수, 𝛽는 평면변형률(plane strain)라고 가정할 경우, $2 - 2 ν$ 값을 갖는 상수이다(𝜈는 포아송 비).

시뮬레이션 초기에 응력확대계수 $K_{I}$ 이 0.8이 될 때까지 변형시킨다. 그리고 0.4 나노초(ns) 동안 $K_{I}$ 을 그 상태로 유지한 채로, 완화(relaxation) 시켜준다. 그 후, 4 피코초(ps)마다 $K_{I}$ 을 0.016만큼씩 증가시킨다.

2.6 파라미터

분자동역학 시뮬레이션을 통해 확인할 수 있는 파라미터 값들 중에 CNA(Common Neighbor Analysis) 파라미터는 특정 원자와 주변 원자들의 위치 관계를 알고리즘을 통해 각각의 격자구조를 수치로 나타내는 파라미터이다(FCC=1, HCP=2, BCC=3, Icosahedral=4, Unknown=5).

비리얼 응력(virial stress)은 원자 스케일에서 각 원자의 위치와 받는 힘으로 각 입자가 받는 응력을 계산할 수 있는 파라미터이며 이를 통해 응력분포를 확인할 수 있고 식 (10)과 같이 나타낼 수 있다.

(10)

σ_{v i r i a l} = \frac{1}{Ω} (- m_{i} {\dot{u}}_{i} \otimes {\dot{u}}_{i} + \frac{1}{2} \sum_{j \neq i} r_{i j} \otimes f_{i j})

여기서, 𝛺는 원자 i의 원자 부피를 의미하고 $r_{i j} = r_{i} - r_{j}$ 은 원자 i와 원자 j 사이의 변위를 의미한다. $f_{i j}$ 은 원자 i와 원자 j 사이에 걸리는 힘을 의미한다.

2.7 시뮬레이션 결과

단위 높이당 수소 원자의 개수를 $N_{H} / L_{z}$ 라 할 때, $N_{H} / L_{z}$ 가 0인 경우, 즉, 수소가 존재하지 않을 때 미소 균열(microcrack)은 Fig. 4와 같이 전위를 방출하거나 흡수하면서 균열의 끝부분이 뭉툭하게 변한다. 이러한 형태의 미소 균열은 진전되지 않고, 구조물의 파괴에 영향을 주지 않는다.

https://cdn.apub.kr/journalsite/sites/jcoseik/2022-035-06/N0040350601/images/Figure_jcoseik_35_06_01_F4.jpg

Fig. 4.

CNA parameter at $N_{H} / L_{z} = 0 / Å$

$N_{H} / L_{z}$ 가 100/Å인 경우(Fig. 5), 초기하중을 가하기 위해 균열을 약간 진전시키면, 균열 선단 끝부분 주변으로 수소 원자가 모여들게 된다. 모여진 수소 원자들이 철 원자 사이의 결합을 약하게 만들고, 그로 인해 벽개가 발생한다.

https://cdn.apub.kr/journalsite/sites/jcoseik/2022-035-06/N0040350601/images/Figure_jcoseik_35_06_01_F5.jpg

Fig. 5.

CNA parameter at $N_{H} / L_{z} = 100 / Å$

$N_{H} / L_{z}$ 를 200/Å으로 늘린 경우(Fig. 6), 균열 끝부분의 아래쪽으로 수소 원자들이 모이는 것을 확인할 수 있다. 수소 원자가 균열 아래쪽에 모일 때 응력 분포도 마찬가지로 점점 균열 아래쪽의 입자들이 큰 값을 갖는 것을 확인할 수 있다. 이러한 응력 집중에 의해 수소 원자가 모이고, 이는 그 주변의 철의 응력을 증가시키고, 다시 증가된 응력이 수소 원자를 모이게 한다. 즉, 이 과정을 반복하면, 수소와 철이 결합한 화합물이 균열의 끝부분을 시작으로 점점 성장해 나간다.

https://cdn.apub.kr/journalsite/sites/jcoseik/2022-035-06/N0040350601/images/Figure_jcoseik_35_06_01_F6.jpg

Fig. 6.

CNA parameter at $N_{H} / L_{z} = 200 / Å$

$N_{H} / L_{z}$ 를 300/Å으로 늘린 경우(Fig. 7), 균열 끝 부분의 위와 아래 모두 수소 원자들이 모이는 것을 알 수 있다. 따라서 전위가 방출되는 것을 모두 막을 수 있고, 수소가 모인 영역을 따라서 균열이 진전한다. 전위가 방출되는 것을 모두 막을 수 있는 수준으로 수소가 모이게 되면, 수소 농도가 증가하더라도 영향을 받지 않게 된다는 것을 알 수 있다.

https://cdn.apub.kr/journalsite/sites/jcoseik/2022-035-06/N0040350601/images/Figure_jcoseik_35_06_01_F7.jpg

Fig. 7.

CNA parameter at $N_{H} / L_{z} = 300 / Å$

$N_{H} / L_{z}$ 를 300/Å으로 유지한 채로 온도를 200~400K으로 변화시킴에 따라 시뮬레이션 상에서 어떤 변화가 있는지 확인해 보았다. 온도가 200K일 경우(Fig. 8), 비교적 낮은 온도이기 때문에 수소의 확산 속도가 느려서 균열 끝부분을 향해 수소가 덜 모이는 것을 확인할 수 있다. 시뮬레이션 최종단계에서 CNA 값을 보면 여전히 균열에서 떨어진 영역에 수소가 남아있는 것을 확인할 수 있다.

https://cdn.apub.kr/journalsite/sites/jcoseik/2022-035-06/N0040350601/images/Figure_jcoseik_35_06_01_F8.jpg

Fig. 8.

CNA parameter at $T = 200 K$

온도를 300K로 증가시킨 경우(Fig. 9), 200K일 때에 비해 상대적으로 확산 속도가 빨라서 수소가 균열의 오른쪽에 많이 모이는 것을 확인할 수 있다. 따라서 균열이 진전될 때, 수소가 모여있는 경로를 따라서 벽개가 발생할 것으로 예상할 수 있다.

https://cdn.apub.kr/journalsite/sites/jcoseik/2022-035-06/N0040350601/images/Figure_jcoseik_35_06_01_F9.jpg

Fig. 9.

CNA parameter at $T = 300 K$

온도가 400K일 경우(Fig. 10), 확산속도가 300K일 때보다 더 빨라서 균열 주변에 모여있는 수소 원자가 부족한 것을 확인할 수 있다. 따라서 균열이 진전될 때, 수소가 부족해서 전위를 완전히 막아주지 못해 방출되는 것을 확인할 수 있다.

https://cdn.apub.kr/journalsite/sites/jcoseik/2022-035-06/N0040350601/images/Figure_jcoseik_35_06_01_F10.jpg

Fig. 10.

CNA parameter at $T = 400 K$

하중 속도(loading rate)의 변경에 따른 시뮬레이션 상의 변화를 살펴보았다. 하중 속도 $\dot{K_{I}}$ 가 기존의 50배인 $2 \times 10^{- 13} M P a \sqrt{m} / s$ 일 경우(Fig. 11), 수소가 균열을 향해 충분히 모여들기 때문에 수소가 모인 영역을 따라서 벽개가 발생할 것으로 예상할 수 있다.

https://cdn.apub.kr/journalsite/sites/jcoseik/2022-035-06/N0040350601/images/Figure_jcoseik_35_06_01_F11.jpg

Fig. 11.

CNA parameter at $\dot{K_{I}} = 2 \times 10^{- 13} M P a m^{1 / 2} s^{- 1}$

하중 속도 $\dot{K_{I}}$ 가 기존의 500배인 $2 \times 10^{- 12} M P a \sqrt{m} / s$ 일 경우(Fig. 12), 수소 원자가 균열 끝 부분으로 모이기도 전에 균열이 벌어지는 것을 확인할 수 있다.

https://cdn.apub.kr/journalsite/sites/jcoseik/2022-035-06/N0040350601/images/Figure_jcoseik_35_06_01_F12.jpg

Fig. 12.

CNA parameter at $\dot{K_{I}} = 2 \times 10^{- 12} M P a m^{1 / 2} s^{- 1}$

기존의 열역학적 관점으로 식을 정리하면, 화학적 퍼텐셜에 의해 추진력이 발생한다고 보고, 화학적 퍼텐셜을 식 (11), 압력장을 식 (12)와 같이 정리하였다(Song and Curtin, 2013). $T$ 는 온도, $k_{B}$ 는 볼츠만 상수, $c$ 는 농도, 𝜇는 화학적 퍼텐셜, $p (r, θ)$ 는 균열 끝부분 주변의 압력장, 𝛺는 철 사이에 투입된 수소의 부분 부피(partial volume)를 의미한다.

(11)

μ = μ_{0} + k_{B} T \ln \frac{c}{1 - c} + p Ω

(12)

p (r, θ) = \frac{2 (1 + ν)}{3 \sqrt{2 π r}} {\dot{K}}_{I} t \cos \frac{θ}{2}

균열 끝으로 향하는 수소 입자의 평균적인 방사상 속도(radial velocity)를 구하면 식 (13)과 같다.

(13)

V_{r} = - \frac{d r}{d t} = - \frac{D}{k_{B} T} (\nabla μ)_{r} = \frac{(1 + ν) D Ω {\dot{K}}_{I} t}{3 \sqrt{2 π}} \frac{\cos (θ / 2)}{r^{3 / 2}}

그리고 반지름 R인 영역 이내에 모든 수소 입자가 균열 끝 $r$ ~0에 도달한다고 하면, 식 (13)을 적분하여 R을 식 (14)와 같이 나타낼 수 있다.

(14)

R = {(\frac{5 (1 + ν) D Ω {\dot{K}}_{I}}{12 \sqrt{2 π} k_{B} T} t^{2} \cos (θ / 2))}^{2 / 5}

BCC 구조에서 단위 높이당 수소원자의 개수를 시간 t의 시점에서 구하면,

(15)

\frac{N_{H}}{L_{z}} = \frac{2 c_{0}}{a_{0}^{3}} \int_{- π}^{π} \frac{R^{2}}{2} d θ = β (\frac{1}{2}, \frac{9}{10}) \frac{2 c_{0}}{a_{0}^{3}} {(\frac{5 (1 + ν) D Ω \dot{K_{I}}}{12 \sqrt{2 π} k_{B} T} t^{2})}^{4 / 5}

𝛽(1/2, 9/10)은 베타 함수이고, 격자 길이는 $a_{0}$ 이다.

식 (16)을 식 (15)에 대입하면, 균열 주변의 수소의 양을 구할 수 있다.

(16)

K_{I} = {\dot{K}}_{I} t

(17)

N_{H} / L_{z} = A_{0} K_{I}^{8 / 5}

(18)

A_{0} = β (\frac{1}{2}, \frac{9}{10}) \frac{2 c_{0}}{a_{0}^{3}} (\frac{5 (1 + ν) D Ω}{12 \sqrt{2 π} k_{B} T {\dot{K}}_{I}})

식 (17), (18)에서 $A_{0}$ 는 하중, 확산, 농도, 온도의 조합으로 균열 주변의 수소 입자들을 제어하는 하나의 파라메터로 정의할 수 있다.

벽개가 발생하기 시작하는 $A_{0}^{*}$ 는 취성화의 유무를 결정하는 척도로 사용될 수 있다. Table 1과 같이 시뮬레이션을 통해 구한 $A_{0}^{*}$ 은 $1.4 \times 10^{11} M P a^{- 8 / 5} m^{- 9 / 5}$ 이다. 기존 연구들에서 수행한 실험결과들을 $A_{0}^{*}$ 을 중심으로 수소 취성화의 여부를 구분할 수 있어야 하기 때문에 취성화가 일어나지 않는 값들 중 최대값인 1.18×10¹¹과 최성화가 일어나는 값들 중 최소값인 1.99×10¹¹ 사이의 값이 되어야 한다(Song and Curtin, 2013).

(19)

1.18 \times 10^{11} M P a^{- 8 / 5} m^{- 9 / 5} < A_{0}^{*} < 1.99 \times 10^{11} M P a^{- 8 / 5} m^{- 9 / 5}

따라서, 본 연구의 시뮬레이션에서 구한 값은 실험 결과와 잘 일치하는 것을 확인할 수 있다.

Table 1.

Parametric study of molecular dynamics simulation

$N_{H} / L_{z}$	$T$	$\dot{K_{I}}$	$A_{0}$	Presence or absence of cleavage
0/Å	300K	$4 \times 10^{- 15} M P a m^{1 / 2} / s$	$0 M P a^{- 8 / 5} m^{- 9 / 5}$	×
100/Å			$1.02 \times 10^{11} M P a^{- 8 / 5} m^{- 9 / 5}$	×
200/Å			$1.41 \times 10^{11} M P a^{- 8 / 5} m^{- 9 / 5}$	○
300/Å			$1.43 \times 10^{11} M P a^{- 8 / 5} m^{- 9 / 5}$	○
300/Å	200K	$4 \times 10^{- 15} M P a m^{1 / 2} / s$	$1.56 \times 10^{11} M P a^{- 8 / 5} m^{- 9 / 5}$	○
	300K		$1.43 \times 10^{11} M P a^{- 8 / 5} m^{- 9 / 5}$	○
	400K		$1.49 \times 10^{11} M P a^{- 8 / 5} m^{- 9 / 5}$	○
300/Å	300K	$2 \times 10^{- 13} M P a m^{1 / 2} / s$	$1.52 \times 10^{11} M P a^{- 8 / 5} m^{- 9 / 5}$	○
300/Å	300K	$2 \times 10^{- 12} M P a m^{1 / 2} / s$	$1.61 \times 10^{11} M P a^{- 8 / 5} m^{- 9 / 5}$	○

3. 결 론

본 연구에서는 수소 취성화 현상에 대해 분자동역학 해석을 수행하였다. 균열 모델을 구성하고, 수정된 EAM/fs 퍼텐셜을 사용하였다. 분자동역학 시뮬레이션을 통해 원자 스케일에서 수소 취성화 현상을 분석하였다. 농도, 온도, 하중 속도에 대해 파라메터 연구를 수행하였다. 농도가 낮은 상태에서 농도를 점차적으로 증가시킴에 따라 취성화 정도가 증가하게 되고, 나중에 특정 수준에 도달하게 되면, 농도에 거의 영향을 받지 않게 되는 것을 확인할 수 있다. 온도가 낮을 때에는 수소가 균열 끝부분까지 완전히 모이지 못해 취성화 정도가 약하고, 너무 높으면 확산 속도가 빨라서 수소가 외부로 나가게 된다. 따라서 온도가 적당한 값을 유지할 때 취성화가 가장 활발하게 일어나는 것을 알 수 있다. 하중 속도는 수소의 확산속도에 비해 상대적으로 빠르면, 수소 원자가 균열 끝부분까지 도달하지 못하기 때문에 이 경우에도 수소 취성화가 잘 일어나지 않는 것을 확인할 수 있다. 또한, 열역학적 계산을 통해 본 연구의 시뮬레이션과 실험 결과를 비교하여 잘 부합되는 것을 확인하였다.

Acknowledgements

본 연구는 해양수산부 재원으로 선박해양플랜트연구소의 주요사업인 “친환경 선박 연료 신뢰성･안전성 평가기술 개발”에 의해 수행되었습니다(1525012984).

References

Cha, S.-H., Ha, S.-H., Cho, S. (2014) Design Sensitivity Analysis of Coupled MD-Continuum Systems Using Bridging Scale Approach, J. Comput. Struct. Eng. Inst. Korea, 27(3), pp.137~145. 10.7734/COSEIK.2014.27.3.137

Matsumoto, R., Taketomi, S., Matsumoto, S., Miyazaki, N. (2009) Atomistic Simulations of Hydrogen Embrittlement, Int. J. Hydrog. Energy, 34(23), pp.9576~9584. 10.1016/j.ijhydene.2009.09.052

Ramasubramaniam, A., Itakura, M., Carter, E.A. (2009) Interatomic Potentials for Hydrogen in α-Iron based on Density Functional Theory, Phys. Rev. B, 79(17), 174101. 10.1103/PhysRevB.79.174101

Song, J., Curtin, W.A. (2013) Atomic Mechanism and Prediction of Hydrogen Embrittlement in Iron, Nat. Mater., 12(2), pp.145~151. 10.1038/nmat347923142843

Song, J., Curtin, W.A. (2011) A Nanoscale Mechanism of Hydrogen Embrittlement in Metals, Acta Mater., 59(4), pp.1557~1569. 10.1016/j.actamat.2010.11.019

Taketomi, S., Matsumoto, R., Miyazaki, N. (2010) Atomistic Study of the Effect of Hydrogen on Dislocation Emission from a Mode II Crack Tip in Alpha Iron, Int. J. Mech. Sci., 52(2), pp.334~338. 10.1016/j.ijmecsci.2009.09.042

Zhang, W. (2016) Evaluation of Susceptibility to Hydrogen Embrittlement - A Rising Step Load Testing Method, Mater. Sci. Appl., 7(08), pp.389~395. 10.4236/msa.2016.78035

Journal of the Computational Structural Engineering Institute of Korea ISSN:1229-3059(Print) 2287-2302(Online) 한국전산구조공학회 논문집

Preview

Parametric Studies on Hydrogen Embrittlement in Liquified Hydrogen Tank using Molecular Dynamics Simulation

ABSTRACT

MAIN

Fig. 1.

Hydrogen embrittlement(Zhang, 2016)

Fig. 2.

Effect of hydrogen in crack propagation(Song and Curtin, 2011)

Fig. 3.

Atomistic configuration of simulation model

(1)

(2)

(3)

(4)

(5)

(6)

(7)

(8)

(9)

(10)

Fig. 4.

CNA parameter at NH/Lz=0/Å

Fig. 5.

CNA parameter at NH/Lz=100/Å

Fig. 6.

CNA parameter at NH/Lz=200/Å

Fig. 7.

CNA parameter at NH/Lz=300/Å

Fig. 8.

CNA parameter at T=200K

Fig. 9.

CNA parameter at T=300K

Fig. 10.

CNA parameter at T=400K

Fig. 11.

CNA parameter at KI˙=2×10-13MPam1/2s-1

Fig. 12.

CNA parameter at KI˙=2×10-12MPam1/2s-1

(11)

(12)

(13)

(14)

(15)

(16)

(17)

(18)

(19)

Table 1.

Parametric study of molecular dynamics simulation

Acknowledgements

References

CNA parameter at $N_{H} / L_{z} = 0 / Å$

CNA parameter at $N_{H} / L_{z} = 100 / Å$

CNA parameter at $N_{H} / L_{z} = 200 / Å$

CNA parameter at $N_{H} / L_{z} = 300 / Å$

CNA parameter at $T = 200 K$

CNA parameter at $T = 300 K$

CNA parameter at $T = 400 K$

CNA parameter at $\dot{K_{I}} = 2 \times 10^{- 13} M P a m^{1 / 2} s^{- 1}$

CNA parameter at $\dot{K_{I}} = 2 \times 10^{- 12} M P a m^{1 / 2} s^{- 1}$