Development of two-dimensional finite volume flow model with edge-based hydrostatic reconstruction and its application to urban flood

Eun Taek Shin; Dong Sop Rhee; Tae Soo Eum; Chang Geun Song

doi:10.3741/JKWRA.2025.58.11.1127

Preview

Special Issue: 지능형 도시홍수 예측

Journal of Korea Water Resources Association. 30 November 2025. 1127-1140
https://doi.org/10.3741/JKWRA.2025.58.11.1127

Development of two-dimensional finite volume flow model with edge-based hydrostatic reconstruction and its application to urban flood

엣지 기반 정수압 재구성 2차원 유한체적 흐름 모형 개발 및 도시 홍수 적용

Eun Taek Shin^a

Dong Sop Rhee^b

Tae Soo Eum^c

Chang Geun Song^d^*

신 은택^a

이 동섭^b

엄 태수^c

송 창근^d^*

^aResearch Professor, Department of Safety Engineering, Incheon National University, Incheon, Korea

^bResearch Fellow, Department of Hydro Science and Engineering Research, Korea Institute of Civil Engineering and Building Technology, Goyang, Korea

^cResearch Professor, Department of Safety Engineering, Incheon National University, Incheon, Korea

^dProfessor, Department of Safety Engineering, Incheon National University, Incheon, Korea

^a인천대학교 안전공학과 학술연구교수

^b한국건설기술연구원 수자원하천연구본부 연구위원

^c인천대학교 안전공학과 학술연구교수

^d인천대학교 안전공학과 교수

^{*Corresponding Author}

License (open-access, https://creativecommons.org/licenses/by/4.0):

©It is identical to the Creative Commons Attribution Non-commercial License (https://creativecommons.org/licenses/by/4.0)

ABSTRACT

This study developed a two-dimensional finite volume flow model with edge-based hydrostatic reconstruction for accurate and stable simulation of flash floods in complex urban environments. The proposed model incorporates MUSCL slope limiting and the HLLC approximate Riemann solver on unstructured triangular cells and employs a well-balanced scheme to maintain the equilibrium between flux and source terms, thereby satisfying the C-property. To enhance numerical stability and computational efficiency for large-scale urban flood simulations, implicit friction treatment and a high-performance parallel computing framework (Taichi) were implemented. The model performance was verified through benchmark tests, including Toro’s Riemann problems and Thacker’s planar solution, demonstrating its capability to capture shock propagation and wetting-drying transitions. The model was further applied to the 2007 Pasha Bulker storm event in Newcastle, Australia, and validated against experimental data as well as results from FLUMEN and HDM-2D models. The analysis confirmed that the proposed model successfully reproduced key urban flood characteristics, such as water depth rise in front of buildings and flow acceleration through narrow passages. Quantitative comparisons at 13 measurement points showed errors comparable to or lower than those of the existing models.

Keywords

Finite volume 2D model

Edge-based hydrostatic reconstruction

MUSCL-HLLC scheme

Urban flood simulation

Wetting-drying

본 연구에서는 복잡한 도시 지형에서 발생하는 돌발 홍수를 정확하고 안정적으로 모의하기 위하여 엣지 기반 정수압 재구성 기법을 적용한 2차원 유한체적 흐름해석 모형을 개발하였다. 개발된 모형은 비구조 삼각 격자에서 MUSCL 경사 제한과 HLLC 근사 리만 솔버를 결합하고, 보존항과 소스항의 균형을 유지하는 well-balanced 기법을 적용하여 C-속성을 만족하도록 구현되었다. 또한, 마찰항의 음해적 처리와 고속 병렬 연산 프레임워크(Taichi)를 도입하여 대규모 도시 홍수 해석에서의 수치적 안정성과 계산 효율성을 확보하였다. 모형의 성능은 Toro의 Riemann 문제와 Thacker의 곡면 회전 문제와 같은 검증 예제를 통해 충격파 전파와 마름-젖음 경계 처리 능력을 확인하였다. 2007년 ‘Pasha Bulker’ 폭풍으로 피해를 입은 호주 뉴캐슬 지역을 대상으로 적용하였으며, 실험 결과 및 기존 모형(FLUMEN, HDM-2D)과 비교·검증하였다. 분석 결과, 건물 전면부 수위 상승, 협소 통로에서의 병목 현상 등 도시 홍수 특성이 정확히 재현되었으며, 13개 검증 지점에서 수심 및 유속 오차는 기존 모형과 유사하거나 일부 구간에서는 더 낮게 나타났다.

키워드

2차원 유한체적 모형

엣지 기반 정수압 재구성

MUSCL-HLLC

도시 홍수

마름-젖음

MAIN

1. 서 론
2. 지배방정식과 수치기법
2.1 천수방정식
2.2 유한체적법
2.3 Runge-Kutta 2차 기법의 적용
2.4 흐름률과 소스항의 계산
2.5 고속 병렬연산 프레임워크를 활용한 향상된 데이터 구조의 활용
3. 모형의 검증
3.1 Toro's Riemann problem
3.2 Thacker’s planar solution
4. 도시 홍수 사례 적용
5. 결론 및 향후 연구

1. 서 론

최근 도시화와 기후 변화로 인해 도시 지역에서 홍수 발생 빈도와 피해가 증가하고 있다. 전 세계의 재난 데이터를 수집 및 통계하여 발표하는 EM-DAT의 재난역학연구센터(Centre for Research on the Epidemiology of Disasters, CRED)에 의하면 2024년 한 해 동안 발생한 재해 중 약 73.5%가 홍수와 폭풍이었다(CRED, 2025). 홍수를 사전에 예측하고 관리하기 위해서는 정밀한 강우예측, 지표면 강우 적체, 지표수 흐름 해석, 관망 흐름 해석 등과 같은 도전적인 과제가 여전히 존재하며, 이는 현재도 많은 공학자와 과학자들의 주요 관심사로 남아있다.

흐름해석 모형은 다양한 자연 현상과 공학적 문제를 이해하고 예측하는데 활용되는 중요한 공학적 도구이다. 흐름해석 모형의 공간 차원은 관심 영역의 크기와 요구되는 해석 정밀도 및 계산 효율성 등에 따라 1차원, 2차원 및 3차원 모형으로 구분된다(Toro, 2009; El Kadi Abderrezzak et al., 2009). 극한 강우나 하천 범람으로 발생하는 도시 홍수는 도시 지형과 인공 구조물의 영향을 받아 다방향으로 확산되므로, 단방향 흐름을 가정하는 1차원 모형으로는 이러한 현상을 적절히 재현하는 데 한계가 있다. 3차원 흐름해석 모형은 구조물 유동 및 정밀 난류 해석에 효과적이지만, 대규모 도시 홍수의 모의에는 적용상의 제약이 크다. 이는 지역별 고해상도 데이터셋의 불균형, 고성능 컴퓨팅 계산 자원 확보의 어려움, 3차원 흐름해석 결과의 검보정 과정의 어려움 때문이다(Ma et al., 2007). 반면에 2차원 데이터로 활용할 수 있는 원격 감지 기술(고해상도 및 고정밀 LiDAR 및 Radar (SAR) 데이터)의 발전과 컴퓨팅 능력의 향상은 2차원 모델의 적용 가능성을 높이고 있다(Liang and Borthwick, 2009). 따라서 도시 홍수의 유량, 수심, 유속 및 침수 면적 예측에는 2차원 흐름해석 모형이 효과적이라 할 수 있다.

최근 수십 년 간 Godunov형 유한체적법은 중력에 의한 자유수면 문제를 다루는 다양한 수치해석과 전산유체역학 분야에서 널리 사용되고 있다. 이 방법은 천이류 흐름, 충격파의 전파, 마름젖음 경계의 이동 등과 같이 복잡한 흐름 현상을 효율적으로 모사할 수 있다는 장점이 있다(Liang and Marche, 2009; Song et al., 2011). 그러나 이 방법은 불규칙한 지형에서 바닥경사 등의 외력항 처리와 관련된 불균형 문제로 인해 수치적 불안정성이 존재한다. 예를 들어, Godunov형 기법에서 공간 정확도를 2차로 높이는 경우 Gibbs 현상과 같은 비물리적 진동이 나타날 수 있으며, 이는 복잡한 지형과 지속적으로 변하는 마름-젖음 경계를 포함한 도시 침수 모의에 적용하는 데 제약으로 작용한다(van Leer, 1979; Toro and Garcia-Navarro, 2007). 이러한 문제를 해결하기 위해 Total Variation Diminishing (TVD) 기법과 같은 고차 정확도 체계가 개발되었으며, 특히 van Leer (1979)가 제안한 MUSCL (Monotone Upstream-centered Scheme for Conservation Laws) 기법은 2차 공간 정확도를 유지하면서도 단조성을 보장하여 널리 활용되고 있다. 한편, Bermudez and Vazquez (1994)가 제시한 C-속성은 정지 흐름에서 보존 법칙 충족 여부를 판단하기 위한 기준으로, 불규칙 지형에서의 안정적 해석 가능성을 검증하는 지표로 활용된다.

상기 연구들이 도시 홍수 수치해석의 기반을 마련하였으나, 비구조 삼각 격자에서 엣지 기반으로 정수압 재구성, 경사 재구성, 소스항 평가를 일관되게 통합하고, 여기에 HLLC-MUSCL 조합과 음해적 마찰 처리를 결합하여 마름/젖음 안정성을 동시에 확보한 사례는 매우 드물다. 특히 도시 규모의 고해상도 격자(수십만~수백만 셀)를 실시간 또는 준실시간 수준에서 계산하기 위해서는 GPU나 다중코어 기반의 고성능 병렬 연산 프레임워크가 필수적임에도, 기존 연구들은 계산 격자 해상도와 연산 효율성 간의 상충 관계(trade-off)를 전제로 하여 두 요소 중 하나를 제한적으로 고려하는 경향이 있었다. 본 연구는 이러한 한계를 극복하기 위해 수치적 안정성과 계산 효율성을 동시에 향상시켜, 고해상도 격자 기반의 도시 홍수 시뮬레이션을 실질적으로 구현 가능한 수준으로 발전시켰다. 이를 위해 엣지 기반 정수압 재구성 기법을 핵심으로 하는 2차원 유한체적 흐름 모형을 개발하였다. 제안된 모형의 핵심은 정수압 재구성-MUSCL/HLLC 체계를 엣지 단위에서 통합하여 재구성하고 리만 상태 보정을 수행함으로써 플럭스와 소스항의 균형을 유지하고 음수심 발생을 방지함으로써 C-속성(well-balanced property)을 안정적으로 보존하는 것이다. 또한 엣지 기반 경사 제한과 HLLC 조합을 통해 젖음-마름 전선과 충격파·희박파를 안정적으로 포착하였으며, 마찰항은 음해적 한계값 접근법으로 처리하여 수치적 경직성을 효과적으로 완화하였다. 계산 효율성 측면에서는 Taichi 기반 JIT/GPU 병렬화 구조를 도입하여 대규모 격자에서도 연산 속도를 크게 향상시켰으며, 이를 통해 대규모 도시 홍수 시뮬레이션을 현실적인 시간 내에 수행할 수 있도록 구현하였다. 개발된 모형은 벤치마크 문제(Toro의 Riemann 문제, Thacker의 곡면 회전 문제)와 2007년 호주 뉴캐슬 Pasha Bulker 폭풍 사례 적용을 통해 수치적 안정성과 재현성을 검증하였으며, 기존 모형(FLUMEN, HDM-2D)과의 정량 비교 결과에서도 높은 신뢰성과 실용성을 입증하였다.

2. 지배방정식과 수치기법

2.1 천수방정식

본 연구에서는 도사 홍수를 수치모의하기 위하여 천수방정식을 사용하였다. 천수방정식은 Navier-Stokes 방정식에서 압력에 의한 가속도가 중력과 균형을 이루는 정수압 상태를 가정하고, 수평 규모가 연직 규모보다 훨씬 크며 연직 속도가 무시될 수 있다고 가정하여 유도된다. 천수방정식은 유동의 깊이와 수평 운동량을 근사화 하는 간소화된 모형으로, 시간에 따라 변하는 비선형 편미분 방정식 시스템을 보존변수로 이루어진 벡터 형태로 나타낸다.

(1a)

\frac{\partial U}{\partial t} + \frac{\partial E}{\partial x} + \frac{\partial G}{\partial y} = S

(1b)

U = [\begin{matrix} h \\ h u \\ h v \end{matrix}], E (U) = [\begin{matrix} h \\ h u^{2} + \frac{1}{2} g h^{2} \\ h v \end{matrix}], G (U) = [\begin{matrix} h \\ h u \\ h v^{2} + \frac{1}{2} g h^{2} \end{matrix}] .

(1c)

S = S_{0} + S_{f} = [\begin{matrix} 0 \\ g h \partial z_{b} / \partial x \\ g h \partial z_{b} / \partial y \end{matrix}] + [\begin{matrix} 0 \\ - C_{f} u \sqrt{u^{2} + v^{2}} \\ - C_{f} v \sqrt{u^{2} + v^{2}} \end{matrix}]

여기서, $h$ 는 수심, $u$ 와 $v$ 는 각각 $x$ 와 $y$ 방향으로의 유속, E(u)는 $x$ 방향으로의 흐름률 벡터, G(u)는 $y$ 방향으로의 흐름률 벡터이고, $S$ 는 소스 벡터로, 경사 소스 $S_{0}$ 항 와 마찰 소스 항 $S_{f}$ 로 세분화된다. $z_{b}$ 는 바닥 고도를 나타내며, $C_{f}$ 는 일반적으로 $g n^{2} / h^{1 / 3}$ 로 계산되는 바닥 거칠기 계수, $n$ 은 Manning 계수이다.

2.2 유한체적법

유한체적법은 유체 역학, 열 전달, 전산 유체 동역학(CFD) 등에서 널리 사용되는 수치 해석 기법으로, 계산 영역을 작은 제어 체적(셀)으로 분할하고 각 체적에 대해 질량, 운동량, 에너지 등의 보존 법칙을 적분 형태로 적용하여 방정식을 이산화하는 방법론이다. 이 기법은 복잡한 기하학적 형상이나 비구조 격자에도 유연하게 적용할 수 있으며, 플럭스 계산을 통해 안정성과 정확성을 확보할 수 있다는 장점이 있다. 본 연구에서는 비구조화된 격자에서 유한 체적법을 적용하였으며, 그 전개는 다음과 같다.

보존변수로 이루어진 벡터형태로 표현된 Eq. (1)을 셀에 대하여 적분하면 Eq. (2)와 같이 표현된다.

(2)

\int_{Ω} \frac{\partial U}{\partial t} d Ω + \int_{Ω} (\frac{\partial E}{\partial x} + \frac{\partial G}{\partial y}) d Ω = \int_{Ω} S d Ω

여기서 는 셀 의 체적을 나타낸다. 발산 정리를 적용하면, Eq. (2)는 다음과 같이 변환된다.

(3)

\int_{Ω} \frac{\partial U}{\partial t} d Ω + \oint_{Ω} F (U) ∙ n d Γ = \int_{Ω} S d Ω

여기서 𝛤는 격자를 구성하는 셀의 표면 경계이고, $n = {(n_{x}, n_{y})}^{T}$ 은 경계면에서의 외향 단위법선 벡터이다. F(U)ㆍn은 외향 단위 플럭스벡터를 의미하며, 아래와 같이 정의된다.

(4)

F (U) ∙ n = E n_{x} + G n_{y} = [\begin{matrix} h u_{⊥} \\ h u u_{⊥} + \frac{1}{2} g h^{2} n_{x} \\ h v u_{⊥} + \frac{1}{2} g h^{2} n_{y} \end{matrix}]

여기서 $u_{⊥} = u n_{x} + v n_{y}$ 와 $u_{∥} = - u n_{y} + v n_{x}$ 는 셀을 구성하는 경계면의 수직 및 접선 방향의 속도를 의미한다. 최종적으로 삼각형 셀의 경우, 경계에 대한 F(U)ㆍn의 원형 적분은 다음과 같이 계산할 수 있다.

(5)

\oint_{Ω} F (U) ∙ n d Γ = \sum_{k = 1}^{3} F_{k} (U) ∙ n k l_{k}

2.3 Runge-Kutta 2차 기법의 적용

본 연구에서는 2단계 양해적 Runge-Kutta 기법(RK2)을 적용하여 시간에 대한 2차 정확도를 확보하였다. RK2 기법은 오일러 방식에 비해 더 높은 시간 정확도를 가지며, 수치적 안정성과 정확성을 개선하기 위해 사용된다. 오일러 방식은 1차 시간 정확도를 가지며, 단일 시간 단계에서 단순히 현재 상태를 기반으로 다음 상태를 예측하기 때문에 오차가 누적될 가능성이 크다. 특히, 복잡한 비선형 문제나 빠르게 변화하는 시스템에서는 이러한 오차가 계산 결과를 크게 왜곡할 수 있다 (Sod, 1985). 반면, RK2 기법은 중간 단계 계산을 통해 시간적 변화율을 더 정밀하게 추정하여 이러한 한계를 극복하고, 2차 정확도를 제공함으로써 보다 안정적이고 정확한 수치 해를 얻을 수 있다. 새로운 시간 단계에서 셀 $i$ 의 상태 는 다음과 같이 갱신된다.

(6)

U_{i}^{n + 1} = U_{i}^{n} + Δ t \cdot Φ_{i}

여기서 $Φ_{i}$ 는 다음과 같이 계산된다.

(7)

Φ_{i} = \frac{1}{2} (K_{1} + K_{2})

K₁과 K₂는 다음과 같이 정의된다.

(8a)

K_{1} = - \frac{1}{A_{i}} \sum_{j} F_{i j} \cdot n_{i j} L_{i j} + S_{i}

(8b)

K_{2} = - \frac{1}{A_{i}} \sum_{j} F_{i j}^{*} \cdot n_{i j} L_{i j} + S_{i}^{*}

위 Eqs. (6), (7), (8)에서 $Δ t$ 는 시간 단계, $A_{i}$ 는 셀 $i$ 의 면적, $S_{i}$ 는 소스 항, $F_{i j}$ 는 플럭스, $n_{i j}$ 는 $j$ 번째 경계의 법선 벡터, $L_{i j}$ 는 $j$ 번째 경계의 길이를 나타낸다. 또한, $F_{i j}^{*}$ 와 $S_{i}^{*}$ 는 중간 단계에서 계산된 플럭스와 소스 항을 의미한다. RK2 기법은 오일러 방법에 비해 계산 비용은 다소 증가하지만, 시간 정확도와 안정성 측면에서 효율성이 높아 많은 연구에서 적용되고 있다(Hubbard, 1999; Liang and Marche, 2009; Wang et al., 2011).

2.4 흐름률과 소스항의 계산

2.4.1 MUSCL 기법의 적용

2차원 삼각 격자 시스템에서 공간 정확도를 향상시키면서 비물리적 진동을 억제하기 위해 MUSCL형 경사 제한 방법이 널리 사용되고 있다. 공간 2차 정확도를 위한 핵심 개념은 경사를 정의하고 그 경사가 단조성을 만족하도록 조정하는 것이다. 다차원 기법은 van Leer (1973)의 선구적인 연구를 기반으로 하며, 경사 제한 방법에 대한 심층적인 논의는 LeVeque (1992)의 연구에 상세히 제시되어 있다. 본 연구에서는 이러한 경사 제한 방법을 바탕으로 Venkatakrishnan (1995)이 제안한 셀 기반 제한 방법을 변형하여 엣지 방식으로 적용하였다.

본 연구에서는 공간 2차 정확도를 확보하기 위해 모형을 다음과 같이 구성하였다. 먼저 각 셀 내 근사값이 균일한 삼각 격자 시스템(Fig. 1(a))에서 1차 정확도의 물리량을 높이로 나타내면 Fig. 1(b)와 같이 표현된다. 이 경우 이웃 셀의 영향을 고려하지 않으므로 흐름률 계산 시 과도한 확산이 발생한다. 반면, 2차 정확도의 경우 인접 셀의 물리량을 고려하여 셀 값을 선형적으로 재구성을 하므로 Fig. 1(c)와 같이 표현된다. 최종적으로 선형적으로 재구성된 경사값에 단조성 원칙을 적용하고 경사 제한자의 최대값을 식별함으로써 비물리적인 진동을 억제할 수 있다. 여기서 단조성 원칙은 재구성된 함수값이 인접 셀의 최대값과 최소값 범위 내에 제한됨을 의미한다.

https://cdn.apub.kr/journalsite/sites/kwra/2025-058-11/N0200581113/images/kwra_58_11_13_F1.jpg

Fig. 1.

Example of cell reconstruction in triangular grid system: (a) Triangular grid system overview, (b) 1st order accuracy physical quantity representation, (c) 2nd order accuracy linear reconstruction

삼각 격자 시스템에 MUSCL 경사 제한 방법을 적용하는 절차는 다음과 같다. 먼저, 인접 셀에서 재구성된 제한되지 않은 변수 ${\bar{q}}_{i, M_{k}}$ 는 아래와 같이 주어진다.

(9)

{\bar{q}}_{i, M_{k}} = q_{i} + \bar{\nabla} q_{i} ∙ r_{i, M_{k}}

여기서 $q_{i}$ 는 $i$ 번째 셀 중심에서의 임의 변수; r는 셀 중심을 기반으로 한 위치 벡터; 첨자 $M_{k}$ 는 $k$ 번째 경계 중심; $\bar{\nabla} q_{i}$ 는 셀 $i$ 내 임의 변수(수면 높이 𝜂, 수심 $h$ , -방향 유량 $q_{x}$ , $y$ -방향 유량 $q_{y}$ )의 제한되지 않은 경사로, 주변 셀 $k$ ( $k$ =1, 2, 3)의 경사 값을 고려하여 Eq. (10)으로 계산된다.

(10)

\bar{\nabla} q_{i} = (- \frac{μ_{2}}{μ_{3}}, - \frac{μ_{1}}{μ_{3}})

여기서 $μ = (μ_{1}, μ_{2}, μ_{3}) = (P_{1} - P_{3}) \times (P_{2} - P_{3}); P_{k} = {(x_{k}, y_{k}, q_{k})}^{T}$ 이다. 위에서 언급한 단조성 원칙을 적용하기 위해, 제한되지 않은 경사는 Eq. (11)과 같이 경사 제한자를 통해 제한된다.

(11)

\nabla q_{i} = Φ q_{i}

여기서 𝛷는 경사 제한 함수로, Eq. (12)와 같이 주어진다.

(12a)

Φ_{i} = \min (Φ_{i, k})

(12b)

Φ_{i, k} = \{\begin{matrix} Φ (\frac{q_{i \max} - q_{i}}{M_{k} - q_{i}}) & if & {\bar{q}}_{M_{k}} - q_{i} > 0, \\ Φ (\frac{q_{i \min} - q_{i}}{M_{k} - q_{i}}) & if & {\bar{q}}_{M_{k}} - q_{i} < 0, 1 \\ 1 & if & {\bar{q}}_{M_{k}} - q_{i} = 0 . \end{matrix}

여기서 $q_{\max} = \max (q_{i}, q_{k 1}, q_{k 2}, q_{k 3})$ 이고 $q_{\min} = \min (q_{i}, q_{k 1}, q_{k 2}, q_{k 3})$ 로 계산한다. 제한 함수 𝛷는 Venkatakrishnan (1995)이 제안한 함수를 사용하였으며, 이는 Barth and Jespersen (1989) 기법을 수정한 형태로 Eq. (13)과 같다.

(13)

Φ (y) = \frac{y^{2} + 2 y}{y^{2} + y + 2}

마지막으로, 제한된 경사를 사용하여 경계 중심에서 보간된 모든 유량 변수를 다음과 같이 계산할 수 있다.

(14)

q_{i, M_{k}} = q_{i} + \nabla q_{i} \cdot r_{i, M_{k}}

이 접근법은 모든 셀에 동일하게 적용되므로 중복 계산이 포함된다. 따라서 제한된 경사를 계산하기 위한 격자 정보를 사전에 계산하고, 이를 엣지 기준으로 계산하는 방식을 적용하여 효율성을 높일 수 있다.

2.4.2 정수압 재구성

앞서 언급한 MUSCL형 재구성에서 경사 제한 방법을 적용하더라도 수심이 음수로 계산될 수 있으며, 이는 수치적 불안정성을 크게 초래한다. 이러한 문제를 해결하기 위해 Audusse et al. (2004)은 셀 인터페이스에서 변수를 수정하는 강력한 양수 유지 방법을 제안하였으며, 그 적용 가능성과 정확성은 여러 선행 연구에서 입증되었다(Marche et al., 2007; Liang and Marche, 2009; Wang et al., 2011). 본 연구에서는 Audusse et al.(2004)의 방법을 확장하여 엣지 기반 압력 재구성 기법을 개발하였으며, 그 절차는 다음과 같다. 첫 번째 단계는 인터페이스의 좌우 바닥 높이를 비교하여 최대값을 Eq. (15)와 같이 선택하는 것이다.

(15)

z_{M} = \max (z_{M}^{L}, z_{M}^{R})

여기서 $z_{M}^{L}$ 과 $z_{M}^{R}$ 은 인터페이스를 기준으로 좌우 방향의 바닥 고도를 의미한다. 두 번째 단계에서는 인터페이스 좌우 수심 $h_{M}^{L}$ 과 $h_{M}^{R}$ 을 조정하여 음수 수심이 발생하지 않도록 한다. 이 조정은 계산된 인터페이스 바닥 고도 $z_{M}$ 을 기반으로 하며, 엣지 좌우 수심은 아래와 같이 재구성된다.

(16)

h_{M}^{L} = \max (0, η_{M}^{L} - z_{M}), h_{M}^{R} = \max (0, η_{M}^{R} - z_{M})

마지막으로, 조정된 수심을 반영하여 인터페이스에서의 유량도 다음과 같이 갱신된다.

(17a)

q_{x M}^{L} = h_{M}^{L} u_{M}^{L}, q_{y M}^{L} = h_{M}^{L} v_{M}^{L}

(17b)

q_{x M}^{R} = h_{M}^{R} u_{M}^{R}, q_{y M}^{R} = h_{M}^{R} v_{M}^{R}

여기서 $u_{M}^{L,} u_{M}^{R} v_{M}^{L,} v_{M}^{R}$ 은 인터페이스를 기준으로 좌우 방향의 속도 벡터를 의미한다. 정수압 재구성 과정을 통해 수정된 변수는 흐름률 계산과 바닥 소스항 계산에 사용된다. 본 연구에서는 흐름률 계산을 위하여 Toro et al. (1994)에서 제시한 HLLC 근사 리만 솔버를 적용하였으며, Hou et al. (2013)에서 제시된 경사 플럭스 방법을 경사 소스 항을 평가하는 데 적용하였다. 흐름률과 바닥 소스항을 구하기 위한 두 방법론은 기존의 연구에서 다양하게 검증되었고 적용되었으므로 본 연구에서는 생략하였다.

2.4.3 마찰 소스항의 적용

Song et al. (2011)은 마찰항을 2단계 양해적 TVD Runge-Kutta 방식으로 처리할 경우 잠재적인 경직 문제가 발생할 수 있으며, 이는 수치 기법의 안정성과 효율성에 큰 영향을 미칠 수 있다고 보고하였다. 특히 복잡한 지형과 마름젖음 현상을 포함하는 도시 홍수를 고려할 때, 마찰항은 수치 모형의 안정성을 좌우하는 핵심 인자이다. 이를 안정적으로 처리하기 위해 유한체적 모형에서는 음해적 기법(Liang and Marche, 2009)과 반음해적 기법(Begnudelli and Sanders, 2006)이 널리 활용된다. 본 연구에서는 Liang and Marche (2009)가 제안한 음해적 마찰 한계값을 도입하여, 마찰항이 최대값에 도달할 경우 유체가 정지하도록 구현하였다. 이 기법에서 마찰항은 다음과 같이 주어진다.

(18a)

u^{n + 1} = \frac{u^{n}}{1 + Δ t {\hat{τ}}^{n}}, v^{n + 1} = \frac{v^{n}}{1 + Δ t {\hat{τ}}^{n}}

(18b)

{\hat{τ}}^{n} = g n^{2} \sqrt{{(u^{n})}^{2} + {(v^{n})}^{2}} {({\hat{h}}^{n})}^{4 / 3}

여기서 ${\hat{h}}^{n}$ , $u^{n}$ 및 $v^{n}$ 은 마찰 소스 항 처리를 위한 초기 상태 변수를 나타낸다. 이 기법은 속도 성분의 방향이 역전되지 않도록 하여 모형의 안정성을 향상시킨다.

2.5 고속 병렬연산 프레임워크를 활용한 향상된 데이터 구조의 활용

본 연구에서 개발된 2차원 흐름해석 모형은 파이썬(Python) 언어를 기반으로 구현되었다. 파이썬은 코드 가독성이 높고, NumPy, SciPy와 같은 다양한 라이브러리를 통해 수치 해석을 용이하게 지원하는 장점이 있다. 하지만 순수 파이썬 환경은 인터프리터 언어의 특성상 반복 연산이 많은 대규모 시뮬레이션에서 성능이 심각하게 낮아진다(van der Walt et al., 2011). 특히 고해상도 2차원 도시 홍수 모형의 경우, 메모리 접근과 병렬 처리의 비효율성으로 인해 계산 시간이 과도하게 증가하며, 이는 실시간 예측이나 대규모 격자 처리를 제약하는 요인으로 작용한다. 이러한 단점을 극복하기 위해, 본 연구에서는 Taichi 프레임워크를 도입하였다. Taichi는 파이썬에 내장된 도메인 특화 언어로, GPU 또는 멀티코어 CPU에서 고성능 병렬 연산을 손쉽게 구현하도록 지원한다(Hu et al., 2019). Taichi는 Just-In-Time (JIT) 컴파일을 통해 파이썬 코드를 하드웨어에 최적화된 코드로 변환하며, 자동으로 벡터화(SIMD)와 병렬화를 수행하여 기존 파이썬 코드의 성능을 수십 배 이상 향상시킬 수 있는 이점이 있다. 본 연구에서는 고속 병렬 연산 프레임워크 적용에 따른 성능 향상을 정량적으로 평가하기 위해, 각각 25,000개, 100,000개, 400,000개, 1,600,000개의 셀을 가진 가상의 비정규 삼각 지형 격자를 생성하였다. 이후, 생성된 격자의 엣지에서 Riemann Solver 계산을 40 s 동안 수행하는 동안의 연산 시간을 측정하였다. Table 1은 격자를 구성하는 셀 수(25,000, 100,000, 400,000, 1,600,000)에 대해 단일 CPU (Single), 멀티 CPU (Multi), GPU를 사용한 실행 시간을 비교한 결과다. 단일 CPU에서는 셀 수가 증가함에 따라 실행 시간이 선형적으로 늘어나며, 25,000 셀에서는 215 s, 1,600,000 셀에서는 39,266 s가 소요되었다. 멀티 CPU를 사용하면 병렬 처리로 인해 실행 시간이 단축되며, 25,000 셀에서 153 s, 1,600,000 셀에서 9,340 s로 개선되었다. GPU를 활용하면 성능이 더욱 크게 향상되며, 25,000 셀에서 45 s, 1,600,000 셀에서 2,644 s로 단축되었다. GPU 대비 단일 CPU의 속도 비율(GPU vs Single)은 4.72배에서 14.85배로, GPU 대비 멀티 CPU의 속도 비율(GPU vs Multi)은 3.37배에서 3.53배로 나타났다. 이는 GPU가 특히 대규모 격자에서 압도적인 성능 우위를 보이며, 셀 수가 많아질수록 GPU의 효율성이 두드러진다는 것을 보여줌과 Taichi의 JIT 컴파일과 병렬화 기술이 고해상도 2차원 흐름 해석에서 대형 그리드 처리를 가능하게 한다는 점을 뒷받침한다.

Table 1.

Comparison of computation time by grid cell count (Single CPU, Multi CPU, GPU)

Cell Count	CPU Single (s)	CPU Multi (s)	GPU (s)	GPU vs Single	GPU vs Multi
25,000	215	153	45	4.72×	3.37×
100,000	835	262	130	6.42×	2.01×
400,000	5,109	1,311	487	10.48×	2.69×
1,600,000	39,266	9,340	2,644	14.85×	3.53×

3. 모형의 검증

도시 홍수에 의한 흐름 특성은 도심지의 지형, 건물 밀도, 하수 시스템 등 다양한 요인에 따라 달라지며, 이로 인해 확산, 충격파 및 고속 흐름과 같은 복잡한 현상이 빈번하게 발생한다. 따라서 본 장에서는 개발된 2차원 흐름해석 모형의 도시 홍수 예측 성능을 검증하기 위해, 해석해가 존재하는 검증 예제를 대상으로 비교하였다.

3.1 Toro's Riemann problem

본 연구에서 개발한 모형은 댐 붕괴로 인한 불연속 흐름을 효율적으로 해석하기 위해 HLLC Riemann 솔버에 MUSCL 기법을 적용하였다. 이 솔버는 복잡한 비선형 문제를 단순화하기 위해 파동의 형태를 가정하지만, 가정이 적절하지 않을 경우 비물리적인 값이 계산된다. Toro의 Riemann problem은 Toro (2001)에 의해 제시된 검증 테스트로, 2차원 흐름 해석 모형에 적용된 근사 Riemann 솔버의 성능을 평가하기 위한 예제로 구성되어 있다. 본 연구에서는 솔버의 성능을 다양한 상황에서 검증하기 위해 Toro의 Riemann problem 중 “Left sonic rarefaction and right shock”와 “Dam break on dry bed”를 검증 예제로 활용하였다. 모든 테스트는 $x$ =0부터 $x$ =50 m까지 길이 50 m인 수평바닥의 직사각형 마찰이 없는 수로에서 수행되었으며, 초기조건은 Table 2와 같다. 여기서, $h_{L}$ 은 좌측 영역의 수심, $u_{L}$ 은 좌측 영역의 유속, $h_{R}$ 은 우측 영역의 수심, $u_{R}$ 은 우측 영역의 유속, $x_{0}$ 는 초기 불연속면이 놓여 있는 좌표이다.

Table 2.

Initial conditions for Toro‘s Riemann problems

Riemann Problem	$h_{L}$	$u_{L}$	$h_{R}$	$u_{R}$	$x_{0}$
Left sonic rarefaction and right shock	1.0	2.5	0.1	0.0	10.0
Dambreak on dry bed	1.0	0.0	0.0	0.0	20.0

개발된 모형의 결과를 Exact Riemann Solver의 해와 비교하여 제시하였다. 또한, 개발 모형에 적용된 2차 정확도 기법의 정확성을 정량적으로 검증하기 위하여 1차 정확도 기법을 적용한 결과와 비교하였으며, 이를 Fig. 2에 나타내었다. Left sonic rarefaction and right shock 문제의 경우, Table 1의 초기조건에 따라 좌측의 1 m 수체가 붕괴하여 우측으로 빠르게 전파된다. 붕괴 후 7 s 시점에서, 우측으로 빠르게 진행하는 흐름과 정지된 흐름이 만나 형성된 충격파(Shock wave)의 선단부(front)가 해석해와 잘 일치하였으며, 중간 영역(Intermediate region)에서의 전단파(Front wave)와 희박파(Rarefaction wave) 또한 잘 일치하였다. Dam break on dry bed 문제는 좌측의 1 m 수체가 우측의 완전 마름 상태의 수로로 붕괴하는 흐름을 모사하는 예제이다. Figs. 2(c) and 2(d)와 같이 본 연구에서 개발된 모형은 충격파가 마름 영역과 접할 때에도 수치적 불안정성이 발생하지 않고 수심과 단위폭당 유량을 안정적으로 계산하였다.

https://cdn.apub.kr/journalsite/sites/kwra/2025-058-11/N0200581113/images/kwra_58_11_13_F2.jpg

Fig. 2.

Numerical results for Toro’s Riemann problems: (a) depth - Left sonic rarefaction and right shock, (b) unit discharge (same case as (a)), (c) depth - dambreak on dry bed, (d) unit discharge (same case as (c))

1차 정확도의 HLLC와 2차 정확도의 MUSCL 결과를 비교하여 Table 3에 제시하였다. 이 표에서 $L_{1}$ 오차는 전체 구간에서의 평균적인 오차를 나타낸다. Table 3에서 알 수 있듯이 HLLC에 비해 MUSCL 기법의 $L_{1} (η)$ 및 $L_{1} (q)$ 가 현저히 감소하였다. Left sonic rarefaction and right shock 문제에서, $L_{1} (η)$ 는 0.008 m에서 0.004 m으로, $L_{1} (q)$ 는 0.022 m²/s에서 0.015 m²/s로 감소하였으며, dam break on dry bed에서는 $L_{1} (η)$ 가 0.019 m에서 0.002 m으로, $L_{1} (q)$ 가 0.045 m²/s에서 0.006 m²/s로 감소하였다. 이는 1차 정확도 기법이 흐름의 선단부에서 더 확산된 결과를 보이는 반면, 2차 정확도 기법인 MUSCL은 날카로운 선단부를 더 정확하게 포착함을 입증한다. 특히, Fig. 2(d)와 같이 마름젖음 선단부에서 유속이 급격히 증가하는 현상을 잘 재현하였다.

Table 3.

$L_{1}$ error comparison of HLLC and MUSCL methods (water depth and unit discharge)

Riemann Problem	$L_{1} (η) (m)$		$L_{1} (q) (m^{2} / s)$
Riemann Problem	HLLC	MUSCL	HLLC	MUSCL
Left sonic rarefaction and right shock	0.008	0.004	0.022	0.015
Dam break on dry bed	0.019	0.002	0.045	0.006

3.2 Thacker’s planar solution

개발된 모형을 Thacker (1981)의 2차원 곡면 위에서 수심이 시간에 따라 변하는 문제에 적용하였다. 이 문제는 수면이 Eq. (19)와 같은 2차원 곡면 영역에서 수면이 시간에 따라 주기적으로 진동하며 이동하는 현상을 모의한다. 이를 통해 매 순간 달라지는 마름젖음 영역의 변화를 수치해와 비교할 수 있어 대표적인 검증 예제로 활용된다.

(19)

z_{b} (x, y) = - h_{0} [1 - \frac{{(x - x_{0})}^{2} + {(y - y_{0})}^{2}}{a^{2}}]

여기서, ( $x_{0}$ , $y_{0}$ )는 격자의 중심; $h_{0}$ 는 중심에서의 수심; $a$ 는 중심에서 지표와 수면이 맞닿는 영역까지의 거리를 나타낸다. Thacker (1981)가 제시한 바와 같이, 바닥 마찰을 무시하면 수위와 속도의 해석해는 Eq. (20)을 이용하여 계산할 수 있다.

(20a)

η (x, y, z) = \frac{σ h_{0}}{a^{2}} (2 x \cos (ω t) + 2 y \sin (ω t) - σ)

(20b)

u (t) = - ω σ \sin (ω t), v (t) = ω σ \cos (ω t)

여기서 𝜎는 운동의 진폭을 결정하는 상수이고, 𝜔는 영역 중심을 기준으로 하는 회전 주파수를 나타낸다. Fig. 3(a)는 Eq. (20)를 3차원 바닥 지형으로 도시한 그림이며, Fig. 3(b)는 Eq. (20a)에서 $t$ =0에서의 초기 수심을 나타낸다.

https://cdn.apub.kr/journalsite/sites/kwra/2025-058-11/N0200581113/images/kwra_58_11_13_F3.jpg

Fig. 3.

Initial conditions for Thacker's planar rotation benchmark test: (a) bed elevation contours, (b) initial depth contours

Eq. (20)의 초기조건에 따라 수체는 $x$ 방향으로 이동하며, 중력과 바닥 형상에 의해 회전 운동을 한다. 마찰이 없으므로 이 운동은 감쇠하지 않고 무한히 반복되며, 이동 과정에서 지속적으로 마름젖음 전선(front)이 형성된다. 특히 마름젖음 경계면에서는 매우 얇은 박막 흐름이 나타나는데, 이때는 수심에 비해 유량이 상대적으로 커서 수치적 불안정성에 의해 비물리적인 속도가 발생할 수 있다. 이러한 비물리적인 속도는 원형 회전 수체의 형상을 왜곡시키고 정상적인 흐름을 방해한다. 따라서 견고한 수치모형이라면 박막 흐름에서 발생하는 수치적 불안정성을 제어하고, 반복적인 회전 운동이 진행되면서 원형 수체의 형상을 유지해야 한다.

본 적용사례에서 $T$ 는 회전 운동 후 수체가 원래 위치로 되돌아오는 주기이며, 수치모의는 3.5 $T$ 까지 수행하였다. 이 중 3 $T$ 에서 3.5 $T$ 구간의 결과를 해석해와 비교하였다. Fig. 4에서 점선은 Eq. (20a)에 따른 해석해의 수심을, 컨투어는 수심 모의 결과를 나타낸다. 이 그림에서 알 수 있듯이, 3.5 $T$ 의 회전 운동 후에도 수체의 형상이 해석해와 잘 일치하였다. Fig. 5에는 $y$ =2m단면에서 3 $T$ 부터 3.5 $T$ 동안의 단면 수위, $x$ 방향 유량, $y$ 방향 유량이 제시되어 있으며, 해석해에서 나타나는 진자 운동과 같은 반복적 거동과 수치해가 잘 일치하였다. 해석해와 수치 모의 결과 간의 오차를 정리한 Table 4에 따르면 그 오차는 10^-3~10^-4 m 수준이었으며, 이는 개발된 수치 모형이 얇은 박막 흐름을 잘 처리하며, 경사 소스 항과 플럭스 간의 균형이 잘 유지됨을 보여준다.

https://cdn.apub.kr/journalsite/sites/kwra/2025-058-11/N0200581113/images/kwra_58_11_13_F4.jpg

Fig. 4.

Depth comparison between analytical solution and numerical result from 3 $T$ to 3.5 $T$ in Thacker's planar rotation problem

https://cdn.apub.kr/journalsite/sites/kwra/2025-058-11/N0200581113/images/kwra_58_11_13_F5.jpg

Fig. 5.

Comparison of the water level, x-direction discharge, and y-direction discharge from 3 $T$ to 3.5 $T$ at the cross-section of $y$ =2 m

Table 4.

Quantitative verification results of Thacker’s planar rotation benchmark

t	$L_{1}$ error
t	𝜂	$q_{x}$	$q_{y}$
3.0T	2.9E-04	7.2E-04	3.7E-04
3.1T	6.0E-04	6.5E-04	3.0E-04
3.2T	6.3E-04	4.4E-04	4.8E-04
3.3T	6.5E-04	5.0E-04	6.7E-04
3.4T	7.4E-04	5.5E-04	1.2E-03
3.5T	4.0E-04	7.7E-04	5.0E-04

4. 도시 홍수 사례 적용

3장에서 제시된 두 가지 검증사례를 통해 개발된 모형이 충격파의 전파와 마름-젖음 경계의 천이 과정과 같은 복잡한 흐름 특성을 안정적으로 재현함을 확인하였다. 이에 따라 본 장에서는 실제 도시 홍수 사례를 대상으로, 복잡한 도시 지형에서 발생하는 구조물 주변부의 침수심 상승, 흐름 병목, 급속한 침수 확산 등 주요 도시 홍수 특성을 재현하여 개발된 모형의 성능을 평가하였다.

도시 홍수 적용 대상지로 호주 뉴캐슬의 Morgan-Selwyn를 선정하였다. 이 지역은 2007년 6월 Pasha Bulker 폭풍으로 큰 피해를 입었으며, Smith et al. (2012)은 뉴사우스웨일스대학교 수자원 연구센터에서 이 홍수를 재현하기 위한 축소 모형을 제작하였다(Figs. 6(a) and 6(b)). 이 모형은 축소된 디지털 지형 데이터를 기반으로 600 mm 간격의 지형 템플릿을 제작하고, 레이저 절단을 통해 수직 오차를 ±1 mm 이내로 유지하였다. 또한 주택 진입로를 측량해 약 20mm 높이로 추가하고, 항공 사진을 이용해 건물 형상을 축소하였다.

https://cdn.apub.kr/journalsite/sites/kwra/2025-058-11/N0200581113/images/kwra_58_11_13_F6.jpg

Fig. 6.

Overview of flash flood scaled model experiment: (a) experimental model photo, (b) detailed structure, (c) watershed location, (d) digital terrain data

본 연구에서는 Shin et al. (2019)의 FLUMEN 및 HDM-2D 모형 결과와 수리 실험값을 활용하여 개발된 모형의 도시 홍수의 적용 가능성을 평가하였다. HDM-2D는 Petrov-Galerkin 안정화 기법을 사용하여 천수방정식을 이산화하는 2차원 흐름 분석 모델이다. 이 모델은 채널 흐름과 도시 홍수를 시뮬레이션할 수 있는 범용 2차원 표면 흐름 모델로, 폭풍 해일, 강 홍수, 내륙 침수와 같은 여러 홍수 원인이 있는 복잡한 홍수 상황을 분석하고 수치적으로 예측하기 위한 수많은 연구에서 사용되고 검증되었다(Park et al., 2020; Eum et al., 2023). FLUMEN (FLUvial Modeling ENgine) 모델은 스위스의 Beffa에 의해 개발되었으며, 현재 스위스, 독일, 오스트리아와 같은 국가에서 홍수 분석 모델로 사용되고 있다(Fluvial.ch, 2010). 수치모의를 위해 구축한 지형은 총 20,791개의 셀과 11,083개의 노드로 구성하였으며, 이는 선행 연구에서 Shin et al. (2019)의 HDM-2D 모델과 FLUMEN에 적용된 격자 구성과 동일하다(Fig. 6(d)). 적용 사례는 상류에서 발생한 홍수파가 빠르게 도시 내부로 유입되어 단시간에 광범위한 침수를 야기하는 돌발홍수로 이 홍수파는 도심을 통과한 뒤 외부 경계를 따라 자유롭게 유출된다. 이는 Smith et al. (2012)의 수치실험과 동일하므로, 상류 경계에는 2007년 홍수 당시 Merewether Street에서 관측된 유출 유량 그래프에서 추출한 최대 유량 19.7 m³/s를 적용하였다. 하류 경계 조건과 주변 경계 조건은 자유 경계로 두어 홍수파가 경계면 외부로 빠져나가도록 설정하였다. 조도계수 또한 Smith et al. (2012)에서 사용한 값과 동일하게 설정하였으며, 초기 조건은 모든 절점에서 유속이 0인 마름 상태로 시작하였다.

수심과 유속 모의 결과를 Fig. 7에 제시하였다. Fig. 7(a) 등수심도에서 건물 전면부와 측면부의 수심이 증가하는 경향을 보이며, 모의 영역 전반에 걸쳐 침수심이 큰 차이를 보였다. 또한, Fig. 7(b)에서 건물의 차폐 효과로 인해 건물 전면 및 후면에서 흐름이 매우 느려지는 반면, 건물 사이의 협소한 통로에서는 병목 현상으로 인해 유속이 급격히 증가하는 특징을 확인할 수 있으며, 이는 실제 도시 홍수에서 흔히 나타나는 수심 집중 및 유속 가속 현상을 잘 재현한 결과라 판단된다.

https://cdn.apub.kr/journalsite/sites/kwra/2025-058-11/N0200581113/images/kwra_58_11_13_F7.jpg

Fig. 7.

Urban flood simulation results: (a) depth distribution, (b) velocity distribution

모의 결과의 정량적 검증은 Shin et al. (2019)이 제시한 FLUMEN 및 HDM-2D 결과와 수리 실험값을 활용하여 수행하였다. Shin et al. (2019)은 Smith et al. (2012)의 수리 실험 결과를 바탕으로 총 13개 지점의 수위와 속도를 분석하였으며, 이 지점은 저지대 또는 건물 사이의 흐름 집중 구간에 위치한다. 본 연구에서는 동일한 검증 지점에 대해 수리 실험값, HDM-2D, FLUMEN, 개발된 모형의 결과를 비교하여 Table 5와 Table 6에 제시하였다. 또한, 각 모형의 오차 지표를 요약하여 Table 7에 나타내었다. Table 7에서 $L_{1}$ , $L_{2}$ , $L_{\infty}$ 오차를 비교한 결과, 개발된 모형의 수심 오차는 FLUMEN ( $L_{1}$ : 0.113 m, $L_{2}$ : 0.125 m, $L_{\infty}$ : 0.182)보다 낮았으며, HDM-2D ( $L_{1}$ : 0.083 m, $L_{2}$ : 0.097 m, $L_{\infty}$ : 0.153)와 유사하거나 약간 높게 나타났다. 속도 오차의 경우 개발된 모형( $L_{1}$ : 0.112 m/s, $L_{2}$ : 0.138 m/s, $L_{\infty}$ : 0.295)은 FLUMEN (MAE: 0.278 m/s, RMSE: 0.309 m/s, $L_{\infty}$ : 0.498)보다 낮았고, HDM-2D (MAE: 0.110 m/s, RMSE: 0.134 m/s, $L_{\infty}$ : 0.256)와 비슷한 수준을 보였다. 이러한 결과는 개발된 모형의 정확도가 기존 선행 연구 모형과 유사하거나 일부 지표에서 우수함을 의미하며, Fig. 8과 같이 중앙의 대각선에 가까울수록 측정값과의 일치도가 높은 산점도에서 우수한 모의 성능을 확인할 수 있다.

Table 5.

Depth comparison in urban flood problem. (unit: m)

point	Measured	HDM-2D	FLUNEM	Present
1	0.721	0.660	0.658	0.661
2	0.442	0.412	0.369	0.427
3	0.100	0.100	0.101	0.111
4	0.406	0.348	0.342	0.345
5	0.442	0.379	0.359	0.383
6	0.784	0.631	0.615	0.624
7	0.523	0.466	0.384	0.403
8	0.442	0.363	0.337	0.311
9	0.766	0.629	0.614	0.632
10	0.712	0.567	0.533	0.560
11	0.604	0.467	0.422	0.606
12	0.397	0.267	0.247	0.257
13	0.388	0.360	0.280	0.347

Table 6.

Velocity comparison in urban flood problem. (unit: m/s)

point	Measured	HDM-2D	FLUNEM	Present
1	0.900	0.901	0.402	0.605
2	0.900	1.098	0.671	0.709
3	0.600	0.673	0.492	0.683
4	0.800	0.703	0.525	0.712
5	0.600	0.555	0.382	0.598
6	0.900	0.644	0.432	0.737
7	0.700	0.764	0.373	0.494
8	0.800	0.707	0.411	0.704
9	1.000	0.913	0.776	0.892
10	0.600	0.711	0.717	0.610
11	1.300	1.128	1.059	1.185
12	0.700	0.707	0.763	0.709
13	0.600	0.380	0.146	0.512

Table 7.

Summary of error metrics by models

Model Error	HDM-2D		FLUMEN		Present
Model Error	h (m)	V (m/s)	h (m)	V (m/s)	h (m)	V (m/s)
$L_{1}$	0.083	0.110	0.113	0.278	0.084	0.112
$L_{2}$	0.097	0.134	0.125	0.309	0.100	0.138
$L_{\infty}$	0.153	0.256	0.182	0.498	0.160	0.295

https://cdn.apub.kr/journalsite/sites/kwra/2025-058-11/N0200581113/images/kwra_58_11_13_F8.jpg

Fig. 8.

Scatter plot comparison of measured and simulated values: (a) depth comparison, (b) velocity comparison

이상의 결과를 종합하면, 개발된 모형은 홍수파가 건물과 같은 장애물에 도달하였을 때 나타나는 수위 상승 및 유속 변화 양상을 정확하게 포착하였다. 특히 건물 사이의 좁은 통로에서 병목 현상으로 인한 수위 상승과 유속 증가를 잘 모의하였다. 13개 측정 지점에서의 비교 결과, 기존 연구와 유사하거나 일부 구간에서는 더 우수한 성능을 보였다. 이러한 점을 고려할 때, 개발된 모형은 도시 홍수와 같은 돌발적 유동 상황에서도 높은 수치적 안정성과 신뢰할 만한 결과를 제공함을 확인할 수 있다.

5. 결론 및 향후 연구

본 연구에서는 복잡한 도시 환경에서 발생하는 돌발 홍수를 정확하고 안정적으로 모의하기 위해 엣지 기반 정수압 재구성 기법을 적용한 2차원 유한체적 흐름 모형을 개발하였다. 비구조 삼각 격자에서 MUSCL 경사 제한과 HLLC 근사 리만 솔버를 결합하고, 플럭스와 소스 항의 균형을 유지하는 well-balanced 기법을 도입함으로써 C-속성을 만족시켰다. 또한, 마찰 항을 음해적으로 처리하고, Taichi 기반 고속 병렬 연산 프레임워크를 활용하여 대규모 도시 홍수 시뮬레이션의 수치 계산 효율성을 크게 향상시켰다.

모형의 성능 검증 결과, Toro의 Riemann 문제와 Thacker의 곡면 회전 문제를 통해 충격파 전파와 마름-젖음 경계 처리 능력을 확인하였다. 또한 2007년 호주 뉴캐슬의 Pasha Bulker 폭풍 사례에 적용한 결과, 건물 전면부의 수위 상승과 협소 통로에서의 병목 현상 등 도시 홍수의 주요 특성을 정확하게 재현하였다. 정량적 비교에서 개발될 모형의 오차는 기존 모형(FLUMEN, HDM-2D)과 유사하거나 일부 지표에서 더 우수한 수준을 보였으며, 이는 본 모형이 복잡한 지형과 동적 흐름 조건에서도 신뢰성 있는 결과를 제공함을 입증한다. 그러나 현재 모형은 지표 흐름만을 해석하기 때문에 도시 홍수에서 중요한 지하 관망이나 인공 배수 시스템의 영향을 포괄적으로 반영하지 못하는 한계가 있다.

향후 연구에서는 본 모형이 지표 흐름을 안정적이고 효율적으로 해석하는 장점을 활용하여, 도시에서의 복잡한 흐름 전달(예: 관망 흐름, 펌프, 저류지, 배수 등)을 반영한 통합 흐름해석 모형으로 확장할 계획이다. 이를 통해 실제 강우 자료와 연계된 홍수 예측 시스템을 구축하고, 기후 변화 시나리오를 반영한 장기 홍수 위험 평가를 수행함으로써 모형의 실용성을 높이고자 한다. 또한 다양한 실제 침수 사례에 대한 추가 적용과 검증을 통해 모형의 범용성과 신뢰성을 지속적으로 강화할 예정이다.

Acknowledgements

본 연구는 환경부의 재원으로 한국환경산업기술원의 기후위기대응 홍수방어능력 기술개발사업의 지원을 받아 연구되었습니다(2480000599).

Conflicts of Interest

The authors declare no conflict of interest.

References

Audusse, E., Bouchut, F., Bristeau, M.-O., Klein, R., and Perthame, B. (2004). “A fast and stable well-balanced scheme with hydrostatic reconstruction for shallow water flows.” SIAM Journal on Scientific Computing, Vol. 25, No. 6, pp. 2050-2065. doi: 10.1137/S1064827503431090.

10.1137/S1064827503431090

Barth, T.J., and Jespersen, D.C. (1989). “The design and application of upwind schemes on unstructured meshes.” 27th AIAA Aerospace Sciences Meeting, Reno, NV, U.S. doi: 10.2514/6.1989-366.

10.2514/6.1989-366

Begnudelli, L., and Sanders, B.F. (2006). “Unstructured grid finite-volume algorithm for shallow-water flow and scalar transport with wetting and drying.” Journal of Hydraulic Engineering, Vol.132, No. 4, pp. 371-384. doi: 10.1061/(ASCE)0733-9429(2006)132:4(371).

10.1061/(ASCE)0733-9429(2006)132:4(371)

Bermudez, A., and Vazquez, M.E. (1994). “Upwind methods for hyperbolic conservation laws with source terms.” Computers & Fluids, Vol. 23, No. 8, pp. 1049-1071. doi: 10.1016/0045-7930(94)90004-3.

10.1016/0045-7930(94)90004-3

Centre for Research on the Epidemiology of Disasters (CRED) (2025). 2024 Disasters in numbers: A hot and stormy year. Brussels, Belgium.

El Kadi Abderrezzak, K., Paquier, A., and Mignot, E. (2009). “Modelling flash flood propagation in urban areas using a two-dimensional numerical model.” Natural Hazards, Vol. 50, pp. 433-460. doi: 10.1007/s11069-008-9300-0.

10.1007/s11069-008-9300-0

Eum, T., Seo, I. W., Shin, E., and Song, C.G. (2023). “Development and application of a user-friendly general-purpose predictive simulation tool for two-dimensional flow analysis.” Environmental Modelling & Software, Vol.163, 105665. doi: 10.1016/j.envsoft.2023.105665.

10.1016/j.envsoft.2023.105665

Fluvial.ch (2010). FLUMEN user manual (Version 07-2010), accessed 10 November, 2025, <https://www.fluvial.ch/m/FlumenMan.pdf.>.

Hou, J., Liang, Q., Simons, F., and Hinkelmann, R. (2013). “A 2D well-balanced shallow flow model for unstructured grids with novel slope source term treatment.” Advances in Water Resources, Vol. 52, pp. 107-131. doi: 10.1016/j.advwatres.2012.08.003.

10.1016/j.advwatres.2012.08.003

Hu, Y., Li, T., Anderson, L., Ragan-Kelley, J., and Durand, F. (2019). “Taichi: A language for high-performance computation on spatially sparse data structures. ACM Transactions on Graphics, Vol.38, No. 6, 201. doi: 10.1145/3355089.3356506.

10.1145/3355089.3356506

Hubbard, M.E. (1999). “Multidimensional slope limiters for MUSCL-type finite volume schemes on unstructured grids.” Journal of Computational Physics, Vol. 155, No. 1, pp. 54-74. doi: 10.1006/jcph.1999.6329.

10.1006/jcph.1999.6329

LeVeque, R.J. (1992). Numerical methods for conservation laws. Birkhäuser Basel, Basel, Switzerland. doi: 10.1007/978-3-0348-8629-1.

10.1007/978-3-0348-8629-1

Liang, Q., and Marche, F. (2009). “Numerical resolution of well-balanced shallow water equations with complex source terms.” Advances in Water Resources, Vol. 32, No. 5, pp. 873-884. doi: 10.1016/j.advwatres.2009.02.010

10.1016/j.advwatres.2009.02.010

Liang, Q., Borthwick, A.G.L. (2009). “Adaptive quadtree simulation of shallow ﬂows with wet-dry fronts over complex topography.” Computers & Fluids , Vol. 38, No. 2, pp. 221-234.

10.1016/j.compfluid.2008.02.008

Ma, D.J., Sun, D.J., and Yin, X.Y. (2007). “Solution of the 2-D shallow water equation with source terms in surface elevation splitting form.” International Journal for Numerical Methods in Fluids, Vol. 55, No. 5, pp. 431-454.

10.1002/fld.1462

Marche, F., Bonneton, P., Fabrie, P., and Seguin, N. (2007). “Evaluation of well-balanced bore-capturing schemes for 2D wetting and drying processes.” International Journal for Numerical Methods in Fluids, Vol. 53, No. 5, pp. 867-894. doi: 10.1002/fld.1311.

10.1002/fld.1311

Park, I., Seo, I. W., Shin, J., and Song, C.G. (2020). “Experimental and numerical investigations of spatially-varying dispersion tensors based on vertical velocity profile and depth-averaged flow field.” Advances in Water Resources, Vol. 142, 103606. doi: 10.1016/j.advwatres.2020.103606.

10.1016/j.advwatres.2020.103606

Shin, E., Shin, J., Rhee, D., Kim, H.-J., and Song, C. (2019). “Integrated inundation modeling of flooded water in coastal cities.” Applied Sciences, Vol. 9, No. 7, 1313. doi: 10.3390/app9071313.

10.3390/app9071313

Smith, G.P., Wasko, C.D., and Miller, B.M. (2012). “Modelling the influence of buildings on flood flow.” In Proceedings of the 52nd Floodplain Management Australia National Conference, Batemans Bay, NSW, Australia.

Sod, G.A. (1985). Numerical methods in fluid dynamics: Initial and initial boundary-value problems. Cambridge University Press, Cambridge, UK.

10.1017/CBO9780511753138

Song, L., Zhou, J., Li, Q., Yang, X., and Zhang, Y. (2011). “An unstructured finite volume model for dam-break floods with wet/dry fronts over complex topography.” International Journal for Numerical Methods in Fluids, Vol. 67, No. 8, pp. 960-980. doi: 10.1002/fld.2397.

10.1002/fld.2397

Thacker, W.C. (1981). “Some exact solutions to the nonlinear shallow-water wave equations.” Journal of Fluid Mechanics, Vol. 107, pp. 499-508. doi: 10.1017/S0022112081001882.

10.1017/S0022112081001882

Toro, E.F. (2001). Shock-capturing methods for free-surface shallow flows. John Wiley & Sons, Hoboken, NJ, U.S.

Toro, E.F. (2009). Riemann solvers and numerical methods for fluid dynamics: A practical introduction (3rd ed.). Springer, Berlin, Heidelberg, Germany. doi: 10.1007/b79761.

10.1007/b79761

Toro, E.F., and Garcia-Navarro, P. (2007). “Godunov-type methods for free-surface shallow flows: A review.” Journal of Hydraulic Research, Vol. 45, No. 6, pp. 736-751. doi: 10.1080/00221686.2007.9521812.

10.1080/00221686.2007.9521812

Toro, E.F., Spruce, M., and Speares, W. (1994). “Restoration of the contact surface in the HLL-Riemann solver.” Shock Waves, Vol. 4, No. 1, pp. 25-34. doi: 10.1007/BF01414629.

10.1007/BF01414629

van der Walt, S., Colbert, S.C., and Varoquaux, G. (2011). “The NumPy array: A structure for efficient numerical computation.” Computing in Science & Engineering, Vol. 13, No. 2, pp. 22-30. doi: 10.1109/MCSE.2011.37.

10.1109/MCSE.2011.37

van Leer, B. (1973). “Towards the ultimate conservative difference scheme I. The quest of monotonicity.” In Lecture notes in physics (Vol. 18), Edited by Cabannes, H., and Temam, R., Springer, Berlin, Heidelberg, Germany, pp. 163-168. doi: 10.1007/BFb0118673.

10.1007/BFb0118673

van Leer, B. (1979). “Towards the ultimate conservative difference scheme. V. A secondorder sequel to Godunov’s method.” Journal of Computational Physics, Vol. 32, No. 1, pp. 101-136. doi: 10.1016/0021-9991(79)90145-1.

10.1016/0021-9991(79)90145-1

Venkatakrishnan, V. (1995). “Convergence to steady state solutions of the Euler equations on unstructured grids with limiters.” Journal of Computational Physics, Vol. 118, No. 1, pp. 120-130. doi: 10.1006/jcph.1995.1084.

10.1006/jcph.1995.1084

Wang, Y., Liang, Q., Kesserwani, G., and Hall, J.W. (2011). “A 2D shallow flow model for practical dam-break simulations.” Journal of Hydraulic Research, Vol. 49, No. 3, pp. 307-316. doi: 10.1080/00221686.2011.566248.

10.1080/00221686.2011.566248

Journal of Korea Water Resources Association ISSN:2799-8746(Print) 2799-8754(Online) 한국수자원학회 논문집

Preview

Development of two-dimensional finite volume flow model with edge-based hydrostatic reconstruction and its application to urban flood

ABSTRACT

MAIN

(1a)

(1b)

(1c)

(2)

(3)

(4)

(5)

(6)

(7)

(8a)

(8b)

Fig. 1.

Example of cell reconstruction in triangular grid system: (a) Triangular grid system overview, (b) 1st order accuracy physical quantity representation, (c) 2nd order accuracy linear reconstruction

(9)

(10)

(11)

(12a)

(12b)

(13)

(14)

(15)

(16)

(17a)

(17b)

(18a)

(18b)

Table 1.

Comparison of computation time by grid cell count (Single CPU, Multi CPU, GPU)

Table 2.

Initial conditions for Toro‘s Riemann problems

Fig. 2.

Numerical results for Toro’s Riemann problems: (a) depth - Left sonic rarefaction and right shock, (b) unit discharge (same case as (a)), (c) depth - dambreak on dry bed, (d) unit discharge (same case as (c))

Table 3.

L1 error comparison of HLLC and MUSCL methods (water depth and unit discharge)

(19)

(20a)

(20b)

Fig. 3.

Initial conditions for Thacker's planar rotation benchmark test: (a) bed elevation contours, (b) initial depth contours

Fig. 4.

Depth comparison between analytical solution and numerical result from 3T to 3.5T in Thacker's planar rotation problem

Fig. 5.

Comparison of the water level, x-direction discharge, and y-direction discharge from 3T to 3.5T at the cross-section of y=2 m

Table 4.

Quantitative verification results of Thacker’s planar rotation benchmark

Fig. 6.

Overview of flash flood scaled model experiment: (a) experimental model photo, (b) detailed structure, (c) watershed location, (d) digital terrain data

Fig. 7.

Urban flood simulation results: (a) depth distribution, (b) velocity distribution

Table 5.

Depth comparison in urban flood problem. (unit: m)

Table 6.

Velocity comparison in urban flood problem. (unit: m/s)

Table 7.

Summary of error metrics by models

Fig. 8.

Scatter plot comparison of measured and simulated values: (a) depth comparison, (b) velocity comparison

Acknowledgements

Conflicts of Interest

References

$L_{1}$ error comparison of HLLC and MUSCL methods (water depth and unit discharge)

Depth comparison between analytical solution and numerical result from 3 $T$ to 3.5 $T$ in Thacker's planar rotation problem

Comparison of the water level, x-direction discharge, and y-direction discharge from 3 $T$ to 3.5 $T$ at the cross-section of $y$ =2 m