1. 서 론
2. 토석류 유동에 관한 지배 방정식의 수식화
2.1 지배방정식의 유도
2.2 지배 방정식의 유한요소 수식화
3.3 가중행렬에 의한 불연속 시험함수()
3. 토석류의 역학적 구성모델
4. 단일경사 사면에 대한 토석류 흐름 해석
4.1 사면경사와 토석류 물성의 영향
4.2 에너지 소산비
4.3 사면 하부 제방의 영향
5. 결 론
1. 서 론
최근 기후변화에 따른 집중강우와 도시 확장에 따른 산지 개발 등으로 인하여 토사재해에 의한 인명과 재산피해가 증가하고 있다. 그 규모도 대형화됨에 따라 토사재해에 대한 잠재적인 위험성은 지속적으로 증가하는 추세이다(Quan Luna et al., 2011). 토사재해의 여러 유형 중 토석류는 많은 토사량, 빠른 속도, 큰 충격하중, 그리고 시-공간적 낮은 예측률 등으로 매우 위협적인 사면파괴로 인식되고 있다.
이동속도가 빠르고 이동거리가 길어 피해 범위가 넓은 토석류의 피해를 최소화하기 위해서는 피해 지역을 정확히 예측하여 이용을 제한하여야 한다. 이를 위해서는 토석류의 유동 특성을 이해하고 이동범위를 평가하는 것이 매우 중요하다. 하지만, 지형학적, 지질학적, 토질공학적 요소 등이 모두 결합되어 나타나는 토석류의 이동 메커니즘을 유동학적 모델만으로 해석하기에는 많은 어려움이 있다 (Iverson, 2003).
토석류 흐름에 대한 수치해석 프로그램은 지배방정식의 유도, 지배방정식의 수식화 그리고 연관된 구성모델로 개발로 구분할 수 있다. 예로, 역학적 응력-변형 해석을 프로그램 개발을 위해서는 지배방정식에 해당하는 Newton’s 2nd law에 의한 힘평형 방정식(
), Galerkin 수식화, 그리고 역학적 구성모델로 Modified- Cam-Clay 모델과 같은 다양한 토질역학 모델이 요구된다.
토석류 흐름에 대한 지배방정식은 토석류를 단일상으로 가정하거나(O’Brien et al., 1993; Iverson and Denlinger, 2001), 토립자와 유체의 2상으로 고려한 지배 방정식의 유도에 대한 연구가 활발히 진행되고 있다(Teisson et al., 1993; Pitman and Le, 2005). 하지만, 토립자와 유체의 상호작용을 고려할 수 있는 통합적인 이론식은 여전히 추가적인 연구가 필요한 실정이다(Pudasaini, 2012).
토석류 흐름 해석을 위한 지배방정식의 수식화는 연속체 해석과 불연속의 개별요소 해석으로 구분할 수 있다. 기존의 수치해석 기법들은 CFD(Computational fluid dynamics, Vassilevski et al., 2012), 유체 유동해석 프로그램(FLO-2D, Kanako 2D 등, Wu et al., 2013), DEM (Discrete element method, Salciarini et al., 2010) 을 이용하고 있다. CFD 기법은 일반적으로 유체에 관한 동역학적 해석법으로 산악지형과 같은 복잡한 해석 영역의 모델링에 매우 어려움이 있고, 증분 형식의 수식화에 의하여 해의 안정성이 저하될 수 있다. 일반적인 토석류 흐름해석을 수행하는 산악지형의 복잡한 운동학적 경계조건을 고려하는 것이 매우 어렵다. 유체 유동해석 프로그램은 대부분 유한차분법/유한체적법을 사용하며, 증분형태의 수식화에 의한 수치적 불안정성을 유발한다. 그리고 토석류를 고점성의 유체의 흐름으로 해석하기 때문에 함수비가 큰 입자흐름(granular flow)에 대한 지반공학적 특성의 고려가 결여되어 있다. DEM 기법은 입자의 거동 해석을 수행하는 미시적 관점의 수치 기법으로 해석영역이 극히 제한적이다. 그리고 단순한 입자간 역학적 모델을 이용하여 입자와 물의 혼합물에 대한 해석의 어려움이 있고, 유체와의 정교한 상호작용 연동해석이 어렵다.
국내의 경우는 상용화된 프로그램을 도입하여 토석류의 흐름특성에 대한 비교 연구가 진행중이며(Lee and Kim, 2013), 수치해석도구 개발은 아직 기초단계에 있다.
본 연구에서는 토석류의 유동 형상을 모사하기 위한 지반공학적 거동 메카니즘(토압, 역학적 구성모델 등)을 고려한 지배방정식을 유도하고, SU/PG(Streamline Upwind/ Petrov-Galerkin) 기법에 의한 수식화와 Newton의 반복과정을 적용하여 해의 수렴성을 향상시키고자 한다. 그리고 토사재해시 발생할 수 있는 다양한 경계조건을 고려할 수 있도록 유한요소 수치해석 기법을 제안하고자 한다.
2. 토석류 유동에 관한 지배 방정식의 수식화
2.1 지배방정식의 유도
토석류의 유동에 대한 지배방정식은 단일상 재료의 연속방정식과 모멘트 보존의 법칙에 의한 힘평형 방방정식으로부터 다음과 정리할 수 있다.
(1)
여기서, 
는 rectangular Cartesian 좌표계에서 각 방향으로 토사의 유속이며, 
는 임의 위치에서 토사의 유입 혹은 유출량이다. 그리고 
는 중력 가속도이며, 
는 토체에 작용하는 단위체적당 외부 하중이다.
제시된 연속방정식과 힘평형 방정식은 토사 유동의 3차원적인 거동을 표현할 수 있으나, 복잡한 산악지형을 해석하기 위해서는 경계조건과 모델링 소요시간에 대한 최적화가 필요하다. 따라서, 지배방정식에 대하여 깊이 방향(
방향)으로 Leibnize integration rule를 이용한 적분을 수행하였다(Vreugdenhil, 1994; Medina et al., 2008).
|
Fig. 1. Main variables and coordinate system to derive governing equations (Delinger & Iverson, 2004) |
토사의 밀도가 일정하다고 가정하고 깊이에 대한 평균 토사속도(
)를 변수로 하여 단위면적에 대한 연속방정식을 정리하면 다음과 같다(Shin, 2014).
(2)
그리고 
방향과 
방향에 대한 힘평형 방정식을 단위 면적에 대하여 정리하면 다음과 같다.
(3)
여기서, 
와 
는 변위율장(strain field)에 의하여 결정되는 토석류 내부의 수평토압계수(lateral earth pressure)를 나타낸다(Delinger and Iverson, 2004). 깊이에 따른 연직응력의 변화(
)에 비하여 
과 
가 상대적으로 작은 값을 나타내므로 무시하였다. 그리고, 연직방향의 유동속도는 
으로 결정되
므로 지배방정식의 주변수로 사용하지 않는다.
2.2 지배 방정식의 유한요소 수식화
토석류 유동에 대한 연속방정식 및 힘평형 방정식들을 정리하면 다음과 같다.
(4)
유한요소 수식화에서 주변수(main variable)의 선정은 수치해석의 안정성에 지대한 영향을 미친다. 동일한 지배방정식이라 할지라도, 적절치 못한 주변수의 선정은 해의 수렴성을 크게 저하시킨다. 특히, 유체 흐름해석에서 사용하는 주변수(
)를 사용하여 토석류의 두께가 없는 영역(
)을 해석하게 되면 수치해결 결과가 수렴하지 못하고 발산하게 된다.
따라서, 본 연구에서는 주변수를 
으로 선정하고 지배방정식을 행렬의 형태로 정리하였다.
(5)
여기서, 
, 
, 
, 
, 
, 
, 
, 
이다.
그리고 지배방정식을 다음과 같은 행렬형태로 표현할수도 있다.
(6)
여기서, 
, 
이다.
제시된 식 (5)와 식 (6)의 두가지 형태의 지배방정식에 대하여 보간함수를 이용한 공간분할(spatial discretization)을 위하여 Petrov-Galerkin 수식화를 수행하였다(Atallah and Hazzab, 2013).
|
(a) |
|
(b) |
Fig. 2. Numerical instability with time in dumping simulation. (a) Galerkin formulation, (b) SU/PG formation (DG scheme) |
(7)
여기서, 
는 일반적인 유한요소법의 형상함수(shape function)에 대한 행렬형태(matrix form)이고, 
(SU-PG 상향가중함수)는 유동방향을 기초로 한 상향 가중인자를 고려한 불연속 시험함수(test function)에 대한 행렬형태이다.
식 (7)의 수식화에 임의의 시간 
과 
사이에서 사다리꼴 법칙으로 시간적분을 수행하고, Newton의 반복과정을 적용할수 있도록 근사화하였다. 그리고 미지의 변수 
에 관하여 정리하여 유한요소 프로그램 코딩을 수행하였다.
3.3 가중행렬에 의한 불연속 시험함수(
)
일반적으로 유한요소법에서 사용되는 Galerkin 수식화는 토석류의 유동과 같이 대류흐름이 지배적인 문제에 있어서는 유사진동(pseudo oscillation)을 보인다(Fig. 2, Brooks and Hughes, 1982). SU-PG(Stramline-Upwind Petrov-Galerkin) 기법을 이용하여 이러한 유사진동과 이로 인한 수치해석의 발산을 억제하였다(Shin, 2014).
SU/PG 기법에서 불연속 시험함수 
는 유동의 흐름방향을 고려한 상향가중행렬을 사용하였다(Brooks and Hughes, 1982).
(8)
여기서, 무차원 계수 
, 
, 
(Berger and Stockstill, 1995)이다. 
와 
는 
와 
흐름방향으로의 가중행렬이며 일반적으로 DG(dissipative Galerkin)과 CDG(characteristic dissipative Galerkin) 기법이 주로 사용되고 있다(Hicks and Steffler, 1992).
방향으로의 가중행렬 
는 다음과 같은 순서로 계산할 수 있다. 우선 
행렬의 eigen-value 
와 eigen-vector 
를 각각 산정한다.
(9)
여기서, 
이다.
새롭게 정의되는 eigen-vector 
를 산정하는 방법은 DG와 CDG 기법이 있다. 토석류 흐름에 대한 SU/PG 수식화에서는 DG 기법이 안정적인 결과를 도출하였으며, 프로그램에 적용한 
는 다음과 같다.
(10)
마지막으로 
방향으로의 가중행렬는 
으로 산정할 수 있다.
방향으로의 가중행렬 
는 
행렬을 이용하여 동일한 방법으로 산정된다. 산정된 가중행렬 
와 
을 이용하여 SU-PG 수식화에서 불연속 시험함수 
(식 (8))를 결정한다.
3. 토석류의 역학적 구성모델
다양한 토석류의 유동에 대한 역해석 결과는 다량의 조립토(coarse particles)를 포함한 토석류의 점성과 항복응력이 토석류내 조립토를 제외한 부분에 대한 점성이나 항복응력보다 수십내지 수백배 큰 것으로 나타났다. 즉, 점도계에 의한 측정결과는 세립토 함유량에 따른 점도와 강도가 증가하는 경향은 확인할 수 있으나, 획득한 자료를 이용하여 실제 토석류의 속도나 퇴적두께를 예측하기는 어렵다(Takahashi, 2007).
Takahashi(2007)은 유동내 평균입자농도에 따라 토립자와 유체의 혼합물의 거동을 구분하였다. 입자농도가 0.2이하인 경우는 입자들간의 충동응력(collision stress)이 흐름의 하부에서 발생하는 미성숙 토석류(immature debris flow)로 구분하였다. 그리고 입자농도가 0.2보다 커지면 동적토석류(dynamic debris flow)로 구분하며, 입자의 충동응력, 난류 혼합응력(turbulent mixing stress), 그리고 점성응력(viscous stress)이 주요한 응력 성분이 된다. 이러한 토석류의 흐름내에서 토립자간의 직접적인 전단마찰은 미비할 수 있으나, 입자간의 충동은 흐름내부의 조립토에 의한 팽창(dilation)을 유발하게 된다.
입자간의 마찰 및 팽창에 의한 전단저항을 고려하지 않고, 유체 유동만을 이용하여 발생할 수 있는 최대 전단 저항은 산정해 보았다. 예로, 사면경사 30°, 토석류의 평균입자직경 d50=4.75mm(#4), 토석류의 두께 1m일 때, 평균흐름속도는 대략 23m/s에 이르게 된다(Julien and Paris, 2010). O’Brien and Julien(1988)의 실험결과에서 유체와 토립자의 혼합체에 대한 최대 수리학적 물성값은 항복응력=1000dyes/cm2(0.01tonf/m2), 점성=100poise (0.001tonf･s/m2)이다. 유체의 유동저항에 의한 저항강도는 0.03tonf/m2로, 정지상태에서 연직으로 자립할수 있는 높이는 5.7cm에 불과하다. 따라서, 토석류 흐름에 대한 현실적인 예측을 위해서는 유동학적 저항 특성 뿐만 아니라, 토질역학에서의 마찰과 팽창에 의한 전단 저항력을 동시에 고려하여야 한다.
본 연구에서 토체와 바닥면의 점탄소성 역학적 구성관계를 
으로 정의하였다. 그리고 토석류의 유동특성에 대한 역학 모델은 유체와 토립자의 특성을 동시에 고려할 수 있도록 다음과 같은 Coulomb-viscous 모델을 사용하였다.
(11)
여기서, 
[tonf/m2]와 
[tonf･s/m2]는 유동계(rheometer) 실험을 통하여 결정되는 토석류의 항복응력과 동점성이며, 
[-]는 토석류의 난류-분산 효과를 고려하기 위한 유사 manning 마찰계수이다(Julien and Lan, 1991).
제안된 역학적 구성모델은 일관접선계수(consistent- tangent modulus)를 이용한 수식화를 수행하였다(Simo and Taylor, 1985). 항복함수와 소성 포텐셜 함수는 다음과 같이 정의하였다.
(12)
여기서, 
과 
는 토석류와 저면의 마찰각(friction anlge)과 팽창각(dilation angle)이다.
일관성조건(consistency condition)에 의하여 소성변형 증분계수 
는 다음과 같이 산정된다.
(13)
여기서, 
는 토석류의 바닥면 전단변형에 대한 탄성 변형계수[tonf/m3]이며, 전단응력의 크기는 
이다.
토석류의 소성변형에 대한 일관접선계수는 다음과 같이 산정할 수 있다.
(14)
여기서, 전단응력 단위 벡터(deviatoric stress unit vector)는 
, 
이다.
4. 단일경사 사면에 대한 토석류 흐름 해석
개발된 유한요소 프로그램의 수치해석상의 안정성을 확인하고, 지형조건과 입력 물성치의 영향에 의한 토석류의 흐름특성을 분석하기 위하여 예제해석을 수행하였다. 경사지와 평지로 이어진 해석 단면(1m×1m, 2,500요소, 7,701절점)의 정상부 부근에 토석류 유발(면적=0.2m×0.3m) 지점이 위치한다. 수치해석에 사용된 사면의 기하학적 조건 및 적용 물성치는 Table 1에 정리하였다.
=0.2 (SU-PG scheme)
=25°, 35°
=2 tonf/m3
=0 tonf/m2, 
=0 tonf･s/m2, n=0 
=5°, 10°, 
=
, 
=10 tonf/m34.1 사면경사와 토석류 물성의 영향
수치해석은 바닥면과 토석류 사이의 마찰각(기본값: 
=5°), 사면경사(기본값 
=35°), 그리고 토석류 유발높이(기본값: h=0.15m)를 변화시키면서 비교 분석을 수행하였다. Fig. 3은 토석류 흐름해석에서 Newton의 반복과정을 통한 해의 수렴성을 나타내고 있다. 본 연구에서는 주변수를 
을 사용하였으며, 토석류의 흐름이 없는 구간에서도 안정적인 수렴성을 보이고 있다. 불평형 하중에 의한 정규화된 error norm보다는 불평형 에너지에 대한 정규화된 error norm이 빠른 수렴성을 보이는 것을 알 수 있다.
Fig. 4은 토석류 유발높이 h=0.15m인 경우에 대하여 사면경사, 저면의 마찰각, 그리고 하부의 제방 설치에 따른 토석류 형상의 변화를 나타내고 있다. Case A(
= 35°, 
=5°, h=0.15m)인 경우, 사면 상부에서 유발된 토석류는 0.63sec 이후에 하부 평지에 도달하고, 평지에서 확산되는 것을 확인할 수 있다. Case B(
=25°, 
=10°, h=0.15m)의 경우, 상대적으로 완만한 사면경사와 높은 마찰각으로 사면내에서 더욱 확산되며 더욱 많은 토석류가 사면내에 지체됨을 알 수 있다(time=4.0sec @ Fig. 4b). 1.55sec 경과후에 하부 평지에 도달한 후, 하부 평지에서 멀리 확산되지 못하고 경사지 단부에 집중되어 있음을 확인할 수 있다.
Fig. 5는 사면의 기하학적 조건(
), 유발 토사량(h), 그리고 저면 마찰저항(
)이 토석류 흐름특성에 미치는 영향을 분석하였다. Fig. 5a-토체높이와 5b-토체높이는 Fig. 3의 “d” 지점(비탈기슭, foot of slope)과 “e”지점(저부 평면의 끝)에서 시간경과에 따른 토석류의 높이 변화를 나타내고 있다. 기본 case(h=0.15m, 
=35°, 
=5°)에 대비하여 흐름특성이 완화될수록(
↓, 
↑, h↓), 사면 저부에 도달하는 시간이 많이 소요되었다. 그리고 저부 평면에서의 평균 확산속도를 산정해 보면 기본 case(h=0.15m, 
=35°, 
=5°)의 경우 0.159m/s인 반면, 나머지 경우들은 0.071m/s(h=0.10m, 
=35°, 
=5°), 0.085m/s(h=0.15m, 
=25°, 
=5°), 0.085m/s(h=0.15m, 
=35°, 
=10°)으로 나타났다. 특히(h=0.15m, 
= 35°, 
=10°)인 경우, 사면내에서의 흐름속도는 기본 case에 매우 근접하게 나타났으나, 평지의 유동속도는 기본 case의 50% 정도에 미치는 것을 알 수 있다. 사면경사가 저면 마찰각보다 크기 때문에, 사면 내부에서의 흐름은 토석류 유발량과 사면경사에 의하여 질량 흐름(Fig. 5c-체적유량)이 결정되고, 하부 평지에서는 유입 질량흐름과 저면 마찰각에 의해서 확산속도가 결정된다.
|
Fig. 3. Numerical convergence in normalized force and energy norm during Newtonian iteration (Case A) |
=35°, 
=5° (default)
=25°, 
=10°
=35°, 
=5° (embank)
4.2 에너지 소산비
|
|
|
|
Fig. 5. Debris flow configuration with time under different geometry and material property (initial debris height h=0.15 m) | |
토석류의 이동과 확산은 위치에너지의 마찰저항에 의한 에너지 감소, 그리고 운동에너지에 의한 이동과 확산으로 이해할 수 있다. 이때 마찰에너지는 토석류와 저면 사이의 마찰 뿐만 아니라, 내부의 입자간의 이동에 의한 에너지 변화도 포함하게 된다. 단일 물체의 사면내의 이동에 의한 에너지 변화는 에너지 소산비 
=(energy dissipation)/(potential energy)=
로 정의할 수 있다(Kokusho and Hiraga, 2012). 
는 토석류의 흐름 특성에 관한 지수로 사용할 수 있는데, 
가 커질수록 토석류가 더욱 먼 거리까지 확산하게 된다.
Fig. 5d는 시간경과에 따른 에너지 소산비의 변화를 정리하였다. 토석류 흐름이 시작되는 초기에는 모든 case의 
값이 거의 동일하며, 질량흐름이 진행함에 따라 
값이 증가함을 알 수 있다. 최종적인 
값은 평지에 퇴적된 토석류의 위치수두에 의하여 결정된다. 기본 case의 
값이 크게 나타난 것은 최종단계에서의 위치수두가 낮은 상태로 확산이 상대적으로 많이 진행된 것으로 판단할 수 있다.
하지만, 본 예제와 같이 마찰저항이 지배적인 입자 흐름에서는 
를 이용한 토석류의 확산 평가는 제한적인 것으로 판단된다. 따라서, 사면내에서의 짧은 진행거리와 좁은 평지 조건 등의 기하학적 경계조건을 개선하여 
의 변화에 대한 추가적인 연구가 필요한 것으로 판단된다.
|
|
(a) | (b) |
|
|
(c) | (d) |
Fig. 6. Effect of embankement at the bottom of the slope. (a) debris height, (b) debris velocity, (c) debris mass flux, (d) reaction force on the embankment | |
4.3 사면 하부 제방의 영향
사면의 하부에 위치한 제방의 토석류 흐름(Case C, Fig. 4)에 대한 영향을 파악하기 위하여 수치해석 결과를 정리하였다. Fig. 6은 제방의 설치에 따른 토석류의 흐름특성의 변화와 제방에 가해지는 하중을 보여주고 있다. 사면에서 진행해오던 토석류는 제방에 막히면서 일시적으로 토석류의 높이가 증가하지만, 토석류의 흐름이 측면으로 진행되면서(time=4.0sec, Case C, Fig. 4) 토석류의 본류가 진행된 후에는 제방이 설치되지 않는 경우보다 낮아지는 것을 확인할 수 있다(점선, Fig. 6a-토체높이). 제방 설치에 의하여 평지의 종점에 토석류가 도달하는 시간은 0.5sec이상 지체되고 토석류의 높이 또한 낮아진 것을 확인할 수 있다. 토석류의 흐름속도와 질량흐름도 제방 설치로 인한 흐름의 측방 확산과 정체로 인하여 낮아짐을 알 수 있다(Fig. 6b-토석류속도 and 6c-체적유량). 하부에 설치된 제방에 가해지는 토압은 제방이 설치되지 않은 경우의 흐름토압과 비교에서 최대 하중이 지속되는 시간이 다소 증가하지만, 제방의 측면으로 토석류가 확산하면서 더 이상 증가하지 않음을 알 수 있다(Fig. 6d-토체압).
5. 결 론
최근 토석류에 의한 피해가 증가하고 있으며, 이를 최소화하기 위하여 토석류의 이동 메커니즘을 이해하고 예측하기 위한 활발한 연구가 진행중이다. 토석류 흐름현상에 관한 수치해석은 주로 유체유동 해석프로그램을 이용하고 있으나, 토석류를 유체이동으로 이해하는 기존의 연구들은 수치적인 안정성과 모델링 그리고 다양한 경계조건의 적용에 대한 제약이 있다.
본 연구에서는 지형 형상을 고려할 수 있는 토석류의 유동현상에 대한 연속방정식과 힘평형 방정식을 유도하였다. 지배방정식에 대한 깊이적분을 수행하고, 토체의 두께 
, 
와 
흐름방향의 평균속도 
, 
를 주변수로 사용하여 흐름이 없는 해석영역에 대한 수치적인 안정성을 확보하였다. 지배방정식의 보간함수를 이용한 공간분할에서 Petrov-Galerkin 수식화 기법을 적용하였다. 이때 불연속 시험함수는 유동방향을 기초로 한 상향 가중인자를 고려하였으며, 가중행렬은 DG기법이 토석류 흐름에 더욱 적합한 것으로 나타났다. 시간적분과 Newton의 반복과정을 통하여 해석결과의 수렴성과 안정성을 개선하였다. 그리고 토석류의 유체 및 토립자의 특성을 동시에 고려할 수 있는 Coulomb-viscous 역학적 구성모델을 제시하였다.
단일경사 사면에서 사면경사, 토사 유발량, 저면 마찰 저항이 토석류 흐름특성에 미치는 영향을 비교 분석하였다. 사면내에서의 질량흐름은 유발량과 사면경사에 의하여 의존하며, 평지의 확산흐름은 유입 질량흐름과 저면 마찰각에 의해서 결정되었다. 마찰저항이 지배적인 흐름에서는 에너지 소산비를 이용한 평가는 토석류의 확산 평가에 제한적이다. 그리고 사면하부에 설치된 제방은 토석류 흐름을 측면으로 확산시키면서 평지에서의 질량흐름을 완화시키는 것으로 나타났다.
개발된 해석프로그램을 활용하여 토석류 발생 예상 지역의 다양한 위험인자에 대한 평가를 수행하고, 토석류 흐름을 조절하여 피해를 최소화하기 위한 시설물의 설계방안에 대하여 제안할 수 있을 것으로 판단된다.
