1. 서 론
2. 수치해석 모델링
2.1 구성방정식(Governing equations)
2.2 모델 검증(Model verification)
2.3 초기조건 및 경계조건(Initial condition and boundary conditions)
3. 메타모델을 이용한 전산 실험 설계(Computer experiment design using meta-model)
4. 결과 및 분석
5. 결 론
1. 서 론
수소는 에너지의 생산과 소비를 효율적으로 연결하는 도구로 활용될 수 있다는 장점이 있다. 이에 호주, 유럽 등을 비롯한 많은 선진국에서는 친환경적인 에너지 시스템으로 전환하기 위한 수소경제 활성화 정책들을 시행·발표하고 있다(FCH JU, 2019; CSIRO, 2018). 국내에서도 2019년 ‘수소경제 활성화 로드맵’을 발표한 이후 수소 경제로의 전환에 박차를 가하고 있으며, 수소의 생산 및 공급 확대는 필연적으로 수소 물류 규모의 확대로 이어져 안정적인 중·대규모의 수소 저장 기술에 대한 개발 이슈가 논의되고 있다. 지상저장의 경우, 동일한 양의 가스를 저장하기 위해 더 많은 공간이 필요하며, 재난재해 및 사고의 위험성에 비교적 쉽게 노출되어 안전상의 우려가 존재하기 때문에 최근에는 지중저장 방식에 주목하고 있다. 지중저장은 대규모 수소 저장이 가능하여 가장 경제적이고 합리적인 저장 방식으로 인식되고 있다(Taylor et al., 1986). 지중저장기술과 관련된 국내 선행연구 동향을 살펴보면, CO2지중저장에 대한 연구들이 주를 이루고 있으며(Kang and Huh, 2008; Park and Cho, 2008; Kim et al., 2009; Kim et al., 2017; Han et al., 2017; Kim et al., 2018), 수소 지중저장에 관한 연구는 아직 미비한 실정이다. 반면 해외에서는 천연가스 및 수소가스를 대상으로 한 암염돔(salt cavern), 대수층(aquifer) 및 폐유전·폐가스전 등 다양한 지중저장방식에 대한 연구가 활발히 진행 중이다(Ozarslan, 2012; Bai et al., 2014; Simon et al., 2015; Heinemann et al., 2018; Heinemann et al., 2021). 암염돔은 두꺼운 지하 암염층(salt rock layer)에 물을 주입하여 만든 원통형 인공공동(cavity)이며 투수성이 낮아 수소 저장에 부합하는 물리화학적 특성을 지니고 있으며, 최대 깊이 2,000m, 부피 1,000,000m3, 높이 300-500m, 직경 50-100m의 규모로 건설이 가능하다(Michalski et al., 2017). 독일 및 폴란드에서는 자국의 지질학적 기준을 제시하고 수소 저장 가능성을 평가하였고(Michalski et al., 2017; Tarkowski and Czapowski, 2018; Lankof and Tarkowski, 2020), 프랑스에는 수소 손실을 야기하는 지중 내 미생물에 의한 영향을 알아보기 위한 지구화학적 모델링 연구를 수행하였다(Ebigbo et al., 2013). 이외에도 수소의 주입에 의한 암반의 유체역학적 거동을 평가하기 위한 연구들이 수행되었는데(Hagemann et al., 2015; Sainz-Garcia et al., 2017; Shi et al., 2020), 일례로 Shi et al.(2020)은 암석과 유체의 상호작용으로 인해 암석의 공극과 투수성이 변화되는 것을 실험적으로 확인하였다.
기존 선행연구들은 대체적으로 암염돔 방식에 대한 내용이 주를 이루고 있는데, 암염돔을 활용한 수소 저장시설은 극히 일부 국가에서만 운영 중이며, 우리나라는 암염돔을 건설할 수 있는 지질학적 조건을 지니고 있지 않다(Fig. 1). 국내의 지질학적인 환경에 부합하는 독자적인 수소저장방식으로서 인공구조물을 활용한 저심도 수소 저장 방식 등을 개발해야 하며(Fig. 2), 이와 관련된 안전기준 확립 및 지반 안정성 평가 연구가 필요하다. 따라서 본 논문에서는 저심도(약 10-20m) 지중 수소저장시설에서의 수소 누출에 따른 지반의 수리역학적 거동을 예측하고자 수리-역학 복합해석을 수행하였다. 이상 유동이 고려되는 벤치마크 실험결과를 통하여 해석 모델의 예측 신뢰성을 검증하였고, 메타모델을 활용한 매개변수연구를 통해 수소가스의 지반 침투에 따른 지표융기현상에 대한 영향 인자들에 대해 분석하였다.
2. 수치해석 모델링
2.1 구성방정식(Governing equations)
본 연구에서는 고압 수소기체의 지반 침투와 이에 따른 지반의 수리·역학적 응답 거동을 예측하기 위한 유한요소모델을 구축하였다. 해석모델은 다상유체흐름을 고려하는 불포화 압밀 방정식을 기반으로 하며 질량보존방정식(mass balance equation)과 운동량균형방정식(momentum balance equation)이 서로 연계되어 해석된다. 이를 위해 Biot(1962)이 제시한 다공성탄성이론(poroelasticity theory)과 불포화 지반의 열-수리-역학적 연계거동에 대해 정리한 선행연구들을 참조하여 구성방정식을 구축하였다(Lewis and Schrefler, 1998; Shin, 2011; Shin and Yoon, 2022).
본 연구에서는 지반을 지하수위 상부인 불포화조건으로 가정하여 흙 입자, 간극수, 공기로 이루어진 3상 구조로 간주하였다. 간극수압과 수소기체압을 주종속변수(main dependent variable)로 가지는 질량보존 방정식은 각각 식 (1), 식 (2)처럼 정리된다.
여기서, α는 공극압과 매질의 응력 사이 상관관계를 정의하는 Biot 계수, n는 지반의 간극률, Sw는 간극수 포화도, Sg는 침투된 수소기체의 포화도, pw는 간극수압(Pa), pg는 수소기체압(Pa), pc는 기체압력(pg)과 간극수압(pw)의 차로 정의되는 모세관 압력(Pa), Ks는 흙 입자의 체적탄성계수(Pa), Kw는 간극수의 체적탄성계수(Pa), g는 중력가속도(m/s2), ρw는 간극수의 밀도(kg/m3), μw는 간극수의 동점성계수(Pa·s), εv는 체적변형률, t는 시간(s)을 의미한다. κ는 고유투수계수(intrinsic permeability, m2)이며, Kg는 수소기체의 체적탄성계수(Pa), μg는 수소기체의 동점성계수(Pa·s)이며, krw 및 krg은 각각 불포화조건에서의 간극수 및 수소기체 흐름에 대한 상대투수계수로서 식 (3)과 식 (4)와 같이 간극수포화도(Sw) 대한 함수로 표현될 수 있다(Brooks and Corey, 1966).
한편, 식 (5) 및 식 (6)은 식 (1)-(2)에서 산정된 유동해석의 결과가 운동량균형방정식과 서로 연계되는 과정을 포함하며, 응력해석의 결과로 산출되는 체적변형률(εv)은 매 time increment 마다 질량보존방정식에 각각 반영되어 유동해석의 결과가 반복적으로 갱신된다.
여기서, DT는 tangent constitutive 텐서이며, u는 변위(m), g는 중력가속도(m/s2)를 의미한다.
2.2 모델 검증(Model verification)
많은 연구자들이 불포화 압밀 해석모델의 예측 신뢰성을 검증하기 위해 Liakopoulos(1964)가 수행한 불포화 압밀실험결과를 Benchmark 예제로 이용하여 왔다(Schrefler and Simoni, 2017; Gawin et al., 1995; Klubertanz et al., 1997; Lewis and Schrefler, 1998). Liakopoulos test는 사전에 시료를 완전히 포화시켜 상부에서 하부로 일정한 유동흐름 조건을 구현하여 초기 조건을 조성한 후, 시험 시작 단계에서 상부의 유량 유입을 중단하고 중력에 의해서 시편 하부로 수분이 서서히 빠져나가도록 유도한 상태에서 시편 내부의 간극수압, 공기압, 침하량 등을 계측한 실험이다. 본 연구에서도 Liakopoulos(1964)의 benchmark 실험 결과와의 비교를 통해 해석모델의 신뢰성을 검증하였다. Fig. 3에 나타난 바와 같이 시료의 관측 깊이에 따른 가스 압력 및 간급수압 변화에 대한 각 시간대 별 예측결과가 Liakopoulos test 실험 결과와 비교적 잘 일치함을 확인하였다. 따라서 본 연구의 사용된 수치해석 모델은 지반의 불포화 압밀거동을 신뢰성 있게 예측하는 것으로 판단된다.
2.3 초기조건 및 경계조건(Initial condition and boundary conditions)
고압 수소기체의 지반 침투과정을 모사하기 위해 고려된 해석 모델은 해석의 효율성을 위해 2차원으로 구현되었다. 유동해석에 대한 초기조건은 전체도메인을 지하수위 상부조건인 불포화상태로 간주하였다. 또한, 도메인 좌측경계의 임의의 지점에 수소가스침투압을 시간에 대한 함수형태(pg=f(t))로 적용하였고, 도메인 하단은 No flow 조건(-ρgu=0)을 부여하였다. 또한, 지반을 통과하여 지표면 밖으로 방출되는 수소가스의 양은 미미할 것으로 판단하여 지표면 상부 또한 No flow 조건(-ρgu=0)을 부여하였다. 유동해석에 대한 초기조건은 지반의 자중에 의해 평형상태에 도달한 static 상태로 설정하였다. 도메인 하단은 x축과 y축 방향에 대해 변위가 구속된 Fixed 경계조건을, 우측경계에는 x축 방향에 대한 변위가 구속된 Roller 경계조건을 설정하였다. 도메인의 좌측경계는 Symmetry 조건을 부여하였다. 상기 기술된 구성방정식, 초기조건 및 경계조건을 적용하기 위해 상용 수치해석 코드인 COMSOL Multiphysics을 이용하였다(Comsol Inc., 2022). 해석에는 총 1,750개의 사각요소망(rectangular mesh)이 사용되었고, 해석시간은 총 100일로 설정하였다. 모델에 사용된 재료물성은 Table 1과 같다.
Table 1.
Material properties used in numerical simulation model
3. 메타모델을 이용한 전산 실험 설계(Computer experiment design using meta-model)
실험 및 모델에서 종속변수 변화에 관여하는 영향 인자 간의 관계식을 정확하게 구현하는 것은 현실적으로 불가능하기 때문에 메타 모델을 이용한 실험설계를 이용하여 매개변수연구의 효율성을 높이는 방안이 제시되고 있다(Booker et al., 1999; Go et al., 2020). 실제 시스템을 대신하여 표현한 모델인 메타모델은 입력 데이터와 출력 데이터 사이의 간단한 관계를 나타내는 수학적 방정식 혹은 알고리즘이다(Booker et al., 1999; Hoffman et al., 2003). 선행연구에서 제시된 다양한 메타모델 방식 중에서 Gaussian Kriging 방식은 적당한 수의 매개변수를 가지는 비선형적인 결정론적 모델에 적용하기에 가장 적합하다고 알려져 있으며, 다음과 같은 식으로 표현될 수 있다(Matheron, 1963).
여기서 βj는 미지 계수; Bj(x) j=1,…,L는 데이터에 값에 근거한 기준값, z(x)는 오차 값이다. Gaussian 상관함수인 r(θ; s, t)는 식 (8)을 통해 메타모델의 오차를 산정한다.
식 (8)에서 θk는 평활 파라미터, n은 디자인 변수의 개수, 그리고 sk와 tk는 샘플의 k번째 요소들이다. y(x)에 대한 Krging 예측값 는 다음과 같이 구할 수 있다.
여기서, R은 상관행렬이고 r'(x)={r(θ; x1, x), …, r(θ; xn, x)}' 이다. 본 연구에서는 다섯 종류의 영향 인자 xid에 대하여 50-300개의 세트를 도출하였고, 이러한 입력변수들을 수치해석 모델에 대입하여 데이터 세트를 구성하였다.
이때, max와 min은 입력변수의 최대 및 최솟값이며(Table 2), rand(0,1)은 Latin Hypercube Sampling을 통하여 0과 1 사이에서 균등하게 추출된 임의의 값이다.
Table 2.
Range of independent variable used in the sensitivity analysis
| Parameter | Minimum | Maximum | Unit |
| Leakage depth | 1 | 6 | m |
| Intrinsic permeability | 1.573E-14 | 1.017E-12 | m2 |
| Young’s Modulus | 15 | 100 | MPa |
| Leakage pressure | 150 | 610 | kPa |
| Leakage duration | 5 | 80 | day |
Fig. 5는 일반적인 Random sampling 함수를 이용하여 임의로 추출된 데이터세트의 분포와 Latin Hypercube Sampling을 통하여 추출된 데이터세트의 분포의 차이를 보여준다. 일반적인 Random 함수를 이용한 샘플링 기법의 경우에는 특정 매개변수 상호간 overlapped points가 나타나지만 Latin Hypercube Sampling 기법의 경우에는 이러한 현상 없이 대체적으로 균등(uniform)한 분포를 만들어냄을 알 수 있다. 특히, 매개변수의 개수가 많이 늘어날수록 Random 함수를 이용한 샘플링에서는 overlapped points 현상이 두드러지게 나타나기 때문에 Latin Hypercube Sampling의 효과는 더욱 극대화된다.
4. 결과 및 분석
Fig. 6은 2.3절에서 정의된 수소누출모델에서 도출된 해석 결과이다. 고압수소가 지반으로 침투되는 시점에서의 지반 내부의 수소압, 간극수압, 수소기체 포화도, 간극수 포화도 값 분포를 보여준다. 누출지점에서 증가된 수소압 및 간극수압이 주변지반으로 전파되는 양상을 나타내었다. 또한, 수소기체 포화도는 누출 전 43%에서 누출 후 최대 55%까지 증가한 반면 간극수 포화도는 누출 전 57%에서 누출 후 45%로 감소하는 양상을 확인하였다. 불포화지반의 수리학적 특성변화에 따른 역학적 거동변화가 예상되어 Fig. 7과 같이 지표면 상부에서의 최대변위를 야기하는 위치에서의 지표변위(displacement)를 평가하였다. 누출 시작시점부터 누출이 중단되는 시점까지 지표면 변위는 로그함수개형으로 증가양상을 나타내었고, 누출 후에는 다시 급격히 처음상태로 회귀하는 변화를 나타내었다. 누출 시 지표면 변위 증가는 고압수소에 의한 상향 침투압에 기인한 것으로 판단되지만, 이 외에도 지반의 탄성계수(E, Young’s modulus), 절대투수계수(intrinsic permeability) 및 누출깊이, 침투압의 크기 등 다양한 요인에 따라 영향을 받을 수 있다.
Fig. 8은 고압수소가스의 누출 깊이에 따른 지반 내부의 수소압 분포 및 수소가스 침투흐름 양상을 보여주며 Fig. 9는 누출 깊이 별 수소가스 상향침투에 의한 지표면의 최대변위를 보여준다. Fig. 9에서 알 수 있듯이 수소가스 누출 깊이가 약 3-4m 일 때 지표면 변위가 가장 크게 나타남을 알 수 있다. 저심도 수소저장소의 방호 구조체의 크기를 감안할 때, 저심도 저장시설의 설치심도(installation cover depth)가 4m 이상으로 깊어질수록 수소침투에 의한 지표면 변위 정도가 줄어들 것으로 판단된다.
한편, 본 연구에서는 메타모델을 활용한 매개변수연구를 통해 수소가스의 지반 침투에 따른 지표융기현상에 대한 영향 인자들의 민감도에 대해 분석하였다. 이를 위해 Table 2과 같이 주요영향인자들(지반의 탄성계수 및 절대투수계수, 누출깊이, 누출압의 크기, 누출지속시간)을 선정하였고 최대 300개의 데이터세트를 구축하여 수치해석을 수행하였다. Fig. 10은 매개변수연구의 결과를 보여준다. 이때, 독립변수와 종속변수는 각각 다른 스케일의 범위를 지니고 있기 때문에 최대 및 최소값을 기준으로 0에서 1사이 값으로 정규화(normalize)하여 표현하였다. 누출깊이의 경우 다른 독립변수들이 통제되어 있을 때에는 특정 깊이에서 지표변위가 최대로 나타나는 양상을 보였지만(Fig. 9) 다른 독립변수들이 통제되어 있지 않은 경우에는 침투깊이와 지표변위 사이에 특별한 상관관계가 관측되지 않았다. 반면, 절대투수계수와 침투압의 크기, 침투지속시간은 종속변수인 지표변위와 양의 상관성이 확인되었다. 지반의 탄성계수(E, Young’s modulus)는 지표변위증가에 대해 오히려 반비례하는 경향을 보였기 때문에 탄성계수의 역수(1/E)와는 상대적으로 높은 양의 상관성을 확인하였다. Fig. 11에서 알 수 있듯이 종속변수인 지표변위는 고려된 독립변수들과 서로 밀접하게 연관되어 영향을 받기 때문에 선형성을 보이는 상위 3개의 인자들에 대해서도 회귀식에 대한 분산정도가 크게 나타나는 것으로 분석되었다. 그럼에도 본 연구에서 수행된 매개변수연구는 지반의 탄성계수가 작을수록, 절대투수계수가 클수록, 수소가스 누출압력이 크고 누출 지속시간이 길어질수록 수소가스의 지반누출에 따른 지표변위가 크게 나타나며 그 영향에 대한 민감도는 지반의 탄성계수가 상대적으로 가장 크다는 것을 확인시켜주었다. 또한, 메타모델 구성에서 Latin Hypercube Sampling 기법을 사용할 시 매개변수연구의 결과에 대한 재현성이 유지하는 최소의 데이터세트의 개수를 확인하고자 데이터세트의 개수를 줄여가며 기준 치 대비 민감도의 상대오차를 평가하였는데, Fig. 12와 같이 데이터세트의 개수가 최소 50개 이상이 되면 어느 정도 일관성 있는 민감도를 산출할 수 있는 것으로 확인되었다. 따라서, Latin Hypercube Sampling 기법을 이용한 메타모델 구성 시 최소 50개 이상의 데이터세트를 구축한다면 전체 경우의 수를 대변하는 결과를 얻을 수 있을 것으로 판단된다.
5. 결 론
본 연구에서는 저심도 지중 수소저장시설에서의 수소 누출에 따른 지반의 수리역학적 거동을 평가하고자 수리-역학 복합해석을 수행하였다. 벤치마크 실험결과와의 비교를 통해 예측 신뢰성이 검증된 해석모델을 통해 고압수소가 지반으로 침투되는 시점에서 발생하는 지반 내부의 수소압, 간극수압, 수소기체 포화도, 간극수 포화도 변화 등을 확인하였다. 또한, 메타모델을 활용한 매개변수연구를 수행하여 수소가스의 지반 침투에 따른 지표융기현상에 대한 영향 인자들의 민감도에 대해 분석하였다. 이를 통해 도출된 결론은 다음과 같다.
(1) 저심도 지반에 설치된 지중 수소저장시설에서 고압수소가스가 지반으로 누출 될 경우, 불포화지반에서는 수리학적 특성변화에 따른 역학적 거동변화가 일어나며, 누출 시작시점부터 누출이 중단되는 시점까지 상향 침투압에 의해 지표면 변위가 로그함수개형으로 증가하는 양상으로 나타났다.
(2) 고압수소가스의 누출깊이는 다른 영향인자들이 통제되어 있을 때에는 특정 깊이(약 3-4m)에서 지표변위가 최대로 나타나는 양상을 보였지만 다른 영향인자들이 통제되어 있지 않은 경우에는 침투깊이와 지표변위 사이에 특별한 상관관계가 나타나지 않았다.
(3) 고압수소가스의 지반 침투 시 지반의 절대투수계수와 침투압의 크기, 침투지속시간은 상향 침투압에 따른 지표변위 변화와 양의 상관성이 확인되었지만, 지반의 탄성계수(E, Young’s modulus)는 지표변위증가에 대해 오히려 반비례하는 경향을 나타내었다.
(4) 고압수소가스 지반 침투에 따른 지표융기현상에 대한 매개변수연구 결과, 지반의 탄성계수가 작을수록, 절대투수계수가 클수록, 수소가스 누출압력이 크고 누출 지속시간이 길어질수록 수소가스의 지반누출에 따른 지표변위가 크게 나타나며 그 영향에 대한 민감도는 지반의 탄성계수가 상대적으로 가장 큰 것으로 확인되었다. 이러한 연구결과는 향후 수소가스 누출뿐만 아니라 수소가스 폭발에 대한 지반 복합해석 평가 시 유용한 기초자료로 활용될 것이다.
