1. 서  론

우리나라의 하천은 여름 홍수기에 대부분의 강우가 집중된다. 이러한 특성 때문에 홍수기에 댐 상류로부터 다량의 부유사가 저수지로 유입되고, 두꺼운 밀도류의 형태를 유지하며 댐까지 전파되는 문제가 발생한다. 우리나라의 대표적인 대형 다목적 댐인 소양강댐이나 대청댐은 매년 홍수기마다 고농도 부유사 밀도류가 발전방류를 통해 하류에 배출되었다(Chung et al., 2009; Ryu et al., 2011). 특히, 임하댐에서는 2003년 태풍 매미에 의해 880 NTU (Nepthelometric Turbidity Unit)를 넘는 농도의 부유사 밀도류가 장기간 발생한 기록이 있다(An and Julien, 2014). 하천에서 고농도의 부유사 밀도류는 부유사의 탁도 및 수반된 오염물질을 통해 하천 환경에 악영향을 미친다. 이런 문제를 해결하기 위해 한국수자원공사에서는 부유사 저감에 대한 계획을 세웠으며, 주요 저감 방법으로 댐에 선택적 취수 시설(selective withdrawal structure)을 설치하는 것을 추진하였다. 선택적 취수 시설의 설계를 위해서는 취수부에서 부유사의 농도 분포에 대한 정보가 필요한데, 이를 위해서 밀도류 해석은 필수적이다(Chung et al., 2009). 위에서 설명한 대규모 댐에서는 하층밀도류(underflow) 형태로 발생한 부유사가 도중에 중층밀도류(interflow) 형태로 전파되기도 하지만, 부유사의 농도가 크거나 댐의 규모가 작으면 하층밀도류 형태를 유지하여 전파되므로 이에 대한 연구는 여전히 중요하다.

많은 해외 연구자들이 실험을 통해 하층밀도류에 대한 연구를 수행하였다. Parker et al. (1987)은 부유사 밀도류 실험을 통해 물 연행(entrainment)과 Richardson 수의 관계를 제시하였고, 부유사 연행과 입자의 레이놀즈 수의 관계를 검토하였으며, 부유사 분포의 형상계수(shape factor)와 하상전단계수(bed shear coefficient)의 값을 제시하였다. Altinakar et al. (1990)은 하상경사가 작은 경우 입자크기와 부력 흐름률(buoyancy flux)이 전파속도 및 선단부의 두께에 미치는 영향을 분석하였다. Garcia (1993)는 염수 밀도류와 부유사 밀도류의 수직구조 및 도수가 유사한 양상을 보이는 것을 관찰하였고, 층적분(layer-averaged)된 밀도류의 운동량보존방정식에서 농도 및 속도의 형상계수가 1인 것을 제시하였다. Bonnecaze et al. (1995)는 2차원적으로 전파되는 밀도류의 퇴적량이 유속과 입자의 침강 특성에 대한 함수임을 제시하였다. Ohey and Schleiss (2007)와 Toniolo et al. (2007)은 수로에 장애물이 있을 경우 밀도류에 의한 부유사 퇴적 양상의 변화에 대해 분석하였다.

실험뿐만 아니라 수치모형을 통해서도 많은 연구자들이 하층밀도류에 대해 연구를 수행하였다. Choi and Garcia (1995)는 dissipative Galerkin FEM을 사용하여 밀도류의 전파를 수치모의하고 선단부의 전파, 내부 도수, 퇴적의 특성을 분석하였다. Choi and Garcia (2002)는 2차원 k-ε 난류모형을 사용하여 밀도류의 수직구조를 수치모의하고 부력에 의한 난류생성과 층적분 모형의 형상계수를 분석하였다. Kostic and Parker (2003)는 수치모의를 통해 밀도류에 의해 모래가 상부에 퇴적되고 점토가 하층에 퇴적되는 양상을 모의하였다. Cao et al. (2015)은 2차원 이중 층적분 모형을 이용하여 중국의 Xiaolangdi 저수지에 적용하고 모형의 적용성 및 배사 효율에 대해 검토하였다. Lai et al. (2015)은 2차원 층적분 모형을 이용하여 대만의 Shihmen 저수지에 적용하고 모형의 적용성을 검토하였다.

위와 같이 하층밀도류에 대하여 많은 실험 및 수치모의를 통한 연구가 진행되어 왔으나 밀도류의 전파과정에 대한 분석은 많지 않은 상황이다. 또한 이를 안정적으로 모의할 수 있는 수치모형을 개발할 필요가 있다. 따라서 본 연구의 목적은 부유사 밀도류의 수치모의를 위하여 유한차분법을 바탕으로 한 ULTIMATE 기법을 이용하여 모형을 구축하고 수치모의를 통하여 사면에서 발달하는 부유사 밀도류의 수리적 특성을 파악하는 것이다.

2. 지배방정식

부유사 밀도류를 모의하기 위한 부정류 1차원 층적분 모형은 다음과 같이 구성된다.

 (1)

 (2)

 (3)

Eqs. (1)~(2)는 각각 유체의 질량보존방정식, 운동량보존방정식, 그리고 Eq. (3)은 부유사의 질량보존방정식이다. 여기서 는 주 흐름방향 위치, 는 시간, 는 밀도류층의 깊이, 는 밀도류의 평균 속도, 는 부유사의 평균 농도, 은 입자의 수중비중(submerged specific gravity), 는 중력가속도, 는 하상고, 그리고 는 입자의 낙하속도(fall velocity)다. 와 는 각각 밀도류에 연행되는 물과 부유사에 대한 연행계수, 는 전단속도, 그리고 는 밀도류의 바닥 근처 농도에 대한 형상계수이다. 부유사 밀도류에 의해서 발생하는 퇴적이나 재연행을 고려하기 위해서는 다음과 같은 Exner 방정식을 사용해야 한다.

 (4)

여기서, 는 하상 입자의 공극률이다. Eqs. (1)~(4)는 다음과 같이 무차원화하여 표현할 수 있다.

 (5)

 (6)

 (7)

 (8)

여기서, 위첨자와 아래첨자는 각각 무차원화된 변수와 유입부의 변수를 의미한다. 모든 변수들은 유입부의 변수들을 이용해 무차원화하였다(, , , , , , ). 는 유입부의 Richardson 수()이다.

3. 폐합 문제

Eqs. (5)~(8)를 통해 , , , 를 계산하기 위해서, 그 전에 계산해야 할 미지수들이 남아있는 폐합 문제가 있다. 이를 해결하기 위해 기존에 다른 연구자들의 연구를 통해 제시된 식들을 사용하였다. 물과 유사의 연행계수인 와 는 각각 Fukushima et al. (1985)과 Garcia and Parker (1993)가 제시한 관계식을 이용하였으며 다음과 같다.

 (9)

 (10)

여기서, 는 이고 이다. 과 는 입자 레이놀즈 수 에 의해 결정되며 다음과 같다. 여기서 는 입자의 대표입경이다.

  (11)

무차원 전단속도 는 무차원 평균 속도를 이용하여 다음과 같이 계산된다.

 (12)

여기서, 는 하상전단계수로 본 연구에서는 Choi and Garcia (1995)가 사용한 0.01을 사용하였다. 형상계수 는 Garica and Parker (1989)가 제시한 범위인 1 - 2에서 입자크기에 따라 다르게 사용하였다.

4. 수치해석 방법

4.1 이동경계

본 모형에서는 밀도류의 전파를 모의하기 위해서 이동경계조건을 사용하였다. 이동경계조건을 적용하기 위해 Landau 변환(Crank, 1984)을 이용하였다. Landau 변환에서는 이동경계조건을 구현하기 위해서 위치와 시간에 대한 계산영역을 각각 새롭게(, ) 정의한다. 여기서 는 시간에 따라 전파되는 밀도류의 유입부에서 선단부까지의 거리()를 유입부의 밀도류의 두께로 무차원화한 값이다. 와 를 계산하기 위해 새로운 변수를 정의하고(, ), 새로운 계산 영역에서 Eqs. (5)~(8)을 적용하면 다음과 같이 변형된다.

 (13)

 (14)

 (15)

 (16)

여기서, 는 밀도류의 무차원 전파 속도이다. 위의 Eqs. (13)~ (16)을 사용하면 절점의 수를 유지하면서 밀도류의 전파 양상을 모의할 수 있다.

4.2 경계조건

본 연구에서 고려해야 할 경계조건으로는 유입부 경계면과 선단부 경계면이 있다. 밀도류의 선단부가 수로의 유출부까지 도달한 경우에는 수로의 유출부를 경계면으로 한다. 모든 변수를 유입부의 값으로 무차원화하였으므로 유입부에서 , , 그리고 는 모두 1이다. 밀도류가 전파중일 때 밀도류의 선단부에서 , , 그리고 의 값은 보외법(extrapolation)을 사용하여 계산하였으며, 는 이동된 선단부에서의 하상자료를 사용하였다. 또한 선단부의 는 로 설정하였다. 밀도류가 수로의 유출부에 도달하는 경우 의 값을 0으로 하고 유출부 경계조건을 적용하였다.

4.3 수치모의 기법

밀도류는 강한 이류(highly-advective flow) 특성을 보이므로 이를 모의할 수 있는 적절한 계산기법이 필요하다. 특히 선단부나 도수영역처럼 급격한 변화를 보이는 구간을 모의할 때는 오버슈트(overshoot)나 언더슈트(undershoot)가 흐름 전체에 미치는 영향을 고려해야 한다. 따라서 본 연구에서는 모형에 QUICKEST 기법과 ULTIMATE 제한자를 사용하여 Eqs. (13)~(15)의 이류항을 계산하고 오버슈트와 언더슈트를 제한하였다.

4.3.1 QUICKEST 기법

QUICKEST (Quadratic Upstream Interpolation for Convective Kinematics with Estimated Upstream Terms) 기법은 Leonard (1979)가 제시한 3차 정밀도의 계산기법이다. 부정류에서 사용되는 본 기법은 이류에 영향을 받는 변수 가 이송률 에 의해 이류될 때, 시간 증분량 동안 가 이류되는 정도를 고려한다. 이를 위해 와 가 동안 하류로 이동하면서 변동되지 않는다고 가정한다. 양의 이송률 로 이송되는 변수 에 대한 1차원 이류 방정식은 다음과 같다.

 (17)

Eq. (17)은 유한차분법을 사용하여 다음과 같이 차분된다.

 (18)

여기서, 는 번째 위치, 번째 시간에서의 를 의미하며 은 동안 값의 평균이다. 아래첨자과 은 각각 번째 절점에서 만큼 양의 방향과 음의 방향에 있는 값임을 의미하며 Fig. 1과 같다. 여기서 우변항을 적분하여 표현하면 다음과 같다.

 (19)

여기서, 로, 는

Courant 수이다. Eq. (19)에서 위의 가정을 추가하여 적절한 QUICK 관계식을 적용하면 다음과 같이 전개된다.

 (20)

이것을 Eq. (18)에 적용하면 에 대해 다음과 같은 식을 얻을 수 있다.

(21)

(22)

아래첨자, , , 은 각 절점의 위치로 Fig. 1과 같다.

4.3.2 ULTIMATE 제한자

ULTIMATE (Universal Limiter for Transient Interpolation Modeling of the Advective Transport Equation) 제한자는 Leonard (1991)에 의해 제시된 제한기법이다. 2차보다 높은 정밀도의 수치기법에서 효과적이며 특히 급격한 변동을 보이는 형상을 모의하는데 매우 우수한 성능을 발휘한다. 제한자를 일반화시키기 위해 NVD (Normalized Variation Diagram) 개념을 이용한다. 각각 정규화된 값들은 다음과 같다.

 (23)

 (24)

오버슈트와 언더슈트에 의한 진동을 억제시키기 위해 단조성 유지 평가를 이용하여 의 값을 다음의 범위로 제한한다.

  (25)

 (26)

Eqs. (25)와 (26)을 거쳐 계산된 값을 Eq. (18)에 대입하여 값을 계산한다.

5. 모형의 적용 및 결과

본 연구에서는 밀도류의 특성을 파악하기 위해 Garcia (1990)의 실내 실험자료를 모형의 적용대상으로 이용하였다. 실험자료에 모형을 적용하여 모의 결과를 실험 결과와 비교하고 분석하였다. Garcia (1990)는 실내 실험을 통해 해저협곡을 타고 흘러가 선상지에 전파되면서 발생하는 도수에 대해 연구하였으며, 이 때 사용된 수로의 재원은 Fig. 2와 같다. 도수의 특성을 분석하기 위해 염수 밀도류와 부유사 밀도류의 부력 흐름률()을 변화시키면서 관찰하였으며 자세한 사항은 Table 1과 같다. 여기서 SAL은 염수 밀도류 실험이고 NOVA, DAPER, GLASS A는 부유사 밀도류 실험이며 각각의 실험에서 사용된 부유사 입자의 종류를 의미한다. 이를 모의하기 위해 종방향 격자의 수는 유입부에서 선단부까지 100개를 사용하였으며 선단부가 하류로 전파되어 계산영역이 증가되면 격자의 수는 유지되고 격자의 크기는 이에 비례하여 증가하게 된다. 계산에 사용 시간 간격은 Courant 수가 1에 가까운 수가 되도록 설정하였다.

5.1 부유사 밀도류의 전파

부유사 밀도류의 전파에 대한 부정류적인 특성은 실험을 통해 파악하기 쉽지 않으므로 수치모의를 통해 자세한 관찰을 수행하는 것이 중요하다. Fig. 3은 모형을 NOVA 2실험에 적용하여 시간에 따른 밀도류의 두께, 속도, 농도, 그리고 부력 흐름률의 변화를 도시하고 측정이 수행된 시간의 모의결과와 측정결과를 비교한 것이다. Fig. 3을 보면 모형이 전반적으로 측정결과와 유사하게 모의하나 도수가 발생하는 구간에서 다르게 모의되는 것을 볼 수 있다. 이것은 도수의 복잡한 형태를 모형과 측정 모두 명확하게 표현하기 어려운 점에 의한 것으로 판단된다. Figs. 3(a)와 (b)를 보면 모형이 경사부에서 발달되어 선단부의 경사가 급한 형태를 보이는 밀도류의 전파 양상을 수치진동 없이 잘 모의하는 것을 알 수 있다. 경사부에서 Richardson 수가 평형상태에 도달하여 유속이 유지되고, 상층 경계면에서 물 연행으로 인해 밀도류의 두께가 일정하게 증가하는 것을 볼 수 있다. 이후 경사 변경점에서 밀도류의 두께가 크게 증가하기 시작하고 유속이 감소되는 것을 볼 수 있다. 이것은 경사에 의해 생성되는 밀도류의 추진력이 사라지기 때문이다. 밀도류가 전파되는 중에는 평탄부에서도 밀도 Froude 수()가 1 이상인 사류(supercritical flow)를 유지한다. 이 때 다시 Richardson 수가 평형에 도달해 유속이 유지되는 현상을 볼 수 있다. 밀도류가 수로 끝에 도달하면 자유낙하로 인해 유속이 국부적으로 증가하며 밀도류의 두께가 고정된다. 수로 끝에 도달한 밀도류는 하류단에서부터 밀도 Froude 수가 1보다 작아져 상류(subcritical flow) 상태가 되며, 하류의 영향을 받아 평탄부에서 밀도류의 두께가 점점 증가하게 된다. 이 영향이 경사 변경점까지 전파되어, 사류를 유지하는 경사부와 만나 도수가 발생하게 된다. Figs. 3(c)와 (d)를 보면 부력 흐름률이 거의 유지되는 상태로 전파되나 농도가 하류로 갈수록 작아지는 것을 볼 수 있다. 이것은 하류로 갈수록 밀도류의 두께가 증가하기 때문이다. 또한 시간이 충분히 지난 후에 부력 흐름률이 수로 전체에 걸쳐서 균일해지는 것을 볼 수 있는데, 이것은 하층경계면에서 부유사의 연행 및 퇴적이 적기 때문인 것으로 보인다.

5.2 전파 속도

밀도류의 전파 속도에 영향을 미치는 요인에 대해 분석하였다. 이를 위해 모형을 부유사 밀도류 및 염수 밀도류 실험에 적용하여 모의 결과를 분석하였다. Fig. 4는 입자의 크기가 다르거나 부력 흐름률이 다른 실험 자료들을 이용하여 밀도류 선단부의 이동 양상을 비교한 것이다. Figs. 4(a)와 (b)를 비교해 보면 입자의 크기는 밀도류의 전파 속도에 큰 영향을 주지 않는 반면, 부력 흐름률은 전파 속도에 뚜렷한 영향을 미치는 것을 알 수 있다. Fig. 4(a)에서 GLASS A2의 유속이 다른 것은 대부분의 입자가 퇴적되어 부력 흐름률이 크게 감소했기 때문으로 볼 수 있다. Britter and Linden (1980)은 부력 흐름률과 전파 속도에 대해 다음과 같은 관계식을 제시하였다.

         for   (27)

Eq. (27)과 염수 밀도류의 전파 속도 모의 결과를 비교하여 Fig. 5에 제시하였다. 비교 결과 모의 결과가 Eq. (27)의 전파 속도와 거의 일치하는 것을 볼 수 있다.

Fig. 6은 밀도류의 전파 속도가 시간에 따라 변화되는 양상을 도시한 것이다. 초기 전파 속도가 평형 Richardson 수에 도달할 때 까지 급격하게 감소된 후, 속도의 변화량이 크게 줄어드는 것을 볼 수 있다. 이후 하상의 경사 변경점을 만나 전파 속도가 다시 급격하게 감소된다. 이러한 양상은 밀도류의 유속 분포 형태와 거의 일치한다.

5.3 내부 도수

Garcia (1990)의 실험은 강한 이류 특성을 가지고 도수가 발생한다. 따라서 유한차분법이나 유한체적법을 통한 계산을 수행하기 위해서는 적절한 수치모의 기법이 필요하다. Figs. 7과 8은 NOVA 2 실험에 대해 각각 Courant 수가 0.1, 0.01일 때 ULTIMATE 제한자의 유무에 따른 수치모의 결과를 비교한 것이다. Fig. 7(a)를 보면 ULTIMATE 제한자를 사용하지 않을 경우 비정상적인 모의 결과가 발생하는 것을 볼 수 있다. 이것은 QUICKEST 기법만을 사용할 경우 경사 변경점에서 오버슈트와 언더슈트가 발생하여 과도한 진동이 나타나기 때문이다. 이러한 현상은 1이 넘지 않는 Courant 수에서 ULTIMATE 제한자를 사용하면 억제할 수 있다(Leonard, 1991). Courant 수를 줄여 오버슈트의 영향을 감소시키면 Fig. 8과 같이 ULTIMATE 제한자 사용 여부와 무관하게 유사한 결과를 모의하는 것을 볼 수 있다. 도수 후의 수심은 0.216 m로, 공액수심(conjugate depth) 공식을 이용하여 계산한 수심 0.20 m와 비슷하다. 약간의 차이는 공액수심 공식이 도수구간에서 물 연행 등의 요소를 고려하지 않기 때문에 발생하는 것으로 보인다. 그러나 Courant 수를 너무 작게 할 경우 과도한 계산량이 필요하다. 이러한 이유로 ULTIMATE 기법을 사용하는 것은 적절한 것으로 판단된다.

5.4 부유사 밀도류에 의한 퇴적

부유사 밀도류는 밀도류의 하층 경계면에서 부유사의 연행과 퇴적이 발생하여 하상을 변동시킨다. Fig. 9는 부력 흐름률과 Richardson 수가 유사한 두 실험에 대해 부유사 밀도류 입자의 퇴적 양상을 실험 결과와 모의 결과를 함께 비교한 것이다. 실험이 시작되고 30분 경과된 후의 결과로, 모형이 밀도류에 의한 하상 변동을 정량적으로 잘 모의하는 것을 볼 수 있다. 부유사 밀도류는 체적 농도에 비례해서 퇴적이 발생하기 때문에 체적 농도 분포와 유사한 형태의 퇴적 양상을 보인다. GLASS A2 실험의 결과를 보면 부유사가 유입부에서 대부분 퇴적되고 유출부 근처에서는 DAPER 4실험보다 적게 퇴적된다. 이것은 GLASS 입자는 유입부에서 대부분 퇴적되어 유출부 근처에서는 부력흐름률이 크게 감소되는 반면, DAPER 입자의 경우 부력 흐름률이 거의 유지하면서 하류로 전파되기 때문이다.

Fig. 10은 Fig. 9에서 경사 변경점 중심으로 확대한 것이다. DAPER 4 실험의 측정 결과와 모의 결과를 보면 평탄부에서 퇴적 두께가 증가하는 모습을 보인다. 이것은 경사 변경점을 지나 도수가 발생하면서 부유사의 연행이 크게 감소하였기 때문으로 보인다. 그러나 GLASS A2 실험에서는 그런 양상을 보이지 않는데, 이는 상대적으로 큰 입자의 크기로 인해 연행이 퇴적에 비해 매우 작기 때문인 것으로 판단할 수 있다. 이러한 현상은 Fig. 11과 같이 부력 흐름률이 크고 입자의 크기가 작은 경우 뚜렷하게 발견되어, 도수하는 영역 부근에서 하상이 불연속적인 형태를 보이는 것을 관찰할 수 있다.

6. 결  론

본 연구에서는 ULTIMATE 기법을 이용하여 밀도류 층적분 모형의 해석을 위한 수치모형을 제시하였다. 제시된 모형을 모형 협곡과 모형 선상지로 구성된 Garcia (1990)의 실내 실험에 적용하여 밀도류의 수중 전파, 전파 속도, 내부 도수 발생, 그리고 부유사의 퇴적에 대해 분석하였으며 그 결과는 다음과 같다.

수치모형이 협곡 및 선상지 모형에서 발달하는 부유사 밀도류를 잘 모의하였다. 밀도류의 선단부에서 수치진동이 발생하지 않았으며 선단부의 경사가 비교적 급하게 잘 모의되었다. 유입부에서 사류로 시작된 밀도류가 평탄부에 도달해서도 사류로 전파되었다. 밀도류가 하류단에 도달한 후 Froude 수가 1 이하로 감소되며 상류를 형성해 영향을 상류 방향으로 전파시켰으며 경사 변경점에서 도수를 발생시켰다. 부유사 입자 크기가 작은 경우 밀도류 전파에 의한 퇴적량이 적어 농도와 부력 흐름률이 시간에 따라 거의 일정하게 유지되는 것을 확인하였다.

부력 흐름률이 클수록 속도가 전파 속도가 증가하는 것을 확인하였다. 부유사 입자가 9 μm 보다 작은 경우 전파 속도에 직접적으로 영향을 미치지는 않으나, 입자 크기가 큰 경우 부유사가 많이 퇴적되어 부력 흐름률이 감소해 전파 속도가 감소하였다. 부력 흐름률과 전파 속도의 관계는 Britter and Linden (1980)이 제시한 식과 거의 일치하였다.

수치모형이 ULTIMATE 제한자를 사용하지 않을 경우 작은 Courant 수에서도 도수가 발생하면 오버슈트와 언더슈트에 의해 모의가 불가능한 것을 확인하였다. Courant 수를 0.01까지 감소시켜 계산할 경우 ULTIMATE 사용에 의한 결과의 차이가 거의 없으며, 도수 높이도 적절하게 계산하였다. ULTIMATE 제한자를 사용할 경우 과도한 수치진동을 제한시켜 안정적인 모의를 수행할 수 있는 것을 확인하였다.

수치모형이 부유사 밀도류의 전파에 의한 하상 변동을 정량적으로 적절히 모의하였다. 입자 크기가 클수록 경사가 급하게 퇴적되고, 부력 흐름률이 클수록 전체 퇴적량이 증가하는 것을 확인하였다. 평탄부에서 부유사 연행이 적기 때문에 입자 크기가 작고 부력 흐름률이 큰 경우 경사 변경점에서 하상이 불연속적인 형태를 보이는 것을 확인하였다. 본 모형은 향후 부유사 밀도류에 의해 발생하는 저수지 퇴적 문제와 국내 저수지에서 부유사 밀도류에 의해 발생하는 문제점들에 대한 분석을 가능하게 할 것으로 기대된다.