1. 서 론
2. DECOVALEX-2019 Task D
3. 수치해석 모델
3.1 구성방정식 및 열-수리-역학적 모델
3.2 초기조건 및 경계조건
3.3 입력물성
4. 수치해석 결과
4.1 히터파워 및 온도
4.2 상대습도
4.3 응력
4.4 해체 후 벤토나이트 블록 특성
5. 결 론
1. 서 론
고준위방사성폐기물을 안전하게 처분하기 위해 제안된 심층처분개념의 대표적인 예로는 스웨덴 SKB(Swedish Nuclear Fuel and Waste Management Company)가 제안한 KBS(Kärnbränslesäkerhet)를 들 수 있다. KBS는 다중방벽 시스템 개념을 기반으로 하여 지하 수 백 미터 깊이에 위치한 균열이 적고 안정한 암반에 처분터널과 처분공을 굴착하고 벤토나이트 완충재와 함께 처분용기를 처분공에 처분한 이후, 처분터널을 뒷채움재로 메워서 안전하게 격리하는 방식으로 수직처분방식과 수평처분방식이 있다(Fig. 1). 이러한 심층처분개념은 그 성능을 직접 평가해 볼 수 있는 장점이 있기 때문에 전 세계적으로 장기간에 걸쳐 고준위방사성폐기물을 안전하게 관리할 수 있는 방법이라는 공감대가 형성되었으며(Berg and Brennecke, 2013), 많은 국가에서 고준위방사성폐기물 처분을 위해 심층처분개념에 근거하여 자국의 현실에 맞는 처분 개념을 도출하고, 그 처분 개념의 안전성과 성능을 평가하고자 처분 환경과 유사한 심부 암반에 지하연구시설(Underground Research Laboratory, URL)을 건설하여 다양한 실험들을 수행하고 있다.
심층처분개념을 기반으로 한 처분시스템에서는 방사성 붕괴열에 의한 온도변화가 예상되며, 처분 시스템으로 유입되는 지하수에 의해 벤토나이트 완충재의 수리적 특성 변화뿐만 아니라 열전도도와 같은 열적 물성의 변화로 인한 열적 거동의 변화 역시 예상된다. 그리고 온도변화에 따른 열응력과 벤토나이트 완충재의 포화도 변화에 따른 팽윤압은 처분 시스템 및 처분장 주변의 응력 상태를 변화시키게 된다. 변화된 응력은 벤토나이트 완충재 및 암반의 투수계수와 흡입력에 영향을 주게 되어 수리적인 거동에 영향을 미치게 된다(Fig. 2). 고준위방사성폐기물 처분시스템에서의 이러한 열-수리-역학적 복합거동(coupled thermo-hydro-mechanical (THM) behavior) 변화는 장·단기 안전성 및 안정성에 영향을 미칠 수 있기 때문에 고준위방사성폐기물 처분장 부지선정, 설계, 그리고 운영을 위해서 THM 복합거동 변화에 대한 이해와 이를 분석하고 예측할 수 있는 수치해석 코드 및 모델이 반드시 필요하다.
처분시스템에서 예상되는 복합거동을 효과적으로 분석하고 예측할 수 있는 수치모델과 해석기법을 개발하고 이를 검증하기 위해 1992년부터 국제공동연구 DECOVALEX(DEvelopment of COupled models and their VALidation against EXperiments, http://www.decovalex.org) 프로젝트가 진행되고 있으며, 최근에는 2016년부터 2019년까지 DECOVALEX-2019가 진행되었다. DECOVALEX-2019에서 수행된 Task D에서는 스위스의 결정질 암반에 위치한 GTS(Grimsel Test Site)에서 수행된 실규모 현장시험(Full-scale Engineered Barriers EXperiment, FEBEX)(ENRESA, 2000)에 대한 수치해석을 수행하고 각 기관별로 제안된 다양한 THM 모델과 수치해석 결과를 상호 비교․검증하고자 하였다. 본 연구에서는 DECOVALEX-2019 Task D에 대해 간단히 소개하고, 한국원자력연구원이 FEBEX 모델링을 위해 사용한 시뮬레이터와 THM 모델을 살펴보고 이를 이용하여 수행된 모델링 결과를 살펴보고자 한다.
2. DECOVALEX-2019 Task D
다중방벽 개념의 처분장에서 완충재는 1)처분장 주변 암반에 존재하는 지하수가 처분장 내로 침투하여 처분용기와 접촉하는 것을 제한하고, 2)처분용기가 기능을 상실하여 침출된 방사성핵종들이 처분장 주변의 암반으로 유출되는 것을 지연시키며, 3)폐기물로부터 발생하는 붕괴열이 주변 암반으로 잘 발산되도록 하여 처분시스템의 온도가 지나치게 상승하는 것을 방지하고, 4)지진과 같은 외부 충격을 완화시켜 처분용기를 보호하도록 설계되기 때문에 심층처분장의 안전성을 확보하는데 매우 중요한 요소이다(Cho, 2017). 따라서 Task D에서는 벤토나이트 완충재 블록에서 예상되는 THM 복합거동 특성을 분석하고 이를 예측할 수 있는 수치모델과 해석기법을 개발하고자 하였으며, 한국원자력연구원, 체코 IGN(Institute of Geonics of Czech Academy of Sciences), 일본 JAEA(Japan Atomic Energy Agency), 그리고 대만 NCU(National Central University of Taiwan)가 참여하였다. 개발된 수치모델과 해석기법의 검증을 위해 Task D에서 선정한 FEBEX(Fig. 3)는 스페인 ENRESA(Empresa National de Residoos Radiactivds, S.A)가 제안한 화강암 기반의 기준 처분 시스템을 실증하고, 처분시스템에서의 THM 복합거동 특성을 살펴보고자 수행된 실험이다. FEBEX 현장시험을 위해 TBM으로 길이 70.39 m, 직경 2.28 m의 터널을 1995년 9월 25일부터 굴착하기 시작하여 1995년 10월 30일에 완공하였다. 그리고 1996년 7월 1일부터 붕괴열을 모사하기 위한 두 개의 히터와 완충재 블록, 그리고 다양한 센서들을 설치하고 2.7 m 두께의 콘크리트 플러그를 설치하여 시험 구간을 폐쇄하였다. 두 히터를 가동한 1997년 2월 27일로부터 61일이 지난 이후에 히터 온도는 100°C에 도달하였으며, 이후 히터의 온도를 일정하게 유지하였다. 2002년 2월 28일에는 히터 1의 전원을 끄고, 콘크리트 플러그, 벤토나이트 블록 및 히터 1을 제거하기 시작하여 2002년 7월 17일에 해체 작업이 완료되었다.
Fig. 3.
General layout of the FEBEX in situ test and selected sections for comparison of modelling results and measured temperature, relative humidity, and total stress (modified afterBárcena et al., 2003)
FEBEX 현장시험은 2000년부터 2003년까지 수행된 DECOVALEX-III의 Task 1로 선정된 바 있다. DECOVALEX-III에서는 1000일간 계측된 현장 데이터를 기반으로 해석모델과 시뮬레이터를 검증하고자 하였지만, DECOVALEX-2019 Task D에서는 18.5년 동안 계측된 벤토나이트 블록에서의 온도 및 상대습도 그리고 벤토나이트 블록과 암반 경계부에서의 응력 데이터를 기반으로 검증을 하였다. 뿐만 아니라 두 번에 걸쳐 진행된 해체 작업에서 확보한 벤토나이트 블록의 포화도, 함수율, 그리고 건조밀도 데이터와 계산 값을 비교하여 사용된 THM 모델과 시뮬레이터를 다양한 방법으로 검증하고자 하였다. 특히, 벤토나이트 완충재의 THM 물성은 건조밀도에 따라 크게 다르기 때문에 벤토나이트 완충재에서의 THM 복합거동 특성을 잘 예측하기 위해서는 건조밀도 변화를 잘 예측해야 하므로, DECOVALEX-2019 Task D에서는 완충재의 건조밀도 변화 예측 모델과 기법을 개발하고 검증하는 것에 중점을 두었다. Task D에서는 해석모델과 기법의 검증을 위해 많은 섹션에서의 현장 데이터를 활용하였지만, 본 연구에서는 열적 거동을 비교하기 위해 두 히터의 끝단인 D1과 I 섹션에서 계측된 온도 데이터를 활용하였으며, 수리적 거동을 비교하기 위해 상대적으로 온도 변화가 크지 않은 C 섹션, 히터 1이 위치하고 있는 E1섹션, 두 히터 사이의 H 섹션, 그리고 히터 2가 위치하고 있는 F2 섹션에서 계측된 상대습도 데이터를 이용하였고, 역학적 거동 비교를 위해 B2와 E2섹션에서 측정된 터널 축방향 및 횡방향 응력 데이터를 활용하였다(Fig. 3). 해체작업이 진행된 이후 현장시험에서 관측된 포화도, 함수율, 그리고 건조밀도 데이터와 계산 값의 비교를 위해 본 연구에서 선정한 지점은 Fig. 4에 나타나있다.
Fig. 4.
Location of sampling sections after dismantling operations and four selected sections, 15, 31, 43, and 56, for comparison of modelling results and measured data (modified afterVillar et al., 2020)
3. 수치해석 모델
3.1 구성방정식 및 열-수리-역학적 모델
FEBEX 현장시험을 모델링하고 벤토나이트 완충재 블록에서의 THM 복합거동 특성을 평가하기 위해 TOUGH2-MP/FLAC3D를 이용하였다. TOUGH2-MP/FLAC3D는 TOUGH-FLAC(Rutqvist et al., 2002)과 동일한 알고리즘으로 개발된 해석 시뮬레이터로써 TOUGH2(Pruess et al., 1999)의 병렬해석 버전인 TOUGH2-MP(Zhang et al., 2008)와 FLAC3D(Itasca, 2012)를 연동하여 THM 복합거동을 계산하도록 개발되었다(Fig. 5). 먼저, 다공성 매질에서의 열-수리 거동을 TOUGH2-MP에서 계산하고, 그 결과로 계산된 온도, 공극수압, 흡입력, 그리고 포화도의 정보를 FLAC3D로 넘겨준다. 이후, TOUGH2-MP에서 전달된 정보들을 이용하여 열응력, 공극수압, 그리고 팽윤압의 변화에 따른 역학 해석을 수행하고 변화된 평균 유효응력, 공극률, 그리고 건조밀도를 TOUGH2-MP에 전달하여 투수계수, 흡입력, 공극률의 정보를 업데이트한 후 다음 스텝으로 넘어간다(Fig. 6). 이러한 일련의 과정을 반복하면서 FEBEX 현장시험에 대한 수치해석을 수행하고 THM 복합거동 특성을 살펴보았다.
Fig. 6.
Numerical procedure of a linked TOUGH2-MP and FLAC3D simulation with explicit sequential solutions (modified afterRutqvist et al., 2002)
3.1.1 TOUGH2-MP
다공성 매질에 대한 다성분(multi-components)의 이상 유동(two-phase flow)을 해석하기 위해 미국 LBNL(Lawrence Berkeley National Laboratory)에서 개발한 TOUGH2-MP에서는 열전달을 위해 전도와 이류가 고려되며 사용되는 에너지 보존 방정식은 다음과 같다.
| $$\frac d{dt}((1-\phi)\rho_RC_RT+\phi\sum_\beta^{}S_\beta\rho_\beta u_\beta)+\nabla(-\lambda\nabla T+\sum_\beta^{}h_\beta F_\beta)=Q$$ | (1) |
여기서, β는 유체의 상, ∅는 공극률, ρR과 CR는 각각 입자의 밀도와 비열이고, T와 t는 각각 온도와 시간이다. Sβ는 포화도(saturation of phase β), ρβ는 밀도(density of phase β), uβ는 내부 에너지(specific internal energy in phase β), λ는 열전도도, hβ는 엔탈피(specific enthalpy in phase β), Fβ는 유체 플럭스(flux of phase β), 그리고 Q는 유입 또는 유출된 에너지(energy source or sink term)이다.
본 연구에서는 에너지 보존 방정식의 열전도도(λ)가 식 (2)와 같이 포화도(Sl)의 함수로 계산되고, 비열(CR)은 식 (3)과 같이 온도의 함수로 계산되도록 TOUGH-MP를 수정하였다.
| $$\lambda=\lambda_{sat}+\frac{\lambda_{dry}-\lambda_{sat}}{1+e^{(S_l\;-x_0\;)/d_x}}$$ | (2) |
여기서, λsat와 λdry는 각각 포화 및 건조 상태에서의 열전도도를 의미하며, x0는 λsat와 λdry의 평균값일 때의 포화도를 의미하고, dx는 상수이다.
| $$C_R=c_1T+c_2$$ | (3) |
여기서, c1은 상수이고 c2는 0°C에서의 비열이다.
TOUGH2-MP에서 사용되는 유체의 질량 보존 방정식은 식 (4)와 같고, Darcy 속도(Darcy velocity in phase β, uβ)는 식 (5)로 계산된다.
여기서, ∅는 공극률, 는 질량비(mass fraction of component κ(air and water) present in phase β), 는 확산 텐서(hydrodynamic dispersion tensor of component κ in phase β), 그리고 Qβ는 유체의 유입 또는 유출(fluid source or the sink term in phase β)을 의미한다.
| $${\boldsymbol u}_\beta=-k\frac{k_{r\beta}}{\mu_\beta}(\nabla P_\beta-\rho_\beta g)$$ | (5) |
여기서, k는 투수계수, krβ, µβ, 그리고 Pβ는 각각 β상에서의 상대투수계수, 점성, 그리고 압력을 의미하고 g는 중력 가속도이다.
벤토나이트 완충재의 포화 수리전도도(Ksat)는 건조밀도가 증가함에 따라 감소하는 경향을 보이기 때문에(Fig. 7), 식 (6)와 (7)을 이용하여 수치해석에서 반영하였다. 건조밀도는 FLAC3D의 역학해석에서 계산된 비체적(specific volume, υ)과 공극률의 상관관계(식 (8))와 식 (9)을 이용하여 계산하였고, 계산된 건조밀도는 다음 스텝의 열-수리해석에서 업데이트 되도록 TOUGH2-MP 코드를 수정하였다(Fig. 6).
| $$k_{sat}=(K_{sat}\times\mu)/(\rho_w\times g)$$ | (6) |
| $$K_{sat}=4.377\times10^{-10}\times\rho_d^{-21.92}$$ | (7) |
여기서, µ와 ρw는 각각 물의 점도와 밀도이며, g는 중력가속도이다.
| $$\phi=(v-1)/v$$ | (8) |
| $$\rho_d=\rho_R(1-\phi)$$ | (9) |
여기서, ρR은 ENRESA(2000)에 보고된 2.7을 이용하였다.
Fig. 7.
Hydraulic conductivity used in numerical model and hydraulic conductivity of trimmed and remoulded samples presented in the previous study (Villar, 2017)
암반에 대한 유체 상의 상대 투수계수(relative permeability)는 식 (10)으로 계산하였으며, 나머지 매질의 유체 상대 투수계수는 식 (11)에 의해 결정되는 것으로 가정하였고, 유효 포화도(S*)는 식 (12)로 계산하였다. 기체 상의 상대 투수계수는 식 (13)으로 계산하였다. 또한 흡입력(s)는 식 (14)에 나타나 있는 van Genuchten 모델을 따르는 것으로 가정하였다.
| $$k_{rl}=\sqrt{S^\ast}\left\{1-(1-\lbrack S^\ast\rbrack^{1/\lambda_{van}})^{\lambda_{van}}\right\}^2$$ | (10) |
여기서, λvan은 van Genuchten 파라미터이다.
| $$k_{rl}=(S^\ast)^n$$ | (11) |
여기서, krl은 유체에 대한 상대 투수계수이고, n은 상수이다.
| $$S ^{*} =(S _{l} -S _{lr} )/(1.0-S _{lr} )$$ | (12) |
여기서, Slr은 잔류 포화도(residual saturation)를 의미한다.
| $$k_{rg}=1-k_{rl}$$ | (13) |
여기서, krg은 기체의 상대 투수계수이다.
| $$s=-P_0(\;\lbrack S^\ast\rbrack^{1/\lambda_{van}})^{1-\lambda_{van}}$$ | (14) |
여기서, λvan과 P0는 van Genuchten 상수이다.
벤토나이트 완충재에서 기체 투수계수(kgas)는 물의 투수계수(kliq)에 비해 월등히 크기 때문에(ENRESA, 2000) 이를 반영하기 위해 본 연구에서는 식 (15)를 이용하였다.
| $$k_{gas}=k_{liq}(1+b/P)$$ | (15) |
여기서, b와 P는 각각 Klinkenberg 파라미터와 압력이다.
확산에 의한 유체 플럭스(diffusive flux of component κ in phase β, )은 식 (16)을 이용하여 계산하였으며, 유선(flow path)의 비틀림도(tortuosity)와 관련된 상수(τβ)는 식 (17)로 정의했고, 확산계수()는 온도와 압력에 따라 달라지기 때문에(Vargaftik, 1975; Walker et al., 1981) 식 (18)을 이용하여 계산하였다.
| $$f_\beta^\kappa=-\phi\tau_0\tau_\beta\rho_\beta d_\beta^\kappa\nabla X_\beta^\kappa$$ | (16) |
여기서, κ와 β는 각각 성분과 상을 의미하며, ∅는 공극률, τ0는 매질의 비틀림도, ρβ은 β상에서의 밀도, 그리고 은 β상에서 κ의 질량분율을 의미한다.
| $$\tau_\beta=k_{r\beta}(S_\beta)$$ | (17) |
여기서, krβ는 β상에서의 상대투수계수를 의미하고 Sβ는 β상의 포화도를 의미한다.
| $$d_\beta^\kappa(P,T)=d_\beta^\kappa(P_0,T_0)\frac{P_0}P\lbrack\frac{T+273.15}{273.15}\rbrack^{1.8}$$ | (18) |
여기서, (P0, T0)는 대기압(1atm)과 0°C에서의 공기와 수증기 혼합물에 대한 확산계수 2.13×10-5 m2/s이다.
현장시험에서 측정된 상대습도(Rh) 데이터와의 비교를 위해 본 연구에서는 TOUGH2-MP에서 계산된 포화도를 Kelvin’s equation(식 (19))을 이용하여 상대습도로 변환하였다.
| $$Rh=\exp(\frac{-sm_vv_w}{R(T+273.16)})$$ | (19) |
여기서, s는 흡입력, mv는 물의 비체적(0.001 m3/kg), vw는 수증기의 몰 질량(18.016×10-3 kg/mol), 그리고 R은 이상 기체 상수(ideal gas constant, 8.31432 J/mol/k)이다.
3.1.2 FLAC3D
역학적 거동 해석을 위해 식 (20)에 나타나 있는 모멘텀 밸런스 방정식이 이용되었다.
| $$\nabla\;\bullet\;\boldsymbol\sigma+\boldsymbol b=0$$ | (20) |
여기서, σ는 응력 텐서 그리고 b는 물체력 벡터를 의미한다.
유효응력은 전응력(σ)과 TOUGH2-MP에서 계산된 공극수압(Pp)과 Biot 상수(α)를 이용하여 식 (21)로 계산된다.
| $$\sigma'=\sigma-\alpha P_p$$ | (21) |
벤토나이트 블록에 대한 역학적 해석을 위해서는 BBM(Barcelona Basic Model)(Alonso et al., 1990)을 사용하였으며, 나머지 물질에 대해서는 탄성 모델을 사용하였다. FLAC3D에서는 BBM을 지원하지 않기 때문에 본 연구에서는 Lee et al.(2019)에 기술되어 있는 BBM을 토대로 개발된 해석모듈(Lee et al., 2020a)을 사용하였다.
3.2 초기조건 및 경계조건
수치해석에서는 높이 50 m, 폭 100 m, 그리고 회전각 10°를 갖는 3차원 축대칭 모델(axis-symmetric model)을 이용하였다(Fig. 8). 암반의 초기 온도는 12°C이고 완전히 포화되었다고 가정하였으며, 초기 공극수압은 0.7 MPa, 초기 응력은 28 MPa의 등방 응력조건으로 설정하였다. Fig. 8에 나타나 있는 붉은 점선에서는 온도와 압력을 고정하였고, 푸른 점선에서는 열과 유체의 이동은 없으며, 붉은 점선과 푸른 점선에서 법선 방향 변위는 모두 고정되었다. 초기조건과 경계조건이 부여된 암반을 이용하여, 시험 터널 굴착, 현장시험 장치 설치, 운영, 그리고 해체를 Table 1과 같이 단순화하여 수치해석을 수행하였다. 본 연구에서는 히터가 발열을 시작하는 시점을 기준일(Day 0)로 설정하였고, 시험 터널은 Day –385에 순간적으로 모두 굴착되는 것으로 가정하였으며 굴착 이후 터널 벽면의 온도는 12°C로, 공극수압은 대기압으로 가정하여 0.1 MPa로 고정하였고, 역학 해석에서는 자유면으로 간주하였다. 시험터널 굴착 이후 250일간의 환기작업을 거쳐 Day –135에 벤토나이트 블록, 히터와 라이너(liner) 그리고 콘크리트 플러그를 동시에 설치하여 시험구간을 폐쇄하였다고 가정하였다. 설치된 시험 장치와 플러그는 모두 12°C의 온도와 대기압 상태로 가정하였다. 그리고 압밀된 벤토나이트 블록은 초기 응력상태를 Gens et al.(2009)와 Sánchez et al.(2012)에서 사용한 값과 동일하게 0.2MPa로 설정하고 나머지는 모두 0.1MPa로 가정하였다. 벤토나이트 블록의 초기 포화도는 0.55, 그리고 히터, 라이너, 그리고 플러그의 초기 포화도는 0.99로 가정하고 해석을 수행하였다.
Table 1.
Modeling procedure for FEBEX
대규모 현장 히터시험에서는 히터와 많이 이격된 곳에서 전력을 공급하기 때문에 열손실이 발생할 수밖에 없으며, 스위스 Mont Terri URL의 HE-D 현장시험과 벨기에 HADES URL의 ATLAS 현장시험에서 추정된 열손실은 약 5%로 보고되었으므로(Garitte et al., 2014), FEBEX 현장시험에서도 두 현장시험과 유사하게 약 5%의 열손실이 있었다고 가정하였다. 시험 구간을 폐쇄하고 Day 0일 때, Fig. 9의 검은 점선으로 표시된 부분에 각각 1,200W의 95%인 1140W의 파워를 입력하였고, Day 20에는 2,000W의 95%인 1900W의 파워를 주입하였으며, Day 53부터 검은 점선의 온도를 100°C로 고정하였다. Day 1,827에서는 히터 1의 온도를 고정하지 않으므로 히터 전원을 끈 것을 모사하였으며, Day 1,966에서 해체 구간에 위치한 벤토나이트 블록의 건조밀도, 함수율, 그리고 포화도를 계산했다. Day 1,974에서는 해체 구간에 있는 모든 요소를 초기온도인 12°C와 대기압으로 설정하고, 역학해석에서는 자유면으로 설정하여 해체를 모사하였으며, 1차 해체 이후에 설치된 플러그는 12°C의 초기 온도, 대기압 그리고 0.1MPa의 초기 응력으로 설정하였다. 이후 Day 6,630에서 히터 2를 끄고, Day 6,758에서 벤토나이트 블록의 건조밀도, 함수율, 그리고 포화도를 계산하였다.
Fig. 9.
Location of the applied power and constant temperature at the heaters in the numerical simulations (Lee et al., 2020b)
3.3 입력물성
수치해석을 위해 사용된 입력물성은 Table 2에 나타나 있으며, 벤토나이트 블록에 대한 BBM 파라미터는 Table 3에 정리하였다. Table 2와 Table 3에 알파벳은 참고문헌을 의미하며 *, **, ***를 포함하고 있는 입력값은 각각 가정 값, 계산 값, 그리고 보정 값을 의미한다.
Table 2.
Input parameters
Note: aTaken from Gens (2017). bTaken from Gens (2018a). cTaken from ENRESA (2000). dTaken from Gens et al. (2009). *Assumed values. **Calculated value by equation(s). ***Calibrated value.
Table 3.
BBM parameters of bentonite block
Note: aTaken from Gens et al. (2009). bTaken from Sánchez et al. (2012).
4. 수치해석 결과
사용된 THM 모델과 해석 시뮬레이터의 검증과 현장 적용성을 검토하고자 현장시험 데이터와 해석 결과를 비교․분석 하였다. 현장시험 데이터는 Papafotiou et al.(2017), Gens(2018b), 그리고 Gens(2019)에서 보고된 값을 이용하였다.
4.1 히터파워 및 온도
4.1.1 히터파워
벤토나이트 블록은 시간이 지남에 따라 주변 암반으로부터 유입되는 지하수로 인해 공극의 공기는 상대적으로 높은 열전도도를 갖는 물로 대체되기 때문에 포화가 진행됨에 따라 벤토나이트 블록의 열전도도는 증가하게 된다. 따라서, 시간이 지남에 따라 증가된 벤토나이트의 열전도도로 인해 히터가 일정한 온도를 유지하기 위해서는 더 큰 파워가 필요하기 때문에 히터 온도를 일정하게 유지하는 구간에서의 히터 파워는 서서히 증가하게 된다. 이러한 경향은 현장시험 데이터뿐만 아니라 수치해석 결과에서도 잘 나타나고 있다(Fig. 10). 하지만, 현장에서 히터 1의 파워가 약 5~15%정도 작은 반면, 수치해석에서 계산된 히터파워는 차이가 없는 것으로 나타났다. 이는 수치해석에서 시험구간 주변 암반은 균질하며 모든 벤토나이트 블록을 동시에 설치되었다고 가정했기 때문으로 판단된다. FEBEX 현장시험이 수행된 터널 주변 암반은 대부분 화강암이지만, 히터 1이 위치한 구역에서는 주변 화강암에 비해 약 15 ~ 20%정도 낮은 열전도도를 갖는 것으로 보고된(Keusen et al., 1989) 황반암(lamprophyre)이 존재한다(Fig. 3). 즉, 히터 1은 히터 2에 비해 상대적으로 낮은 열전도도를 갖는 암반에 위치하고 있기 때문에 히터의 온도를 100°C로 유지하는데 상대적으로 적은 히터 파워가 필요했을 것으로 보인다. 그리고 1996년 7월 1일부터 1996년 10월 15일까지 수행된 일련의 설치 과정에서 터널 안쪽부터 벤토나이트 블록이 설치되었기 때문에 히터 가동시점에서 터널 안쪽의 벤토나이트 블록의 포화도가 상대적으로 높았을 것이다. 즉, 히터 2 주변 벤토나이트 블록의 열전도도가 상대적으로 높기 때문에 동일하게 히터 온도를 100°C로 유지하기 위해서는 히터 2에 더 많은 파워를 주입하여야 했을 것이다. 이러한 두 가지의 이유로 현장에서는 히터 파워가 차이를 보였지만 모델링에서는 차이를 보이지 않았으며, 추후 황반암과 설치 프로세스를 반영한 모델링을 모두 수행하고 어떠한 영향이 더 큰지도 조사해 볼 필요가 있을 것으로 판단된다.
4.1.2 온도
히터에서 상대적으로 가까운 지점에서 계측된 온도와 계산된 온도는 전반적으로 유사한 거동을 보였지만(Fig. 11(a)), 계측된 온도가 다소 편차를 보인 반면, 수치해석에서 계산된 D1섹션과 I섹션에서의 온도는 큰 차이를 보이지 않았다. 벤토나이트 중간지점과 암반 근처에서의 계측된 온도를 살펴보면, 섹션 I가 두 히터 사이에 위치하기 때문에 히터로부터 동일한 이격거리에 위치했음에도 불구하고 섹션 I에서의 온도가 높게 나타났는데, 수치해석에서도 열 중첩의 영향이 잘 모사된 것으로 판단된다(Fig. 11(b)와 (c)). 또한, FEBEX 시험터널 중심으로부터 이격된 거리에 따른 온도분포 역시 계측값과 계산값이 유사하게 나타난 것을 알 수 있으며, 섹션 I에서 전반적으로 온도가 높게 나타나고 있음을 알 수 있었다(Fig. 11(d)). 시간에 따른 히터 파워와 벤토나이트에서 온도 변화 그리고 온도분포가 전반적으로 유사한 경향을 보였기 때문에, 본 연구에서 사용된 TOUGH2-MP의 열전도도 모델(식 (2))이 벤토나이트 블록에서의 열적 거동을 모사하기에 적절한 것으로 판단된다.
4.2 상대습도
섹션 C, E1, H, 그리고 F2에서 상대습도 변화와 시험터널 중심으로부터 이격된 거리에 따른 섹션 F2에서의 상대습도 분포를 Fig. 12에 나타내었다. 암반 근처에서 벤토나이트 블록의 상대습도는 암반으로부터의 지하수 유입으로 빠르게 증가하는 것으로 나타났고, 지하수가 서서히 안쪽으로 이동하여 벤토나이트 블록의 중간 지점에서도 상대습도가 증가하는 것으로 나타났다(Fig. 12(a)와 (b)). 수치해석에서는 섹션별로 큰 차이를 보이지 않는 것으로 나타났고, 전반적으로 지하수 유입에 따른 상대습도의 변화를 수치모델이 잘 재현한 것으로 나타났다.
히터와 인접한 벤토나이트 완충재에서는 Fig. 13에 나타난 것처럼 온도 증가로 인해 물이 수증기로 변하고 가스의 압력구배에 의해 수증기는 바깥쪽으로 이동하게 된다. 이동된 수증기는 온도 감소로 인해 물로 변화되어 히터와 가까운 지점에서의 상대습도는 순간적으로 증가하는 것으로 나타났고, 이후 지속적인 온도의 상승으로 히터와 가까운 지점에서의 상대습도는 감소하다가 지하수의 유입으로 인해 상대습도는 서서히 증가하는 것으로 나타났다. 수치해석에서도 히터 가동 이후의 순간적인 상대습도 증가가 나타났고 상대적으로 온도가 더 높은 섹션 E1에서 더 크게 발생한 건조현상도 잘 모사되었다(Fig. 12(c)). 그리고 섹션 F2에서의 상대습도 분포 변화에서도 히터 부근에서의 건조현상과 시간이 지남에 따라 전반적인 상대습도 증가 경향이 잘 나타났다(Fig. 12(d)). 하지만, 섹션 H에서는 상대적으로 낮은 증가를 보였으며, 섹션 E1에서는 히터 가동 초기에 상대습도 증가 이후 감소되는 속도가 느린 것으로 나타났고, 500일 이후 상대습도 증가가 계측된 값에 비해 빠르게 진행되는 것으로 계산되었다. 히터 부근에서의 상대습도 변화를 조금 더 현실적으로 재현하기 위해서는 사용된 상대투수계수 모델 또는 확산 모델의 수정이 필요할 것으로 판단되며, 다양한 모델을 이용하여 최적의 수리모델을 도출할 필요가 있을 것으로 보인다.
Fig. 13.
Scheme of the coupled thermo-hydro-mechanical and chemical (THMC) processes in the high-level radioactive waste disposal system (Lee et al., 2020c)
4.3 응력
사용된 BBM이 벤토나이트 블록의 역학적인 거동을 모사하기에 적합한지를 판단하기 위해 섹션 B2의 벤토나이트 블록 중간지점에서 계측된 터널 축방향 응력 데이터와 섹션 E2의 암반과 인접한 곳에서 계측된 반경방향 응력 데이터를 이용하여 계산된 값과 비교하였다. 벤토나이트 블록은 현장에 설치된 이후, 지하수 유입에 따른 팽윤압의 발생과 온도 변화에 따른 열응력의 발생으로 축방향과 반경방향의 응력이 증가하는 것으로 나타났으며, 수치해석 상에서도 이러한 경향이 잘 재현되는 것으로 나타났다(Fig. 14). 하지만 수치해석에서 축방향과 반경방향의 응력 모두 상대적으로 빠르게 증가하는 경향을 보였는데, 이는 현장 조건과 달리 수치해석에서는 암반과 벤토나이트가 완벽하게 일치하는 조건을 가정하였기 때문인 것으로 판단된다. 실제 현장에서는 Fig. 15에서처럼 암반과 벤토나이트 블록 사이에 작은 틈이 존재하기 때문에 지하수의 유입으로 인한 벤토나이트의 팽윤으로 그 틈을 메우는데 시간이 걸리지만, 수치해석에서는 완벽하게 접촉하고 있는 조건을 가정하였기에 상대적으로 응력 증가가 빠르게 나타난 것으로 판단된다.
4.4 해체 후 벤토나이트 블록 특성
FEBEX 현장시험은 약 5년간의 운영 이후 히터 1이 위치한 구간을 해체하였으며, 이후 약 13.5년간 더 운영을 하고 18.5년 이후 FEBEX 현장시험을 완전히 해체하였다. 두 차례의 해체 작업을 수행하면서 완충재 블록에서 샘플을 채취하여 포화도, 함수율, 그리고 건조 밀도를 측정하여 분포 특성을 살펴보았다. 본 연구에서는 사용된 TOUGH2-MP/FLAC3D 시뮬레이터와 수치해석 모델이 최종 단계에서의 포화도, 함수율, 그리고 건조밀도의 분포특성을 잘 모사할 수 있는지를 살펴보고자 Fig. 4에 나타나 있는 1차 해체 작업 구역의 섹션 15와 31 그리고 2차 해체 작업 구역의 43과 56에서 획득한 샘플의 결과와 비교하였다.
4.4.1 포화도 및 함수율
해체 이후 측정된 포화도와 함수율, 그리고 수치해석으로부터 계산된 값을 FEBEX 시험터널의 중심으로부터 이격된 거리와 함께 Fig. 16과 Fig. 17에 나타내었다. 첫 번째 해체 이후 측정된 벤토나이트 블록의 포화도와 함수율의 변화를 살펴보면, 섹션 15에서는 암반으로부터의 지하수 유입으로 인해 전 구간에서 초기값보다 증가한 것으로 나타났고, 히터가 위치한 섹션 31 주변 암반부근에서는 섹션 15와 유사하게 증가하는 것으로 나타났지만 히터 부근에서는 온도증가에 따른 기화로 인해 섹션 15에 비해 상대적으로 크게 변화지 않는 것으로 조사되었다(Fig. 16(a)와 Fig. 17(a)). 수치해석으로 계산된 포화도와 함수율은 암반부근 뿐만 아니라, 섹션 31의 히터부근에서 역시 유사한 것으로 나타났지만, 섹션 15의 터널 중심부에서 상대적으로 낮게 나타났다.
2차 해체 작업에서 측정된 포화도와 함수율은 히터가 위치한 섹션 43과 위치하지 않은 섹션 56 모두 첫 번째 해체작업에서 측정된 포화도와 함수율 보다 높게 나타났으며, 히터 부근의 벤토나이트의 경우 포화상태로 변화된 것으로 조사되었다(Fig. 16(b)와 Fig. 17(b)). 수치해석에서도 첫 번째 해체 때보다 포화도와 함수율이 증가했지만 전반적으로 계측값보다 낮게 나타났다. 첫 번째 해체와 최종 해체에서 계측된 값이 계산값과 차이를 보이는 것은 수치해석에서 사용된 투수계수가 작았거나, 상대투수계수 모델 또는 흡입력 모델에 사용된 상수들이 벤토나이트 블록의 수리적 거동을 모사하기에 부족했을 것으로 판단되며, 추후 모델에 사용된 파라미터에 대한 보정이 필요할 것으로 보인다.
4.4.2 건조밀도
해체 이후 측정된 건조밀도와 수치해석으로 계산한 값을 터널중심으로부터 이격된 거리에 따라 초기 건조밀도와 함께 나타내었다(Fig. 18). 암반과 인접한 벤토나이트 블록은 설치 이후 암반으로부터의 지하수 유입으로 인해 상대습도와 포화도가 증가하는 것으로 나타났다(Fig. 12(a), Fig. 17). 이러한 지하수의 유입으로 벤토나이트 완충재는 팽윤하게 되어 부피가 증가하여 건조밀도는 감소하게 된다. 반면, 암반과 인접한 벤토나이트 블록의 팽윤으로 인해 터널 중심부와 히터부근에서는 벤토나이트 블록이 압밀되어 건조밀도가 증가하는 것으로 조사되었다. 그리고 포화도와 함수율의 경향과는 달리 첫 번째 해체 작업이 이루어지는 동안 측정된 건조밀도는 18년 이후의 최종 해체 작업에서 측정된 건조밀도와 큰 차이를 보이지 않았다.
현장에서 관찰된 건조밀도의 분포뿐만 아니라, 두 해체 시점에서 건조밀도 분포가 유사한 점 역시 BBM을 이용한 수치해석에서 잘 재현되었지만, 계측된 건조밀도보다 조금 크게 계산이 된 것으로 보인다. 터널 중심과 히터부근에서의 계산된 포화도와 함수율이 상대적으로 작게 나타난 반면, 건조밀도는 크게 나타난 것으로 보아 본 연구에서 사용한 건조밀도에 따른 투수계수 변화(식 (5))가 실제 완충재의 특성과 차이를 보이기 때문으로 판단된다. 다시 말해, 수치해석에서 건조밀도가 증가함에 따라 투수계수의 감소정도가 실제 완충재 블록에서의 감소보다 더 크게 감소하여 포화도와 함수율이 작게 나타난 것으로 판단된다. 또 내부에서의 포화도와 함수율이 작게 나타나 팽윤이 상대적으로 덜 발생하여 건조밀도가 크게 나타나게 된 것으로 판단된다. 식 (5)를 수정하여 건조밀도 증가에 따른 투수계수의 감소를 완화한다면, 지하수의 유입이 상대적으로 빠르게 진행되어 계측된 포화도와 함수율이 좀 더 유사하게 나타날 것이고 터널 중심과 히터 부근에서의 팽윤이 조금 더 빠르게 진행되어 건조밀도가 조금 감소할 것으로 판단된다.
5. 결 론
본 연구에서는 TOUGH2-MP/FLAC3D 시뮬레이터와 BBM을 이용하여 스위스 GTS에서 18년 이상 수행된 FEBEX 현장시험에 대한 수치해석을 수행하고 벤토나이트 완충재 블록에서의 THM 복합거동을 살펴보았다. 그리고 사용된 THM 모델의 적합성과 BBM의 현장 적용성을 살펴보기 위해 계측된 히터 파워, 온도, 상대습도, 응력 데이터와 해체 이후 조사된 완충재 블록의 포화도, 함수율, 그리고 건조밀도 분포 데이터를 수치해석에서 계산된 값과 비교·분석하였다. 수치해석 결과, 시간에 따른 히터 파워와 온도 변화는 전반적으로 잘 재현되었지만, 히터 1과 히터 2에서의 파워 차이를 계산할 수는 없었으며, 이를 개선하기 위해서는 FEBEX 터널 주변에 분포하는 황반암과 3개월 이상 진행된 시험장치 및 벤토나이트 블록의 설치 공정을 고려할 필요가 있을 것으로 판단된다. 상대습도 변화 특성과 분포 역시 전반적으로 잘 모사되었으나, 수치해석에서 히터 부근에서의 포화과정이 상대적으로 빠르게 진행된 것으로 보아 상대투수계수 모델 또는 확산 모델에 일부 수정이 필요할 것으로 보인다. 열응력 및 팽윤압의 발생으로 인한 응력변화는 현장조건과 수치해석에서의 가정조건 차이로 운영 초기에 다소 차이를 보였지만 전반적으로 잘 계산되는 것으로 나타났다. 해체 이후 측정한 포화도, 함수율, 그리고 건조밀도의 분포 역시 전반적으로 잘 재현되었지만, 건조밀도가 터널 중심과 히터부근에서 조금 크게 계산되었기 때문에, 벤토나이트 블록의 투수계수가 상대적으로 작은 값으로 반영되어 포화도와 함수율이 작게 계산된 것으로 보이며, 이를 개선하기 위해서는 건조밀도에 따른 투수계수 모델에 수정이 필요할 것으로 판단된다. 현장 데이터와 계산값을 비교·분석하여 도출된 THM 모델의 수정 및 보완사항을 토대로 개선하고, 현장시험 조건을 최대한 반영하여 수치해석을 수행한다면 더 좋은 해석결과를 도출할 수 있을 것으로 판단된다. 뿐만 아니라, THM 모델에서 사용된 입력 파라미터의 온도 의존성을 검토하여 추가적으로 수치해석 모델을 수정할 필요가 있을 것으로 판단된다. 본 연구에서 사용한 THM 모델과 TOUGH2-MP/FLAC3D 시뮬레이터는 향후, 실험실 시험과 현장시험 설계 및 데이터 분석에 활용될 수 있을 뿐만 아니라 지금까지 개발된 한국형기준처분시스템(Korean reference disposal system)의 성능평가에도 활용될 수 있을 것으로 판단된다.
