1. 서 론

고준위방사성폐기물(high-level radioactive waste, HLW)을 안전하게 격리하기 위하여 스웨덴의 처분 전담기관(Swedish Nuclear Fuel and Waste Management Company, SKB)은 Fig. 1과 같은 심층처분방식(deep geological disposal concept)의 Kärnbränslesäkerhet(KBS)를 제안하였다. 지하 수백 미터 깊이의 암반에 처분 터널을 굴착하고 고준위방사성폐기물을 공학적방벽물질인 벤토나이트 완충재와 함께 처분공에 처분한 후, 뒷채움재(backfill materials)로 처분 터널을 완전히 메워 안전하게 격리하는 방식이다. 한국의 경우, 스웨덴의 심층처분 개념을 토대로 한국원자력연구원이 지하 500 m 심도에 고준위방사성폐기물을 직접 처분하는 한국형 기준 처분시스템(Korean Reference HLW disposal System, KRS)을 Fig. 2와 같이 제안하였다(Lee et al., 2007, Lee et al., 2008).

Fig. 1.

Principles of final disposal of spent nuclear fuel according to the KBS-3V method (SKB, 2010)

Fig. 2.

Concept of Korean reference HLW vertical disposal system, KRS-V1 (9, Lee et al., 2008)

이러한 심층처분 시스템에서는 고준위방사성폐기물로부터 발생하는 붕괴열에 의한 온도 변화와 주변 암반으로부터의 지하수 유입에 따른 완충재의 포화도 변화가 예상된다. 그 뿐만 아니라, 온도 및 포화도 변화에 따른 열응력 및 팽윤압의 발생으로 공학적방벽 및 주변 암반의 응력 상태가 변화될 것으로 예상된다. 또한, 변화된 응력은 벤토나이트와 암반의 수리적 물성을 변화시켜 처분시스템 전반의 수리적인 거동에 영향을 주며, 포화도에 따라 열전도도가 변화하는 벤토나이트 완충재의 특성으로 인해 열적 거동에도 영향을 미치게 되어 Fig. 3과 같은 복잡한 열-수리-역학적 복합거동(coupled thermo-hydro-mechanical(THM) behavior)이 처분시스템 전반에서 일어날 것으로 예상된다. 이러한 THM 복합거동의 변화는 고준위방사성폐기물 처분장의 안전성 및 안정성에 영향을 미칠 수 있기 때문에, 처분장 부지선정 및 처분시스템 설계를 위해서는 THM 복합 거동 변화에 대한 이해뿐만 아니라 이를 분석･평가하고 예측할 수 있는 THM 모델 및 수치해석 기법이 반드시 필요하다.

Fig. 3.

Mutual relationships between THM processes in a porous medium (CODE BRIGHT, 2004)

고준위방사성폐기물 처분장 설계 및 성능평가 관점에서 볼 때, 벤토나이트 완충재가 100°C 이상의 온도에 장기간 노출될 경우 완충재의 열적 변성(thermal alteration)으로 인해 완충재의 성능이 저하될 수 있기 때문에, 심층처분방식을 고려하고 있는 여러 나라들은 처분용기와 완충재의 접촉면에서의 최고 온도가 100°C를 넘지 않도록 규정하고 있다(JNC, 2000, SKB, 2009, Lee et al., 2007). 이러한 온도 기준을 만족하는 처분장을 설계하기 위해서는 벤토나이트 완충재의 열전도도가 중요한 설계인자이다. 또한, 벤토나이트 완충재의 열전도도는 포화도에 따라 크게 변화하는 특성이 있기 때문에, 벤토나이트 완충재의 열적 특성 뿐만 아니라, 수리적 특성인 투수계수와 수분보유곡선(water retention curve) 역시 중요한 설계 인자가 된다. 처분장 설계 및 처분장의 성능평가에 있어 중요한 인자인 열전도도, 투수계수, 그리고 흡입력(suction, s)은 벤토나이트 완충재의 건조밀도에 크게 영향을 받기 때문에(ENRESA, 2000, Tang et al., 2008), 벤토나이트 완충재에서의 복잡한 THM 복합거동을 이해하고 장기거동을 정확하게 예측하기 위해서는 열-수리적 복합거동에 따른 건조밀도 변화를 예측할 수 있는 역학적 모델이 반드시 필요하다.

이에, 본 연구에서는 불포화토의 역학적인 거동을 잘 모사하는 것으로 알려져 있는 Barcelona Basic Model(BBM)을 간단하게 소개하고, TOUGH2-MP/FLAC3D(Lee et al., 2016) 시뮬레이터의 BBM 해석 모듈 개발 방법에 대해 자세히 기술한 다음, 개발된 BBM 해석 모듈의 검증 과정을 상세히 기술하고자 한다.

2. BARCELONA BASIC MODEL

BBM에 대한 자세한 내용은 Lee et al.(2019)에 기술되어 있기 때문에 본 논문에서는 간단하게 언급하고자 한다. BBM은 불포화토의 역학적인 거동 해석을 위해 Modified Cam Clay Model(MCCM)을 확장하여 탄-소성 모델(elasto-plastic model)로 처음 제안되었다(Alonso et al., 1990). 포화 상태(s =0)에서는 MCCM과 동일하지만, 불포화 상태(s ≠0)에서는 항복면이 Fig. 4와 같이 증가하는 모델이다. BBM은 정규 압밀 압축지수(compressibility parameter in virgin soil states, λ(0))가 흡입력에 대한 함수의 형태(λ(s))로 표현되는 특징이 있으며, 불포화토의 포화도(saturation)가 증가하여 흡입력이 감소할 경우 Fig. 5에서처럼 팽윤(swelling) 거동을 모사할 수 있다. Gens(1995)는 탄소성 모델인 BBM을 온도(T) 변화에 따른 항복면의 변화를 고려하여 열･탄소성 모델(thermo-elasto-plastic model)로 새롭게 제안하였다. 새롭게 제안된 BBM은 온도가 증가함에 따라 항복면이 작아지는 특징을 갖는 모델로써, 항복면을 삼차원 p―q―T 공간에 도시하면 Fig. 6과 같이 표현이 된다.

Fig. 4.

Three-dimensional representation of the yield surface in the p ‐ q ‐ s space (modified afterAlonso et al., 1990)

Fig. 5.

Compression curves for saturated and unsaturated soil (modified afterAlonso et al., 1990)

Fig. 6.

Three-dimensional representation of the yield surface in the p ‐ q ‐ T space (modified afterGens, 1995)

3. BARCELONA BASIC MODEL 해석 모듈 개발

본 연구에서는 Alonso et al.(1990)에 기술되어 있는 흡입력 변화에 따른 항복면의 변화와 Gens(1995)에 기술되어 있는 온도 변화에 따른 항복면 변화, 그리고 과압밀 압축지수를 Rutqvist et al.(2011)에서와 마찬가지로 흡입력의 함수(κ(s)))로 고려한 열･탄소성 버전의 BBM 해석 모듈을 개발하였다. 기본적으로 Rutqvist et al.(2011)에 기술된 TOUGH2-FLAC 시뮬레이터에서의 BBM 해석 모듈 개발과 동일하게 FLAC에서 제공하고 있는 User Defined Model(UDM)과 MCCM을 이용하였지만, 본 연구에서는 평균유효응력의 변화뿐만 아니라 흡입력의 변화에 따른 소성변형률을 모두 고려하였으며, 소성변형에 따른 항복면 변화를 명확하게 기술하고 이를 모두 반영하여 BBM 해석 모듈을 개발하였다.

3.1 응력상태

BBM에서 압축응력은 양의 값으로 정의되고, 평균유효응력(p')과 von Mises 응력(q)은 각각 식 (1)과 식 (2)로 표현된다.

여기서, p', p, 그리고 p^∅는 각각 평균유효응력(mean effective stress), 전응력(total stress), 그리고 공극수압(pore pressure)을 의미한다. 또한 σ₁, σ₂, 그리고 σ₃는 각각 주응력(principal stress)을 의미하며, σ'₁, σ'₂, 그리고 σ'₃는 각각 유효 주응력(effective principal stress)을 의미한다.

TOUGH2-MP/FLAC3D를 이용하여 THM 복합거동을 해석하는 경우, 공극수압과 흡입력은 TOUGH2-MP로부터 계산되어 FLAC3D로 전달된다. TOUGH2-MP에서 단상 유동 현상(single phase fluid condition)을 해석하는 경우, 단상 유동(single fluid phase)이 기상(gas phase)인 경우는 p^∅는 p^g로 계산되며, 액상(liquid phase)에 대해서는 p^∅는 p^l로 계산된다. 그리고 이상 유동 현상(two phase fluid condition)에 대해서는 기상 압력(gas pressure)이 액상 압력(liquid pressure)에 비해 월등히 크기 때문에 식 (1)은 식 (3)에 의해 평균유효응력이 계산되며, 흡입력은 식 (4)에 의해 계산된다(Pruess, 1999).

여기서, p^g와 p^l는 각각 기상과 액상의 압력(gas and liquid phase pressure)을 의미한다.

3.2 응력, 온도, 그리고 흡입력 변화에 따른 변형률

체적 변형률(volumetric strain)은 식 (5)와 같이 계산되며, 등가 편차 변형률(equivalent deviatoric strain)은 식 (6)에 의해 계산된다.

여기서, ε₁, ε₂, 그리고 ε₃는 주응력 방향으로의 변형률(principal strain)을 의미한다.

BBM에서 변온조건에서의 변형률 변화는 식 (7)과 같이 계산된다.

여기서, dε^e, dε^p, dε^s, 그리고 dε^T는 각각 탄성 변형률(elastic strain) 변화, 소성 변형률(plastic strain) 변화, 흡입력에 의한 변형률(suction strain) 변화 그리고 열적 변형률(thermal strain) 변화를 의미한다.

3.2.1 탄성 변형률

개발된 BBM 해석 모듈에서 식 (8)의 체적 탄성변형률의 변화$\left(d{\epsilon }_v^e\right)$는 식 (9)에 나타나 있는 체적변형계수(bulk modulus, K)와 평균유효응력의 변화로 계산되며, 과압밀 압축지수(κ(s))는 식 (10)과 같이 흡입력의 함수로 계산된다. 그리고 편차 탄성 변형률(deviatoric elastic strain) 변화$\left(d{\epsilon }_q^e\right)$은 전단변형계수(shear modulus, G)를 이용하여 식 (11)과 같이 계산된다.

여기서, κ_ps0는 포화 상태의 과압밀 압축지수를 의미하고, α_ps는 물질상수, 그리고 e는 간극비(void ratio)이다.

3.2.2 소성 변형률

Fig. 4와 Fig. 6에서 나타나 있듯이, 탄성영역의 경계는 식 (12)에 나타나 있는 온도와 흡입력의 변화에 따른 항복면 함수(loading-collapse(LC) yield surface)에 따라 결정된다. 흡입력과 온도변화에 따른 점착력(cohesion, p_s(s,T))변화는 식 (13)으로 계산되고, 흡입력과 온도에 따른 선행압밀응력(preconsolidation stress, p₀(s,T))은 식 (14)로 계산된다. 식 (14)에서 온도 변화에 따른 포화 상태에서의 선행압밀응력($p_0^{*}\left(T\right)$)은 식 (15)에 의해 계산되며, 흡입력의 함수인 정규압밀 압축지수(λ_ps (s))는 식 (16)에 의해 계산된다.

여기서, M은 p―q공간에서 한계상태선(critical state line)의 기울기이며, k_s는 경험적 물질 상수이고, p_s0는 포화상태에서의 인장강도를 의미한다. p^c는 선행압밀응력에 대한 기준압력(reference stress)이고 $p_0^{*}$는 완전 포화 조건하에서의 기준 온도에서의 선행압밀응력을 의미하며, r_λ와 β_λ는 물질상수이다.

응력상태가 항복면에 도달하게 되면, 평균유효응력 변화에 따른 소성 변형률$\left(d{\epsilon }_q^p\right)$과 von Mises 응력 변화에 따른 소성 변형률$\left(d{\epsilon }_q^p\right)$은 식 (17)에 나타나 있는 소성 포텐셜(plastic potential) 함수를 각각 p'와 q로 미분하여 식 (18)과 식 (19)로 나타낼 수 있다. 또한 식 (8)과 식 (11)에 식 (18)과 식 (19)를 각각 다시 대입하여 정리하면, 식 (20)과 식 (21)과 같이 정리된다.

여기서, c_a는 M²(2p'-p₀(s,T) + p_s (s,T))이고, c_b는 α_a2q이다.

FLAC3D를 이용한 BBM 해석 모듈에서 새로운 응력들(p'_new, q_new)은 MCCM에서와 동일하게 식 (22)와 식 (23)으로 계산되는데, 이를 식 (17)에 대입하여 정리하면 식 (24)와 같이 정리된다. 이때, 각각의 상수 a, b, 그리고 c는 식 (25) ~ 식 (27)과 같다. FLAC3D에서 제공하는 MCCM과 마찬가지로 식 (24)의 해인 소성 계수(plastic multiplier)의 변화량(dΛ)은 상수 a, b, 그리고 c를 이용하여 도출할 수 있다(Itasca, 2012).

Rutqvist et al.(2011)에서는 식 (21)의 마지막 항에 있는 3G가 체적변형계수(K)로 표기되어 있고, 식 (23)에 나타나 있는 3G 모두 체적변형계수로 표기되어 있는 오류가 있어 본 연구에서는 이를 정정하여 표기하였다. 또한, BBM에서 소성변형은 온도와 흡입력의 변화에 따른 항복면 함수(loading-collapse(LC) yield surface)에 의해 결정될 뿐만 아니라, Fig. 4에 나타나 있는 흡입력 증가(suction increase, SI)가 기존의 불포화토 최대 흡입력(s₀))을 넘어서게 되도 소성변형은 일어난다(Alonso et al., 1990). 선행연구인 Rutqvist et al.(2011)에서는 흡입력 증가에 따른 소성변형(($d{\epsilon }_s^p$)을 고려하지 않았지만, 본 연구에서는 Alonso et al.(1990)에 기술되어 있는 식 (28)을 이용하여 흡입력 증가에 따른 소성변형$\left(d{\epsilon }_s^p\right)$을 추가적으로 고려하였다.

여기서, υ, λ_s 그리고 p_atm는 각각 비체적, 초기 상태(virgin state)에서 흡입력 변화에 대한 강성 파라미터, 그리고 대기압이며, κ_sp (p',s)는 불포화 상태에서의 수리적 팽윤지수를 의미하고 식 (29)와 같이 계산된다.

여기서, p_ref는 기준 압력(reference stress)을 의미하고, κ_sp0는 포화 상태와 기준 압력 p_ref에서의 수리적 팽윤지수를 의미한다. 또한 α_sp와 α_ss는 물질상수이다.

본 연구에서는 평균유효응력뿐만 아니라 흡입력의 변화에 따른 소성변형률을 모두 고려하여 소성변형률 변화$\left(d{\epsilon }_v^p\right)$를 식 (30)으로 계산되도록 BBM 해석 모듈을 개발하였으며, 식 (31)과 식 (32)를 이용하여 경화법칙(hardening law)에 따른 두 항복면 LC와 SI의 변화를 모두 반영할 수 있도록 BBM 해석 모듈을 개발하였다.

3.2.3 열적 변형률

BBM에서 온도 변화에 따른 열응력으로 야기되는 체적 변형률은 식 (33)로 계산될 수 있다(Alonso et al., 1990). 기본적으로 FLAC3D에서는 열적 변형률이 식 (34)와 같이 자동으로 계산되기 때문에 본 연구에서는 선열팽창계수(α_T)를 FISH 함수를 이용하여 수정하여 열적 변형률을 식 (33)과 같이 계산되도록 하였다.

여기서, α₀와 α₂는 온도변화에 따른 열팽창 관련 물질 상수, α_T는 선열팽창계수(linear thermal expansion coefficient), 그리고 ΔT는 초기 온도와 현재 온도의 차이를 의미한다.

3.2.4 흡입력에 따른 변형률

흡입력 변화로 야기되는 체적 변형률$\left(d{\epsilon }_v^s\right)$은 식 (35)를 이용하여 계산되며, 흡입력 강성 상수(suction bulk modulus, K^s)는 식 (36)으로 계산되고 그에 따른 팽윤압의 변화(dp^s)는 FLAC3D에서 열응력 계산 과정과 유사하게 식 (37)과 같이 계산된다. 계산된 팽윤압의 변화는 선행연구(Rutqvist et al., 2011)와 마찬가지로 식 (38)과 같이 등방으로 작용하는 것으로 가정하였으며 FISH 함수를 이용하여 이를 반영하였다.

여기서, K^s는 흡입력 강성 상수를 의미하며, 수리적 팽윤지수 κ_sp (p',s)는 식 (29)에서 정의된 값이다.

3.3 건조밀도

비체적은 식 (39)와 같이 정의되며, TOUGH2-MP/FLAC3D를 이용한 수치해석에서 (i+1)번째 해석 스텝의 비체적(υⁱ⁺¹), 간극비(eⁱ⁺¹) 그리고 공극률(∅ⁱ⁺¹)은 각각 식 (40) ~ 식 (42)로 계산된다. 또한 입자의 부피(V^s)는 변하지 않고 공극의 부피(V^∅)만 변하는 것으로 가정하여 요소의 건조밀도(${\rho }_d^{i+1}$)는 FLAC3D에서 식 (43)과 같이 계산된다.

여기서, V^s는 토질의 부피 V에 포함된 입자의 부피, V^∅는 공극의 부피이고, υⁱ⁺¹, ∅ⁱ⁺¹, eⁱ⁺¹, 그리고 ${\rho }_d^{i+1}$는 새로운 스텝(i+1)에서 계산된 비체적, 공극률, 간극비, 그리고 건조밀도를 의미한다. ρ_grain은 입자의 밀도를 의미한다.

4. BARCELONA BASIC MODEL 해석 모듈 검증

본 연구에서는 개발된 BBM 해석 모듈을 다음의 세 가지 방법으로 검증하였다. 먼저, FLAC3D에서 제공하고 있는 MCCM의 등방 압밀시험 예제(Itasca, 2012)를 활용하여 포화 상태에서의 BBM 해석 모듈을 검증하였고, 이후, COMSOL Multiphysics를 이용하여 Alonso et al.(1990)에서 제시한 예제를 수치시험한 결과와 개발된 BBM 해석 모듈을 이용한 해석결과를 비교하여 포화 상태뿐만 아니라 불포화 상태에서의 BBM 해석 모듈을 검증하였다. 마지막으로, Quick Tools의 검증을 위해 Kristensson and Ảkesson(2008)에서 수행된 일련의 검증 방법을 동일하게 적용하여 개발된 BBM 해석 모듈을 최종적으로 검증하였다.

4.1 MCCM 예제를 이용한 포화 상태에서의 해석 모듈 검증

BBM은 포화 상태에서 MCCM과 동일한 모델이므로 불포화 상태에서 모델을 검증하기 에 앞서 개발된 BBM 해석 모듈이 포화 상태에서 MCCM과 동일한 거동을 보이는 지를 먼저 검증할 필요가 있다. 이를 위해 사용된 예제는 Itasca(2012)에 기술되어 있는 MCCM의 예제로 단일 요소에 대한 등방 압밀 시험이다(Fig. 7). 사용된 등방 압밀 시험은 x=0 평면, y=0 평면, 그리고 z=0 평면의 법선방향 변위를 구속하고 반대쪽 세 방향에서 0.5×10^-4m/s로 재하(loading)와 제하(unloading)를 반복하도록 설계되었다. 검증을 위해 사용된 BBM 물성은 Table 1에 정리되어 있으며, MCCM에서 사용된 물성과 동일하다.

MCCM과 BBM을 이용하여 수행된 두 가지의 수치시험에서 일련의 재하와 제하 과정을 거치는 동안 평균유효응력과 변위의 변화, 로그스케일의 평균유효응력과 비체적 변화, 그리고 변위와 체적 탄성계수 및 전단 탄성계수 변화는 Fig. 8에 나타나 있듯이 정확히 일치하는 것으로 나타났다. 이러한 결과를 토대로 볼 때, 개발된 BBM 해석 모듈이 포화 상태에서 MCCM과 동일하게 계산되고 있음을 알 수 있다.

Fig. 7.

Configuration of simulation test to verify the BBM implementation into FLAC3D under fully saturated condition using MCCM example in FLAC3D manual

Table 1. Input parameters for the modeling of isotropic compression test to verify BBM module using MCCM example in FLAC3D manual

Fig. 8.

Result of simulation tests to verify the implementation BBM into FLAC3D within the UDM capability of FLAC3D

4.2 Alonso 예제와 COMSOL Multiphysics 해석결과를 이용한 해석 모듈 검증

두 번째 검증에서는 Alonso et al.(1990)의 예제를 이용하여 흡입력과 응력변화에 따른 비체적 변화를 계산한 후 그 결과를 개발된 BBM 해석 모듈을 이용한 결과와 비교하였다. BBM 해석 모듈을 이용한 결과와 Alonso et al.(1990)에 나타나 있는 결과가 일치하지만, Alonso et al.(1990)에 나타나 있는 결과가 낮은 해상도이기 때문에 본 논문에 수치시험의 결과와 함께 나타내는 것이 적절하지 않을 것으로 판단되어, 본 연구에서는 추가적으로 COMSOL을 이용한 해석결과와 개발된 BBM 해석 모듈을 이용한 수치시험의 결과를 비교함으로써 간접적인 방법으로 다시 검증하고자 하였다. 사용된 물성은 Table 2에 나타난 것처럼 Alonso et al.(1990)에서 나타나 있는 것과 동일하다.

Table 2. Input parameters for numerical simulations of examples presented in Alonso et al. (1990)

4.2.1 Alonso 예제 1

Alonso et al.(1990)의 예제는 흡입력과 평균유효응력을 서로 다른 경로로 변화시키면서 그에 따른 비체적의 변화를 계산하는 것이다. 첫 번째 예제의 첫 번째 경로는 $\overrightarrow{ABDF}$로, 평균유효응력이 0.15 MPa이고 초기 흡입력이 0.2 MPa인 탄성영역 내의 A지점에서 포화도의 증가로 완전 포화상태(s = 0)에 이르는 $\overrightarrow{AB}$ 구간, 이후 평균유효응력이 증가하는 $\overrightarrow{BD}$ 그리고 $\overrightarrow{DF}$ 구간으로 이루어져 있다(Fig. 9(a)). 흡입력이 감소하는 $\overrightarrow{AB}$ 구간에서는 팽윤에 의해 비체적은 증가하는 것으로 계산되었으며, 평균유효응력이 증가하는 $\overrightarrow{BD}$ 구간에서는 Fig. 9(a)에 나타나 있는 초기 항복면(LC_i)을 만나기 직전까지(p < 0.2 MPa) 탄성 거동을 보이며 감소하고, 항복면을 만나는 순간부터는(p ≥ 0.2 MPa) 소성 거동을 보이며 비체적이 급격하게 감소하는 것으로 나타났다(Fig. 9(b)).

두 번째 경로는 $\overrightarrow{ACDF}$로, 초기 평균유효응력이 0.15 MPa에서 0.35 MPa로 증가하는 $\overrightarrow{AC}$ 구간, 이후 소성영역에서 흡입력이 0.0 MPa로 감소하는 $\overrightarrow{CD}$ 구간과 평균유효응력이 0.6 MPa로 증가하는 $\overrightarrow{DF}$ 구간으로 이루어진다(Fig. 9(a)). 흡입력의 변화는 없고 평균유효응력의 증가하는 $\overrightarrow{AC}$ 구간에서는 항복면에 도달하기 전까지는 탄성 거동을 보이다가 항복면을 만난 이후에는 소성 거동을 보이며 비체적이 급격하게 줄어드는 것으로 계산되었다(Fig. 9(b)). 이후, $\overrightarrow{CD}$ 구간에서는 흡입력의 감소로 인해 급작스러운 소성변형이 일어나면서 붕괴(collapse)가 일어나게 되는데, 이러한 소성변형 역시 잘 계산되는 것으로 나타났다. 첫 번째 경로 $\overrightarrow{ABDF}$구간에서의 $\overrightarrow{AB}$ 구간에서와 마찬가지로 $\overrightarrow{CD}$ 구간에서 팽윤에 의한 비체적의 증가도 있지만, 소성영역에서 LC 곡선이 지속적으로 변화함에 따라 압축변형이 더욱 크게 발생하기 때문에 전체적으로 비체적이 크게 줄어드는 것으로 계산되었다. 이후, 또 다시 평균유효응력이 증가하는 $\overrightarrow{DF}$ 구간에서는 소성변형이 발생하면서 비체적은 더욱 감소하는 것으로 계산되었다.

마지막 세 번째 경로 $\overrightarrow{ACEF}$는 초기 평균유효응력 0.15 MPa에서 0.35 MPa을 지나 0.6 MPa로 증가하는 $\overrightarrow{AC}$ 및 $\overrightarrow{CE}$ 구간과 이후, 소성영역에서 흡입력이 0.0 MPa로 감소하는 $\overrightarrow{EF}$ 구간으로 이루어진다(Fig. 9(a)). 평균유효응력이 증가되는 $\overrightarrow{AC}$ 및 $\overrightarrow{CE}$ 구간에서는 두 번째 경로와 마찬가지로 항복면에 도달하기 전까지는 탄성 거동을 보이다가 항복면을 만난 이후에는 소성 거동을 보이며 비체적이 급격하게 줄어드는 것으로 계산되었다. 이후, 흡입력이 0.0 MPa로 변화되는 $\overrightarrow{EF}$ 구간에서는 탄성영역을 벗어난 영역에서 지속적으로 흡입력이 변화하기 때문에 두 번째 경로의 $\overrightarrow{CD}$ 경로와 유사하게 포화 상태의 선행압밀 응력($p_0^{*}$)이 지속적으로 변화하게 되어 큰 소성변형이 일어나 급격한 붕괴가 일어나는 것으로 계산되었다(Fig. 9(b)).

Fig. 9.

Response of model to decreasing suction (wetting) at increasing net mean stress

4.2.2 Alonso 예제 2

두 번째 검증 예제의 첫 번째 $\overrightarrow{ABDF}$ 경로와 세 번째 $\overrightarrow{ACEF}$ 경로는 각각 첫 번째 검증 예제의 세 번째 $\overrightarrow{ACEF}$ 경로와 첫 번째 $\overrightarrow{ABDF}$ 경로와 동일하기 때문에 두 경로에 대한 설명과 수치시험 결과에 대한 언급은 생략하기로 한다. Fig. 10(a)의 푸른색 파선으로 나타나 있는 두 번째 $\overrightarrow{ACDF}$ 경로는 초기 평균유효응력이 0.15 MPa로 유지되면서 흡입력이 0.2 MPa에서 0.1 MPa로 감소하는 $\overrightarrow{AC}$ 구간, 이후 평균유효응력이 0.15 MPa에서 0.6 MPa로 증가되는 $\overrightarrow{CD}$ 구간, 그리고 다시 흡입력이 0.0 MPa로 감소하는 $\overrightarrow{DF}$ 구간으로 이루어진다. Fig. 10(b)를 살펴보면, 흡입력이 변화되는 $\overrightarrow{AC}$ 구간에서의 비체적은 팽윤에 의해 증가했지만, 흡입력의 변화가 $\overrightarrow{AE}$ 구간에 비해 상대적으로 작기 때문에 비체적 변화가 $\overrightarrow{AE}$ 구간보다 상대적으로 작게 발생하는 것으로 계산되었다. 이후 평균유효응력이 0.15 MPa에서 0.6 MPa로 증가하는 $\overrightarrow{CD}$ 구간에서는 탄성 거동을 보이다가 항복면을 만난 이후 소성 거동으로 비체적은 감소하는 것으로 계산되었다. 이후 항복면을 지난 구간에서 흡입력이 감소하는 $\overrightarrow{DF}$ 구간에서는 첫 번째 검증 예제와 마찬가지로 급작스러운 붕괴로 인해 비체적이 감소하는 것으로 계산되었다.

첫 번째와 두 번째 검증 예제에서는 불포화 상태인 A지점에서 서로 다른 경로를 따라서 포화상태인 F지점으로 이동하는 과정에서의 비체적 변화를 살펴보았다. 경로에 따라 서로 다른 비체적 변화이력을 갖지만 최종 비체적은 경로에 따라 상관없이 동일한 값으로 계산되는 것으로 나타났다(Fig. 9(b), Fig. 10(b)). 이는 첫 번째와 두 번째 수치시험이 경로에 상관없이 초기 항복면 곡선(LC_i)을 지나 소성변형이 일어난 이후, 최종적으로 F지점에서 동일한 항복면 곡선(LC_F)을 가지기 때문이다. 즉, 경로는 다르지만, 최종 위치(F)에 도달할 때까지 변화되는 흡입력 변화와 평균유효응력이 같고, 포화 상태의 선행압밀 응력($p_0^{*}$)의 최대값이 같기 때문에 변화된 비체적은 동일하게 계산된 것이다.

Fig. 10.

Response of model to decreasing suction (wetting) at increasing net mean stress

4.2.3 Alonso 예제 3

앞의 두 예제와는 달리 세 번째 검증 예제에서는 포화 상태에서 불포화 상태로 변화하면서 비체적의 변화를 계산하는 것이다(Fig. 11(a)). 붉은색 실선으로 나타나 있는 첫 번째 $\overrightarrow{ABDF}$ 경로에 대한 수치시험 결과를 보면, 평균유효응력이 0.15 MPa에서 0.60 MPa로 증가하는 구간에서는 앞의 두 예제와 유사하게 탄성 거동 이후, 소성 변형이 발생하면서 비체적이 감소하게 되고, 포화 상태의 선행압밀 응력($p_0^{*}$)이 지점에 위치하게 되는 LC 곡선(LC_B)을 최종적으로 갖게 된다. 따라서 포화 상태에서 불포화 상태로 변하면서 흡입력이 증가하게 되는 $\overrightarrow{BD}$ 구간 및 $\overrightarrow{DF}$ 구간은 탄성영역에 위치하게 되고, 흡입력의 증가에 따른 수축으로 비체적이 감소하는 것으로 계산되었다(Fig. 11(b)).

Fig. 11.

Response of model to increments of net mean stress (p) and increasing suction (drying)

푸른색의 파선으로 나타나 있는 두 번째 $\overrightarrow{ACDF}$ 경로에 대한 수치시험 결과를 보면, 탄성영역에서 변화되는 $\overrightarrow{AC}$ 구간에서는 포화 상태(s = 0.0 MPa)에서 불포화 상태(s = 0.1 MPa)로의 변화함에 따라 비체적은 감소하고, 평균유효응력이 증가하는 $\overrightarrow{CD}$ 구간에서는 탄성변형과 소성변형이 일어난다. 이때, D지점을 지나는 LC 곡선(LC_D)이 형성되어 흡입력이 다시 증가하는 $\overrightarrow{DF}$ 구간은 탄성영역에 위치하게 되고(Fig. 11(a)), 흡입력의 증가에 따른 수축으로 비체적이 감소하는 것으로 계산되었다.

검은색 쇄선으로 나타나 있는 세 번째 $\overrightarrow{ACEF}$ 경로에 대한 수치시험 결과를 보면, 포화 상태에서 불포화 상태(s = 0.2 MPa)로 변화하는 $\overrightarrow{AC}$구간 및 $\overrightarrow{CE}$ 구간에서 비체적의 감소가 나타났다. 이후 $\overrightarrow{EF}$ 구간에서는 탄성 거동을 보이다가 항복면을 만난 이후 소성 거동으로 비체적은 감소하는 것으로 계산되었고, 세 번째 $\overrightarrow{ACEF}$ 경로에서는 최종적으로 F지점을 지나는 LC 곡선(LC_F)이 형성되었다(Fig. 11(b)).

포화 상태에서 불포화 상태로 변화하는 검증 예제에 대한 수치시험 결과를 살펴보면, 첫 번째와 두 번째의 예제에 대한 수치시험 결과와는 달리 경로에 따라 F지점에서의 최종 비체적은 서로 다르게 계산된다. 이는 첫 번째와 두 번째 예제에서는 경로에 상관없이 최종지점(F)에서 갖는 LC 곡선이 동일하지만, 세 번째 예제에서는 경로에 따라 최종지점(F)에서 갖는 LC 곡선의 크기가 LC_B < LC_D < LC_F 순으로 서로 다르기 때문에 소성 변형이 일어난 정도가 달라 최종 지점 F에서의 비체적이 포화 상태의 선행압밀 응력($p_0^{*}$)의 크기($p_{0,F}^{*}$ < $p_{0,D}^{*}$ < $p_{0,B}^{*}$)와는 반대로 계산되었다(Fig. 11(b)).

4.2.4 Alonso 예제 4

지금까지의 Alonso 예제들에서는 초기 불포화토 최대흡입력(s₀)이 0.3 MPa이지만, Fig. 12(a)에 나타나 있듯이 네 번째 예제에서의 초기 불포화토 최대흡입력(s₀)은 0.025 MPa이다. 먼저, Fig. 12(b)에 나타나 있는 첫 번째 경로인 $\overrightarrow{AB}$ 구간의 결과를 살펴보면, 앞서 언급한 예제들과 유사하게 최대 흡입력 변화는 없이 응력만 증가하기 때문에 흡입력 증가(SI)에 따른 항복면 변화는 없으며 탄성 변형 이후 초기 항복면 곡선(LC_i)을 지나 소성변형이 일어난다.

반면, 두 번째 경로인 $\overrightarrow{ACDE}$ 구간은 초기 최대 흡입력(s₀)을 넘어 0.3 MPa까지 흡입력이 증가하는 $\overrightarrow{AC}$ 구간과 이후 다시 초기 상태인 포화 상태로 돌아오는 $\overrightarrow{CD}$ 구간, 그리고 마지막으로 응력이 증가하는 $\overrightarrow{DE}$ 구간으로 이루어져 있다. $\overrightarrow{AC}$ 구간에서는 Fig. 11(a)에 나타나있는 세 번째 예제의 $\overrightarrow{ACE}$ 경로에서처럼 흡입력의 증가에 따른 수축으로 비체적은 감소한다(Fig. 12(b)). 이때, 초기 최대 흡입력(s₀)을 넘어서게 되면 식 (28)에 의해 소성변형이 일어나게 되어 Fig. 12(b)에 나타나 있듯이 세 번째 예제의 $\overrightarrow{ACE}$ 경로보다 더 큰 변형이 일어나게 되고 식 (31)에 의해 새로운 포화 상태의 선행압밀 응력($p_0^{*}$)을 갖게 되어 새로운 LC_C를 형성하게 된다. 이후, $\overrightarrow{CD}$ 구간은 새로운 LC_C에 의해 탄성영역에 위치하게 되며 흡입력의 감소에 의한 팽윤 거동으로 비체적은 증가하게 된다. $\overrightarrow{CD}$ 구간에서는 소성변형이 일어나지 않기 때문에, 항복면은 C지점에서와 동일하다(LC_C = LC_D). 이후, 평균유효응력이 증가하는 $\overrightarrow{DE}$ 구간에서는 탄성변형이 일어나다가 LC_D와 만나게 되면서 소성변형이 일어나는 것으로 나타났다(Fig. 12(b)).

Fig. 12.

Response of model to change of SI and increments of net mean stress (p)

4.3 실험실 시험 결과 및 선행연구들 이용한 해석 모듈 검증

Kristensson and Åkesson(2008)은 BBM 파라미터를 효과적으로 도출하기 위해 MathCad를 이용하여 프로그램(Quick tools)을 개발하였고 이를 검증하기 위해 Dueck and Nilsson(2010)에 보고된 MX80 벤토나이트에 대한 실험실 시험 결과와 Code_Bright를 이용한 수치시험 결과를 이용하였다. 또한, Rutqvist et al.(2011) 역시 TOUGH-FLAC을 이용한 BBM 해석 모듈의 검증을 위해 Kristensson and Åkesson(2008)에서 수행한 방법과 동일한 방법을 이용하였기에 본 연구에서도 동일한 방법으로 개발된 BBM 해석 모듈을 검증하고자 한다. 검증을 위한 수치시험에 사용된 각각의 BBM 입력 물성값은 Table 3에 나타나 있으며 Kristensson and Åkesson(2008)에서 사용한 값과 정확히 일치한다. 본 연구에서 사용된 초기 비체적(υ_ini)은 Kristensson and Åkesson(2008)에 나타나 있는 값과 동일한 값을 사용하였지만, Rutqvist et al.(2011)에서 사용된 값은 초기 비체적(υ_ini)은 Kristensson and Åkesson(2008)와 본 연구에서 사용한 값과는 차이를 보인다(Table 3).

Table 3. Input parameters of numerical simulations for the laboratory tests

4.3.1 불포화 조건에서의 압축시험

개발된 BBM 해석 모듈의 검증을 위해 사용된 첫 번째 실험실 시험은 Dueck and Nilsson(2010)에 언급된 O4_1105 시험으로 직경 35 mm, 높이 14 mm의 MX80 벤토나이트 블록에 대한 압축시험이다(Fig. 13). O4_1105 시험의 초기 축방향 압축응력과 횡방향 압축응력은 각각 0.18 MPa과 2.97 MPa이고 시험이 진행되는 동안 초기 흡입력 28 MPa는 일정하게 유지되었다. 이후 횡방향으로 구속된 상태에서 축방향으로 세 번의 재하(0.18 MPa → 1.00 MPa → 8.13 MPa → 19.77 MPa)와 두 번의 제하(19.77 MPa → 9.12 MPa → 1.00 MPa)를 하면서 간극비의 변화와 횡방향의 응력을 측정하였다.

개발된 BBM 해석 모듈을 이용한 O4_1105 시험의 수치해석 결과는 선행연구의 결과들과 함께 Fig. 14에 나타내었다. 실선과 파선은 각각 α_a = 1인 연관 소성 모델(associative model)과 α_a ≠ 1인 비연관 소성 모델(non-associative model)의 결과를 의미한다. 계산된 간극비와 횡방향 응력은 각각 붉은색과 푸른색으로 나타내었고(Fig 14(a)), 계산된 평균유효응력에 대한 von Mises 응력 값은 붉은색으로 나타내었다(Fig. 14(b)). 실험실 시험 결과는 ‘×’로 나타냈으며, Kristensson and Åkesson(2008)에 나타나 있는 Code_Bright 해석 결과는 ‘○’로 나타내고 Quick tools의 결과는 검은색 실선으로 나타내었다. 실험실 시험결과는 세 번의 재하와 두 번의 제하를 수행하면서 총 5회 측정하였으므로 연속적인 결과를 알 수 없기 때문에 각각의 실험데이터를 기준으로 하여 점선으로 연결하여 표기하였다.

Fig. 13.

Device used for the oedometer tests (modified from Dueck and Nilsson, 2010)

개발된 BBM 해석 모듈을 이용하여 계산된 간극비, 횡방향 응력, 평균응력, 그리고 von Mises 응력 모두 Kristensson and Ảkesson(2008)에 나타나 있는 Code_Bright의 결과(○)와 비교하면 일치할 뿐만 아니라, Kristensson and Ảkesson(2008)에 나타나 있는 Quick tools의 결과(검은색 실선)와는 완벽하게 일치하여 검은색 실선이 보이지 않는다(Fig. 14). 이러한 사실로 볼 때, 개발된 불포화 조건에서 BBM 해석 모듈의 해석 결과에 대한 신뢰성이 충분히 확보된 것으로 판단된다. 또한, 연관 소성 모델의 경우, 축방향 응력에 대한 간극비와 횡방향 응력 변화 곡선은 큰 차이를 보이고 있지 않지만(Fig. 14(a)), 평균유효응력에 대한 von Mises 응력의 변화 이력 곡선은 큰 차이를 보이고 있다(Fig. 14(b)). 비연관 소성 거동을 보이는 매질에 대한 수치해석을 수행할 때에는 비연관 소성 파라미터(α_a)의 영향을 반드시 고려하여 입력 물성 산정에 주의를 기울여야 할 것으로 판단된다.

Fig. 14.

Numerical simulation results of compression test under constant suction with experimental and modeling results presented in Kristensson and Ảkesson (2008) for MX80 bentonite

4.3.2 팽윤시험

검증을 위한 두 번째 실험실 시험은 Dueck and Nilsson(2010)에 언급된 O3_1105 및 O4_1005 시험으로 앞에서 언급한 O4_1105 시험과 같은 장비를 사용하였으며, 동일한 크기의 MX80 벤토나이트 블록을 이용한 팽윤시험이다(Fig. 13). 두 시험의 초기 축방향 및 횡방향 압축응력은 각각 8.9 MPa과 2.54 MPa이고 시험이 진행되는 동안 축방향 응력은 일정하게 유지되었다. 이후, 횡방향으로 구속된 시편에 물을 주입하여 초기 벤토나이트의 흡입력을 101.5MPa에서 12.6 MPa까지 감소시키며 간극비와 횡방향 응력 변화를 계측하였다.

개발된 BBM 해석 모듈을 이용하여 O3_1105 및 O4_1005에 대한 수치시험을 수행한 결과를 Kristensson and Ảkesson(2008)에 나타나 있는 결과와 함께 Fig. 15에 나타내었다. 예제로 선정이 된 두 시험은 횡방향으로 변위가 구속된 조건에서 물 주입이 이루어지기 때문에 벤토나이트의 팽윤으로 인한 간극비가 증가와 횡방향 응력이 동시에 증가하게 된다(Fig. 15(a)). 지속적인 물 주입으로 흡입력이 약 80 MPa로 감소하게 되면 벤토나이트는 초기 항복면에 도달하게 되고(Fig. 15(b)), 이후 항복구간에서는 흡입력 감소로 인한 팽윤과 소성변형으로 인한 압축변형이 서로 상쇄되어 간극비의 변화가 크게 발생하지 않게 된다(Fig. 15(a)). 전체적으로 개발된 BBM 해석 모듈의 경우, 벤토나이트 블록의 팽윤과 항복면 이후의 소성 거동을 잘 모사할 수 있는 것으로 나타났고, Kristensson and Ảkesson(2008)의 결과와 개발된 BBM 해석 모듈을 이용한 수치시험의 결과는 앞에서의 경우와 마찬가지로 정확히 일치하기 때문에 Fig. 15(a)에서 구분이 되지 않고 있으며, Code_Bright의 결과와도 일치하는 것으로 계산되었다.

4.3.3 포화 조건에서의 압축시험

검증을 위해 사용된 세 번째 실험실 시험은 Dueck and Nilsson(2010)에 언급된 T39 시험으로 직경 152 mm, 높이 30 mm의 MX80 벤토나이트 블록에 대한 압축시험이다(Fig. 16). T39시험은 먼저 벤토나이트 블록을 완전히 포화시키고 난 이후, 초기 축방향 및 횡방향 압축응력을 모두 1.1 MPa로 설정한 이후, 횡방향 응력을 일정하게 유지하면서 축방향 변형률을 10% 까지 증가시키면서 축방향 응력과 축방향 변위를 측정한 압축시험이다.

Fig. 15.

Numerical simulation results of swelling test under constant axial load with experimental and modeling results presented in Kristensson and Ảkesson (2008) for MX80 bentonite

Fig. 16.

Rowe oedometer with MX80 sample having a diameter of 152 mm and an initial height of 30 mm (modified from Dueck and Nilsson, 2010)

수치시험에서는 수직 변형률을 증가시키기 위해 상부에서 5.0×10^-7m/s의 속도로 가압을 하였으며, 축방향 가압에 따른 평균유효응력과 von Mises 응력의 변화를 실험실 데이터와 함께 Fig. 17에 도시하였다. 실험실 시험 및 수치시험 모두, 축방향의 응력이 서서히 증가하여 초기 항복면을(initial yield surface) 만날 때까지(q≤300 kPa) 탄성 거동을 보이고 있으며, 항복면을 지난 이후부터 소성변형을 나타내는 것으로 나타났다(Fig. 18(a)). von Mises 응력의 변화 뿐만 아니라, 축방향 압축변형에 의한 탄성영역 및 소성영역에서의 계산된 부피 변화 역시 실험실 시험의 결과를 잘 모사하고 있으며, 기존의 선행연구의 결과와 일치하는 것을 알 수 있다(Fig. 18(b)). 팽윤시험과 포화 조건에서의 압축시험에 대한 모델링 결과를 살펴볼 때, 개발된 BBM 해석 모듈의 포화조건에서의 수리-역학적 해석 결과에 대한 신뢰성 역시 확보된 것으로 판단된다.

Fig. 17.

Evolution of deviatoric stress during the numerical simulation in this study and laboratory test (T39) presented in Dueck and Nilsson (2010)

Fig. 18.

Numerical simulation results of swelling test under fully saturated condition with experimental and modeling results presented in Kristensson and Ảkesson (2008)

5. 결 론

본 연구에서는 흡입력과 온도변화를 고려하여 불포화토의 역학적 거동 해석을 수행할 수 있는 BBM에 대해 살펴보고, TOUGH2-MP/FLAC3D를 이용하여 해석을 수행할 수 있도록 BBM 해석 모듈을 개발하였다. 개발된 BBM 해석 모듈은 1990대 초반에 처음 제안된 탄소성 버전의 BBM과 온도 변화에 따른 항복면 변화를 고려할 수 있도록 제안된 열-탄소성 버전의 BBM을 토대로 개발되었으며, 과압밀 압축지수를 흡입력의 함수 형태로 반영할 수 있도록 개발하였다. 본 연구에서 개발된 BBM 해석 모듈은 TOUGH-FLAC을 이용하여 BBM 해석 모듈을 개발한 선행연구와 동일한 방식으로 개발되었으나, 소성변형에 따른 항복면 변화를 추가적으로 고려하여 열-탄소성 버전의 BBM 해석 모듈을 개발하였다. 개발된 BBM 해석 모듈의 검증을 위해 각종 다양한 예제에 대한 수치시험을 수행하였으며, 다양한 실험실 시험에 대한 수치시험을 추가적으로 수행하여 그 결과를 Quick Tools 및 Code_Bright의 결과와 비교하였다. 개발된 BBM 해석 모듈을 이용한 수치시험 결과와 검증 예제의 결과 그리고 실험실 시험 결과 및 선행연구의 결과가 모두 동일하게 계산된 것을 토대로 개발된 BBM 해석 모듈을 신뢰할 수 있을 것으로 판단되지만, 사용된 검증방법 모두 수리-역학적인 거동에 대한 검증에 국한되어 있는 한계점이 있다. 따라서 온도 변화에 따른 불포화토의 역학적 거동 변화에 대해서는 추가적으로 검증해야만 할 것으로 판단된다. 또한 개발된 BBM 해석 모듈에 대한 현장 적용성도 파악해 볼 필요가 있을 것으로 판단된다. 개발된 BBM 해석 모듈의 검증과 적용성 검토가 완료된 이후, 개발된 BBM 해석 모듈과 TOUGH2- MP/FLAC3D를 이용하여 한국원자력연구원에서 수행중인 In-DEBS(In-situ Demonstration of Engineered Barrier System) 현장시험 결과 대한 예측 및 시험 결과 분석 그리고 한국형 기준 처분시스템에서의 THM 복합거동 특성 평가 및 성능평가에 이용될 수 있을 것으로 판단된다.