1. 서 론

호수나 저수지 등의 사면에서 발생하는 산사태(landslides) 및 토석류(debris flows)는 수면과 접촉하면서 수면충격파(impulse wave)를 발생시킬 수 있으며, 이 파의 전파는 제체(堤體) 및 반대편 사면에서의 파 쳐오름(wave run-up), 댐 제체에서의 월파(over-topping) 등으로 인해 심각한 피해를 유발할 수 있다(Fritz et al., 2003; Panizzo et al., 2005). 이러한 현상은 저수지 사면에서 산사태 및 토석류의 발생, 사면에서의 전파뿐만 아니라 수면 부근 가까운 영역(near field)에서 수중으로 토석류 파의 진입 및 수면충격파 발생, 저수지 바닥에서의 수중 토석류 전파, 수면충격파의 먼 영역(far field)으로의 전파, 먼 거리 영역에서 파의 쳐오름 및 월파 등의 종합적인 사건으로 구성된다. 이 현상은 복잡한 지형학적 조건에서 연직방향으로 토석류 층과 물 층이 거동하는 3차원 자연현상으로서, 저수지 수중 바닥을 따라서 이동하는 토석류와 수면에서 발생한 충격파를 동시에 고려해야만 정확한 현상과 이로 인한 피해를 예측할 수 있다.

산사태와 토석류로 인한 충격파 현상을 해석하는 연구는 실험적 연구와 수치해석 연구로 나눌 수 있다. 실험 연구는 특정 지형에서의 현상 규명으로 목적이 국한되며, 축적효과(scale effect)로 인해서 대부분 단순한 실험수조에서의 현상 규명을 중점적으로 수행되고 있다. Slingerland and Voight (1979) 그리고 Huber and Hager (1997)는 물리 모형실험을 통해서 산사태로 인한 충격파의 파고를 추정할 수 있는 관계식을 제시하였다. Synolakis (1987)와 Huber and Hager (1997)는 산사태 발생 사면의 반대편에 있는 사면이나 저수지 댐 제체에 도달한 후 발생하는 파의 쳐오름을 예측할 수 있는 관계식을 개발하였다. Akgün (2011)은 이들 경험공식을 이용하여 터키의 Kurtun 댐의 저수지 유역에서 발생한 산사태에 의한 수면충격파로 인한 피해를 평가한 바 있다. 하지만, 이러한 경험공식에 근거한 충격파 전파 현상의 예측은 단순 지형에만 적용이 가능한 한계가 있다. 주목해야 할 것은 속도는 다르지만, 수면충격파뿐만 아니라 수중 토석류도 전파되어 궁극적으로 댐 제체 또는 산사태/토석류가 발생한 반대편 사면에서 쳐오름 및 월파가 발생할 수 있다는 것이다. 따라서 이러한 충격파와 토석류의 전파 현상을 지배하는 토석류와 물의 경계면에서의 혼합(mixing), 토석과 물의 상호작용, 토석류 및 물과 댐 구조물과의 상호작용 등을 높은 정확도로 해석하면서 실제 문제에 적용성이 높은 3차원 수치모델의 개발 및 적용이 요구된다.

유체와 고체 혼합물의 큰 규모 운동(large scale movements)은 산사태와 토석류 흐름, 충격파 전파, 댐 붕괴 등을 포함하는 자연과학 분야의 많은 문제에서 중요한 역할을 한다. 전 세계적으로 저수지, 호수, 해안, 해양 그리고 하천 및 댐 공학자들에게 특히 주목을 받는 것 중 하나는 산사태와 토석류에 의해서 유발된 높은 파고와 빠른 파속을 갖는 수면충격파의 전파 현상이다. 소규모 쯔나미라고도 할 수 있는 이들 수면충격파는 상대적으로 짧은 파장(wave lengths)과 높은 진폭(amplitude)이 특징이며, 수역(water basin)의 경계에서 가공할 높이의 파 쳐오름을 발생시킬 수도 있다(Akgün, 2011). 이러한 현상들은 침식, 지진에 의한 진동, 폭풍, 집중호우 또는 인공저수지에서의 수문 조작에 따른 수위 변동 등과 관련된 자연재해에 의해서 시작될 수 있다. 결과적으로, 복잡한 저수지의 지형을 고려해야 하고 저수지 수중 바닥을 따라서 이동하는 토석류와 수면에서 발생한 충격파를 서로 다른 시간 및 길이 규모로 동시에 3차원적으로 고려해야만 정확한 현상과 이로 인한 피해를 예측할 수 있다. 국내 토석류 전파에 대한 수치모의 연구는 대부분 실용적인 2차원 수치해석 연구에 국한되어 있으며, 바닥물질의 연행을 고려하지 못하는 상용 프로그램을 이용한 연구(Kim et al., 2013)을 시작으로 이후 바닥물질 연행에 따른 토석류 거동 특성을 분석하는 연구들(Kim et al., 2015; Kim et al., 2020; Lee et al., 2020)이 주를 이루고 있다.

3차원 수치해석 연구는 대상 영역의 지형학적 특성을 합리적으로 고려할 수 있으며, 수리모형 실험에서 발생하는 축척효과 문제가 없는 실제 수면충격파 해석에 적용되고 있다. 하지만, 비뉴튼 유체인 토석류와 뉴튼 유체인 물 및 공기의 유동 특성 그리고 이들 서로 다른 유체 흐름의 상호작용을 해석할 수 있는 수치모형를 개발하기는 쉽지 않다. Pastor et al. (2009)는 a) 토석류 흐름은 유체화된(fluidized) 토양 유동모형을 병합한 수심평균 모형으로 해석하고, b) 가까운 영역에서 토석류와 수면파와의 상호작용은 다상흐름 해석을 위한 레벨셑(level set) 알고리즘을 이용한 후, c) 마지막으로 먼 영역으로 충격파의 전파는 수심평균 모형을 이용하였다. 이 연구 결과에 의하면 산사태 및 토석류 입자들이 물로 침투하는 현상을 정확하게 재현하는 것이 초기 수면충격파 발생에 매우 중요하게 영향을 미치며, 토석 입자와 물과의 상호작용을 재현하기가 쉽지 않다. Pastor et al. (2009)의 수치해석 결과는 수면충격파의 월파를 크게 과소평가하고 있으며, 토석류는 과도하게 반대편 사면으로 전파됨을 보여준다. Abadie et al. (2010)은 물, 공기, 토석류를 3개의 유체로 고려하고 VOF (volume of fluid) 알고리즘을 이용하여 Navier-Stokes 방정식을 해석하는 진보된 모델을 개발하여 저수지에서 산사태에 의한 충격파의 전파 현상을 모의하였다. 이들은 VOF 기법을 이용하여 토석류로 인한 충격파의 영향을 예측할 수 있음을 보여주었지만, 토석류가 물로 침입하는 현상을 정확하게 모의하는 것이 전체 해석의 정확도를 확보하는 데 매우 중요하며, 이에 대한 심도 있는 추가 연구가 필요함을 보여주었다. 이들 연구는 토석류를 뉴튼 유체로 단순화하여 고려함으로써 전체 충격파 해석 결과에 한계가 있다. Ataie-Ashtiani and Yavari-Ramshe (2011)는 4차 정도(order) Boussinesq 근사화(approximation)를 이용한 수심평균 2차원 모형을 개발하여 이란(Iran)의 Maku와 Shafa-Roud 두 개 댐에서 발생한 산사태로 인한 저수지 내 파의 전파를 모의하고 발생 파고, 파 쳐오름, 댐 정상부에서의 최대 파고, 월파 유량 등을 평가하였다. 하지만 이들은 2차원 수심평균 모형을 이용하여 깊은 수심에서 바닥면을 따라 흐르는 토석류의 전파 영향과 수면 부근에서 발생하는 파의 전파를 각각 반영하여 모의할 수 없으므로 계산의 정확도에 한계가 있다. Vacondio et al. (2013)은 계산격자를 사용하지 않는 라그랑지 기법을 이용한 3차원 SPH (smoothed particle hydrodynamics) 기법을 적용하여 1963년 약 2500명의 사망자를 발생시킨 이탈리아 Vajont 저수지 사면에서 발생한 산사태/토석류가 야기한 수면충격파가 콘크리트 댐 제체를 월류하는 현상을 수치모의였다. 시뮬레이션 결과는 파의 최대 쳐오름과 저수지에 남아 있는 물의 최종 수위로 적용성을 평가하였다. 이 연구는 Navier-Stokes 방정식을 해석하여 인공저수지에서의 토석류와 물의 흐름 거동을 재현한 3차원 수치해석 연구이다. 이들 연구에 따르면 저수지 내 충격파와 물의 흐름 거동에 가장 큰 영향을 미치는 영향인자는 산사태/토석류 흐름의 지속시간인 것으로 나타났다. 이 결과는 충격파 모의 시 선단부의 유속, 선단부 높이, 입수 시 토석류와 물의 상호작용, 수중에서 토석류의 전파 유속 및 이동 거리 등 산사태/토석류의 흐름 특성을 물리적이고 합리적으로 재현하는 것이 전체 충격파의 전파 및 이로 인한 피해를 예측하는 데 매우 중요함을 보여준다. 하지만, SPH기법에서 경계값 문제(boundary value problem)의 강한 공식화(strong formation)는 큰 오차를 발생시키는 문제가 있으며, 이에 따라 일반적으로 SPH기법은 수치기법의 안정성(stability)과 견고성(robustness)을 증가시키기 위해서 약한 공식화된(weak formulated) 평형방정식을 이용하게 되고 계산시간도 더 증가하는 단점이 된다. 따라서 비뉴턴(non-Newtonian) 유체인 토석류 그리고 부분적으로 과농도류와 뉴튼유체의 특성을 동시에 보일 수 있는 물 흐름의 상호작용을 물리적으로 합리적인 근거를 갖고 유도된 지배방정식을 이용하여 해석할 수 있는 3차원 다상흐름(multi-phase flow) 수치해석 모형의 적용이 필요하다.

이 연구에서는 토석, 물, 공기를 고려하는 3차원 다상 난류 흐름 해석을 위한 수치모형을 개발하고, 토석류 유발 수면충격파에 대한 Fritz et al. (2001, 2003)의 수리실험 자료를 재현함으로써 개발된 수치모형의 적용성을 평가한다. 토석류의 유변학적 거동(reological behavior)을 모의하기 위해서 수정된 Herschel-Bulkley 모형을 적용한다.

2. 수치해석 방법

2.1 지배방정식

수치해석의 기본 개념은 서로 다른 유체 간의 경계면 거동 및 상호작용을 모의하기 위해서 VOF 기법을 이용한다. 흐름 운동을 나타내는 지배방정식은 Navier-Stokes 운동량방정식과 연속방정식이다. VOF 기법은 각 유체 상(phase)에 대한 체적분율(volume fraction) 변수의 이송방정식을 해석하고 각 유체의 체적분율에 근거한 가중평균을 이용하여 각 유체의 물리적 특성을 해석하게 된다. 이 연구에서는 3상(three phases)의 유체, 즉 공기, 물 그리고 토석 혼합물을 고려한다. 이 논문에서는 지면의 제약으로 인해서 지배방정식을 가능한 간략한 벡터 형식으로 표현하고, 이미 잘 알려진 난류방정식은 지배방정식의 기술 없이 간략한 설명만을 포함하기로 한다.

운동량방정식은 다음과 같이 표현된다.

여기서, U는 유체의 유속, ρ는 유체의 밀도, μ와 μt는 유체의 점성계수(dynamic viscosity)와 난류점성계수이다. 그리고 F_s는 외력항(external force term)이다. 이 식과 일반적인 Navier-Stokes 방정식과 차이는 총 압력(total pressure) 대신에 피에조 압력(piezometric pressure), prgh=p-ρg를 이용한다는 것이다.

이 연구에서 물은 비압축성 유체로 고려(ρ = 상수)되며, 연속방정식은 다음과 같이 간략화된다.

난류 해석을 위해서는 큰와모의(large eddy simulation, LES)를 수행하였으며, 하부-격자 규모(sub-grid scale, sgs) 모형으로 표준 Smagorinsky 모형을 이용하였다. 비뉴튼 유체인 토석류의 흐름에 대한 전단응력은 Papanastasiou (1987)의 방법을 이용하여 수정된 Herschel-Bulkley (HB)모형을 적용하였다.

여기서 τ는 전단응력(shear stress), τy는 항복응력(yield stress), K는 비뉴튼 유체의 컨시스턴시(consistency) 지표, γ˙는 흐름 변형율, m은 Papanastasiou 방법에서 이용하는 매개변수, n는 멱지수이다. 일반적으로, 멱지수 n = 1로 설정하면 Eq. (3)은 Bingham 관계식에 해당하고, Herschel-Bulkley 관계식의 경우 진흙 흐름(mudflow)에 대해서는 n = 0.33 그리고 토석류 흐름에 대해서는 n = 0.5의 값이 이용된다.

하나의 계산 셀에서 VOF은 Fvol=αiVcell로 계산되며, 여기서 Vcell은 계산격자에서 단위 셀의 체적이고, αi는 계산셀에서 유체 i의 유체분수이다. 각각의 분수 α는 0과 1사이의 값을 가지며, 이 연구에서는 α1을 공기, α2를 물 그리고 α3를 산사태/토석류 구성 물질로 설정하였다. 예를 들어 하나의 셀이 완전히 토석류로만 채워져 있으면 α3은 1이며, 완전히 공기로만 채워져 있으면, α2와 α3는 0의 값을 갖는다. 이 스칼라 함수인 각 α 유체에 대한 이송방정식으로 계산된다.

여기서, Cα는 분리흐름(segregated flow)의 정도를 결정하는 계수이다. 실제 계산에서는 두 개의 유체에 대한 운동량방정식과 연속방정식을 풀고 나머지 유체에 대해서는 다음과 같이 산정한다.

운동량방정식에 포함된 표면장력은 다음과 같이 계산된다.

여기서 σ는 표면장력계 계수이며, n은 계산셀 경계면에 대한 수직 단위벡터로서 다음과 같이 정의된다.

그리고 θx=∇n는 셀 경계면의 곡률이다.

2.2 수치해석 기법

지배방정식의 수치해석을 위한 기본 도구로는 오픈소스 CFD (computational fluid dynamics) Toolbox인 OpenFOAM (OpenFOAM, 2021)에 포함되어 있는 ‘multiphaseInterFoam’을 이용하였다. 이 연구에서 기존 프로그램을 수정하여 개발한 부분은 비뉴튼 유체인 토석류의 유변학적 거동 해석을 위해 Papanastasiou (1987)의 수정된 Herschel-Bulkley 관계식 모듈을 추가로 개발하여 적용한 것이다. 지배방정식은 2차 정확도의 유한체적법(finite volume method)을 이용하여 이산화하였다. 운동량방정식에서의 이송항은 유계중앙차분(bounded central differencing)의 일종인 Gamma 기법(Jasak et al., 1999)을 이용하였으며 정확도 조절계수 ψ = 0.1을 적용하여 수치분산(numerical dissipation)을 최소화하였다. 시간항은 2차 정도의 후방차분(backward difference)기법을 이용하였다. VOF 지배방정식의 이송항은 van Leer의 TVD (total variation diminishing) 기법을 이용하여 이산화하였으며, 기타의 항들은 중앙차분기법을 이용하여 이산화하였다. 비압축성 흐름의 시간 종속(time-dependent) 계산에서 압력-유속 조합(pressure-velocity coupling)을 위해서는 PISO 알고리즘을 이용하였으며, 각 계산시간 단계에서 수정 반복 계산 횟수는 5회로 설정하였다.

3. 수치해석 대상

토석류와 수면충격파의 상호작용을 해석하기 위한 3차원 다상흐름 수치모형을 검증하기 위해서 미국 알라스카 주 남쪽 연안에 위치한 Lituya Bay의 상류에 있는 Gilbert Inlet에서 1958년 7월 10일 발생한 산사태로 인한 수면충격파를 수치모의하고 결과를 실험자료와 비교하여 분석하였다. 사고 당시 산사태는 매우 빠른 유속으로 수면에 충격을 가해서 기록적인 파의 쳐오름 높이를 기록하였다. Fig. 1에서 보인 오른쪽 사면에서 발생한 토석류가 수면에 충격을 가하여 대규모 수면충격파를 발생시켰고 이 수면충격파는 반대편에 있는 산을 EL. 524 m 높이까지 넘어가면서 사면을 침식시키는 사건이 발생하였다.

이 사건을 재현하기 위해서 Fritz et al. (2001, 2003)은 스위스연방공대에 특수 제작된 장치를 이용하여 수리모형실험을 수행하였다. 이 실험은 1:675 축척으로 만들어진 단면에서 Froude 상사법칙에 근거해서 실시하였다. 폭이 일정한 2차원 형상의 모형수로의 길이, 폭 그리고 수심은 각각 11 m, 0.5 m 그리고 1.0 m이다. 사면길이는 총 3 m이며 사면경사는 45°를 적용하였다. 실제 Gilbert Inlet에서 발생한 산사태의 폭은 823 m이고 단위폭당 총 토석량은 37,200 m³/m 이다. 실험에서 사용한 수로의 제원과 단면 모식도는 Fig. 2와 같다. 사건 발생 당시 Gilbert Inlet에서의 수심 h은 약 122 m이다. 산사태의 동적 활동 충격(dynamic slide impact) 특성은 공기압력 산사태 발생기(pneumatic landslide generator)를 이용하여 제어하였다.

수리모형실험에서는 레이져 거리센서(LDS), 입자 영상 유속계(PIV), 파고계(CWG)를 이용하여 각각 토석류의 두께, 단면 유속벡터, 수면충격파의 파고 변화를 관측하였다. 실험에서는 비중 2.64, 직경 4 mm의 입자들을 이용하여 전체밀도 ρs = 1.62 t/m³, 유효 내부 마찰각 ϕs = 43°, 동적 바닥 마찰각 δ = 24°의 토석류를 이용하였다. 실험에서 발생시킨 최대 토석류의 두께는 134 m이고 수면을 충격하기 전의 토석류의 평균 길이는 748 m이다. 토석류의 평균 유속은 110 m/s이며, 토석류가 입수하는데 소요되는 시간은 약 6.8초이다. 결과적으로 토석류의 유속과 Gilbert Inlet에서의 수심에 근거한 Froude 수는 3.18이다.

Fig. 1.

Illustration of landslide dimensions in Gilbert Inlet and landslide-induced-tsunami runup to 524 m on spur ridge directly opposite to landslide impact [adapted from Fritz et al., 2001]

Fig. 2.

Schematic of cross section of Gilbert Inlet along slide axis in NE to SW: Geometry correnponsindg to simplified physical model configurations [adapted from Fritz et al., 2001]

4. 적용 및 결과

4.1 적용 조건

Fritz et al. (2001, 2003)가 이용한 폭이 일정한 2차원 형상의 수리모형실험 수로를 동일하게 수치적으로 재현하기 위해 구성한 계산격자의 총 계산셀 수는 약 2.8 × 10⁶개이며, 격자의 형상은 Fig. 3과 같다. 수치해석에서 수로 바닥과 측벽에 대해서는 비활조건(no-slip condition) 그리고 상부 경계에서는 대기압 조건을 설정하였다. 수로의 하류단에서는 실제 표고 524 m를 넘어 월파한 물이 자연스럽게 배출되도록 무경사(zero-Gradient) 조건을 설정하였다.

Fig. 3.

Computational mesh

Fritz et al. (2001, 2003)은 모형실험에서 산사태로 토석류 흐름을 재현하기 위해서 직육면체의 상자(box)에 토석 입자를 넣고 압축공기로 상자를 가속시킨 후 사면에 배출하는 공기압력 산사태 발생기를 이용하였다. 본 수치해석에서는 가속된 상자에서 사면으로 배출된 토석 물질의 초기 형상을 Fig. 4와 같이 수치적으로 적용하였다. Fritz et al. (2001)이 실험에서 재현한 거동 전(pre-motion) 산사태의 두께는 92 m이었으며, 사면을 흘러내리면서 최대 두께가 134 m까지 점진적으로 증가하는 것으로 나타났다. 그리고 수면과 충돌 전 산사태의 길이는 약 748 m이고 속도는 약 110 m/s이였다. 이 연구에서는 Fritz et al. (2001)이 재현한 길이 748 m와 두께 92 m의 산사태를 ‘Case A’로 명명하여 수치해석을 수행하고, 사면에서 두께가 최대로 증가한 134 m 두께와 이에 상응하는 길이 514 m의 산사태를 ‘Cace B’로 명명하여 수치모의를 수행하여 산사태의 두께가 수면충격파의 발달 및 전파에 미치는 영향을 분석한다. 두 가지 수치모의 조건을 정리하면 Table 1과 같으며, 그림으로 표현하면 Fig. 4와 같다.

Fig. 4.

Initial conditions of relatively thick (Case A) and thin (Case B) landslides on the slope and water in the bay

산사태로 인해 발생한 토석류의 유변학적 거동을 재현하기 위해 이 연구에서 적용하는 수정된 Herschel-Bulkley 모형에서 전단응력을 지수분포형으로 증가하도록 이용되는 매개변수 값은 m = 1,000으로 일반적인 값(Mitsoulis, 2007)을 적용하였으며, 거동을 제어하는 멱지수 n의 값은 이류(mud flow)보다는 사석을 포함하는 토석류(stony debris flow)에 가까움으로 이에 해당하는 0.5값을 적용하였다.

4.2 결과 분석

두 가지 산사태 발생 조건을 이용한 수치모의를 통해 구한 일련의 토석류파와 수면파의 전파 순간들을 이에 상응하는 Fritz et al. (2001)가 촬영한 실험영상과 비교하면 Fig. 5와 같다. 이 그림에서 실험 영상들은 산사태로 인해 발생한 높은 유속의 토석류에 의한 수면 충격, 이로 인한 수면충격파의 생성과 전파 그리고 맞은편 사면에서의 쳐오름과 관련된 큰 비정상성(unsteadiness)를 갖는 일련의 과정 중에 복잡한 현상들이 나타남을 보여주고 있다. 그리고 두 가지 초기조건을 이용한 3차원 수치해석은 재현된 수면충격파의 생성, 전파 그리고 쳐오름의 형태를 부분적으로 차이는 있지만, 전반적으로 실험에서 나타난 현상을 양호하게 예측하는 것으로 나타났다.

Figs. 5(a) and 5(b)의 실험자료에서 토석류는 수면에서의 충격 그리고 바닥에서의 굴절로 인해 팽창(bulking)되는 것을 볼 수 있다. 아울러 토석류의 어깨(shoulder)에 해당하는 선단부 후면에서는 토석류가 수체(water body)로 빠르게 침입해 들어가면서 발생한 흐름분리로 인해 토석류 후면에 공동(air cavity)이 발생함을 보여준다. Figs. 5(a) and 5(b)의 수치해석 자료를 보면 Case A의 경우 토석류 선단부의 전파는 실험자료와 잘 일치하지만, 상당량의 토석류가 아직 사면에 남아 있으며, 수면충격파의 높이를 상대적으로 작게 모의한다. 상대적으로 두꺼운 토석류 층을 초기조건으로 설정한 Case B는 대부분 토석류가 바닥면에 도달하였으며, 수면충격파의 높이가 상대적으로 높게 그리고 실험자료를 더 유사하게 모의하는 것으로 나타났다. 아울러 Figs. 5(a) and 5(b)를 비교해보면 두꺼운 토석류가 수중으로 유입되는 경우 수면충격파의 수면이 상대적으로 오랫동안 높이 유지되는 것으로 나타났다. 이러한 결과는 토석류의 유입 두께가 수면충격파의 초기 생성 형태와 높이에 큰 영향을 미침을 보여준다.

Figs. 5(c) and 5(d)에서 실험자료는 토석류 층과 물 층 사이에 흐름분리가 발생하는 동안 물, 토석류 그리고 공기의 3개 상은 분명한 경계선으로 구분이 되며, 토석류가 바닥을 따라 전파하는 동안 공동이 붕괴하면서 물과 공기 사이에 강한 혼합이 발생함을 보여준다. 아울러 토석류가 수평 바닥의 중간지점을 넘어서면서부터는 수면충격파의 선단부가 토석류 파의 선단부보다 빠르게 이동하는 것을 수치적으로 잘 모의하고 있다. 이 경우도, 초기 토석류의 두께를 높게 설정한 Case B가 Case A보다 실험결과를 상대적으로 잘 재현하는 것으로 나타났다. 한편, 초기 토석류의 두께가 낮게 설정된 Case B의 경우 토석류 진행과 반대의 방향으로 전파되는 수면파가 수치해석 Case A 및 실험보다 빠르게 산사태가 발생한 사면쪽으로 진행하는 것으로 나타났다.

Fig. 5(e)에서 실험자료는 하나의 에너지 맥파(pulse of energy)로 구성된 비선형 파가 반대편 사면에 거꾸러지는 것을 보여준다. 그리고 경사면을 치고 올라가면서 파의 붕괴(breaking)가 시작되지만 가파른 경사 때문에 완전히 발달하지는 않음을 보여준다. 그리고 얇게 깔린 물이 계속해서 맞은편 사면을 따라 치고 올라가는 쳐오름현상을 보여준다(Fig. 5(f)). 두 가지 수치해석 결과도 이러한 수면충격파의 전파 및 쳐오름 형상을 조금 차이는 있지만 양호하게 실험결과를 예측하는 것으로 나타났다. Figs. 5(e) and 5(f)에서의 수치해석 결과를 비교해보면, Case A의 결과는 토석류 파의 선단부가 실험과는 달리 바닥으로부터 떨어지는 현상을 보여주고, 이로 인한 흐름 저항으로 바닥을 따라 전파되는 토석류의 선단부가 맞은편 사면에 Case B보다 상대적으로 늦게 도달하는 것으로 나타났다. 이것은 실험에서 점토 성분이 없는 입자들로만 구성된 토석류의 거동을 Herschel-Bulkley 모델이 적절히 고려하지 못하는 것도 원인 중의 하나일 것으로 판단된다. 따라서, Fritz et al. (2001)의 실험자료를 더 정확하게 모의하기 위해서는 토석류의 구성물질 간 점성력의 영향을 고려하는 Herschel-Bulkley나 Bingham 모형보다는 순수 입자흐름(granular flow)를 모의할 수 있는 유변학적 모형의 적용이 필요한 것으로 판단된다.

Fig. 5.

Comparison of [lower] computed propagations of debris flow induced impact water wave and underwater debris flow with [upper] experimental snapshots of Fritz et al. (2001)

Fritz et al. (2003)은 총 137개의 수리실험을 수행하여, 이들 중 선정된 3개의 연속자료 모음에 대해서 PIV를 이용해서 계산한 일련의 순간 유속 벡터 자료를 제시하였다. 이 3개의 자료는 1) 상대적으로 낮은 충격 속도로 인해 토석류 상부에서 흐름분리가 발생하지 않는 경우, 2) 높은 충격 속도에 의해서 흐름분리가 발생하는 경우, 그리고 3) 반대 방향 즉 산사태 발생 사면 방향으로 진행하는 2차 파에 의해서 형성되는 단파(bore)의 경우를 포함한다. 이들 자료 중 Fritz et al. (2001)의 실험조건과 유사한 두 번째 PIV 실험 결과와 이 연구에서 계산된 수치해석 자료와 비교하여 분석하였다.

수로 중앙 단면에서 PIV 실험을 통해 계산된 유속벡터장과 두 가지 수치해석을 통해 수로 중앙단면에서 모의된 유속벡터장을 비교하면 Fig. 6과 같다. Fig. 6(a)에서 실험자료는 토석류가 수면에 충격을 가하면서 물 흐름이 분리되기 시작하는 순간을 보여준다. 물은 초기에 토석류의 침입으로 수면 부근에서는 우상향으로 분출되고 그 아래에서는 토석류 침입 방향으로 밀리게 된다. 두 수치해석 모두 이 순간에는 이러한 물의 분출과 밀림을 모두 실험자료와 유사하게 모의하는 것으로 나타났다. Fig. 6(b)에서 실험자료는 토석류의 계속되는 강한 충격력에 의해서 수심의 약 1.5배만큼 더욱 솟구쳐 오른 물의 표면이 약한 S 형상의 배후면을 형성함을 보인다. 여기까지는 두 가지 수치해석 모두 수리실험에서 보인 물이 솟구쳐 오르는 형상과 수면충력파의 유속성분을 매우 양호하게 재현하는 것으로 나타났다.

Fig. 6(c)에서 실험자료는 솟구쳐 오른 수면충격파의 배후면이 거의 90°로 연직선의 형상을 보이고 이 연직선 상의 물은 약 45°의 우상향 유속벡터 성분으로 전파되는 것을 보여준다. 수치해석 Case A의 모의 결과는 수면충격파의 상부가 더 이상 솟구쳐 오르지 못하고 정점(peak)이 완만하게 감쇄되는 것으로 나타났다. Case B의 수치모의 결과는 수면충격파 정점부의 위치는 실험결과와 유사하게 높이 위치하며 수면충격파의 배후면 형상은 실험 결과보다는 다소 왜곡되지만, Case A보다 우수하게 실험결과를 재현하는 것으로 나타났다. 이러한 해석 결과는 수면충격파가 발생할 때의 초기 형상과 솟구쳐 오르는 정점부의 높이가 침입하는 토석류의 유속뿐만 아니라 토석류의 두께에 큰 영향을 받음을 보여준다.

Fig. 6(d)에서 실험자료의 유속벡터를 보면 솟구쳐 올랐던 수면충격파가 약 -45° 방향으로 낮아지고 수면충격파의 대부분은 우측으로 이동하며, 후면은 여전히 약 70° 정도의 직선 형태임을 알 수 있다. 이에 상응하는 Case A 수치모의 자료는 수면의 각도가 이미 약 30° 정도까지 낮아졌으며, 수면충격파 내의 물입자 일부는 아래 방향으로 그리고 미약하지만 꼬리부분 일부는 토석류 진행과 반대방향으로 흐르기 시작함을 보여준다. Case B 수치모의 결과는 수면의 각도는 완만하지만, 계산된 유속벡터장을 보면 Case A보다 상대적으로 양호하게 실험결과를 재현하는 것으로 나타났다.

Fig. 6(e)에서는 실험에서 관측된 수면충격파의 배후면 각도는 더욱 낮아지고 수면충격파의 정점 부근에서 유속벡터가 빠르게 하강하는 한편, 수면충격파의 꼬리부분이 토석류 진행과 반대 방향으로 흐르기 시작하는 것을 볼 수 있다. 이 순간에 수면충격파의 형상과 바닥을 따라 전파되는 토석류의 형태를 모의한 두 가지 수치해석 결과는 모두 실험결과와는 다소 차이가 있는 것으로 나타났다. 즉, Case A는 수면충격파의 꼬리부분이 역방향으로 전파되는 것을 과대 산정하고 있으며, Case B는 볼록한 형태의 수면을 모의하는 반면에 실험결과는 오목한 형태의 수면을 보여주고 있다. 이러한 결과는 두 수치해석 모두 계산 후반부로 갈수록 수면충격파의 형상을 정확하게 예측하기에는 한계가 있음을 의미한다.

Figs. 5 and 6에서 제시된 실험결과와 수치해석 결과를 종합해 보면, 산사태 발생으로 형성된 토석류가 수면으로 침입하면서 발생시키는 수면충격파의 초기 솟구쳐 오르는 높이와 수면 형태 그리고 유속분포는 토석류 형상에 대한 수치해석 초기 조건을 적절히 설정한다면 양호한 정확도로 수치모의가 가능함을 알 수 있다. 하지만, 연직상향으로 솟구친 수면충격파의 정점이 중력에 의해 하강하면서 감쇄되는 단계에서는 유입되는 토석류의 초기 두께 조건에 따라 모의 결과가 상이하게 나타났다. 따라서 수로 바닥을 따라 전파되는 토석류의 특성을 재현하기 위해서는 토석류 특성에 맞는 유변학적 모형의 적용이 필요해 보인다. 한편, 실험 및 수치해석의 주요 목적 중 하나인 수면충격파가 반대편 사면을 따라 전파되는 쳐오름 높이는 두 수치해석 모두 양호하게 예측하는 것으로 나타났다. 결국, 수면충격파의 발생 형상과 전파 양상뿐만 아니라 수로 바닥을 이동하는 토석류의 전파 양상까지 보다 정확하게 모의하기 위해서는 입자로만 구성된 토석류 흐름을 정확하게 모의할 수 있는 유변학적 모형 또는 토석 입자들을 직접 추적하는 기법을 병행하는 오일러-라그랑지 수치기법의 적용이 필요하다고 판단된다.

Fig. 6.

Comparison of simulated velocity vector fields (left: Case A and right: Case B) with experimental PIV data of Fritz et al. (2003)

5. 결 론

토석류의 비뉴튼 유체 거동 특성을 수정된 Herschel-Bulkley 모형으로 제어한 3차원 다상흐름 수치모형을 구축하여 산사태/토석류 유발 수면충격파의 전파를 수치모의하였으며, 해석 결과를 타 연구자의 수리실험 자료와 비교 분석하여 적용한 수치모형의 적용성을 평가한 결과를 요약하면 다음과 같다.

산사태 발생으로 형성된 토석류가 수면으로 침입하면서 발생시키는 수면충격파가 솟구쳐 오르는 높이와 수면 형태 그리고 수중 유속분포를 정확하게 모의하기 위해서는 사면에서 내려와 수면으로 침입하는 토석류 거동의 정확한 모의가 중요한 것으로 나타났다. 즉, 수면으로 유입되는 토석류 선단부의 두께와 유속이 적절히 모의 된다면, 수면충격파의 정점부가 솟구치는 높이와 수면형은 매우 우수한 정확도로 예측이 가능한 것으로 나타났다.

토석류의 초기 형상을 다르게 설정한 두 가지 수치해석의 결과는 연직상향으로 솟구친 수면충격파가 최고점에 도달한 후 중력에 의해 하강하면서 감쇄되는 단계에서부터 상이해지는 것으로 나타났다. 솟구쳐 오른 수면충격파의 높이가 실험값보다 낮은 경우 수면충격파의 전파가 상대적으로 빠르게 진행되는 한편, 토석류가 전파되는 방향과 반대방향으로 단파 형태의 2차 파가 조기에 전파되는 것으로 나타났다. 토석류 초기 두께를 상대적으로 크게 설정한 수치모의 결과는 만을 가로지르는 수면형과 함께 반대편 사면에서의 쳐오름 현상까지 양호하게 실험자료를 재현할 수 있는 것으로 나타났다. 반대편 사면에 도달한 수면충격파가 사면을 거슬러 흐르는 최고 쳐오름 높이는 토석류 총량이 같은 경우 수면으로 유입되는 토석류의 초기 두께에 민감하지 않은 것으로 나타났다.

발생된 수면충격파가 만을 가로질러 전파되는 속도와 수면형은 초기에 솟구쳐 오른 높이와도 관련이 있지만, 바닥을 따라 전파되는 토석류 선단부의 형태와 유속에도 영향을 받는 것으로 나타났다. 따라서, 수로 바닥을 따라 전파되는 토석류의 전파 특성을 보다 정확하게 재현하기 위해서는 실험에서 점토 성분이 없는 입자만을 이용하여 재현한 토석류의 물질 특성에 맞는 유변학적 모형의 적용 또는 토석 입자들을 직접 추적하는 기법을 병행하는 오일러-라그랑지 수치기법의 적용도 고려할 수 있다고 판단된다.