Development of a Numerical Simulator for Swelling Behavior of the Buffer Material in the Engineered Barrier Induced by the Alteration under Hydro-Mechanical-Chemical Conditions

Taehyun Kim; Jeonghwan Hwang; Changsoo Lee; Jin-Seop Kim

doi:10.7474/TUS.2025.35.5.584

Preview

Original Article

Tunnel and Underground Space. 31 October 2025. 584-601
https://doi.org/10.7474/TUS.2025.35.5.584

Development of a Numerical Simulator for Swelling Behavior of the Buffer Material in the Engineered Barrier Induced by the Alteration under Hydro-Mechanical-Chemical Conditions

수리-역학-화학적 조건을 고려한 변질에 따른 공학적방벽 내 완충재 팽윤 거동 모사를 위한 해석 시뮬레이터 개발

Taehyun Kim¹^*

Jeonghwan Hwang¹

Changsoo Lee²

Jin-Seop Kim²

김 태현¹^*

황 정환¹

이 창수²

김 진섭²

¹Senior Researcher, Korea Atomic Energy Research Institute

²Principal Researcher, Korea Atomic Energy Research Institute

¹한국원자력연구원 선임연구원

²한국원자력연구원 책임연구원

^{*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 non-commercial use, distribution, and reproduction in any medium, provided the original work is properly cited.

ABSTRACT

A buffer material, a key component of the engineered barrier in a geological disposal system, exhibits swelling behavior when saturated. The chemical composition of groundwater can significantly influence its swelling capacity. Therefore, numerical analysis of the buffer’s swelling behavior, incorporating chemical effects, is essential to ensure the long-term stability of the disposal system over hundreds of thousands of years. In this study, we developed the OGSPHREEQC-FLAC3D simulator, which enables coupled Hydro-Mechanical-Chemical (HMC) analysis. OGSPHREEQC-FLAC3D integrates OGS-FLAC3D for thermal-hydro-mechanical (THM) simulation with PHREEQC for chemical analysis using a sequential coupling approach. The Barcelona Basic Model was employed to simulate swelling behavior induced by suction pressure. Chemical effects were incorporated by considering changes in ion concentration in the solute, and a coefficient reflecting montmorillonite concentration variation was applied. A series of numerical simulations were conducted to evaluate the swelling behavior of Na-type bentonite under one-directional flow conditions. To investigate the effect of potassium ion (K⁺) concentration in the injection water, K⁺ concentration was increased in the sample through advection and diffusion. As K⁺ concentration rose, the amount of illite produced by chemical alteration also increased, leading to a reduction in swelling pressure. These results demonstrate the capability of the developed simulator and provide a foundation for future research aimed at advancing fully coupled thermal-hydro-mechanical-chemical (THMC) numerical simulation for geological disposal systems.

Keywords

OGSPHREEQC-FLAC3D

Buffer material

Illitization

Swelling behavior

Hydro- Mechanical-Chemical coupled behavior

처분시스템 내 공학적방벽의 핵심 요소인 완충재는 포화 시 팽윤하는 특성을 지니고 있으며, 지하수의 화학 조성에 따라 팽윤하는 특성 또한 영향을 받게 된다. 따라서 수십만 년에 이르는 기간 동안 처분시스템의 장기 안정성 확보를 위해서는 화학적 영향을 고려한 완충재 팽윤 거동에 대한 수치해석적 연구가 필요하다. 본 연구에서는 수리-역학-화학 해석이 가능한 OGSPHREEQC-FLAC3D 해석 시뮬레이터를 개발하였다. OGSPHREEQC-FLAC3D 해석 시뮬레이터는 순차적 연동해석 기법을 적용하였으며, 열-수리-역학 복합거동 해석이 가능한 OGS-FLAC3D와 화학 거동 해석이 가능한 PHREEQC를 연동하였다. 흡입력에 따른 팽윤은 Barcelona Basic Model을 적용하였으며, 용질 내 이온 농도 변화에 따른 영향을 반영하였고, 몬모릴로나이트 함량 변화에 따른 영향 계수를 적용하였다. 개발된 해석 시뮬레이터의 적용성을 확인하기 위해 Na-벤토나이트를 대상으로 일차원 유동 시 팽윤 해석을 수행하였다. 주입수의 K⁺ 이온 농도에 따른 변질 영향을 평가하고자 하였으며, 이류 및 확산에 따라 시료 내 K⁺ 이온 농도가 증가함을 확인하였다. K⁺ 이온 농도 증가 시 변질에 따라 생성되는 일라이트 양이 증가하였으며, 일라이트화에 따라 팽윤압이 감소함을 확인하였다. 해석 결과는 향후 열-수리-역학-화학 복합거동을 고려한 심층 처분 해석 시뮬레이터 개발의 기반으로 활용할 수 있을 것으로 기대된다.

키워드

OGSPHREEQC-FLAC3D

완충재

일라이트화

팽윤 거동

수리-역학-화학 복합거동

MAIN

1. 서 론
2. 벤토나이트 완충재 변질 특성
3. OGSPHREEQC-FLAC3D 해석 시뮬레이터
3.1 해석 알고리즘 및 시뮬레이터 개발
3.2 완충재 변질에 따른 팽윤압 변화 모델
4. 완충재 변질을 고려한 HMC 복합거동 해석
4.1 해석 모델
4.2 입력변수
4.3 해석 결과
5. 결 론

1. 서 론

고준위방사성폐기물 처분장 내 공학적방벽의 핵심 구성요소인 완충재의 대표 재료는 팽윤 특성을 지닌 벤토나이트가 활용되고 있다. 벤토나이트 완충재 내에서는 방사성폐기물에서 발생하는 고온의 열로 인한 유체의 상변화, 가스 이동 등이 유발되며, 동시에 암반으로부터 지하수 유입으로 인해 수분 함량이 시간에 따라 크게 변화하게 된다. 이와 같은 수분 함량 변화는 벤토나이트 완충재의 팽윤 특성에 직접적인 영향을 미치게 되고, 팽윤 특성의 변화는 투수 특성의 변화로 이어져 향후 발생할 수 있는 장기적인 핵종 이동 특성에도 변화를 일으킬 수 있다. 따라서 공학적방벽 내에서 발생하는 복잡한 열-수리-역학-화학적(Thermal-Hydro- Mechanical-Chemical, THMC) 복합거동에 대한 이해를 높이기 위해 전 세계적으로 다양한 연구가 이뤄져 왔으며, 수치해석은 장기 거동 예측, 효율성 및 예측성, 반복성에 장점이 있으므로 그 필요성이 점차 높아지고 있다. 수치해석은 해석대상 영역의 특성을 모사하는 방법에 따라 크게 연속체 해석과 불연속체 해석으로 나누어 볼 수 있으며, 연속체 해석에서는 대상 영역을 유한개의 요소로 이산화(discretization)한 후 질량 보존, 운동량 보존, 에너지 보존 방정식을 미분방정식 형태로 구성하여 풀이하게 된다. 유한요소법, 유한차분법 등이 이에 해당하며, 고준위방사성폐기물 처분장의 안정성 평가를 위하여 다양한 수치해석 연구가 진행되어 왔다.

연속체 해석을 적용하고 있는 THM 복합거동 관련 해석 코드로는 TOUGH-FLAC (Rutqvist et al., 2002, 2011, Lee et al., 2020, 2021a), CODE_BRIGHT (Guimarães et al., 2006, Gens et al., 2010), COMSOL (Nardi et al., 2014, Kim et al., 2018), OpenGeoSys (Kolditz et al., 2012, Benisch et al., 2013), OGS-FLAC (Park et al., 2019, Kim et al., 2021, 2025) 등이 있다. 해당 코드들에서는 벤토나이트 완충재의 팽윤 거동을 모사하기 위해 다양한 구성 모델을 활용할 수 있으며 흡입력에 따라 선형적 팽윤 거동을 가정한 선형 팽윤 모델, 반대로 비선형적 팽윤 거동을 모사할 수 있는 Barcelona Basic Model (Alonso et al., 1990) 및 완충재 미시 구조 특성을 고려한 Barcelona EXpansive Model (BExM)(Gens et al., 2006) 등이 활용되고 있다.

이외에도 처분시스템 내 화학적 복합거동을 고려하기 위해 고준위방사성폐기물 처분장 안정성 검토를 위한 THC 복합거동 해석을 위해 TH 해석을 위한 COMSOL과 C 해석을 위한 PHREEQC를 연동한 연구가 소개된 바 있으며(Nardi et al., 2014), Rutqvist et al. (2014)는 화학적 특성 해석을 수행할 수 있는 TOUGHREACT와 역학적 해석을 담당하는 FLAC3D를 연동하여 벤토나이트로 충전된 점토암 기반의 처분장을 대상으로 THMC 해석을 수행하였다. Zheng et al. (2015, 2017)는 Rutqvist et al. (2014)와 같은 선형 탄성 모델에 기반한 TOUGHREACT-FLAC3D 해석 시뮬레이터를 활용하여 Kunigel 벤토나이트에 대한 팽윤압을 모사한 바 있다.

그러나 언급한 연구들은 처분시스템 해석 시 화학적 특성 혹은 역학적 특성에 의한 영향을 반영하지 않은 THM 및 THC 관련 연구만을 수행하였거나 THMC 해석을 수행하였더라도 완충재의 팽윤 거동 특성을 선형 거동으로 가정하였다는 한계점을 지니고 있다. 흡입력에 따른 벤토나이트의 비선형적 팽윤 거동 특성과 변질에 따른 팽윤 특성 변화를 함께 반영할 수 있는 해석 시뮬레이터는 연구가 미비한 상황이며, 따라서 본 연구에서는 완충재 변질에 따른 비선형적 팽윤 거동 특성 변화를 모사할 수 있는 해석 시뮬레이터 개발을 목표로 연구를 수행하였다. 단계적 개발을 통해 최종적으로 THMC 복합거동 해석이 가능한 해석 시뮬레이터 개발을 수행하고자 하며, 이번 논문에서는 선행 연구로 개발된 HMC 복합거동 해석 시뮬레이터에 대한 개발 내용 및 적용성 검증을 위해 수행한 해석 결과를 소개하고자 한다.

2. 벤토나이트 완충재 변질 특성

처분시스템 성능평가에 있어서 핵심 요소인 완충재의 대표 재료인 벤토나이트는 미시적으로 층상 구조(layer lattice) 및 확산 이중층 구조(diffused double layer structure)를 지니고 있으며, 수분과 접촉하면 팽윤하는 특성을 갖고 있다. 이에 따라 처분시스템 내로 지하수 유입 시 팽윤 과정을 통해 처분공 혹은 처분터널과 완충재 사이의 틈(gap)을 메우고 이를 통해 기밀성을 유지할 수 있다. 또한 낮은 투수계수를 지니고 양이온 교환 성능이 높기 때문에 수착을 통해 핵종 유출을 억제하는 역할을 할 수 있다. 따라서 벤토나이트의 팽윤 특성을 발생시키는 광물 및 해당 광물의 구조적인 특성을 이해할 필요가 있다.

벤토나이트의 주요 구성 광물로는 스멕타이트(몬모릴로나이트), 석영, 장석, 운모, 방해석, 제올라이트 등이 있으며 벤토나이트 종류별로 함량의 차이가 있으나 고준위방사성폐기물 처분장 완충재 후보군에 해당하는 벤토나이트의 경우 스멕타이트 함량이 60% 이상, 기타 광물이 나머지를 차지하는 것으로 보고된 바 있다(Lee et al., 2011, Lee et al., 2021b). 이중 팽윤 특성을 발현하는 것은 스멕타이트 그룹으로써 특히 벤토나이트 내에는 몬모릴로나이트가 주요 구성 광물로 자리하고 있다. 몬모릴로나이트의 화학 조성은 $R_{0.33}^{+} (A l_{1.67} M g_{0.33}) S i_{4} O_{10} (O H)_{2}$ 로 표현할 수 있고, 여기서 $R_{0.33}^{+}$ 은 층간 양이온으로써 Ca 혹은 Na 이 자리하여 벤토나이트 유형을 결정하게 된다.

Ca 혹은 Na가 들어가는지 여부에 따라 층간 두께가 달라지게 되며, Na 이온의 경우 Ca 이온보다 더 큰 양이온 크기 및 더 작은 수화 반경을 갖고 있다. 따라서 Na 이온이 보다 쉽게 층간 교환이 이뤄질 수 있으며, 일반적으로 Ca 벤토나이트보다 더 큰 팽윤압을 갖게 된다(Park et al., 2023).

이때 몬모릴로나이트는 지하수와 접촉 시 층간 양이온의 교환에 따라 팽윤 정도가 떨어지는 변질 현상이 발생할 수 있다. Fig. 1에서 볼 수 있는 바와 같이 사면체판 구조의 $S i^{4 +}$ 양이온을 $A l^{3 +}$ 양이온이 대체하게 되는 과정은 일라이트화의 주요 단계 중 하나이다(Lee and Choi, 2016). $A l^{3 +}$ 로 치환되면 입자 표면에 음전하가 증가하고, 그 결과로써 공극수 중의 $K^{+}$ 이온이 입자 층간에 끼어들게 되면서 일라이트화가 진행될 수 있다(Lee and Cho, 2007). 식 (1)은 용해-침전에 따른 완충재 내 스멕타이트의 일라이트화 과정을 표현한 식에 해당한다(Hwang et al., 2021).

(1)

S m e c t i t e + 0.52 H^{+} + 0.63 A l O_{2}^{-} + 0.6 K \to I l l i t e + 0.26 H_{2} O + 0.08 M g^{2 +} + 0.33 N a^{+} + 0.5 S i O_{2}

https://cdn.apub.kr/journalsite/sites/ksrm/2025-035-05/N0120350507/images/ksrm_2025_355_584_F1.jpg

Fig. 1.

A schematic diagram for illitization from smectite (Lee, 2024)

일라이트화가 진행될 경우, 층 간 간격이 감소하게 되고, 층간 양이온의 인력이 강해지게 되며 벤토나이트의 팽윤능을 점차적으로 상실하게 된다. 따라서 일라이트화 진행 정도를 예측하는 것은 처분시스템의 핵심 구성요소인 벤토나이트의 성능을 평가하고 유지하기 위한 중요한 요소라고 할 수 있다.

$K^{+}$ 이온 및 $A l^{3 +}$ 이온의 경우 주로 정장석( $K A l S i_{3} O_{8}$ ) 및 사장석( $C a A l_{2} S i_{2} O_{8}$ , $N a A l S i_{3} O_{8}$ )에 분포하고 있으며, 지하수에 용해/침전되어 벤토나이트 내부로 유입될 경우 일라이트화 반응을 가속화 할 수 있다. 이외에도 벤토나이트로부터 용해된 지하수 중의 실리카가 완충재에 침전되어 굳어지는 실리키화 반응 또한 완충재의 팽윤성을 크게 저하시킬 수 있으나 고준위폐기물 처분 환경에서 발생하는 열수반응은 대부분 일라이트화에 의한 것으로 알려져 있다(Lee and Cho, 2007, Lee and Choi, 2016, Park et al., 2023).

따라서 처분 환경을 고려한 고온/고압 조건에서 일라이트화 반응을 평가하기 위한 다양한 연구들이 수행되어 왔으며(Allen et al., 1984, Wood, 1983, Yau et al., 1987, Eberl and Hower, 1976, Anderson, 1983), $K^{+}$ 이온 증가에 따라 일라이트화 속도가 진전됨을 밝힌 바 있다. 다만 100°C 조건에서 일라이트화에 걸리는 시간이 104~105년이 소요됨을 보고한 바 있고(Oscarson and Cheung, 1983), 정상 조건의 지하수 조건에서는 일라이트화가 느리거나 전환이 거의 없을 것으로 예측한 바 있다(Oscarson and Cheung, 1983, Anderson, 1983). 따라서 본 연구에서는 벤토나이트 종류 및 다양한 처분 환경 조건에 따라 벤토나이트의 일라이트화 진행 정도를 모사할 수 있는 연구를 수행할 수 있는 해석 도구를 개발하고자 하였다.

3. OGSPHREEQC-FLAC3D 해석 시뮬레이터

3.1 해석 알고리즘 및 시뮬레이터 개발

열-수리-역학-화학적 복합거동 해석을 위한 연계해석 방식은 먼저 지배방정식을 풀이하는 방식에 따라 크게 일체식(monolithic)과 순차적(sequential) 연계 해석 방식으로 나누어 볼 수 있다. 일체식 연계해석 방법은 다중 물리 거동을 풀이하기 위해 구성된 미분방정식을 하나의 행렬로 구성하여 해석하는 방식(Park and Park, 2022)이며, 순차적 연계해석 방식은 다중 물리거동에 대한 각 미분방정식을 순차적으로 풀면서 먼저 풀이된 미분방정식의 결과를 다음 방정식의 초기 조건으로 넘겨주는 형태로 해석하게 된다. 각 연계해석 방식은 장단점이 있으며, 본 연구에서는 기존에 개발된 순차적 연계해석 방식을 적용한 OGS-FLAC3D 해석 시뮬레이터를 기반 해석 시뮬레이터로 선정하였다.

OGS-FLAC3D는 열-수리 해석이 가능한 OpenGeoSys와 역학 해석이 가능한 FLAC3D를 연계한 해석 시뮬레이터이다. 대표적 THM 복합거동 해석 시뮬레이터인 TOUGH-FLAC3D와 같은 연계해석 방식을 적용하고 있으나 메쉬 생성의 효율성 및 다상유체 모델의 다양성 측면에서 장점을 지니고 있으며(Park et al., 2019), 따라서 본 연구에서는 OGS-FLAC3D를 확장하여 개발을 진행하고자 하였다.

화학 해석과 관련하여 OpenGeoSys에서는 이류(advection)에 의한 유체 이동 및 핵종 이동 해석은 가능하지만, 광물 또는 핵종 간 평형 반응(equilibrium reaction)이나 동역학 반응(kinetic reaction)과 같은 핵심적인 화학 반응을 해석하는 기능은 탑재하고 있지 않다. 따라서 Xie et al.(2006)은 화학 반응 해석 수행이 가능한 PHREEQC를 OpenGeoSys와 순차적 연계해석 방식 기반으로 연동하여 THC 해석을 진행한 바 있다. PHREEQC (pH-REdox-EQuilibrium)는 다양한 수성 환경(Aqueous Environment)에서 지구화학적 계산을 수행하도록 설계된 C언어 기반의 무료 소프트웨어이다(Parkhurst et al., 1999). PHREEQC는 LLNL, Sit, Wateq4f, Sit, Pitzer 등 다양한 열역학 데이터베이스를 기본적으로 제공하며, 데이터베이스에 입력된 화학반응식 및 평형상수를 기반으로 (1) 이온 화학종(Aqueous Speciation) 분포 및 광물의 포화 지수(Saturation Index) 계산, (2) 광물(Mineral), 기체(Gas), 표면 착화(Surface Complexation), 양이온 교환(Cation Exchange) 등 복잡한 지화학 반응을 시뮬레이션하고 정량적으로 평가할 수 있다.

해당 기능은 OpenGeoSys에 내장되어 있으며, 사용자가 PHREEQC 관련 기능을 활성화하고 입력 파일을 설정하면 OGS PHREEQC 연동 해석을 수행할 수 있다. 따라서 본 연구에서는 OGSPHREEQC 해석을 수행할 수 있도록 화학 반응 관련 입력 파일을 설정한 후 OGS-FLAC3D 연동 해석을 위해 구성된 FLAC3D 내 연계 모듈을 수정하여 OGSPHREEQC를 통해 해석된 화학 반응 관련 결과가 입력되어 HMC 해석이 순차적으로 수행될 수 있도록 새로운 연계 모듈을 개발하였다.

다음은 연계해석 과정을 순서대로 설명하였으며, Fig. 2는 해당 과정을 도식화한 것이다.

1)먼저 OpenGeoSys에서는 수리 해석을 진행하기 위해 Darcy 법칙에 기반하여 다상 유체 및 광물의 이류-분산에 의한 이동(Advection-Dispersion Equation, ADE)을 풀게 된다.

2)OpenGeoSys에서 결정한 시간 단계(timestep) 내에서 수리 해석이 수행되고, 다종 요소의 질량 이동 방정식(multicomponent mass transport balance equation)에 근거하여 ADE에 의한 광물 혹은 핵종 이동이 계산된다.

3)PHREEQC에서는 요소별로 변경된 이온 농도를 업데이트한 후, 사용자가 설정한 평형 반응식(equilibrium equation) 및 동역학 반응식(kinetic reaction equation)에 따라 화학 반응 해석을 수행하게 된다.

4)해석 완료 후 압력, 흡입력, 매질 내 이온 농도 결과가 도출되면 FLAC3D 해석을 위한 연동 모듈로 해당 정보가 전달된다.

5)OGSPHREEQC 해석 결과는 절점(node)에서 얻어진 정보이며, FLAC3D 내 응력-변형률 해석은 요소(element)에서 진행되므로 연동 모듈에서는 요소를 구성하는 절점과 요소의 중심점 사이의 거리를 통해 계산된 가중산술평균(Fig. 3)을 활용하여(Kwon et al., 2021) FLAC3D 각 요소에 OGSPHREEQC 해석 결과를 전달한다.

6)FLAC3D에서는 전달된 결과를 바탕으로 역학 해석을 수행하게 된다.

7)역학 해석이 완료 후 응력, 변형률 및 공극률 해석 결과는 다음 단계의 수리 해석을 위해 OpenGeoSys로 전달되어 해석이 수행된다.

https://cdn.apub.kr/journalsite/sites/ksrm/2025-035-05/N0120350507/images/ksrm_2025_355_584_F2.jpg

Fig. 2.

A coupling sequence for THMC coupled simulation in OGSPHREEQC-FLAC3D

https://cdn.apub.kr/journalsite/sites/ksrm/2025-035-05/N0120350507/images/ksrm_2025_355_584_F3.jpg

Fig. 3.

Calculation process based on weighted arithmetic mean for node to element value from OGSPHREEQC to FLAC3D

3.2 완충재 변질에 따른 팽윤압 변화 모델

Zheng et al.(2015)은 완충재의 팽윤 특성에 영향을 미칠 수 있는 열-수리-화학적 요소들을 고려하여 완충재의 팽윤압 변화를 다음 식 (2)와 같이 표현한 바 있다.

(2)

d σ_{s} = 3 K β_{s w} d S l + A_{n} d C + A_{s c} d M_{s}

위 식에서 $d σ_{s}$ 는 팽윤압 변화량, 첫 번째 항인 $3 K β_{s w} d S l$ 는 흡입력에 따른 완충재 변형을 고려한 팽윤압을 표현한 것으로 $K$ 는 체적변형계수, $β_{s w}$ 는 수분 함량에 따른 팽윤 정도를 나타내는 계수(moisture swelling coefficient)이며, $S l$ 은 함수비(water saturation)를 의미한다. 두 번째 항( $A_{n} d C$ )은 이온 농도 변화 시 삼투압에 따른 팽윤압 변화를 모사하기 위한 항이며, $A_{n}$ 은 이온 농도 변화( $d C$ )와 관련된 상수이다. 세 번째 항은 몬모릴로나이트의 일라이트화에 따른 팽윤압 변화를 고려하기 위한 항이며, $A_{s c}$ 는 팽윤 특성을 지닌 점토 광물의 몰 분율 변화( $d M_{s}$ )와 관련된 상수에 해당한다.

위 관계식에서 완충재 팽윤 거동은 실제로는 포화도 변화에 따라 비선형적 거동을 보이지만 Zheng et al.(2015)의 연구에서는 흡입력에 따른 팽윤압 변화를 선형적 거동으로 가정하였다는 한계가 있다. 따라서 본 연구에서는 포화도 변화에 따라 비선형적 팽윤 거동을 보이는 BBM 모델을 기반으로 완충재 변질에 따른 팽윤압 변화 모델을 구축하고자 하였다. 따라서 본 연구에서는 식 (2)의 첫 번째 항( $3 K β_{s w} d S l$ )을 BBM에 따른 팽윤압 계산 항( $(K / K^{s}) d S l$ )으로 바꿔 적용하였다. 따라서 화학적 변질을 고려한 팽윤압 변화는 개발된 연동 모듈 내에서 다음 식 (3)과 같이 계산된다.

(3)

d σ_{s} = \frac{K}{K^{s}} d S l + A_{n} d C + A_{s c} d M_{s}

여기서 K는 완충재 체적변형계수(bulk modlus), K^s는 흡입력에 따른 강성 계수(suction bulk modulus)를 의미하며, K^s를 통해 흡입력에 따른 비선형적 거동을 표현할 수 있게 된다(Rutqvist et al., 2011, Lee et al., 2020, Kim et al., 2025).

4. 완충재 변질을 고려한 HMC 복합거동 해석

화학 조건을 정밀하게 제어하고 실험 시 실시간으로 몬모릴로나이트 감소량을 측정하는 것이 현실적으로 어렵기 때문에 개발된 해석 시뮬레이터의 정확한 비교 검증을 수행할 수 있는 실험 데이터 확보가 불가능하였다. 따라서 본 연구에서는 벤토나이트 완충재의 팽윤압 평가를 위해 실내 시험 시 일반적으로 활용되는 일방향 주입 시험을 모사하였고, 관련 물성들을 가정하여 일라이트화에 따른 팽윤압 변화 해석이 가능한지를 정성적으로 평가하고자 하였다.

4.1 해석 모델

실내 팽윤 실험은 벤토나이트 완충재의 수화에 따른 팽윤압 측정을 목적으로 수행된다. 일반적으로 바닥면에서 일정 압력으로 공극수를 주입하며 시간에 따라 시료 최상단에서 압력 변화를 측정하게 된다. 이때 공극수압과 팽윤압을 구분하여 측정하기 어렵기 때문에 시료 내 전응력을 측정하는 것이 일반적이다. 또한 수화에 걸리는 시간을 단축하기 위해 시료는 짧은 디스크 형태의 시편을 활용하며, 본 연구에서도 이와 같은 방식의 개념실험모델을 구성하고자 하였다. Fig. 4는 생성된 수치모델을 나타내고 있으며, 직경 5cm, 높이 3cm의 디스크 형태 시료를 활용하는 것으로 가정하고 축대칭 모델을 생성하였다. 수치 모델의 도메인은 직육면체(hexahedral) 요소로 구성하였으며, 시료의 하부, 중심, 상부에서 결과 계측을 수행하였다. FLAC3D의 경우 응력 및 변형률 등의 결과는 사용자가 지정한 지점과 가장 가까운 거리에 위치한 요소의 중심점에서 계측되지만(FLAC_O1, FLAC_O2, FLAC_O3) OGSPHREEQC의 압력, 흡입력, 이온 농도 등의 결과는 지정한 계측 지점과 가장 가까운 절점(O1, O2, O3)에서 계측이 수행된다(Fig. 4). 따라서 보다 정밀한 해석을 원할 시에는 계측 지점이 위치한 영역의 그리드는 조밀하게 구성할 필요가 있다.

https://cdn.apub.kr/journalsite/sites/ksrm/2025-035-05/N0120350507/images/ksrm_2025_355_584_F4.jpg

Fig. 4.

Axisymmetic numerical model with the monitoring points

실제 실험은 금속 셀 안에서 수행되며, 따라서 시료의 상하부는 변위 고정 경계 조건을 적용하였고, 시료의 측면부는 팽윤 발생이 가능하도록 측면에 수직한 변위만을 고정하였다. 실제 실험의 경우 금속 셀에 의한 마찰력이 팽윤압에 영향을 미칠 수 있으나 본 연구에서는 팽윤압을 정확하게 모사하는 것이 목적이 아니므로 마찰에 의한 영향은 고려하지 않았다. 주입압은 0.1 MPa로 시료의 하부에서 일정 압력으로 주입되는 것을 가정하였으며, 시료 내 초기 포화도는 0.39, 초기 유효응력은 0.1MPa이 되도록 응력 조건을 설정하였다.

4.2 입력변수

4.2.1 역학적 물성

BBM 활용을 위한 입력변수는 Table 1과 같다. BBM 모델의 입력변수는 증류수 사용 시 최대 팽윤압이 약 1.4 MPa이 될 수 있도록 가정한 값을 활용하였으며, 향후 한국원자력연구원에서 수행하고 있는 Ca-벤토나이트인 Bentonil-WRK에 대한 BBM 입력값이 측정되면 이를 활용하여 해석 및 비교 검증을 수행할 계획이다.

Table 1.

BBM input parameters for swelling test

Variable (unit)	Value	Variable (unit)	Value
$v_{c}$ (0) (-)	1.9	$r_{λ}$ (-)	0.74
$κ_{p s 0}$ (-)	0.0102	$β_{λ}$ (MPa-1)	0.07
$κ_{s p 0}$ (-)	0.108	$κ_{s}$ (-)	0.1
𝜈 (-)	0.4	$p_{s 0}$ (MPa)	0.0
$α_{s s}$ (-)	-0.010	$p^{c}$ (MPa)	1.325
$α_{p s}$ (-)	-0.002	$M$ (-)	1.261
$α_{s p}$ (-)	-0.166	$p_{a t m}$ (MPa)	0.1
$α_{α}$ (-)	0.53	$p_{r e f}$ (MPa)	0.007
$λ_{p s 0}$ (-)	1.50	$p_{0}^{*}$ (MPa)	1.45

화학적 변질에 따른 팽윤압 변화를 반영하기 위한 $A_{n}$ 과 $A_{s c}$ 는 기존 Laredj et al.(2010) 및 Zheng et al.(2015, 2017)이 사용한 관계식을 활용하였다. Laredj et al.(2010)은 팽윤성 점토를 대상으로 실험실에서 측정된 팽윤량과 흡입력으로 계산되는 변형이 다른 것을 확인했으며, 이러한 차이를 줄이기 위해 용질 내 이온 농도 변화로 인한 삼투압 변화에 따라 발생하는 변형을 추가함으로써 팽윤압 변화를 모사하고자 했다. 이를 위해 이온 농도 변화( $d C$ )에 따른 응력 증가분을 계산하기 위한 계수인 $A_{n}$ 을 도입하였으며, 다음 식 (4)와 같다.

(4)

A_{n} = - \frac{(5.312 \ln C - 23.569)}{\sqrt{C}} - \frac{7.252 \times 10^{- 4}}{C^{2}}

Zheng et al.(2015)은 기존에 보고된 바 있는 세계 각국 벤토나이트의 최대 팽윤압과 스멕타이트 함량을 조사하여 도시한 후 선형적 관계를 통해 $A_{s c}$ 의 범위를 2.0~9.3MPa로 제시한 바 있으며(Fig. 5), 본 연구에서는 이를 참고하여 $A_{s c}$ 값을 2MPa로 가정하였다.

https://cdn.apub.kr/journalsite/sites/ksrm/2025-035-05/N0120350507/images/ksrm_2025_355_584_F5.jpg

Fig. 5.

Swelling pressure versus smectite mass fraction based on the previous experimental observation (Zheng et al., 2015)

4.2.2 수리적 물성

Table 2는 해석에 사용된 수리해석과 관련된 입력변수들을 나타내고 있으며, 불포화토에서의 다상유동 해석을 위해 Richards 유동 모델을 확장한 Philip and de Vries의 모델을 기반으로 하였다. Richards 유동 모델에서는 액체 상태 물의 이동만을 고려하지만 Philip and de Vries 모델에서는 열에 의한 수증기 확산을 추가하여 TH 해석을 수행하게 된다. 따라서 고준위폐기물 처분장과 같이 온도 변화에 의한 영향이 발생하는 곳에서는 Philip and de Vries 모델의 활용이 필요하다고 판단된다. 해당 모델의 경우 Richards 유동 모델에서 사용된 가정을 그대로 적용하였으며, 매질 내에서 기체 이동은 충분히 빠르게 발생하므로 기체압은 소산되는 것으로 가정한다. 이에 따라 기체압은 0.1MPa로 고정된 상태가 되며, 확산 계수의 적용을 통해 온도 변화에 따른 기체의 분산(diffusion) 거동과 포화도, 압력 변화에 따른 유동 변화를 모사하게 된다(Kim et al., 2024).

Table 2.

Hydraulic input parameters for swelling test

Variable	Value
Initial saturation [-]	0.39
Permeability [m2]	2.1e-20
Porosity [-]	0.44
Storage coefficient [-]	3.0e-11
Entry pressure [MPa]	22.4
van Genuchten paramter, m [-]	0.359
Residual saturation [-]	0.0
Maximum saturation [-]	1.0

모세관압(capillary pressure)과 따른 포화도 간 관계는 다음 식 (5)와 같이 van Genuchten 경험식을 따르는 것으로 적용하였으며, 상대투수 계수는 식 (6)과 같다.

(5)

S_{e} = \frac{S_{l} - S_{r l}}{S_{s l} - S_{r l}} = {[1 + {(\frac{P_{c}}{P_{o}})}^{\frac{1}{1 - m}}]}^{- m}

(6)

k_{r l} = \sqrt{S_{e}} {[1 - {(1 - S_{e}^{\frac{1}{m}})}^{m}]}^{2}

식 (5)에서 $S_{e}$ 는 유효 포화도, $S_{l}$ 은 현재 포화도, $S_{s l}$ 은 최대 포화도, $S_{r l}$ 은 잔류 포화도를 의미한다. $P_{c}$ 는 모세관압, $P_{o}$ 는 진입압력(entry pressure), $m$ 은 물질상수에 해당하며, 식 (6)에서 $k_{r l}$ 은 상대투수계수를 의미한다. 따라서 다상유동 해석 시 모세관압과 포화도 간 관계에 따라 유효 포화도가 결정되며, 유효 포화도 변화에 따라 상대투수계수가 바뀌면서 유체 흐름의 변화가 발생하게 된다.

4.2.3 화학적 물성

화학 해석을 위해서는 먼저 벤토나이트 내의 몬모릴로나이트와 변질 시 일라이트에 대한 화학식을 구성할 필요가 있으며, 앞서 언급한 바와 같이 본 연구에서는 Na-벤토나이트를 가정하여 각각 다음 식 (7) 및 식 (8)과 같이 구성하였다.

(7)

N a_{0.33} M g_{0.33} A l_{1.67} S i_{4} O_{10} (O H)_{2} = 0.33 M g^{2 +} + 0.33 N a^{+} + 1.67 A l O^{2 -} + 0.66 H_{2} O + 4 S i O_{2} + 0.68 H^{+}

(8)

K_{0.6} M g_{0.25} A l_{1.8} (A l_{0.5} S i_{3.5} O_{10}) (O H)_{2} = 0.6 K^{+} + 0.25 M g^{2 +} + 2.3 A l O^{2 -} + 0.4 H_{2} O + 3.5 S i O_{2} + 1.2 H^{+}

또한 동역학 반응을 모사하기 위해서는 반응과 관련된 상수 및 온도 변화 시 반응 변화를 반영할 수 있는 상수를 설정해야 하며, 이에 대한 입력변수는 Table 3에 나타난 바와 같다. 용이한 분석을 위하여 $K^{+}$ 이온만 포함된 증류수를 주입하는 것으로 가정하였고, 주입부의 $K^{+}$ 농도는 0.1 mol, 0.5 mol, 5.0 mol로 증가시켜 가며 농도 증가에 따른 변화를 살펴보았다. 일반적으로 지하수 내 $K^{+}$ 농도는 채취된 지역에 따라 차이를 보일 수 있으며, 본 연구에서는 그 영향을 살펴보기 위해 큰 편차를 두고 설정하였다. 향후 팽윤 실험에 대한 실험값이 확보될 시에는 지하수의 조성도 함께 측정하여 해석에 고려할 필요가 있다.

Table 3.

Chemical input parameters and initial conditions for swelling test

	Variable	Initial value
Chemical species with initial concentration	K	1.0e-8
	Ca	1.0e-8
	Mg	1.0e-8
	Na	1.0e-8
	Cl	1.0e-8
	Si	1.0e-8
	C	1.0e-8
	Al	1.0e-8
Chemical conditions	pH	1.0e-8
	log k at 25°C of Na-Montmorillonite (Temperature dependent reaction coefficient)	-34.622
	delta H of Na-Montmorillonite [kJ/mol] (Calculated enthalpy of reaction)	133.98
	log k at 25°C of Illite (Temperature dependent reaction coefficient)	-42.6981
	delta H of Illite [kJ/mol] (Calculated enthalpy of reaction)	165.83

매질 내 초기 이온 농도 역시 Table 3에 기록하였으며, 해석 모델 내 초기 일라이트는 존재하지 않는 것으로 하였고 해석 편의성을 위해 벤토나이트는 모두 몬모릴로나이트로 이루어진 것으로 가정하였다. 따라서 해석 결과 시 계측되는 일라이트는 모두 몬모릴로나이트의 변질로 인해 발생하며, 변질은 초기 포화도에 의한 평형 반응(equilibrium reaction) 및 유체 흐름에 의해 전달되는 $K^{+}$ 에 의한 동역학 반응(kinetic reaction)의 합으로 발생하게 된다.

4.3 해석 결과

4.3.1 포화도 변화

해석 시간은 포화 및 일라이트화 반응이 발현되는 시간을 고려하여 20,000일까지 수행하였다. Fig. 6(a)는 $K^{+}$ 0.1 mol을 포함한 이온수를 주입하였을 때 계측 위치별 시간에 따른 포화도 변화를 나타내고 있으며, 주입수 내 $K^{+}$ 농도가 0.5 mol 및 5.0 mol 일 때에도 유사한 양상을 나타냈다. 시료 전체 도메인의 초기 포화도는 0.39이며, O1 지점은 시료의 최하부 및 공극수 주입부이기 때문에 거의 주입 즉시 포화 상태에 도달하는 것을 확인할 수 있다. 시료 위치에 따라 O2, O3의 순서로 순차적으로 포화가 진행됨을 확인할 수 있으며, 최종적으로 포화가 완료되는 시점은 주어진 조건에서 약 20일 경과 후인 것으로 나타났다. O2, O3 지점에서 포화도가 감소 후 증가하게 되는데, 이는 해석 요소 크기 등의 영향으로 수렴도가 불안정하여 일시적으로 발생하는 현상으로 판단된다. 투수계수는 일정한 값을 활용하였으며, 공극률 변화는 변형률을 반영하여 FLAC3D 내 BBM을 통해 계산된 후 연동 모듈을 통해 OGSPHREEQC로 전달된다. 화학적 조성 변화에 따른 공극률의 변화는 현 단계에서는 고려하지 않았기 때문에 주입수 이온 농도에 따른 포화도 변화는 나타나지 않았다.

https://cdn.apub.kr/journalsite/sites/ksrm/2025-035-05/N0120350507/images/ksrm_2025_355_584_F6.jpg

Fig. 6.

Saturation and $K^{+}$ change depending on time at each monitoring point when injecting $K^{+}$ 0.1 mol

4.3.2 $K^{+}$ 이온 농도 변화

Fig. 6(b)는 시간 경과 시 공극수 내 $K^{+}$ 농도에 따라 각 계측 지점에서 관측된 $K^{+}$ 농도를 보여주고 있다. O1, O2, O3는 OGSPHREEQC 해석을 통해 출력된 결과로 절점에서 출력된 결과이며, FLAC_O1, FLAC_O2, FLAC_O3는 OGSPHREEQC에서 FLAC3D로 전달된 결과를 앞서 상술한 가중산술평균을 통해 계측 지점과 가장 가까운 요소의 중심부에서 계산된 결과에 해당한다.

O1 지점의 경우 시료의 최하부로써 OGSPHREEQC에서는 주입부에 해당하는 위치의 절점이기 때문에(Fig. 4), $K^{+}$ 이온 농도가 주입 조건과 동일하게 계속 유지되는 것을 확인할 수 있다. FLAC_O1의 경우 주입부 절점을 포함한 요소를 구성하는 전체 절점들의 가중산술평균 값에 해당하기 때문에 FLAC_O1 요소가 완전 포화되기 전까지는 OGSPHREEQC 값과 차이를 보이게 된다. 하지만 10^-1일 경과 후 포화도가 높아짐에 따라 점차적으로 OGSPHREEQC 결과 값과 일치해 감을 볼 수 있다. O2와 O3 지점은 주입부에서 떨어진 지점이기 때문에 OGSPHREEQC 해석값과 FLAC3D 출력값이 비교적 잘 일치하는 것을 확인할 수 있다.

주입수 내 이온의 이동은 이류(advection)가 가장 큰 영향을 미치게 되므로 포화도 변화와 밀접한 관련이 있으며, 앞서 포화도 해석 결과와 마찬가지로 약 20일 경과 후 모든 지점에서의 $K^{+}$ 이온 농도가 초기 10^-8 mol에서 주입수 이온 농도까지 상승함을 확인할 수 있다.

4.3.3 몬모릴로나이트 및 일라이트 함량 변화

Fig. 7의 (a), (c) 및 (e)는 시간 경과 시 주입수 내 $K^{+}$ 이온 농도에 따라 변화하는 시료 내 몬모릴로나이트의 농도 변화를 나타내고 있다. $K^{+}$ 0.1 mol 주입 시 O1, O2, O3 지점에서 최종적으로 계측되는 몬모릴로나이트 농도는 (0.619, 0.637, 0.650) mol이며, $K^{+}$ 0.5 mol 주입 시에는 (0.617, 0.625, 0.634) mol, $K^{+}$ 5.0 mol 주입 시에는 (0.591, 0.593, 0.600) mol이 남은 것으로 나타났다.

https://cdn.apub.kr/journalsite/sites/ksrm/2025-035-05/N0120350507/images/ksrm_2025_355_584_F7.jpg

Fig. 7.

Montmorillonite and illite change depending on time and injection $K^{+}$ ion concentration at each monitoring point

다음으로 Fig. 7의 (b), (d), (f)는 일라이트화 반응에 따라 새롭게 생성된 시료 내 일라이트의 농도를 나타낸 결과이다. $K^{+}$ 0.1 mol 주입 시 O1, O2, O3 지점에서 최종적으로 계측되는 일라이트 농도는 (0.276, 0.264, 0.254) mol이며, $K^{+}$ 0.5 mol 주입 시에는 (0.278, 0.273, 0.266) mol, $K^{+}$ 5.0 mol 주입 시에는 (0.297, 0.296, 0.291) mol이 생성된 것으로 나타났다. 결과를 살펴보면 주입부에 가까울수록 $K^{+}$ 이온에 접촉된 시간이 늘어나기 때문에 일라이트화가 많이 진행되어 생성된 일라이트 양은 많고, 이에 따라 남아있는 몬모릴로나이트의 양은 적은 것을 알 수 있다.

계측 지점 O3에서 20,000일 경과 후, 주입수 내 $K^{+}$ 농도 변화의 영향을 살펴보면 주입수 $K^{+}$ 이온 농도가 증가함에 따라 최종적으로 시료 내 생성된 일라이트의 양은 증가하고 몬모릴로나이트의 양은 감소한 것을 확인할 수 있다(Table 4). 다만 주입수 내 $K^{+}$ 농도 변화에 따라 일라이트 및 몬모릴로나이트의 변화량이 확연한 차이를 보이지는 않았는데 이는 화학 해석 시 적용한 반응 상수 등에 의한 영향이 있을 수 있는 것으로 판단된다. 따라서 보다 정밀한 해석을 위해서는 사용 대상이 되는 벤토나이트에 대한 화학 반응 상수를 실험적으로 측정하여 적용하는 것이 필요할 것으로 판단된다.

Table 4.

Summary of the simulation results after 20,000 days at the monitoring point O3 (top surface)

$K^{+}$ injection (mol)	0.0	0.1	0.5	5.0
Saturation (-)	1.0	1.0	1.0	1.0
$K^{+}$ concentration (mol)	0.0	0.1	0.5	5.0
Montmorillonite concentration (mol)	1.0	0.650	0.634	0.600
Illite concentration (mol)	0.0	0.254	0.266	0.291
Total axial stress (MPa)	1.410	0.698	0.679	0.618
Max. total axial stress (MPa)	1.410	1.411	1.411	1.408

또한 실제 실험 시에는 일라이트화 반응 시 중간 단계인 몬모릴로나이트/일라이트 혼합층(mixed layer)을 형성하게 되며, 이러한 혼합층 또한 팽윤압 감소에 큰 영향을 미치는 것으로 알려진 바 있다(Lee et al., 2023). 그러나 본 해석에서는 이러한 중간 단계를 모사할 수 없으므로 화학 반응이 완료되어 몬모릴로나이트에서 일라이트로 완전히 전환되었을 때에만 농도 변화가 나타나며, 이에 따라 반응에 걸리는 시간이 크게 나타나는 것으로 생각된다. 따라서 이러한 중간 단계 혼합층 형성에 대한 연구 역시 향후 수행될 필요가 있을 것으로 판단된다.

4.3.4 팽윤압 변화

Fig. 8은 주입수 이온 농도 변화에 따른 축방향 전응력 해석 결과를 나타내고 있다. 앞서 상술한 바와 같이 실제 실험 시 팽윤압과 공극압을 구분하여 측정하기 어렵기 때문에 일반적으로 전응력 변화를 측정하며, 본 해석에서도 향후 실험 결과와 비교 검증을 위해 전응력 해석 결과를 비교하였다. 또한 일라이트화가 없을 때의 결과와 비교를 위하여 증류수 주입을 가정한 해석 결과(검은색 실선)을 함께 도시하였다.

증류수 주입 조건 시 최대 전응력은 약 1.41MPa였으며, 최대 전응력은 시료가 완전 포화에 도달하는 시점인 대략 20일 경과 후 발현되는 것으로 나타났다(Fig. 8(a)). 증류수 주입 시 최대 전응력 도달 후에는 큰 변화 없이 일정하게 유지되는 것으로 나타났으며, 이는 시료가 완전 포화됨에 따라 흡입력이 일정하게 유지되었기 때문이다. BBM은 흡입력에 따라 강성이 변화하는 함수를 갖고 있으며, 따라서 초기 포화도가 증가하는 구간에서는 흡입력 감소에 따라 강성이 증가하여 팽윤압이 증가하게 되지만 완전 포화 이후에는 흡입력이 0으로 일정하게 유지됨에 따라 일정한 강성 값을 갖게 되고 전응력 또한 변화하지 않는 것을 확인할 수 있다.

https://cdn.apub.kr/journalsite/sites/ksrm/2025-035-05/N0120350507/images/ksrm_2025_355_584_F8.jpg

Fig. 8.

Total axial stress change and enlarged view depending on time

하지만 $K^{+}$ 이온 주입에 따라 일라이트화가 진행될 경우 다른 결과가 나타나는 것을 확인할 수 있다. 먼저 최대 전응력 도달 시점인 약 20일 경과 시까지는 주입수 이온 농도 변화에 따라 전응력이 큰 차이를 보이지 않는 것을 알 수 있다. 그러나 약 20일 경과 이후부터 전응력의 감소가 급격하게 발생함을 확인할 수 있다(Fig. 8(b)). 이는 앞서 포화시점 및 $K^{+}$ 이온 농도가 최대값에 도달하는 시점과 거의 유사함을 확인할 수 있다(Fig. 6). 전체 도메인에 걸쳐 $K^{+}$ 농도가 최대값에 도달하게 되면 일라이트화 반응이 가속화되어 몬모릴로나이트의 감소량 또한 빠르게 증가하는 것을 확인할 수 있으며(Fig. 7(a, c, e)), 반대급부로 일라이트의 양 또한 20일을 기준으로 증가 속도가 빨라지는 것을 확인할 수 있다(Fig. 7(b, d, f)). 따라서 팽윤압의 감소는 현재 적용한 관계식에서 몬모릴로나이트의 양에 지배적인 영향을 받는 것을 알 수 있으며, 20,000일 결과 후 주입 수 농도에 따른 최종 전응력은 $K^{+}$ 0.1 mol 주입 시 0.698 MPa, 0.5 mol 주입 시 0.679 MPa, 5.0 mol 주입 시 0.618 MPa인 것으로 나타났다(Table 4). 이를 통해 주입수의 $K^{+}$ 이온 농도 증가 시 일라이트화 반응이 가속화됨을 알 수 있으며, 더 정확한 정량적 해석을 위해서는 관련 상수들에 대한 후속 실험 연구가 추가될 필요가 있음을 확인할 수 있다.

5. 결 론

본 연구에서는 고준위방사성 폐기물 처분 시스템의 핵심 구성요소인 완충재의 팽윤 거동을 수치해석을 통해 모사하고자 하였다. 특히 지하수의 화학 조성에 따라 발생할 수 있는 완충재 내 일라이트화에 따른 팽윤 거동을 구현하는 것을 목표로 하여 연구를 진행하였다.

수치해석 코드로는 순차적 연동해석 방식을 활용하여 처분 시스템 내 열-수리-역학 복합거동 해석을 위해 개발된 OGS- FLAC3D 해석 시뮬레이터를 기반으로 활용하였다. 일라이트화 반응 모사를 위해서는 화학 반응 해석이 수행되어야 하므로 화학 반응 해석을 수행할 수 있는 PHREEQC를 추가하여 연동 해석을 할 수 있도록 해석 알고리즘 및 시뮬레이터 개발을 수행하였으며, 최종적으로 순차적 연계 해석을 기반으로 한 OGSPHREEQC-FLAC3D 해석 시뮬레이터를 개발하였다.

화학적 변질에 따른 팽윤 거동 변화 모사를 위해서 흡입력에 따른 팽윤압 변화, 용질 내 화학 조성에 따른 삼투압 변화, 몬모릴로나이트 함량 변화를 고려한 항들을 산술적으로 합산하여 최종적으로 팽윤압 변화를 모사할 수 있는 모델을 구성하였다. 이때 흡입력에 따른 팽윤압 변화는 비선형적 팽윤 거동 모사가 가능한 Barcelona Basic Model을 활용하여 현실적 팽윤 거동을 모사하고자 하였다.

구축한 모듈의 검증을 위해서 팽윤압 확인을 위해 대표적으로 활용되고 있는 실내 팽윤실험을 모사하여 해석을 수행하였다. 디스크 형태의 시편을 적용하였으며 주입수 내의 $K^{+}$ 이온 농도를 변화시켜 가며 팽윤압 변화를 관측하였다. 해석 결과 $K^{+}$ 이온 농도 증가 시 몬모릴로나이트 함량의 감소량 및 일라이트의 증가량이 늘어나는 것을 확인할 수 있었고, 이를 기반으로 팽윤압 감소 또한 모사됨을 확인하였다. 따라서 개발된 해석 시뮬레이터가 수리-화학-역학적 연계해석을 고려하여 완충재의 팽윤압 거동 특성을 모사하는 것에 활용될 수 있을 것으로 판단된다.

다만 해석에 활용된 입력 변수들은 개발된 해석 시뮬레이터의 적용성 확인을 위해 기존의 참고문헌들을 통해 가정한 값이며, 향후 한국원자력연구원에서 실험 중인 Bentonil-WRK 벤토나이트에 대한 실험적 계측을 수행한 후 실험 결과와 해석 결과의 비교 검증을 수행할 계획이다. 또한 실제 실험 시 관측되는 몬모릴로나이트/일라이트의 혼합층은 본 해석에서는 구현되지 않았으며, 혼합층이 팽윤압 변화에 미치는 영향이 클 경우, 관련 항의 추가적인 고려가 필요할 수 있다. 추가로 본 연구는 온도에 대한 영향을 고려하지 않았으며 온도의 변화를 고려할 경우, 화학 반응 속도 및 불포화 매질에서의 유체 거동에 영향을 미칠 수 있다. 따라서 추후 온도를 고려한 열-수리-역학-화학 복합거동 해석 시뮬레이터의 단계적 개발 완료 시 실험값 및 타 해석 시뮬레이터와의 단계적 비교 및 검증을 통해 다양한 고준위방사성폐기물 처분 관련 연구에 적용할 수 있을 것으로 기대된다.

Acknowledgements

본 연구는 과학기술정보통신부의 재원으로 사용후핵연료관리핵심기술개발사업단(RS-2021-NR056198)의 지원을 받아 수행되었습니다.

References

Allen, C.C., Lane, D.L., Palmer, R.A., and Jonston, R.C., 1984, Experimental studies of packing material stability, In: G.L. McVay (Editor), The Scientific Basis for Nuclear Waste Management, Material Research Society, pp. 105-112.

10.1557/PROC-26-105

Alonso, E.E., Gens, A., and Josa, A., 1990, A constitutive model for partially saturated soils, Geotechnique, 40(3), 405-430.

10.1680/geot.1990.40.3.405

Anderson, D.M., 1983, Smectite alteration, Proceedings of a Colloquium, State Uni. of New York at Buffalo, Karn Bransle Sakerhet Teknist Rapport, TR-83-03.

Benisch, K., Graupner, B., and Bauer, S., 2013, The coupled OpenGeoSys-Eclipse simulator for simulation of CO2 storage-code comparison for fluid flow and geomechanical processes, Energy Procedia, 37, 3663-3671.

10.1016/j.egypro.2013.06.260

Eberl, D. and Hower, J., 1976, Kinetics of illite formation, Geol. Soc,. Am. Bull. 87, 1326-1330.

10.1130/0016-7606(1976)87<1326:KOIF>2.0.CO;2

Gens, A., Guimarães, L.do N., Olivella, S., and Sánchez, M., 2010, Modelling thermo-hydro-mechano-chemical interactions for nuclear waste disposal, Journal of Rock Mechanics and Geotechnical Engineering, 2(2), 97-102.

10.3724/SP.J.1235.2010.00097

Gens, A., Sanchez, M., and Sheng, D., 2006, On constitutive modelling of unsaturated soils, Acta Geotechnica, 1, 137-147.

10.1007/s11440-006-0013-9

Guimarães, L.do N., Gens, A., Sánchez M., and Olivella, S., 2006, THM and reacvitve transport analysis of expansive clay barrier in radioactive waste isolation, Communications in Numerical Methods in Engineering, 22, 849-859.

10.1002/cnm.852

Hwang, J., Choung, S., Han, W.S., and Yoon, W., 2021, Effects of bentonite illitization on Cesium sorption, Journal of Soil Groundwater Environment, 26(5), 29-38.

Kim, M.J., Lee, S.R., Yoon, S., Jeon, J.S., and Kim, M.S., 2018, Effect of thermal properties of bentonite buffer on temperature variation, Journal of the Korean Geotechnical Society, 34(1), 17-24.

10.7843/KGS.2018.34.1.17

Kim, T., Park, C.H., Lee, C., Kim, J.S., Park, E.S., and Bastian, G., 2024, A numerical analysis of Thermo-Hydro-Mechanical behavior in the FE experiment at Mont Terri URL: Investigating capillary effects in bentonite on the disposal system, Geomechanics for Energy and the Environment, 40, 100597.

10.1016/j.gete.2024.100597

Kim, T., Kim, G.I., Lee, C., and Kim, J.S., 2025, Implementation of Barcelona Basic Model into OGS-FLAC3D for simulating elasto-plastic behavior of unsaturated soils, Tunnel and Underground Space, 35(2), 112-124.

Kim, T., Park, C.H., Lee, C., and Kim, J.S., 2021, A numerical study on the Step 0 benchmark test in Task C of DECOVALEX-2023: Simulation for thermo-hydro-mechanical coupled behavior by using OGS-FLAC, Tunnel and Underground Space, 31(6), 610-622.

Kolditz, O, Bauer, S., Bilke, L. Böttcher, N., Delfs, J.O., Fischer, T., Görke, U.J., Kalbacher, T., Kosakowski, G., McDermott, C. et al., 2012, Opengeosys: an open-source initiative for numerical simulation of thermo-hydro-mechanical/chemical (thm/c) processes in porous media, Environmental Earth Sciences, 67(2), 589-599.

10.1007/s12665-012-1546-x

Kwon, S., Kim, G.I., Lee, C., Lee, J., and Kim, J.S., 2021, Development and Verification of TOUGH-3DEC Linking Algorithms for Coupled Hydro-Mechanical Analyses, Korea Atomic Energy Research Institute Technical Report, KAERI/TR-8833/2021, p. 52.

Laredj, N., Missoum, H., and Bendani, K., 2010, Modeling the effect of osmotic potential changes on deformation behavior of swelling clays, Journal of Porous Media, 13(8), 743-748.

10.1615/JPorMedia.v13.i8.50

Lee, C., Lee, J., and Kim, G.Y., 2021a, Numerical analysis of coupled hydro-mechanical and thermo-hydro-mechanical behavior in buffer materials at a geological repository for nuclear waste: Simulation of EB experiment at Mont Terri URL and FEBEX at Grimsel test site using Barcelona basic model, International Journal of Rock Mechanics and Mining Sciences, 139, 1-22.

10.1016/j.ijrmms.2021.104663

Lee, G.J., Yoon, S., and Cho, W.J., 2021b, Effect of bentonite type on thermal conductivity in a HLW repository, Journal of the Nuclear Fuel Cycle and Waste Technology, 19(3), 331-338.

10.7733/jnfcwt.2021.19.3.331

Lee, C., Lee, J., Kim M., and Kim, G.Y., 2020, Implementation of Barcelona Basic Model into TOUGH2-MP/FLAC3D, Tunnel & Underground Space, 30(1), 39-62.

10.7474/TUS.2020.30.1.039

Lee, J.O., and Cho, W.J., 2007, Hydrothermal behaviors and long-term stability of bentonite buffer material, Journal of the Korean Radioactive Waste Society, 5(2), 145-154.

Lee, J.O., Cho, W.J., and Kwon, S., 2011, Thermal-hydro-mechanical properties of reference bentonite buffer for a Korean HLW repository, Tunnel and Underground Space, 21(4), 264-273.

10.7474/TUS.2011.21.4.264

Lee, J.O., and Choi, H.J., 2016, Illitization and long-term stability prediction of domestic Ca-bentonite, Korea Atomic Energy Research Institute Technical Report, KAERI/TR-6694/2016, p. 58.

Lee, K.J., Ohazuruike, L., and Xiong, X., 2023, Effect of smectite illitization on swelling behavior of Gyeongju bentonite, Progress in Nuclear Energy, 166, 104938.

10.1016/j.pnucene.2023.104938

Lee, S.Y., 2024, Bentonite as an engineered barrier material: Understanding of Illitization and the application to an experiment, Special seminar in Korea Atomic Energy Research Institute, Powerpoint presentation.

Nardi, A., Idiart, A., Trinchero, P., de Vries, L.M., and Molinero, J., 2014, Interface COMSOL-PHREEQC (iCP), an efficient numerical framework for the solution of coupled multiphysics and geochemistry, Computers and Geosciences, 69, 10-21.

10.1016/j.cageo.2014.04.011

Oscarson, D.W., and Cheung, S.C.H., 1983, Evauation of phyllosilicates as a buffer component in the dispersal of nuclear fuel waste, AECL-7812, AECL, Canada.

Park, C.H., Kim, T., Park, E.S., Jung, Y.B., and Bang, E.S., 2019, Development and verification of OGSFLAC simulator for hydro-mechanical coupled analysis: Single-phase fluid flow analysis, Tunnel and Underground Space, 29(6), 468-479.

Park, D.H. and Park, C.H., 2022, Performance evaluation of OGS-FLAC simulator for coupled thermal-hydrological-mechanical analysis, Tunnel and Underground Space, 32(2), 144-159.

10.7474/TUS.2022.32.2.144

Park, T.H., Cha, Y.H., Kim, M.S., Hong, C.H., Kim, J.W., Lee, C., Kim, J.S., Lee, K.J., Xiong, X., and You, J., 2023, Experimental Report on the Swelling Behavior of Engineered Barrier in response to Salinity and Temperature Variations, KAERI/ iKSNF-SD2/1-23-SCR-05, SCR report.

Parkhurst, D.L., and Appelo, C.A.J., 1999, User’s Guide to PHREEQC (Version 2) A Computer Program for Speciation Batch- reaction, One-dimensional Transport, and Inverse Geochemical Calculations, USGS Water-Resources Investigations Report 99-4259.

Rutqvist, J., Ijiri, Y., and Yamamoto, H., 2011, Implementation of the Barcelona Basic Model into TOUGH-FLAC for simulations of the geomechanical behavior of unsaturated soils, Computers & Geosciences, 37, 751-762.

10.1016/j.cageo.2010.10.011

Rutqvist, J., Wu, Y.S., Tsang, C.F., and Bodvarsson, G., 2002, A modeling approach for analysis of coupled multiphase fluid flow, heat transfer, and deformation in fractured porous rock, International Journal of Rock Mechanics and Mining Sciences, 39(4), 429-442.

10.1016/S1365-1609(02)00022-9

Rutqvist, J., Zheng, L., Chen, F., Liu, H.H., and Birckholzer, J., 2014, Modeling of coupled thermo-hyrdro-mechanical processes with links to geochemistry associated with bentonite-backfilled repository tunnels in clay formations, Rock Mechanics and Rock Engineering, 47, 167-186.

10.1007/s00603-013-0375-x

Wood, M.L., 1983, Experimental investigation of sodium bentonite stability in Hanford basalt, In: D.G. Brookins (Editor), The Scientific Basis for Nuclear Waste Management, Material Research Society, pp. 727-743.

10.1557/PROC-15-727

Xie, M., Bauer, S., Kolditz, O., Nowak, T., and Shao, H., 2006, Numerical simulation of reactive processes in an experiment with partially saturated bentonite, Journal of Contaminant Hydrology, 83, 122-147.

10.1016/j.jconhyd.2005.11.003

Yau, Y.C., Peacor, D.R., Essene, E.J., Lee, J.H., Kuo, L.C., and Cosca, M.A., 1987, Hydrothermal treatment of smectite, illite, and basalt to 460: Comparison of natural with hydrothermally formed clay minerals, Clays and Clay Minerals, 35, 241-250.

10.1346/CCMN.1987.0350401

Zheng, L, Rutqvist, J., Xu, H., and Birkholzer, J.T., 2017, Coupled THMC models for bentonite in an argillite repository for nuclear waste: Illitization and its effect on swelling stress under high temperature, Engineering Geology, 230, 118-129.

10.1016/j.enggeo.2017.10.002

Zheng, L., Rutqvist, J., Birkholzer, J., and Liu, H.H., 2015, On the impact of temperatures up to 200°C in clay repositories with bentonite engineer barrier systems: A study with coupled thermal, hydrological, chemical, and mechanical modeling, Engineering Geology, 197, 278-295.

10.1016/j.enggeo.2015.08.026

Tunnel and Underground SpaceISSN:1225-1275(Print) 2287-1748(Online)한국암반공학회

Preview

Development of a Numerical Simulator for Swelling Behavior of the Buffer Material in the Engineered Barrier Induced by the Alteration under Hydro-Mechanical-Chemical Conditions

ABSTRACT

MAIN

(1)

Fig. 1.

A schematic diagram for illitization from smectite (Lee, 2024)

Fig. 2.

A coupling sequence for THMC coupled simulation in OGSPHREEQC-FLAC3D

Fig. 3.

Calculation process based on weighted arithmetic mean for node to element value from OGSPHREEQC to FLAC3D

(2)

(3)

Fig. 4.

Axisymmetic numerical model with the monitoring points

Table 1.

BBM input parameters for swelling test

(4)

Fig. 5.

Swelling pressure versus smectite mass fraction based on the previous experimental observation (Zheng et al., 2015)

Table 2.

Hydraulic input parameters for swelling test

(5)

(6)

(7)

(8)

Table 3.

Chemical input parameters and initial conditions for swelling test

Fig. 6.

Saturation and K+ change depending on time at each monitoring point when injecting K+0.1 mol

Fig. 7.

Montmorillonite and illite change depending on time and injection K+ ion concentration at each monitoring point

Table 4.

Summary of the simulation results after 20,000 days at the monitoring point O3 (top surface)

Fig. 8.

Total axial stress change and enlarged view depending on time

Acknowledgements

References

Saturation and $K^{+}$ change depending on time at each monitoring point when injecting $K^{+}$ 0.1 mol

Montmorillonite and illite change depending on time and injection $K^{+}$ ion concentration at each monitoring point