Advanced
Effective Analysis for Rapidly Varying Flows through Improvement in Spatial Discretization of Horizontal Advection Terms
Effective Analysis for Rapidly Varying Flows through Improvement in Spatial Discretization of Horizontal Advection Terms
Journal of Ocean Engineering and Technology. 2014. Aug, 28(4): 324-330
Copyright © 2014, Korean Society of Ocean Engineers
  • Received : February 17, 2014
  • Accepted : August 19, 2014
  • Published : August 30, 2014
Download
PDF
e-PUB
PubReader
PPT
Export by style
Article
Author
Metrics
Cited by
TagCloud
About the Authors
남식 홍

Abstract
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.
Keywords
1. 서 론
조석과 같은 유동수치모델에서 수평이류항의 공간이산화 방식은 해의 정확성, 수치기법의 단순화 및 효율성에 큰 영향을 미친다. 예를 들어, 중앙차분(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 )을 적용하였는데 이는 수평이류항중에 u ∂u/∂ξ 와 같이 동일한 방향의 곱의 성분은 2차 중앙차분(Second order central difference)으로 v ∂u/∂η 와 같이 서로 다른 방향의 곱의 성분의 1계 미분항에 대하여 3차 역방향 유한차분 (Upwind finite difference)화를 2개의 2차항으로 분리하여 이산화 하는 것으로 유속의 공간적인 요동(Oscillation)을 급격히 감소시켰다.
상기의 기법 등을 통해 인위적 점성항을 제거하면서 유속의 공간적인 요동은 감소시켰으나, 육상지역 침수나 경사가 급변하는 수심(또는 장애물)영역에서의 흐름전이(Flow transition)와 같은 유동장이 급변하는 유동을 해석하는 경우, 격자나 시간간격을 상당히 줄이야 하므로 비효율적일 뿐만 아니라 유속경사가 상당히 클 때는 해의 요동현상도 나타나 정해(Correct solution)를 얻지 못할 수도 있다. 그러므로, 본 연구의 목적은 기존의 수치모형( Hong et al., 2008 )에 수평방향 이류항의 공간이산화 방식 개선을 통해 상기의 비효율성을 개선하고 요동현상을 제어하여 급변하는 유동의 효율적인 해석을 수행하기 위함이다. 흐름이 팽창하는 경우는 운동량 보존방정식에, 축소하는 경우는 Bernoulli 방정식에 2차 중앙차분을 기본으로 공간이산화를 적용하였으며, 수치요동을 줄여 해의 정확성을 높이기 위해 이산화 차수를 1차로 줄이는 Minmod 제한자(Limiter)를 사용하였다. 수립된 개선 수치모형을 검증하기 위해 넓은 마루웨어(Broad crested weir)에 대하여 Bernoulli 방정식을 적용하여 유도한 웨어 한계수위와 본 연구의 개선수치모형을 적용하여 얻은 수치실험 결과를 비교하여, 거의 일치함을 확인하였고, 정확성을 평가하기 위하여 흐름이 충분히 Smooth하고 세밀한 격자간격을 가진 조건에 대하여 운동량 방정식과 Bernoulli 방정식을 별도로 각각 이산화 하여 얻은 해가 서로 일치함을 확인하였다. 수심경사가 급변하는 유동장에 기존의 기법( Hong et al., 2008 )과 비교 적용하여 개선 수립된 수치모형이 효율적이라는 것도 나타내었으며, 개선된 모델을 적용하여 육상지역 침수유동도 재현하였다.
2. 수학 모델
본 연구의 수학모델은 Hong et al.(2008) 에서 적용한 것으로 σ 좌표계를 사용하였으며, z 는 실제 연직좌표계, ζ 는 기준면( z = 0)상의 자유수면 변위, d 는 기준면하의 수심, 그리고 H (= d + ζ )는 총 수심으로 나타낼 때 σ 좌표계에서 수심평균 연속방정식은 식 (1)과 같으며 수평방향 ξ, η 방향의 운동량방정식은 식 (2)및 식 (3)과 같다.
PPT Slide
Lager Image
여기서, Jξξ Jηη 는 수평좌표계( ξ, η )에서의 Jacobian을 U, V 는 수평방향의 수심평균유속을 나타내며 Q 는 단위면적당 유입, 용출량이다.
PPT Slide
Lager Image
PPT Slide
Lager Image
여기서, f 는 Coriolis 상수, 𝜈는 와동점성계수를 나타내며 Sξ Sη 는 각각 ξ, η 방향의 운동량 Source나 Sink이다. 또한, u v 는 각각 수평방향인 ξ, η 방향의 유속이고,
PPT Slide
Lager Image
는 연직 ζ 방향의 이동 σ -plane에 대한 상대유속으로 식 (4)의 연속방정식으로 부터 계산되며 고정좌표계에서의 실제 수심방향의 유속 w 는 식 (5)로부터 얻어진다.
PPT Slide
Lager Image
PPT Slide
Lager Image
천해역의 경우 수심방향의 운동량방정식은 식 (6)에서와 같이 정수압으로 가정되어 표현된다. 즉, 부력영향과 수심의 급격한 변화로 인한 수심방향의 가속도는 무시된다고 보았다.
PPT Slide
Lager Image
여기서, ρ g 는 각각 단위밀도와 중력가속도를 나타낸다.
경계조건은 불투과성(Impermeability) 가정에 의하여 자유수면( σ = 0)과 바닥( σ = −1)에서의 연직 ζ 방향의 이동 σ -plane에 대한 상대유속
PPT Slide
Lager Image
는 식 (7)과 같은 운동조건(Kinematic condition)으로 주어진다.
PPT Slide
Lager Image
그리고, 해저바닥과 자유수면에서의 운동량방정식에 대한 경계조건은 각각 식 (8), (9)와 같이 주어진다.
PPT Slide
Lager Image
PPT Slide
Lager Image
여기서 𝜈는 와동점성계수, τ , τ 는 각각 ξ, η 방향의 해저전단응력으로 해저바닥 바로 상층(Above layer)의 수평유속벡터
PPT Slide
Lager Image
와 절대치
PPT Slide
Lager Image
의 항으로 식 (10)과 같이 벡터로 표현된다.
PPT Slide
Lager Image
여기서 C 3D 는 Chezy계수이다.
난류모델은 k ε 모델을 적용하였으며 혼합길이(Mixing length) 는 Kinetic energy k 와 Energy dissipation ε 로부터 다음 식 (11)로 부터 결정된다.
PPT Slide
Lager Image
3. 수치 모형
- 3.1 수치 모형 개요
본 연구의 수치모형은 기본적으로 Hong et al.(2008) 에 의해 개발된 모형을 기본으로 수평이류항의 공간이산화 방식을 개선 수정한 것으로서 사용한 격자체계는 Fig. 1 에서 보듯이 꼭지좌표에 대한 Mapping에 의해 Physical cell을 Computational cell로 변환(식 (1) ~ 식 (5)의 Jξξ , Jηη 참조)시킨 후 Computational cell의 중앙에서는 조위(Water level) η 또는 압력(pressure) p 를, cell 각 면의 직각방향으로는 유속 u , v , w 를 위치하는 Staggered grid ( Fig. 1 , Fig. 2 참조)로서 3차원 지배방정식을 시간과 공간에 대하여 유한 차분화 한다. 연직방향으로는 수심의 변화에 관계없이 수평방향의 전 영역에 걸쳐 동일한 층(Layer)의 개수를 갖는 σ 좌표격자를 적용한다. σ 격자에서는 각 층의 두께를 서로 다르게 줄 수 있어 국부적인 관심영역(예로 취송류 해석에서 수표면 근처 층)의 해석에 유용하다.
PPT Slide
Lager Image
Transformation from physical cell to computational cell
PPT Slide
Lager Image
View of staggered grid system (3D and top view)
ADI scheme ( Leendertse and Gritton, 1971 ; Leendertse et al., 1973 )이 시간이력 적분을 위하여 기본적으로 적용되며, 수직방향의 각 층간의 수평유속은 수직방향의 이류항과 점성항에 의해 서로 결합(coupling)되어 해의 불안전성을 없애기 위해 수직방향의 상호 연결항에 대하여 Fully implicit time integration을 적용한다. 공간적으로 이류항은 아래 식 (12)와 같이 2차 중앙차분을 사용하고 점성항은 1차 중앙차분을 사용하여 이산화 한다.
PPT Slide
Lager Image
여기서, hm,n,k 는 수직방향층격자두께로서 hm,n,k = ΔσkHm,n 로 정의된다. 또한, 이산화 방정식이 상호 결합된 비선형항(Nonlinear term)을 가지고 있으므로 반복계산(Iteration)이 필요하며 본 연구에서는 이산화 운동량방정식에 수렴하는 인접수위( H + ζ )경사를 곱하여 비선형 항을 제거한 후, 각 반복단계에서 질량보존(Mass conservation)을 검토하였다.
난류모델에서 Eddy viscosity(𝜈 v )는 앞 반 시간단계(Previous half time step)의 값을 사용하며, k, ε 과 함께 해저바닥 및 수면에서의 경계조건을 정식화하기 위해 각 층의 경계에서 계산하며 수평방향으로는 계산격자의 중앙에서 계산하는 Staggered 격자체계로 구성하였다. 이류항에 대하여 1차 역방향 유한차분을 적용하여 항상 양해(Positive solution)를 갖도록 하였다.
- 3.2 수평이류항의 공간이산화 방식 개선
서론과 3.1절에서 기술한 바와 같이 육상지역 침수나 경사가 급한 수심(또는 장애물)영역에서의 흐름전이와 같은 유동장이 급변하는 유동을 효율적으로 해석하고 요동현상을 제어하여 위하여 수평이류항의 공간이산화 방식을 개선하였다.
Fig. 3 은 질량 및 운동량보존의 Control volume을 나타낸 것으로 해저바닥은 Fig. 3(c) 에서 보듯이 계단모양의 타일로 근사화하며, 격자 중앙에 조위를, 격자와 격자가 맞닿은 변의 중앙지점에는 유속을 부여하였다. 국부적으로 압력항을 정수압으로 가정하는 것이 의미가 없으므로 질량 및 운동량보존이 중요할 수밖에 없다. 그러므로, 흐름이 팽창하는 경우는 운동량 보존방정식에, 축소하는 경우는 Bernoulli 방정식을 사용하여 격자중앙 주위의 Control volume( Fig. 3 참조)에 대하여 2차 중앙차분을 기본으로 공간이산화를 적용하여 근사(Approximation)화 하되 수치요동을 줄여 해의 정확성을 높이기 위해 이산화 차수를 1차로 줄이는 Minmod 제한자를 사용하였다. 즉, 흐름이 팽창하는 경우 식 (2), (3)의 운동량 방정식의 수평이류항 중에
PPT Slide
Lager Image
와 같이 동일한 방향의 곱의 성분은 q = Hu 관계를 이용하여 식 (13)과 같이 나타내었다.
PPT Slide
Lager Image
Control volume (a) for mass, (b) for momentum in horizontal direction, (c) for momentum in vertical direction
PPT Slide
Lager Image
PPT Slide
Lager Image
식 (13)에서
PPT Slide
Lager Image
은 격자 중앙의 ξ 방향 평균유량( Fig. 3 참조)으로 식 (14)와 같으며, 여기서
PPT Slide
Lager Image
은 식 (15)와 같이 산정하였다.
PPT Slide
Lager Image
Fig.3 (b) 에서 보듯이 격자 중앙에서의 유속으로 식 (16)으로 표현되고, 식 (16)에서 l ( r )는 앞서 기술하였듯이 경계영역에서 발생하는 수치요동을 줄여 해의 정확성을 높이기 위해 이산화 차수를 1차로 줄이는 Minmod 제한자로서 식 (17)과 같다.
PPT Slide
Lager Image
PPT Slide
Lager Image
PPT Slide
Lager Image
PPT Slide
Lager Image
여기서 r 은 유속 증분 비율로 식 (18)과 같다.
PPT Slide
Lager Image
또한,
PPT Slide
Lager Image
와 같이 서로 다른 방향의 곱의 성분은 p = Hv 관계를 이용하여 식 (19)와 같이 나타내었다.
PPT Slide
Lager Image
식 (19)에서
PPT Slide
Lager Image
는 격자간의 ξ 방향 평균유량( Fig. 3 참조)으로 식 (20)과 같으며, 여기서
PPT Slide
Lager Image
은 식 (21)과 같이 산정하였다.
PPT Slide
Lager Image
Fig.3 (b) 에서 보듯이 격자간의 유속으로 식 (22)로 표현되고, 식 (22)에서 l ( r )는 앞서 기술하였듯이 이산화 차수를 1차로 줄이는 Minmod 제한자로 식 (23)과 같다.
PPT Slide
Lager Image
PPT Slide
Lager Image
PPT Slide
Lager Image
PPT Slide
Lager Image
여기서 r 은 유속 증분 비율로 식 (24)와 같다.
PPT Slide
Lager Image
흐름이 축소하는 경우에는 식 (2), (3)을 적분한 에너지보존방정식을 사용하여 일정한 유량 q 를 가진 정상상태 조건하에서 격자중앙 주위의 Control volume( Fig. 3 참조)에 대하여 공간이 산화 및 근사화하여 식 (25) 및 (26)에 나타내었다.
PPT Slide
Lager Image
PPT Slide
Lager Image
그리고, 경계조건 근처에서 발생하는 유동의 불안정을 제거하기 위하여 Stelling(1984) 에 의하여 제안된 기법을 적용하였다.
4. 수치모형 평가
- 4.1 수치모형 검증
수립된 수치모형을 검증하기 위하여 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를 사용하였으며, 수조 상류와 하류는 각각 수위 및 유량조건을 적용하였다.
PPT Slide
Lager Image
Flow over broad-crested weir
PPT Slide
Lager Image
상기의 조건을 가지고 한계흐름(Critical flow) 상태가 된 후의웨어 수위 yc 를 계산하여 식 (27)에 의한 결과와 비교하여 Fig. 5 에 나타내었다. Fig. 5 에서 보듯이 상류유속 V 1 이 0.1m/sec 및 0.5m/sec일 경우에 3%정도 차이가 날 뿐 나머지 경우는 거의 일치하는 것으로 나타났다.
PPT Slide
Lager Image
Comparison between numerical experiment result and results by equation (27)
- 4.2 수치모형 정확성 및 효율성
3.2절에서 이미 기술하였듯이 수평이류항의 공간이산화 방식을 개선할 때 흐름이 팽창하는 경우는 운동량 보존방정식에, 축소하는 경우는 Bernoulli식에 공간이산화를 적용하였다. 그러나 흐름이 충분히 Smooth하고 세밀한 격자간격을 가진 조건에 대하여는 운동량 방정식과 Bernoulli식을 각각 별도로 사용하여 이산화한 결과가 일치할 수밖에 없을 것이다. 그러므로, 본 연구에서 수립된 수치모형의 정확성을 평가하기 위하여, 폭이 2m이며 시점 및 종점부 수심이 각각 1m 및 2m인 바닥경사 1°(충분히 smooth한 흐름 유도)를 가진 2차원 수조 (길이 58m)에 대하여, 시점에서 종점 방향의 확대 흐름(Flow expansion) 및 종점에서 시점 방향의 수축 흐름(Flow contraction)의 경우에 운동량 방정식과 Bernoulli 식을 별도로 적용 이산화한 후 수치모형실험을 수행하였다. 수조 양단은 0.4m 3 /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
PPT Slide
Lager Image
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인 경우에 적용, 웨어 yc 를 격지크기를 1m × 1m부터 가로 및 세로방향으로 반씩 줄여나가 0.03125m × 0.03125m까지 줄여나가면서 비교하여 Table 2 에 나타내었다.
Comparison of water level (m) at weir according to the change of grid size
PPT Slide
Lager Image
Comparison of water level (m) at weir according to the change of grid size
Table 2 로 부터 본 연구의 개선모형은 12.5cm × 12.5cm 격자를 사용하더라도 결과가 0.2616m로 식 (27)에 의한 웨어 수위 yc 0.2697m에 비해 3.6%정도 차이가 나나, 개선 전 모형은 3.125cm × 3.125cm 격자를 사용하여도 0.2242로 17.1% 차이가 났다.
즉, Minmod 제한자의 적용으로 이산화 차수를 1차로 줄여 상대적으로 큰 격자간격을 사용하더라도 정해에 가까운 결과를 얻어 효율적이었다.
- 4.3 수치모형 적용
본 연구에서 수립된 개선모델을 적용하여 Fig. 6 의 육상지역 침수유동을 재현하였다. 10m폭을 가진 수조에 대하여 격자간격 1m × 1m, 조위진폭 0.2m, 주기 30deg/hour의 개방경계조건을 적용하였으며, 적분시간간격 3sec로 2일간 수치시뮬레이션을 수행하였다. Fig. 7 은 시간대별 유속분포를 나타낸 것으로 육상지역( Fig. 6 의 우측영역)이 침수되는 상황을 잘 표현하였다. 그러므로, 본 연구에서 수립된 개선모형은 육상지역 침수유동과 같이 경사가 급변하는 수심영역에서의 흐름전이를 잘 표현하는 것으로 나타났다. 즉, 이를 위해 적용한 Minmod 제한자가 적절히 작동하는 것을 확인하였다.
PPT Slide
Lager Image
Schematic figure for numerical simulation of the inundation of dry land
PPT Slide
Lager Image
Velocity vector distribution through the numerical simulation for the inundation of dry land; (a) flood, full dry (b) flood, partial wet (c) flood, full wet (d) neap, full wet (e) neap, partial dry (f) neap, full dry
5. 결 론
본 연구에서 기존의 수치모형( Hong et al., 2008 )이 유동장이 급변하는 유동의 효율적인 해석을 수행하기 위해 수정 개선되었다. 흐름이 팽창하는 경우는 운동량 보존방정식에, 축소하는 경우는 Bernoulli 식에 2차 중앙차분을 기본으로 공간이산화를 적용 개선하였으며, 수치요동을 줄여 해의 정확성을 높이기 위해 이산화 차수를 1차로 줄이는 Minmod 제한자를 사용하였다.
수립된 수정 수치모형은 넓은 마루 웨어에 적용, 웨어 상부수위(한계수위)의 이론결과와 비교하여 거의 일치함을 확인하였고, 흐름이 충분히 Smooth하고 세밀한 격자간격을 가진 조건에 대하여 운동량 방정식과 Bernoulli 식을 별도로 각각 이산화 하여 얻은 해가 서로 일치함도 확인하였다. 수심경사가 급변하는 유동장에 기존의 기법( Hong et al., 2008 )과 비교 적용하여 개선 수립된 수치모형이 효율적이라는 것도 나타내었으며, 개선된 모델을 적용하여 육상지역 침수유동도 잘 재현하였다.
향후, 수정 수립된 수치모형의 CFL 안정성을 조사하고, 해안에서의 급경사의 해저지형에 대한 유동해석이나 해수범람 모델링과 같은 해안 및 해양 공간에서의 유동해석에 확대 적용함으로서 수정 수립된 모형의 적용범위, 유용성 및 효율성을 조사할 것이다.
References
Hong N.S. , Kim G.Y. , Kang Y.G. 2008 Three-Dimensional Numerical Model for Flow with Silt Protector Journal of Ocean Engineering and Technology http://koix.ksci.re.kr/KISTI1.1003/JNL.JAKO200827448605589 22 (3) 1 - 7
Gresho P.M. , Lee R.L 1981 Don’t suppress the wiggles they are telling you something Jounal of Computer and Fluids http://dx.doi.org/10.1016/0045-7930(81)90026-8 9 223 - 253
Leendertse J.J. 1967 Aspects of a Computational Model for Long-period Water-wave Propagatin. RM-5294-RR Rand Corporation Santa Monica
Leendertse J.J. , Gritton E.C. 1971 A Three-Dimensional Model for Estuaries and Coastal Seas: Vol II, Computation Procedures. Report R-708-NYC The Rand Corporation Santa Monica
Leendertse J.J. , Alexander R.C. , Liu S.K. 1973 A Three-Dimensional Model for Estuaries and Coastal Seas: I~IV, Report R-1417, 1764, 1884, 2187-OWRT The Rand Corporation Santa Monica
Stelling G.S. 1984 On the Construction of Computational Methods for Shallow Water Problems, Technical Report, 35 Rijkswaterstaat
Stelling G.S. , Leendertse J.J. 1991 Approximation of Convective Processes by Cyclic AOI Methods Proceedings of 2nd ASCE Conference on Estuarine and Coastal Modelling
Stelling G.S. , Leendertse J.J. , Spaulding M. L. , Bedford K. , Blumberg A. 1992 Approximation of Convective Processes by Cyclic AOI Methods Estuarine and coastal modeling, Proceedings 2nd ASCE Conference on Estuarine and Coastal Modelling
Vreugdenhil C.B. 1994 Numerical Methods for Shallow Water Flow. Water Science and Technology Kluwer Academic Publisher