In this study, the numerical model developed by Hong et al.(2008) was improved to be applied to rapidly varying flows such as the inundation of dry land or flow transitions due to large gradients of the bathymetry. A numerical approximation was applied that was consistent with the conservation of momentum in flow expansions and with the Bernoulli equation in flow contractions. The approximation was second order, but the accuracy reduced to first order near extreme values by the use of a minmod limiter. The modified model was verified by acomparison with the theoretical critical depth of weir, and for sufficiently smooth conditions and a fine grid size, both approximations converged to the same solution. In terms of the grid size, it was more effective at obtaining solutions than the previous model and reproduced the inundation of dry land.
조석과 같은 유동수치모델에서 수평이류항의 공간이산화 방식은 해의 정확성, 수치기법의 단순화 및 효율성에 큰 영향을 미친다. 예를 들어, 중앙차분(Central differences)의 적용은 2차 정밀도를 가지나 해에서 불필요한 요동현상(Wiggles)을 나타낸다(Gresho and Lee, 1981). 요동현상은 특히 Closed boundary나 Thin dams와 같이 급격한 유동장의 변화를 가진 영역근처에서 나타난다. 반면에, 1차 Upwinding scheme은 절대적으로 요동현상이 없어 해의 안정성은 있지만 해를 얻는 과정에 2차 인위적 점성항(Artificial viscosity term)의 형태가 생겨 단절오차(Truncation error)를 발생 시킨다(Vreugdenhill, 1994). 그러므로, 이류항이 지배적인 유동장에서는 이런 인위적 점성항은 실제 물리적 점성을 압도하여 계산된 해가 실제 해보다 훨씬 더 매끄럽게(Smoothing) 나타난다. 다차 Upwinding scheme은 수치적 요동은 나타나지 않으나 1차 Upwinding scheme의 경우와 유사한 4차 인위적 점성항이 생겨 해의 Smoothing도 없이 바로 수치요동을 제거하여 실제 해를 얻는데 장애가 된다.
이런 인위적 점성항을 제거하면서 수치 요동현상을 최소화하기 위하여 여러 수치기법이 도입되어 왔다. Leendertse(1967)에 의해 수평이류항만을 별도로 이산화 시켜 ADI(alternating direction integration)기법을 수정 적용하였다. 이후, Stelling(1984)은 수평이류항의 동일한 방향의 곱의 성분은 중앙차분을 적용하고 서로 다른 방향의 곱의 성분은 위상차 오류(Phase error)를 점차적으로 줄여 요동을 소산시켰으며, 유속의 공간적인 요동은 4차 소산(Dissipation)방정식을 적용하였다(Stelling, 1984; Stelling and Leendertse, 1992). 이 밖에 Stelling and Leendertse(1992)은 수평방향 이류항의 동일한 방향 및 서로 다른 방향의 곱의 성분에 동시에 위상차 오류를 줄여 요동을 소산시키는 기법을 적용하였다. 최근에는 Hong et al.(2008)은 오탁방지막을 포함하는 3차원 유동모델 구축에서 수평방향의 이류항에 대한 공간이산화(Spatial discretization)를 위하여 ADI scheme을 개선한 방법(Stelling and Leendertse, 1991)을 적용하였는데 이는 수평이류항중에
상기의 기법 등을 통해 인위적 점성항을 제거하면서 유속의 공간적인 요동은 감소시켰으나, 육상지역 침수나 경사가 급변하는 수심(또는 장애물)영역에서의 흐름전이(Flow transition)와 같은 유동장이 급변하는 유동을 해석하는 경우, 격자나 시간간격을 상당히 줄이야 하므로 비효율적일 뿐만 아니라 유속경사가 상당히 클 때는 해의 요동현상도 나타나 정해(Correct solution)를 얻지 못할 수도 있다. 그러므로, 본 연구의 목적은 기존의 수치모형(Hong et al., 2008)에 수평방향 이류항의 공간이산화 방식 개선을 통해 상기의 비효율성을 개선하고 요동현상을 제어하여 급변하는 유동의 효율적인 해석을 수행하기 위함이다. 흐름이 팽창하는 경우는 운동량 보존방정식에, 축소하는 경우는 Bernoulli 방정식에 2차 중앙차분을 기본으로 공간이산화를 적용하였으며, 수치요동을 줄여 해의 정확성을 높이기 위해 이산화 차수를 1차로 줄이는 Minmod 제한자(Limiter)를 사용하였다. 수립된 개선 수치모형을 검증하기 위해 넓은 마루웨어(Broad crested weir)에 대하여 Bernoulli 방정식을 적용하여 유도한 웨어 한계수위와 본 연구의 개선수치모형을 적용하여 얻은 수치실험 결과를 비교하여, 거의 일치함을 확인하였고, 정확성을 평가하기 위하여 흐름이 충분히 Smooth하고 세밀한 격자간격을 가진 조건에 대하여 운동량 방정식과 Bernoulli 방정식을 별도로 각각 이산화 하여 얻은 해가 서로 일치함을 확인하였다. 수심경사가 급변하는 유동장에 기존의 기법(Hong et al., 2008)과 비교 적용하여 개선 수립된 수치모형이 효율적이라는 것도 나타내었으며, 개선된 모델을 적용하여 육상지역 침수유동도 재현하였다.
본 연구의 수학모델은 Hong et al.(2008)에서 적용한 것으로
여기서,
여기서,
천해역의 경우 수심방향의 운동량방정식은 식 (6)에서와 같이 정수압으로 가정되어 표현된다. 즉, 부력영향과 수심의 급격한 변화로 인한 수심방향의 가속도는 무시된다고 보았다.
여기서,
경계조건은 불투과성(Impermeability) 가정에 의하여 자유수면(
그리고, 해저바닥과 자유수면에서의 운동량방정식에 대한 경계조건은 각각 식 (8), (9)와 같이 주어진다.
여기서 𝜈는 와동점성계수,
여기서
난류모델은
본 연구의 수치모형은 기본적으로 Hong et al.(2008)에 의해 개발된 모형을 기본으로 수평이류항의 공간이산화 방식을 개선 수정한 것으로서 사용한 격자체계는 Fig. 1에서 보듯이 꼭지좌표에 대한 Mapping에 의해 Physical cell을 Computational cell로 변환(식 (1) ~ 식 (5)의
ADI scheme (Leendertse and Gritton, 1971; Leendertse et al., 1973)이 시간이력 적분을 위하여 기본적으로 적용되며, 수직방향의 각 층간의 수평유속은 수직방향의 이류항과 점성항에 의해 서로 결합(coupling)되어 해의 불안전성을 없애기 위해 수직방향의 상호 연결항에 대하여 Fully implicit time integration을 적용한다. 공간적으로 이류항은 아래 식 (12)와 같이 2차 중앙차분을 사용하고 점성항은 1차 중앙차분을 사용하여 이산화 한다.
여기서,
난류모델에서 Eddy viscosity(𝜈
서론과 3.1절에서 기술한 바와 같이 육상지역 침수나 경사가 급한 수심(또는 장애물)영역에서의 흐름전이와 같은 유동장이 급변하는 유동을 효율적으로 해석하고 요동현상을 제어하여 위하여 수평이류항의 공간이산화 방식을 개선하였다.
Fig. 3은 질량 및 운동량보존의 Control volume을 나타낸 것으로 해저바닥은 Fig. 3(c)에서 보듯이 계단모양의 타일로 근사화하며, 격자 중앙에 조위를, 격자와 격자가 맞닿은 변의 중앙지점에는 유속을 부여하였다. 국부적으로 압력항을 정수압으로 가정하는 것이 의미가 없으므로 질량 및 운동량보존이 중요할 수밖에 없다. 그러므로, 흐름이 팽창하는 경우는 운동량 보존방정식에, 축소하는 경우는 Bernoulli 방정식을 사용하여 격자중앙 주위의 Control volume(Fig. 3참조)에 대하여 2차 중앙차분을 기본으로 공간이산화를 적용하여 근사(Approximation)화 하되 수치요동을 줄여 해의 정확성을 높이기 위해 이산화 차수를 1차로 줄이는 Minmod 제한자를 사용하였다. 즉, 흐름이 팽창하는 경우 식 (2), (3)의 운동량 방정식의 수평이류항 중에 와 같이 동일한 방향의 곱의 성분은
식 (13)에서 은 격자 중앙의
여기서
또한, 와 같이 서로 다른 방향의 곱의 성분은
식 (19)에서 는 격자간의
여기서
흐름이 축소하는 경우에는 식 (2), (3)을 적분한 에너지보존방정식을 사용하여 일정한 유량
그리고, 경계조건 근처에서 발생하는 유동의 불안정을 제거하기 위하여 Stelling(1984)에 의하여 제안된 기법을 적용하였다.
수립된 수치모형을 검증하기 위하여 Fig. 4와 같은 넓은 마루웨어에 대하여 Bernoulli식을 적용하여 유도한 웨어 한계수위(식 (27))과 본 연구의 개선수치모형을 적용하여 얻은 수치실험 결과를 비교하였다. 수치모형 실험은 수조 중앙에 길이 5.0m, 폭 1.0m 및 높이 0.5m의 웨어 설치 후 상류유속은 0.1부터 0.5m/sec까지 0.1씩 증가시키고, 상류수위는 0.3m부터 0.1m씩 증가시켜 0.5m까지 변화시키면서 웨어상부(한계)수위를 계산하였다. 시간간격은 0.5sec를 격자는 0.03125 × 0.03125m를 사용하였으며, 수조 상류와 하류는 각각 수위 및 유량조건을 적용하였다.
상기의 조건을 가지고 한계흐름(Critical flow) 상태가 된 후의웨어 수위
3.2절에서 이미 기술하였듯이 수평이류항의 공간이산화 방식을 개선할 때 흐름이 팽창하는 경우는 운동량 보존방정식에, 축소하는 경우는 Bernoulli식에 공간이산화를 적용하였다. 그러나 흐름이 충분히 Smooth하고 세밀한 격자간격을 가진 조건에 대하여는 운동량 방정식과 Bernoulli식을 각각 별도로 사용하여 이산화한 결과가 일치할 수밖에 없을 것이다. 그러므로, 본 연구에서 수립된 수치모형의 정확성을 평가하기 위하여, 폭이 2m이며 시점 및 종점부 수심이 각각 1m 및 2m인 바닥경사 1°(충분히 smooth한 흐름 유도)를 가진 2차원 수조 (길이 58m)에 대하여, 시점에서 종점 방향의 확대 흐름(Flow expansion) 및 종점에서 시점 방향의 수축 흐름(Flow contraction)의 경우에 운동량 방정식과 Bernoulli 식을 별도로 적용 이산화한 후 수치모형실험을 수행하였다. 수조 양단은 0.4m3/sec 유량 조건을 부여하였으며, 격자크기는 1m × 1m부터 가로 및 세로방향으로 반씩 줄여나가 0.03125m × 0.03125m까지 줄여나가면서 중앙지점에서의 유속을 Table 1에 비교하였다.
Comparison of velocity(m/sec) at mid of channel according to horizontal discretization method for flow expasion and flow contraction respectively
Table 1에서 보듯이 흐름이 팽창하는 경우는 운동량 보존방정식에, 축소하는 경우는 Bernoulli식에 공간이산화를 적용하는 경우가 상대적으로 큰 격자간격에 대하여서도 정해에 가까운 값을 나타내었으며, 격자간격이 세밀해질수록 운동량 방정식과 Bernoulli 식을 각각 별도로 사용하여 이산화한 결과가 거의 일치함을 알 수 있었다.
본 연구에서 수립된 개선모형의 효율성을 조사하기 위하여, 개선수치모형과 개선 전 수치모형(Hong et al., 2008)을 Fig. 4와 같은 넓은 마루 웨어의 상류 유속 0.3m/sec 수위 0.4m인 경우에 적용, 웨어
[Table 2.] Comparison of water level (m) at weir according to the change of grid size
Comparison of water level (m) at weir according to the change of grid size
Table 2로 부터 본 연구의 개선모형은 12.5cm × 12.5cm 격자를 사용하더라도 결과가 0.2616m로 식 (27)에 의한 웨어 수위
즉, Minmod 제한자의 적용으로 이산화 차수를 1차로 줄여 상대적으로 큰 격자간격을 사용하더라도 정해에 가까운 결과를 얻어 효율적이었다.
본 연구에서 수립된 개선모델을 적용하여 Fig. 6의 육상지역 침수유동을 재현하였다. 10m폭을 가진 수조에 대하여 격자간격 1m × 1m, 조위진폭 0.2m, 주기 30deg/hour의 개방경계조건을 적용하였으며, 적분시간간격 3sec로 2일간 수치시뮬레이션을 수행하였다. Fig. 7은 시간대별 유속분포를 나타낸 것으로 육상지역(Fig. 6의 우측영역)이 침수되는 상황을 잘 표현하였다. 그러므로, 본 연구에서 수립된 개선모형은 육상지역 침수유동과 같이 경사가 급변하는 수심영역에서의 흐름전이를 잘 표현하는 것으로 나타났다. 즉, 이를 위해 적용한 Minmod 제한자가 적절히 작동하는 것을 확인하였다.
본 연구에서 기존의 수치모형(Hong et al., 2008)이 유동장이 급변하는 유동의 효율적인 해석을 수행하기 위해 수정 개선되었다. 흐름이 팽창하는 경우는 운동량 보존방정식에, 축소하는 경우는 Bernoulli 식에 2차 중앙차분을 기본으로 공간이산화를 적용 개선하였으며, 수치요동을 줄여 해의 정확성을 높이기 위해 이산화 차수를 1차로 줄이는 Minmod 제한자를 사용하였다.
수립된 수정 수치모형은 넓은 마루 웨어에 적용, 웨어 상부수위(한계수위)의 이론결과와 비교하여 거의 일치함을 확인하였고, 흐름이 충분히 Smooth하고 세밀한 격자간격을 가진 조건에 대하여 운동량 방정식과 Bernoulli 식을 별도로 각각 이산화 하여 얻은 해가 서로 일치함도 확인하였다. 수심경사가 급변하는 유동장에 기존의 기법(Hong et al., 2008)과 비교 적용하여 개선 수립된 수치모형이 효율적이라는 것도 나타내었으며, 개선된 모델을 적용하여 육상지역 침수유동도 잘 재현하였다.
향후, 수정 수립된 수치모형의 CFL 안정성을 조사하고, 해안에서의 급경사의 해저지형에 대한 유동해석이나 해수범람 모델링과 같은 해안 및 해양 공간에서의 유동해석에 확대 적용함으로서 수정 수립된 모형의 적용범위, 유용성 및 효율성을 조사할 것이다.