Abstract

For communicating urban flood risk to authorities and the public, a realistic three-dimensional visual display is frequently more suitable than detailed flood maps. Virtual reality could also serve to plan short-term flooding interventions. We introduce here an alternative approach for simulating three-dimensional flooding dynamics in large- and small-scale urban scenes by reaching out to computer graphics. This approach, denoted ‘particle in cell’, is a particle-based CFD method that is used to predict physically plausible results instead of accurate flow dynamics. We exemplify the approach for the real flooding event in July 2016 in Innsbruck.

INTRODUCTION

Flooding in urban areas is recognized as a significant risk to public infrastructure, welfare and even as a danger to life (Jha et al. 2012). The reason for the increased number of incidents and detrimental effects has been discussed at length, see e.g. Ashley et al. 2005. For the planning of mitigation measures (both long-term infrastructure adaption and short-term disaster prevention) tools are needed to forecast the magnitude of the flooding event. Over the last few years, flood risk modelling has progressed significantly towards two-dimensional shallow wave simulation of surface flows (Hunter et al. 2008). However, for communicating flood risks to authorities and the public, map-based results are not always suitable and there is a need for more visual display of the action. But an accurate three-dimensional computational fluid dynamics (CFD) simulation of large-scale urban flooding dynamics is out of question – even on modern supercomputers the calculation time would be enormous and thus prohibitive for any practical application.

In this paper, we develop an alternative approach for introducing virtual reality in flood-risk communication by reaching out to computer graphics (CG). We apply particle-based CFD, more precisely the particle in cell (PIC) method, for simulating urban flooding dynamics in a three-dimensional mesh representation of a city. The key here is to balance uncertainty in the model predictions against visual display. The PIC method is based on physical laws for flow calculation but predicts plausible results instead of accurate flow dynamics. On the other hand, PIC allows for numerically stable CFD simulations of whole cities in real time. In the following we introduce the numerical background of the approach and show two example applications, with one addressing the real flooding event in July 2016 in Innsbruck.

METHODS

Computational fluid dynamics

Simulating transient fluid flow generally requires solving the underlying physical model, which is expressed by means of the Navier-Stokes equations (NSE). These equations for incompressible flow consist of a momentum term expressed as  
formula
and the continuity equation , where t denotes time, the velocity vector, the external force (gravity), p the pressure, the density and the kinematic viscosity of the fluid. The momentum equation governs the effect of the individual forces on the fluid advection while the divergence free condition of the continuity equation constrains the incompressibility. Since an analytical solution of the NSE exists only for selected flows, solving the equations for arbitrary problems requires a discretization to a numerical method. The discretization to a finite number of interpolation points generally follows either the Eulerian approach, where points are fixed in space, or the Lagrangian approach, where the interpolation points follow the deformation of the continuum (Bridson 2015).

Lagrangian methods, where smoothed particle hydrodynamics (SPH) is the most prominent mesh-free representative, are advantageous for the simulation of distorted and rapidly moving free surface flows, multi-phase or multi-physics problems and fluid structure interactions (Violeau & Rogers 2016). Despite some inherent drawbacks due to its Lagrangian nature, SPH is validated against a vast number of acknowledged test cases and has been applied to a variety of fields. But a severe limitation for the application of SPH in real-time computer simulations is the high computational effort. For flows and scales relevant in environmental engineering applications, the time step is approximately in the order seconds (Winkler et al. 2017), where in every step each particle has to compute the interaction with several hundred neighbouring particles (Crespo et al. 2015). The simulation of 3 seconds of a relatively small validation test case containing particles requires the execution of highly optimized DualSPHysics SPH code (Crespo et al. 2015) for 11.7 hours on an Intel Xeon E5-1620 3.50 GHz quad core CPU – see the example in Figure 1.

Figure 1

Visual comparison of SPH (left) and PIC/FLIP (right) flow simulations – the impact of water on a rigid object is rendered as a translucent fluid using Mitsuba renderer (Jakob 2010). While the flow dynamics are comparable, the detailed effect of the rigid obstacle is different. By visual comparison to the video and validation of the test case, the flow profile of the PIC/FLIP simulation is even more realistic (Kleefsman et al. 2005).

Figure 1

Visual comparison of SPH (left) and PIC/FLIP (right) flow simulations – the impact of water on a rigid object is rendered as a translucent fluid using Mitsuba renderer (Jakob 2010). While the flow dynamics are comparable, the detailed effect of the rigid obstacle is different. By visual comparison to the video and validation of the test case, the flow profile of the PIC/FLIP simulation is even more realistic (Kleefsman et al. 2005).

To overcome this issue of excessive computational effort, fluid simulation in CG often applies hybrid methods, most notably the PIC method. PIC is a combination of Lagrangian particles with a Eulerian mesh, and was developed in the 1950s for solving hydrodynamic problems. Figure 2 outlines the basic steps in a PIC iteration, where first the quantities of the fluid particles are transferred to the grid, then the NSE are solved on the grid and last the updated velocities are transferred back to the particles. The back and forth interpolation of quantities from particles to the grid creates considerable numerical dissipation, mainly caused by the averaging process. Brackbill & Ruppel (1986) improved the method to reduce this dissipation to its derivative fluid implicit particles (FLIP), which is achieved by transferring back the velocity change instead of the cell velocity to update the particle velocities.

Figure 2

Principles of the PIC/FLIP method, visualizing particles as circles, arrows as velocities, white background as fluid cells and hatched background as boundary cells. (a) shows the current state of the system, (b) interpolates the particle velocities to the grid , (c) adjusts the previous cell velocity (light grey) by adding the velocity change (red) after solving the NSE to , (d) updates the particle velocities by transferring , (e) advects the particle positions according to , (f) shows the state of the system after the PIC update. Please refer to the online version of this paper to see this figure in colour: http://dx.doi.org/10.2166/wst.2017.567.

Figure 2

Principles of the PIC/FLIP method, visualizing particles as circles, arrows as velocities, white background as fluid cells and hatched background as boundary cells. (a) shows the current state of the system, (b) interpolates the particle velocities to the grid , (c) adjusts the previous cell velocity (light grey) by adding the velocity change (red) after solving the NSE to , (d) updates the particle velocities by transferring , (e) advects the particle positions according to , (f) shows the state of the system after the PIC update. Please refer to the online version of this paper to see this figure in colour: http://dx.doi.org/10.2166/wst.2017.567.

Solving the NSE is done with the grid-based finite differences method. The fluid solver is implemented on a marker and cell (MAC) grid, which is a staggered grid that stores different fluid quantities at different grid positions. Thus the pressure is stored at the cell centre while the normal components of the velocity are stored at cell faces. For example, at the face between two adjacent cells in x direction only the x velocity component of the fluid is stored. This method has the advantage that central differences can be solved accurately for the velocity divergence and the pressure gradient (Bridson 2015). Using the Helmholtz-Hodge decomposition theorem allows the NSE to be rewritten in a form that relates the divergent velocity field of the advected particles to the pressure gradient in the form . This equation is called the pressure poison equation, which is a system of linear equations that we solve using a simple Jacobi iteration. There exist more sophisticated algorithms to solve this equation but due to its simplicity it is a reasonable choice for a parallel (computing) implementation in a graphics engine (Pharr & Fernando 2005). Subtracting the final pressure gradient from the divergent velocity field results in the divergence-free velocity that is used to update the particle velocities.

In comparison to SPH the numerical effort is drastically reduced by replacing the computational stencil from hundreds of particles to only the adjacent cells used in central differences. Since the method is very stable (Zhu & Bridson 2005) a significantly higher time step can be chosen than for SPH, such that the simulation of the above-mentioned validation case is executed in 5.5 min on the same system. Comparing the execution times shows that the PIC/FLIP method is more than two orders of magnitude faster than SPH. Um et al. (2016) examined the dambreak setup of Kleefsman et al. (2005) and experienced a speedup of 400 when comparing FLIP to SPH. They also showed that the surface heights measured at various positions match the SPH simulation and the experiment well. Optimizing the code further for graphics processing units (GPUs) allows almost real-time performance for PIC/FLIP execution (NVIDIA Corporation 2008).

Boundary conditions

Unreal engine 4 (UE4) (Epic 2017) provides a suite of development tools for the creation of real-time graphics software such as video games or physics-based simulations. State-of-the-art implementations for photorealistic rendering, particle systems, terrain modelling and virtual (VR) and augmented reality (AR) are provided by the framework. The software is published as open source software (OSS) and the licensing model is solely based on a 5% royalty of the worldwide gross revenue. Thus, UE4 was chosen as CG framework for the development of environmental VR applications, as it can be extended to user requirements without any restriction and is of zero cost for non-commercial products.

At their core, real-time graphics engines are highly optimized software implementations that interpret the geometric representation of an environment in the form of polygon meshes. These meshes model objects in terms of vertices, edges and faces, which allows the engine to render 3D objects to 2D images in screen space by projecting the scene depending on the desired perspective. By adding physical quantities to the objects, engines such as UE4 that contain a physics solver are able to simulate and visualize physical phenomena. We are building on top of these simulation capabilities using the cataclysm research project, which provides a highly optimized PIC/FLIP fluid solver that is integrated in the physics and rendering pipeline of UE4 (NVIDIA 2015).

Creating scenes for environmental fluid simulations requires the incorporation of the visual representations as well as the definition of the physical boundary conditions. Both terms are reflected by the mesh representation of the scene objects, which is constrained by the required level of detail and the available data. For example, in an urban context the correct shape of the buildings is required to predict an accurate effect on the fluid flow, while detailed modelling of the windows or roofs is negligible. The missing visual detail is effectively compensated by properly texturing the object representation with images.

Computational and memory-efficient boundary conditions are defined to prevent particles from penetrating rigid obstacles, which is achieved by introducing an artificial pressure constraint at the object boundaries (Guay 2014). For efficiency reasons the voxelization of object meshes to boundary cells in the solver grid is avoided by computing signed distance field (SDF) volumes around each rigid object. This functionality was originally implemented in UE4 to improve shading with ambient occlusion or ray tracing, but can be utilized for the definition of fluid boundaries. The SDF is based on the distance function defined as  
formula
that returns the closest distance to a closed set of points S for any point , which is extended to the signed distance function as  
formula

The discretization of to SDFs, in graphics applications often called level set (LS), provides a distance function that is smooth everywhere except on the medial axis, consisting of points that are equidistant to different points on the surface (Bridson 2015). Since this axis is by definition disjointed from the surface, the gradient is defined near the surface and represents the outward-pointing normal. These properties are used in the next step to define a pressure gradient for the velocity update computation of the fluid solver (Guay 2014).

Based on this implicit representation of rigid obstacles it is sufficient to compute only the fluid domain for solving the NSE on the grid, without needing to explicitly discretise boundaries (visualized by the cubes in Figure 4 right). Similar to the SDF, the particle density is projected on a sparse grid based on a sharp kernel definition for each particle (Ando et al. 2012). The zero iso-contour of the resulting particle density field describes the fluid surface, which is provided to the solver to define the computational domain.

Visualization

Apart from the efficient computation of particle advection, a fast visualization technique is necessary to attain real-time fluid simulation performance. A straight forward visualization technique for particle fluid simulations is the representation of particles as solid spheres. The most plausible choice of particle volume arises from the initial discretization relationship of fluid volume to the number of particles. This option is a valid choice for the analysis of the fluid flow but not sufficient for realistic visualization of a translucent continuum. Applying partially transparent materials to the particle spheres adds transparency to the rendered fluid, but the high number of meshes cannot be rendered efficiently and the shading of the individual spheres creates artefacts.

To overcome this issue, point clouds have traditionally been triangulated to polygon meshes, e.g. using the marching cubes (MC) algorithm (Lorensen & Cline 1987). Equal to the fluid domain definition above, MC transfers the fluid density to the nodes of a uniform grid. While for the definition of the computational domain it is sufficient to determine the grid cells that have a fluid density greater than 0, MC requires an iso-value in addition to the scalar field. The iso-surface is extracted from the scalar field for each individual cube by determining edges that connect vertices higher and lower than the iso-value. The algorithm is efficient in creating a surface that passes through the interpolated positions of the iso-value on those edges. This algorithm is applied to fluid simulations in Figure 1 but as it handles all possible intersections of the iso-surface it is computationally too expensive for the execution in a real-time visualization environment.

For real-time computation of rendered surfaces, an intermediate approach is implemented that uses an LS uniform grid with a resolution twice as large as the solver grid. The increased resolution allows for a times higher level of detail in the surface generation for the visualization. Extraction of the zero iso-contour is interpolated on the LS grid between vertices similar to the MC algorithm. Since the fluid rendering only requires this from the rendering perspective, the surface extraction is performed after the projection onto 2D screen space (Pharr & Fernando 2005). As discussed in the previous section, the LS of the fluid surface furthermore defines the surface gradient, which is necessary for realistic shading of the fluid simulation. Shading defines the visual properties for physically based rendering such as colour, reflection and refraction.

RESULTS AND DISCUSSION

The strength of the approach presented herein is the rapid integration of various geometric data sources for the creation of physically plausible drainage simulations. Automatic reconstruction of urban scenes receives a lot of attention and has made much progress due to the wide availability of imagery from the air (Musialski et al. 2013). Recent advances in 3D indoor scene reconstruction enables real-time mixed reality applications (Dai et al. 2016). The combination of such approaches with the UE4 graphics engine allows for real-time VR and AR simulations. Integration of the fluid solver discussed provides a framework to augment such visualizations with the simulation of fluid flow for millions of particles.

The users of such an approach would be, among others, offices for disaster control management, emergency control and education, public utilities and politics. They could employ virtual reality for exemplifying the magnitude of the event as well as mitigating measures. Especially for communicating to non-professional audiences, a realistic display of the scene is easier to comprehend than the correct interpretation of conventional tools such as flood maps. The interactive capabilities of the software enable users to improve their understanding of fluid flow in general and allows them to develop skills with respect to the assessment and management of flooding events.

Case study 1: village flooding

In line with the intended purpose of this approach, we simulate an extreme rain event (return period 500 years) that happened in the area of Innsbruck in July 2016 (Kleidorfer et al. 2017). The short-term intensity of the event caused flooding throughout the region with high media coverage featuring videos from an upland village (see Figure 3 left). In the simulation the prototyped village is flooded not by rain runoff but by artificially creating surface water runoff uphill of the depicted scene. The GPU-based implementation of the fluid solver allows the simulation to run interactively with one million particles on a single graphics card, including scene rendering.

Figure 3

Comparison of the actual footage (Life Radio Tirol 2016) and the simulation of an extreme rain event in a village near the city of Innsbruck. The perspective has been altered for better visualization of the scene. DTM for the simulation provided by (c) Land Tirol.

Figure 3

Comparison of the actual footage (Life Radio Tirol 2016) and the simulation of an extreme rain event in a village near the city of Innsbruck. The perspective has been altered for better visualization of the scene. DTM for the simulation provided by (c) Land Tirol.

For the simulation of surface flows, we consider a triangle mesh, based on a digital terrain model (DTM) for this area. The DTM has a resolution of 1 × 1 m and is imported to UE4 after conversion to the supported 16-bit heightmap format. For visualization the landscape is textured with an orthophoto of 0.2 × 0.2 m resolution. This information is augmented with geospatial vector data, which is used to define the volumetric obstacles such as buildings in the region of interest. While there exists no technical challenge to create the 3D mesh representation from the shape data that are available through OpenStreetMap, we rely on the modelling capabilities of ESRI CityEngine (CE). CE enables the automatic retrieval of the required data that serve as the basis for the procedural generation of an urban environment, so that a realistic-looking prototype based on publicly available data can be created within minutes (see Figure 3 right).

In Figure 3 a snapshot of the simulated scene is displayed next to a video frame that was captured during the event. Although the DTM resolution is not capable of capturing fine details such as sidewalks or drainage, the simulation is able to predict the intersection of the two flows at the street junction in a physically plausible way, thus demonstrating the validity of the approach. The movement of the car is captured by the fluid dynamics and diffuse particles have been added for better visualization of the flow. Since the hedge is represented in the DTM as elevation the simulation predicts the flow of the fluid around it.

Case study 2: underground inundation

In general, pluvial flooding of public areas constitutes a rare event due to urban drainage measures. Nonetheless the infrastructure is designed with a certain failure return period such that in case of stronger events, underground inundation cannot be prevented. As such incidents pose a severe threat to public safety, Ishigaki et al. (2008) modelled the effect of underground inundation on evacuation planning using a life-size model of a stair and corridor.

We demonstrate a complementary approach to visualize such an event by applying the solver to a small-scale urban scene of a subway station. The geometric model is provided by UE4 as a showcase example for realistic lighting and rendering capabilities. Adding a fluid source upstairs of the environment shows the capabilities of interactive fluid computation and rendering (see Figure 4). Depending on the simulated inflow discharge, the fluid attains a plausible water depth. Visualization of fluid quantity in the context of such events is a useful tool for the training of rescue staff on the effects of inundation.

Figure 4

Subway scene showing the combination of detailed CG with real-time fluid simulation (left). The right image illustrates the numerical foundation of simulation, visualizing the individual particles. The cubes show the fluid domain that is calculated for every frame.

Figure 4

Subway scene showing the combination of detailed CG with real-time fluid simulation (left). The right image illustrates the numerical foundation of simulation, visualizing the individual particles. The cubes show the fluid domain that is calculated for every frame.

CONCLUSIONS

For the purpose of communicating urban flood risk to authorities and the public, virtual reality is a novel but sound approach. This could also serve to plan short-term disaster prevention measures or rescue interventions. The discussed particle-based CFD method allows for extremely efficient and fast calculation of 3D flow dynamics in large urban scenes. While not aiming for an accurate computation of all features of the fluid dynamics, the method simulates physically plausible results within an acceptable uncertainty range. The results for a highly dynamic test case show accurate prediction of the flow, such that the accuracy of applications in urban drainage is mainly dominated by the quality of the boundary conditions regarding inflow and outflow. The validity of the approach has been exemplified for a real flooding case in an urban settlement.

ACKNOWLEDGEMENT

This research is part of projects that are funded by the Austrian Science Fund (FWF): P26768-N28 and the Austrian Research Promotion Agency (FFG): 850738. We gratefully acknowledge the support of NVIDIA Corporation with the donation of the GTX Titan X GPU used for this research. We thank Doyub Kim for his assistance in creating the renderings for Figure 1 and Land Tirol for providing the elevation data and imagery for the simulation of the rain event.

REFERENCES

REFERENCES
Ando
,
R.
,
Thurey
,
N.
&
Tsuruno
,
R.
2012
Preserving fluid sheets with adaptively sampled anisotropic particles
.
IEEE Transactions on Visualization and Computer Graphics
18
(
8
),
1202
1214
.
Ashley
,
R. M.
,
Balmforth
,
D. J.
,
Saul
,
A. J.
&
Blanskby
,
J. D.
2005
Flooding in the future–predicting climate change, risks and responses in urban areas
.
Water Science and Technology
52
(
5
),
265
273
.
Brackbill
,
J. U.
&
Ruppel
,
H. M.
1986
FLIP: A method for adaptively zoned, particle-in-cell calculations of fluid flows in two dimensions
.
Journal of Computational Physics
65
(
2
),
314
343
.
Bridson
,
R.
2015
Fluid Simulation for Computer Graphics
.
CRC Press
,
Boca Raton, FL, USA
.
Crespo
,
A. J. C.
,
Domínguez
,
J. M.
,
Rogers
,
B. D.
,
Gómez-Gesteira
,
M.
,
Longshaw
,
S.
,
Canelas
,
R.
,
Vacondio
,
R.
,
Barreiro
,
A.
&
García-Feal
,
O.
2015
DualSPHysics: open-source parallel CFD solver based on smoothed particle hydrodynamics (SPH)
.
Computer Physics Communications
187
,
204
216
. .
Dai
,
A.
,
Nießner
,
M.
,
Zollhöfer
,
M.
,
Izadi
,
S.
&
Theobalt
,
C.
2017
BundleFusion: Real-Time Globally Consistent 3D Reconstruction Using On-the-Fly Surface Reintegration
.
ACM Trans. Graph. 36, 3, Article 24 (May 2017), 18 pages. DOI: https://doi.org/10.1145/3054739
.
Epic
,
G.
2017
Unreal Engine 4
. .
Guay
,
M.
2014
Simple rasterization-based liquids
. In:
GPU Pro 5: Advanced Rendering Techniques
(
Engel
,
W.
, ed.).
CRC Press
,
Boca Raton, FL, USA
, pp.
55
64
. .
Hunter
,
N. M.
,
Bates
,
P. D.
,
Neelz
,
S.
,
Pender
,
G.
,
Villanueva
,
I.
,
Wright
,
N. G.
,
Liang
,
D.
,
Falconer
,
R. A.
,
Lin
,
B.
&
Waller
,
S.
, &
others
2008
Benchmarking 2D hydraulic models for urban flood simulations
. In:
Proceedings of the Institution of Civil Engineers: Water Management
,
ICE Publishing, Scotland, UK
, pp.
13
30
.
Ishigaki
,
T.
,
Onishi
,
Y.
,
Asai
,
Y.
,
Toda
,
K.
&
Shimada
,
H.
2008
Evacuation criteria during urban flooding in underground space
. In:
Conference Proceedings 11th International Conference on Urban Drainage
,
Edinburgh International Conference Centre, Scotland, UK. 11 ICUD, 31st August–5th September 2008. International Association for Hydro-Environment Engineering and Research (IAHR), Madrid
.
Jakob
,
W.
2010
Mitsuba renderer
.
Jha
,
A. K.
,
Bloch
,
R.
&
Lamond
,
J.
2012
Cities and Flooding: a Guide to Integrated Urban Flood Risk Management for the 21st Century
.
World Bank Publications
,
Washington, DC, USA
.
Kleefsman
,
K. M. T.
,
Fekken
,
G.
,
Veldman
,
A. E. P.
,
Iwanowski
,
B.
&
Buchner
,
B.
2005
A volume-of-fluid based simulation method for wave impact problems
.
Journal of Computational Physics
206
(
1
),
363
393
.
Kleidorfer
,
M.
,
Tscheikner-Gratl
,
F.
,
Vonach
,
T.
&
Rauch
,
W.
2017
What can we learn from a 500-year event? Experiences from urban drainage in Austria
. In:
Proceedings of 14th International Conference on Urban Drainage (ICUD 2017)
,
Prague
.
Life Radio Tirol
2016
Schwere Unwetter in Tirol: Wenn innerhalb von 5 Minuten mitten in Aldrans die Welt untergeht
.
(accessed 20 June 2017)
.
Lorensen
,
W. E.
&
Cline
,
H. E.
1987
Marching cubes: a high resolution 3D surface construction algorithm
. In:
ACM Siggraph Computer Graphics
.
ACM, New York, NY, USA
, pp.
163
169
.
Musialski
,
P.
,
Wonka
,
P.
,
Aliaga
,
D. G.
,
Wimmer
,
M.
,
Gool
,
L. v.
&
Purgathofer
,
W.
2013
A survey of urban reconstruction
.
Computer Graphics Forum
32
(
6
),
146
177
.
NVIDIA Corporation
2008
Programming guide
.
Pharr
,
M.
&
Fernando
,
R.
2005
Gpu gems 2: Programming Techniques for High-Performance Graphics and General-Purpose Computation
.
Addison-Wesley Professional
,
Boston, MA, USA
.
Um
,
K.
,
Hu
,
X.
&
Thuerey
,
N.
2016
Analysis of free surface simulation using breaking-dam benchmark
. In:
Proceedings of the 11th International SPHERIC Workshop
,
Munich
, pp.
410
415
.
Violeau
,
D.
&
Rogers
,
B. D.
2016
Smoothed particle hydrodynamics (SPH) for free-surface flows: past, present and future
.
Journal of Hydraulic Research
54
(
1
),
1
26
.
http://dx.doi.org/10.1080/00221686.2015.1119209
.
Winkler
,
D.
,
Meister
,
M.
,
Rezavand
,
M.
&
Rauch
,
W.
2017
gpuSPHASE – A shared memory caching implementation for 2D SPH using CUDA
.
Computer Physics Communications
213
,
165
180
.
Zhu
,
Y.
&
Bridson
,
R.
2005
Animating sand as a fluid
. In:
ACM Transactions on Graphics (TOG)
,
ACM Trans. Graph. New York, NY, USA
, pp.
965
972
.