1. 서 론
고준위방사성폐기물(High-Level radioactive Waste, HLW)로부터 발생하는 붕괴열(decay heat)은 고준위방사성폐기물 처분시스템 전반에 걸쳐 온도 변화를 야기한다. 그리고 붕괴열에 의한 온도변화는 처분장 주변 암반의 공극수압을 변화시키고 암반의 유효응력 변화를 야기하는 열-수리-역학적(Thermo-Hydro-Mechanical, THM) 복합거동에 영향을 준다(Fig. 1). 이러한 THM 복합거동의 변화는 고준위방사성폐기물 처분장 주변 암반의 장･단기 안전성에 영향을 미칠 수 있기 때문에 처분장 부지선정, 설계와 운영을 위해서는 처분장 주변 암반에서의 THM 복합 거동 변화에 대한 이해와 이를 분석하고 예측할 수 있는 해석 모델 및 해석 기법은 반드시 필요하다. 하지만, 이러한 THM 복합 거동은 비선형적인 거동을 보이며 구성방정식에 필요한 인자 수가 많기 때문에 해석모델과 해석기법을 통해 처분장 주변에서 예상되는 복합 거동을 정확하게 예측하는 것은 매우 어렵다. 세계 각국에서는 THM 복합 거동을 이해하고 예측하기 위해서 THM 관련 실험실 실험 및 현장시험을 바탕으로 다양한 모델링 기법을 제안하였으나 복잡한 THM 복합 거동 특성 때문에 이를 제대로 예측하는 모델링 기법은 아직 개발되지 못한 실정이다. 이러한 이유로 THM 복합거동 해석 기술 개발을 효과적으로 하기 위해 세계 각국의 모델링 전문가들이 모여 국제 공동연구 DECOVALEX(DEvelopment of COupled models and their VALidation against EXperiments) 프로젝트를 추진하면서 이러한 문제점을 해결하고 있다. DECOVALEX는 1992년 DECOVALEX-I(1992-1994)을 시작으로 DECOVALEX-II(1995-1999), DECOVALEX-III(2000-2003), DECOVALEX-THMC(2004-2007), DECOVALEX-2011(2008-2011), DECOVALEX-2015(2012-2015), 그리고 DECOVALEX-2019 (2016-2019)를 거쳐 현재 DECOVALEX-2023(2020-20203)이 진행 중에 있다(Birkholzer et al., 2019).
본 연구에서는 한국원자력연구원이 DECOVALEX-2015 Task B1을 참여하면서 수행한 1단계 연구내용과 결과를 소개하고자 한다. Task B1의 1단계는 스위스 Mont Terri Rock Laboratory에서 수행된 현장시험 HE-D(Heater Experiment–D)에 대한 수치해석을 수행하고 강한 이방성을 갖는 점토질 암반(Opalinus Clay)에서 일어나는 THM 복합거동 특성을 파악하는 것이다. HE-D 현장시험 모델링을 위해 FLAC3D(Fast Lagrangian Analysis of Continua 3D) ver. 4.0(Itasca, 2009)을 사용하였으며, 수치해석 결과를 현장시험에서 계측된 암반의 온도, 공극수압, 그리고 변형률과 비교함으로써 이방성이 강한 Opalinus Clay의 THM 복합거동에 대한 FLAC3D의 현장 적용 가능성을 검토하였다.
2. Heater Experiment-D(HE-D)
HE-D 현장 시험은 2003년 초에 프랑스 ANDRA에 의해 TER이라 불리는 히터 시험의 mock-up 시험으로 제안된 것으로 prototype의 히터 및 계측장비들의 성능을 평가하고 가열과 냉각에 의한 Opalinus Clay에서의 THM 복합거동 특성을 평가한 이후, 현장 시험 결과를 수치해석 결과와 비교함으로써 수치해석 모델을 검증하고자 하였다. 높이 4.6 m 그리고 폭 4.5 m의 단면을 갖는 Mont Terri Rock Laboratory는 알프스 산의 Jura산을 관통하는 N16 국립고속도로에 있는 Mont Terri 도로 터널의 인근에 위치하고 있으며(Fig. 2), 심도가 약 300 m인 지점의 Opalinus Clay층에 위치하고 있고(Fig. 3), SE 방향(140-150°)으로 층리가 발달해 있다. 층리의 영향으로 Opalinus Clay의 열, 수리, 역학적인 물성은 이방성을 크게 보이는 것으로 보고되었다(Wileveau, 2005, Millard et al., 2013, Garitte et al., 2014). 현장 암반의 온도는 약 15°C로 유지되는 것으로 보고 되었으며(Wileveau and Rothfuchs, 2007), 대상 지역에서의 현지응력은 최대 주응력이 수직방향으로 6 - 7 MPa가 작용하고 있으며, 중간 주응력은 Fig. 2에 나타나 있는 피난대피터널(security gallery)의 축 방향으로 4 - 5 MPa로 나타나고 있으며, 최소 주응력은 중간 주응력 방향에 수직한 방향으로 2 - 3 MPa이 작용하고 있다(Bossart and Wermeolle, 2003).
Fig. 3.
Geological profile along the motorway tunnel showing the location of the Mont Terri Rock Laboratory (Wileveau and Rothfuchs, 2007)
HE-D 시험은 갤러리(Gallery) 98에 위치하고 있는 HE-D 횡갱(niche)에서부터 굴착된 BHE-D0공에 히터를 설치하여 수행이 되었다. BHE-D0공은 MI niche로부터 약 8 m 이격된 지점에 위치하고 있다(Fig. 4). 갤러리 98과 MI 횡갱의 굴착으로 인해 야기된 굴착손상영역(excavation damaged zone, EDZ)은 1-2 m 정도인 것으로 조사되었기 때문에 HE-D 시험지역은 굴착손상영역의 영향을 받지 않은 것으로 볼 수 있다.
Fig. 4.
Layout of the HE-D experiment with heater and measuring boreholes (Wileveau and Rothfuchs, 2007)
HE-D 현장 시험은 2003년 1월부터 9월 까지 시험 설계와 준비를 하였다. 이후, 2003년 10월부터 HE-D 횡갱을 굴착하고 시험 지역에서 지질조사를 수행하였다. 2004년 1월까지 계측장비를 설치하고 2004년 3월에 히터공을 굴착하고 히터를 설치하였다. 2004년 4월 6일에 650 W로 히터를 가동 하였고 히터 패커에서는 약 1 MPa 압력이 히터공벽에 작용하도록 하였다. 첫 번째 가열을 3개월 이상 유지한 이후 2004년 7월 7일에 1950 W로 두 번째 가열을 수행하였다. 8개월 이후 2005년 3월 16일에 히터를 끄고 냉각을 시작하였다(Fig. 5).
Fig. 5.
HE-D test procedure in terms of heat power supply and packer pressure (Wileveau and Rothfuchs, 2007)
3. 수치해석 모델
본 연구에서 암반 모델 크기는 16 m × 16 m × 22 m로 결정하였으며, 암반의 열, 수리, 역학적인 이방성을 고려하여 FLAC3D에서 제공하는 이방적 열전도 모델(anisotropic heat conduction model), 이방적 유동 모델(anisotropic fluid flow model), 그리고 횡등방 탄성 모델(transversely isotropic elastic model)을 이용하였으며, 총 43,561개의 요소를 가지는 해석 모델을 생성하였다. 본 연구에서는 히터를 하나의 발열체로 모사하지 않고 Fig. 6(a)에 나타나 있는 히터의 다양한 구성 물질 반영하여 히터를 자세히 모사하고자 하였다(Fig. 6(b)). FLAC3D 모델에서 온도, 공극수압, 그리고 변위는 격자점에서 계산이 된다. 따라서, 수치해석을 수행한 이후, 계산지점 주변의 값들로부터 후처리 과정에 소요되는 시간을 줄이기 위해서는 계산지점을 격자점에 위치하는 것이 효율적이다. 이에, 본 연구에서는 해석 대상 암반에 대해 육면체 메쉬를 이용하여 전체적으로 성글게 모델을 생성시켜 전체 요소의 수를 줄인 후, 계측지점을 포함하고 있는 하나의 육면체 메쉬를 Fig. 7과 같이 6개의 사각뿔의 메쉬로 변환하여 특정 계측지점이 사각뿔의 꼭지점에 위치하도록 만들어 후처리에 소비되는 시간과 전체 해석기간을 모두 단축하였다.
수치해석에서 사용된 암반 물성과 히터 구성물질의 물성 값은 선행연구(Wileveau, 2005; Millard et al., 2013; Garitte et al., 2014)를 참고하였으며 입력 물성 값은 Table 1과 Table 2에 정리되어 있다. 수치해석에서 히터는 불투수층 매질로 가정하였으며, 역학적 물성은 암반과 동일하다고 가정하였다.
암반의 초기온도는 15°C로 설정하고 6면체의 경계면에서의 온도를 15°C로 고정하였다. 시험 구간의 초기 공극수압은 BHE-D16공에서 계측된 초기 공극수압을 근거로 1.2 MPa로 정하고(Wileveau and Rothfuchs, 2007), 6면체의 경계면에서 공극수압을 1.2 MPa로 일정하게 유지하였다. 현지응력은 최대주응력, 중간주응력, 최소주응력을 각각 6.0 MPa, 4.5 MPa, 2.5 MPa로, 경계면에서의 수직 변위는 발생하지 않는 것으로 설정하였다. 현장 시험을 위해 히터공을 굴착한 이후에 히터공 벽면에서의 초기 온도는 15°C로, 공극수압은 대기압 상태의 0.1 MPa로 고정하였고 역학적 경계조건은 히터 벽면을 자유 경계조건으로 설정하였다.
Table 1.
Properties of Opalinus clay (revised afterWileveau, 2005)
**: ⊥ means that the measurement is made according to the normal axis of the bedding plane.
Table 2.
Properties of the multi-component heating system (Wileveau and Rothfuchs, 2007)
현장 히터시험에서 히터에 가해지는 전력은 시험장치와 상당 거리 이격되어 있는 외부에서 공급되기 때문에 전력의 손실이 불가피하게 발생한다. HE-D 현장시험의 경우, 약 5 %의 전력 손실이 발생한 것으로 보고되었기 때문에(Garitte et al., 2014), 수치해석에서는 실제 공급된 전력의 95 %만 히터에 직접 입력하였으며 2004년 12월 17일에 있었던 전력 차단도 모사하였다. 반면, 가열과 동시에 히터가 암반에 가한 응력은 단순하게 1.0 MPa로 가정하고 시험이 진행하는 동안 지속적으로 가하도록 설정하였다.
4. 해석결과
4.1 온도
본 연구에서는 히터를 가동하기 시작한 2004년 4월 6일을 기준 일(reference day, Day 0)로 결정하였으며, 계측된 온도과 수치해석을 이용하여 계산된 온도를 비교하고자 두 히터에 설치된 각각 세 개의 센서와 다양한 방향으로 분포되어 설치되어 있는 암반 센서들을 선정하였다. 센서들의 자세한 위치는 Wileveau(2005)에 나타나 있으며, 본 연구에서 선정한 센서들은 히터홀의 입구를 바라보는 단면에 표시하였으며, 히터 표면에 설치된 온도 센서뿐만 아니라 암반에 설치된 온도 및 공극수압계의 위치와 이방성이 뚜렷이 나타나는 경사 45도의 층리면(bedding plane)을 함께 표기하였다(Fig. 8).
히터에서 계측된 온도를 살펴보면, 두 번째 히터에서 계측된 온도가 첫 번째 히터에서 계측된 값보다 조금 크게 나타났지만, 수치해석에서는 차이를 보이고 있지 않는 것으로 나타났다(Fig. 9). 이는 수치해석에서는 암반이 균질하지만 현장 암반은 그렇지 않기 때문에 계측된 온도가 서로 차이를 보이는 것으로써, 상대적으로 첫 번째 히터가 설치된 부근의 암반의 열전도도가 높기 때문에 두 번째 히터에서의 온도보다 첫 번째 히터에서의 온도가 낮게 계측되었을 것으로 판단된다. 히터 가동 이후 시간에 따라 온도가 증가하며 계측된 히터의 온도와 계산된 히터 온도는 유사한 패턴을 보이고 있지만, 계측 값과 계산 값은 최대 8°C 정도 차이를 보이고 있다. 이러한 온도 차이는 수치해석에서 사용된 암반의 열전도도가 실제 암반의 열전도도 보다 낮았거나, 수치해석에서 가정한 히터의 열적 물성 값에서 차이를 보이기 때문으로 판단된다.
이방성이 강한 Opalinus Clay 암반에서의 온도변화를 살펴보기 위해 층리면과 수직한 방향으로 설치된 온도 센서에서 계측된 값(Fig. 10)과 평행한 방향으로 설치된 온도 센서에서 계측된 값(Fig. 11)을 살펴보았다. 모델링에서 계산된 히터의 온도가 계측된 온도보다 높게 나타났기 때문에(Fig. 9) 암반에서 계측된 값들 중, 상대적으로 히터 가까이 설치된 센서의 계산값이 계측값보다 최대 약 3 ~ 4°C 정도의 크게 나타났지만, 히터의 가동이후 암반의 온도는 계측 값과 유사하게 증가하였으며 2004년 12월 17일에 발생한 전력의 중단으로 인한 온도 강하도 모사되었고 히터를 끄고 난 이후의 냉각기간동안의 온도 변화 역시 매우 유사하게 나타났다(Figs. 10과 11). 수치해석의 결과가 전반적으로 조금 높게 나타나고 있는 것은 히터의 온도가 전반적으로 높게 계산되었기 때문으로 판단된다. 층리면과 평행한 방향과 수직한 방향에서 모두 시간에 따른 온도변화를 잘 모사할 수 있었던 것은 사용된 암반의 열전도도가 층리와 평행한 방향뿐만 아니라 수직한 방향 역시 적절했던 것으로 판단되며, 평행한 방향과 수직한 방향으로의 비(λ// / λ⊥) 역시 적절한 것으로 판단된다.
4.2 온도변화에 따른 공극수압 변화
층리면과 수직한 방향과 평행한 방향에 설치된 공극수압계에서 계측된 데이터와 수치해석으로 계산된 공극수압 변화를 각각 Figs. 12와 13에 나타내었다. 히터가 가동됨에 따라 증가한 온도 변화(Figs. 10과 11)로 인해 발생하는 공극수압(thermally induced pore pressure)의 변화는 전반적으로 잘 재현되는 것으로 나타났다. 2004년 4월 6일에 650 W로 히터를 가동함에 따라 온도뿐만 아니라 공극수압도 증가하는 것으로 나타났고, 2004년 7월 7일에 1950 W로 두 번째 가열을 수행함으로써 증가하는 온도에 따른 공극수압 역시 적절하게 잘 모사되고 있는 것으로 나타났다. 뿐만 아니라, 히터를 끄고 냉각을 시작하면서 감소되는 공극수압의 변화 역시 수치해석적으로 재현이 가능하였다.
수치해석에서는 암반이 균질하여 동일한 입력물성을 갖는 반면, 실제 암반에서는 불연속면 또는 손상대가 존재하여 암반이 균질하지 못하기 때문에 암반은 국부적으로 서로 다른 열, 수리, 역학적인 물성을 갖는다. 일반적으로 암반에 불연속면이 존재하거나 손상대에서는 암반 투수계수는 증가하며(Bossart et al., 2002, Cho et al., 2013), 열전도도는 손상대에서 감소하는 것으로 알려져 있으며(Lee et al., 2011), 이러한 물성 차이로 인해 열-수리-역학적 복합거동은 영향을 받게 된다(NEA, 2002, Bäckblom and Martin, 1999, Bäckblom, 2008). HEDB10번 센서에서 계측된 온도와 공극수압의 변화를 살펴보면, 온도는 계측값이 상대적으로 작게 나타났으며, 공극수압은 큰 변화를 보이지 않는 것으로 나타났다(Fig. 12). 이는 암반의 불균질성으로 인해 국부적으로 암반의 열전도도가 낮아서 온도가 낮게 계측되었으며, 이로 인해 공극수압의 변화 역시 크지 않았을 뿐만 아니라, 높은 투수계수로 인해 지하수의 유동이 용이하기 때문에 압력의 생성(pore pressure buildup)이 상대적으로 어려웠을 것으로 판단된다. 추후, HEDB10번 센서가 위치한 곳에 불연속면 또는 손상영역이 존재하는지 암반에 대해 면밀한 조사할 필요가 있겠지만, 암반의 전반적인 거동과 다른 거동을 보이는 것으로 보아 HEDB10번 센서가 설치된 영역의 암반에는 불균질성이 존재할 것으로 판단된다.
4.3 변형률
암반 변형률은 현장에서 직접 계측을 할 수 없다. 따라서 현장에서 계측된 변위변화를 이용하여 암반의 변형률을 아래의 식 (1)과 같이 정의하고자 한다.
| $$\in_{i,i-1}=\frac{\delta PT_i-\delta PT_{i-1}}L$$ | (1) |
여기서 δPTi는 각각의 계측지점 i에서의 시추공 축 방향으로의 변위(mm)를, L은 각각의 계측지점 사이의 이격거리(약 1 m)를 의미한다.
Figs. 14와 15에는 보어홀 BHE-D4와 BHE-D5에서 계측된 변위를 식 (1)로 환산한 변형률의 변화와 수치해석에서 계산된 변형률이 나타나 있다. 보어홀 D04에 설치되어있는 PT01, PT02, PT03, 그리고 PT07에서 계산된 변형률은 현장에서 계측된 값과 상이하게 나타나고 있지만, PT04, PT05, PT06, 그리고 PT08에서 계산된 변형률은 어느 정도 유사하게 나타나고 있음을 알 수 있다(Fig 14). 그리고 BHE-D5 보어홀의 PT07과 PT09는 유사하지만 나머지 계측값은 계산값과 차이를 보였다(Fig. 15). 수치해석에서 암반의 이방성을 고려하였음에도 불구하고 결과 차이가 나타난 이유는 암반이 가지고 있는 불연속면과 손상대와 같은 불균질성 때문으로 판단된다. 균질한 암반으로 가정한 수치해석과는 다르게 시험구간 암반에는 불연속면과 국부적으로 다른 물성을 갖는 영역이 존재하여 암반 전체의 역학적 거동이 다르게 나타났을 것으로 보인다. 또한, 현장의 경계조건과 수치해석상에서의 경계조건이 다르기 때문에 역학적인 거동이 서로 차이를 보였을 수도 있을 것이다. 즉, 현장에서는 MI 횡갱과 갤러리 98의 영향이 존재하지만 수치해석에서는 이를 단순화하고 경계부에서의 변위가 없다는 것으로 가정하였기 때문에 결과값이 차이를 보이는 것으로 판단된다. 현장 시험을 계획하고 모델링 할 경우 이러한 경계조건에 유의해야 할 것으로 보이며, 현장 시험 설계 시 경계부에서의 변위를 항상 계측하여 정확한 경계조건을 제시하여야만 합리적인 수치해석 결과가 도출될 수 있을 것으로 판단된다.
5. 결 론
암반에서의 열-수리-역학적 복합거동을 예측하는 수치해석 기법으로 FLAC3D가 타당한지 조사하기 위해 스위스 Mont Terri Rock Laboratory에서 수행한 HE-D 현장시험에 대해 수치해석을 실시하고 현장에서 계측된 온도, 공극수압, 그리고 변형률 값과 비교하였다. 히팅 시스템과 암반에서 계측된 온도와 계산된 온도는 전체적으로 큰 차이를 보이지 않고 유사한 패턴을 보였으며, 온도 변화에 따른 공극수압 역시 수치해석적으로 잘 재현되었다. 전반적으로 온도와 공극수압의 변화가 유사하지만, 수치해석에서 계산된 온도와 공극수압이 큰 차이를 보이는 특정구간은 수치해석 모델과 기법의 문제가 아니라 암반의 불연속면 또는 물성이 크게 차이가 나는 구간으로 예상이 된다. 하지만 이러한 사실을 확인하기 위해서는 계측지점 암반에 대한 면밀한 조사가 필요할 것이다. 본 연구에서는 암반의 열-수리-역학적 이방성을 모두 고려하였음에도 불구하고, 암반에서 계측된 변형률은 계산된 변형률과 일부분에서만 유사하게 나타났다. 이러한 사실로 보아 대상 암반을 탄성적으로 예측하는 것은 무리가 있는 것으로 보이며, 추후 Opalinus Clay 암반에 대한 수치해석을 수행할 경우 이방성을 고려한 모델과 Cam Clay 모델과 같이 점토질 암반에 적합한 모델을 비교하여 최적의 모델을 먼저 도출하고 심도있는 해석을 수행하는 것이 바람직할 것으로 판단된다.
