Strojniški vestnik - Journal of Mechanical Engineering 59(2013)10, 575-584 © 2013 Journal of Mechanical Engineering. All rights reserved. D0l:10.5545/sv-jme.2013.944 Original Scientific Paper Received for review: 2013-01-02 Received revised form: 2013-04-24 Accepted for publication: 2013-06-17 Reducing the Computational Time of the Smoothed Particle Hydrodynamics Method with a Coupled 2-D/3-D Approach Elvira Džebo1* - Dušan Žagar1 - Matjaž Četina1 - Gregor Petkovšek2 1 University of Ljubljana, Faculty of Civil and Geodetic Engineering, Slovenia 2 HR Wallingford, United Kingdom Smoothed particle hydrodynamics (SPH) is a particle-based Lagrangian method that can be adapted for simulating free surface flows. The method is particularly suitable for phenomena in which the flow changes rapidly. A weakness of the SPH method is the very long computational time. However, this can be significantly reduced using new techniques. In this research, the two-dimensional (2-D) and three-dimensional (3D) models Tis Isat, developed at the University of Ljubljana, were used with an appropriate coupling procedure. The new model was validated and calibrated against the results of laboratory experiments. The simulation results of the 2-D/3-D coupled model were compared to the measurements, to the results of the fully 3-D Tis Isat model and to the results of a one-dimensional (1-D) finite difference (FD) model. The performed SPH simulations showed good agreement with measurements and the 1-D model results in the symmetry axis of the channel. The two greatest advantages of the coupled model are a more realistic description of the water-level below the expansion and the significantly shorter computational time as a result of the adopted coupling procedure. Keywords: SPH method, 2-D/3-D coupling, free-surface flow 0 INTRODUCTION Despite the fact that the Eulerian grid-based approach is commonly used for fluid flow simulations, new Lagrangian particle-based techniques have recently been developed and applied. One of the most widely used particle-based approaches to solving fluid motion equations is smoothed particle hydrodynamics (SPH). The SPH method was developed by Gingold & Monaghan [1] and Lucy [2] for solving astrophysical problems. The method has been modified to simulate the dynamics of solids [3] and [4], liquids [5] and [6], explosions [7] etc. Although in SPH water can be considered incompressible [8], a weakly compressible approach is mostly used. In this study we employed the Tis Isat model, developed at the University of Ljubljana. Tis Isat is a typical weakly compressible SPH model that introduces a non-discrete boundary condition with friction into the SPH method [9]. It has been verified against two benchmark cases [10] and [11]. The results were in a good agreement with the measurements and the results of other available models, while behaviour at the boundaries was better [12], [13] and [9]. The Tis Isat model was also verified on a real case study [14] and the results of the Tis Isat model were compared to measurements and results of a mesh-based model [15]. Probably one of the biggest challenges in the development of the SPH method is reducing the computational time. SPH models are usually coupled with other, computationally less demanding models. So far, different coupling approaches have been proposed: a) coupling the SPH model with mesh methods [16], b) coupling SPH models of different resolutions [17]to [21] and c) coupling SPH models of different dimensionalities [22]. In this work we propose a 2-D/3-D coupling approach in order to significantly reduce the computational time. The coupled 2-D/3-D model can be used in all cases where flow transits from a 2-D to a 3-D state and vice versa. As a good example of such transient flow we simulated two channels, one with an expansion (channel A) and the other with a contraction and an expansion (channel B). 1 METHOD 1.1 SPH Method Smoothed particle hydrodynamics is a meshless particle-based computational method for simulating fluid flows. The interpolating kernel function solves Eq. (1) for any quantity A in an arbitrary point rt: )~jA(r)V(r - r )dr. (1) In discrete form, it can be described as: A (r ) = E jp-AjW ( - r), (2) where mj is the mass of the particle j, Pj is its density and Aj the value of the quantity A of the particle j. Many forms of smoothing kernels W are known. In this research the following form was used: *Corr. Author's Address: University of Ljubljana, Faculty of Civil and Geodetic Engineering, Jamova 2, 1000 Ljubljana, Slovenia, elvira.dzebo@fgg.uni-lj.si 575 W (q ) = - 1 - — q1 + — q3 0 < q < 1 2 4 4 (2 - q )3 1 < q < 2 2 < q (3) Eq. (3) was proposed by Monaghan and Lattanzio [23], where q = l/h (h is the smoothing length and l is the distance between the particles) and w is the normalization factor (w = 7nh2/10 for 2-D simulations and w = nh3 for 3-D simulations). The general equations of fluid motion, the dynamic equation (Eq. (4)) and the equation of state (Eq. (5)) were derived by integration by parts, where gradients are replaced by kernel derivatives: (r i u vr 3 j = v„e„mW„ dt " " 1 " mWj + g, (4) (5) where: (6) (7) vi is the velocity of particle i, ri is the position of particle i, pi is the pressure of particle i, v' is the quotient between kinematic viscosity u and the average density of particles i and j, eij is a directional vector, Wj is a derivative of W with respect to l (distance between particles i and j) and g is the gravity acceleration. The treatment of boundary conditions is one of the most demanding parts of the SPH method. In particle methods, wall boundaries are usually converted into particles [24]. The Tis Isat treats solid boundaries differently; the integration approaches are described in Petkovsek et al. [9]. Furthermore, two eddy viscosity coefficients were defined in the Tis Isat model, as described in [9] and [13]: one is used for the wall-particle (bvis) and another for the particle-particle interactions (avis). The turbulent viscosity between particles is then calculated as: ua = avis ■ d2 \v\jla (8) and the turbulent viscosity between the wall and particles is calculated as: where d is the size of the particle, v is the velocity of the particle, la is the distance between particles and lb is the distance between the boundary and the particle. 1.2 Computation of Water Depth The SPH method is particle-based and therefore the accuracy of the computed water depth depends on the particle size d. The real water level can be located anywhere between the computed value ±d/2 and therefore the particle size should be as small as possible. The Tis Isat model calculates water depth hdep using the equation: (10) hdep = d2 for 2-D simulations and hdep = d3 (11) ub = bvis ■ d (9) for 3-D simulations, as described in [13], where W is the smoothing kernel and d is the particle size. Eqs. (10) and (11) are used to calculate water depth at any chosen point. Particles within a 2h distance from the observed point in the horizontal plane influence the water depth at the observed point (red/darker particles, Figs. 1a and b). The influences of the smoothing kernel (Fig. 1c) for any particle within the 2h distance from the observed point are calculated using Eq. (3). As seen from Fig. 1c, particles closer to the observed point have a greater impact on water depth. The graph curve shows the value of the smoothing kernel for a chosen particle and its influence on water depth. 1.3 2-D/3-D coupling For simulating flow in complex channels, such as a channel with a contraction and/or an expansion (Fig. 5), 3-D models should be applied. While flow phenomena near contractions or expansions cannot be satisfactorily described with any width-averaged 2-D model, such channels usually consist of relatively long straight sections, where width-average modelling is sufficient. In such cases, using the extremely time consuming fully 3-D SPH modelling is not necessary over the entire channel. The key idea of this research was to combine the 2-D and 3-D Tis Isat models and use them simultaneously in order to save on computational time. Both versions of the model can be coupled in two ways: a) coupling a 2-D model upstream with a 3-D model downstream and b) coupling a 3-D model upstream with a 2-D model downstream. When water flows into an expansion, a 2-D model is used in the narrow part of the channel. Near u v.... r 2 h fl Observed point Fig. 1. Calculation of water depth with the 2-D Tis Isat model: a) ground plane of the channel with the observation point, b) longitudinal profile of the channel with the observation point, c) smoothing kernel influence (Eq. (3)). Fig. 3. A surge wave is formed by reflexion from the contracted section, travelling in the upstream direction; this phenomenon usually occurs in channels with a contraction, such as channel B; the arrow denotes flow direction [11] Fig. 2. a) 2-D/3-D and b) 3-D/2-D coupling Fig. 4. Flow area shortly downstream of the expansion with locally supercritical flow, which terminates with an oblique hydraulic jump; this phenomenon usually occurs in channels with an expansion, such as channel A; the arrow denotes flow direction [11]. the expansion the coupling procedure begins. When the 2-D simulation transits into a 3-D simulation, the particles at the distance csdt above the expansion are taken into account (c5 being the speed of sound and dt the time step). This is the distance at which the influence of any flow disruption travelling with the speed of sound fades in time step dt and therefore the smallest coupling length recommended in the Tis Isat model. In this area particles from the 2-D model are multiplied in rows across the channel width to be ready for the flow simulation using the 3-D model (Fig. 2a). Initially the particle velocities and pressures in the 3-D model are the same as at the last cross-section of the 2-D model. In the following time step, all particle parameters are calculated using the general 3-D equations for fluid motion (Eqs. (4) and (5)). This approach was used in channel A (Fig. 5a). The same approach can be applied for coupling an upstream 3-D model with a downstream 2-D model. The particles from the 3-D model, located within |0.5 d| from the symmetry axis of the channel (d being the particle size) and within the distance cs dt above the expansion, are transformed into particles for the 2-D model, maintaining their parameters. Tis Isat moves the centres of these particles into the symmetry axis of the channel and the transformed particles are ready for the 2-D flow simulation using the 2-D Tis Isat model (Fig. 2b). In the 2-D SPH model, solid boundaries are located at a distance ±d/2 from the centre of the channel and in the 3-D SPH model they are placed at the real distance from the symmetry axis. Critical flow occurs at the site of sudden expansion of the channel in all the experiments on the physical model [11]. Thus, at the downstream end of the 2-D SPH model an open boundary condition was prescribed where particles can freely leave the computational domain. However, this solution is not entirely general, as there can be cases where the flow at the expansion is not critical during the whole computational time. For example, if there is an obstacle downstream from the expansion, the flow will be critical in the first moments and later a surge wave can propagate upstream even into the narrow part of the channel. In such cases the methodology should be complemented. One of the possible solutions would be to make an iterative procedure between the first section upstream of the 2-D/3-D limit and the first section downstream, calculating each region alternatively with the 2-D and the 3-D model and iteratively approaching the solution. We plan to complement the model in this manner. We are aware, however, that with such a procedure the computational time will probably increase significantly. In the case of contraction/expansion (channel B, Fig. 5b), the 2-D model was used in the part of the channel above the contraction. However, in this particular case the coupling procedure began sooner (9 m downstream from the dam), as the hydraulic phenomenon occurring upstream from the contraction (described in detail in the following section) is much better simulated using a 3-D model. 2 CASE STUDY In channels with expansion and/or contraction, two typical phenomena occur (shown in Figs. 3 and 4). The first is the reflection of a surge-wave from the constricted cross-section, which propagates in the Table 2. Parameters of the numerical experiments Parameter Abbreviation Value Case 1 Case 2 Case 3 Case 4 0.039 m Particle size d or 0.02 m 0.044 m 0.029 m 0.031 m Number of particles (t = 0 s) np Coupled 2-D/3-D model Tis Isat 8,350 2-D particles or 33,403 2-D particles 9,252 2-D particles + 46,326 3-D particles 10,968 2-D particles 7,911 2-D particles + 37,773 3-D particles Fully 3-D 83,663 3-D particles or 668,4873-D particles 129,710 187,806 116,969 model Tis Isat 3-D particles 3-D particles 3-D particles Coefficient used to calculate turbulent Coupled 2-D/3-D model Tis Isat 0.001 or 0.0011 0.001 0.001 0.0013 viscosity between the wall and particles bvis Fully 3-D model Tis Isat 0.004 or 0.0041 0.003 0.001 0.0015 upstream direction (Fig. 3). This phenomenon occurs in channels with a contraction. The second phenomenon occurs in channels with an expansion. Water depth changes rapidly from a high to a low stage. In the low-stage area the flow is supercritical and velocities are highly non-uniform. At the end of this area an oblique hydraulic jump is formed (Fig. 4). Phenomena shown in Figs. 3 and 4 also occur in channels A and B. Simulations were performed in the two channels shown in Fig. 5. In both cases the width of the accumulation was 0.4 m and the channel slope was 0.087%. The physical model and the 1-D finite difference (FD) model are described in more detail in [11] and [25]. Fig. 5. a) Channel A with an expansion; and b) Channel B with a contraction and an expansion [11] In both channels, simulations over dry and wet channel bed were performed. The four cases and the initial water depths are described in Table 1. Water depths in the accumulation were measured just above the dam and water depths in the channel just below it. Fig. 6. Wave front propagation for a) Case 1; b) Case 2; Case 3 and d) Case 4 Fig. 7. Water surface profile in channel axis at a) t = 15.5 s; b) t = 27 s and c) t = 45 s after the dam break for Case 1 Table 1. Case studies and initial parameters Case Channel Water depth in the accumulation just above the dam [cm] Water depth in the channel just below the dam [cm] Case 1 Channel A -dry channel bed 38.9 0.0 Case 2 Channel A-wet channel bed 43.7 7.0 Case 3 Channel B-dry channel bed 29.2 0.0 Case 4 Channel B-wet channel bed 30.7 5.8 Parameters used in the Tis Isat model are given in Table 2. The particles were 10 times smaller than the water depth in the accumulation. Good results were obtained using this particle size and the computational time was still acceptable. Case 1 was also simulated with smaller particles (d = 2 cm; 20 times smaller than the water depth in the accumulation). The number of particles in each simulation depends on the particle size and the model: with the 3-D model the number of particles is invariable, while in the coupled model, the number of particles varies throughout the simulation. The following parameters were calibrated as described in [22]. The proposed value of the smoothing length ratio is between 1.0 and 2.0 [26] and the Tis Isat model has given good results for the value 1.0 [22]. The value of the coefficient used in the calculation of turbulent viscosity between particles avis was set to 0.01, which corresponds to the value proposed by Gesteira et al. [27]. The coefficient for calculating wall-particle turbulent viscosity was Fig. 8. Axonometric view of the part of the channel with an expansion and water flow simulated with a) a fully 3-D model with larger particles (d = 0.039 m); b) a fully 3-D model with smaller particles (d = 0.02 m); c) a coupled 2-D/3-D model with larger particles (d = 0.39 m); and d) a coupled 2-D/3-D model with smaller particles (d = 0.02 m); Tis Isat model (t = 45 s) for Case 1 (pink and blue particles show areas with higher and lower water depths, respectively) calibrated by comparing the simulated results with the chosen parameter bvis to experimental measurements. Some differences between the used values bvis in the coupled 2-D/3-D and the fully 3-D model are to be expected. In the 2-D Tis Isat model the channel sidewalls are located next to the particles and this gives rise to greater friction between water particles and the wall in comparison with the 3-D model. For this reason lower values of the parameter bvis are needed in the coupled model. The best simulation results are presented herein. a) b) c) 0.07 0.05 0.03 0.01 -0.01 -0.03 0.07 0.05 0.03 0.01 -0.01 -0.03 0.07 0.05 0.03 0.01 -0.01 -0.03 • ■ % .c 5% • • ..©. o feu P**««** *« % • y [ml 0.7 -0.5 -0.3 -0.1 0.1 0.3 0.5 0.7 IQ, •* Q. 8 •• P.' \ % O ■ • Os ./Br * jS** «•• y [m] 0.7 -0.5 -0.3 -0.1 0.1 0.3 0.5 0.7 ■ q • • * H »»HOM °s> r#o •• ¿5 • • < I-*, y [mj -0.7 -0.5 -0.3 -0.1 0.1 0.3 0.5 0.7 ■ measurement —1-D finite difference model 03-D model Tis Isat (d = 0.039 m) • 3-D model Tis Isat (d = 0.02 m) 02-D/3-D model Tis Isat (rf= 0.039 m) • 2-D/3-D model Tis Isat {d = 0.02 m) Fig. 9. Cross-sections at a) x = 32.22 m, b) x = 33.11 m and c) x = 34.00 m at t = 45 s simulated with the 3-D and 2-D/3-D Tis Isat models with both larger and smaller particles for Case 1; the thickness of the slice is equal to the particle size, d 3 RESULTS AND DISCUSSION The results of both models (coupled 2-D/3-D model and the fully 3-D model) were compared to the measurements. Although the 1-D mesh-based models give poorer results, the available 1-D FD model results are also shown [11]. Fig. 6 shows the wave-front propagation for all cases. 3.1 Case 1 Fig. 7 shows a comparison of the measured and simulated data 15.5, 27 and 45 s after dam failure. Compared to measurements in the channel axis, neither of the SPH models gave significantly better results than the 1-D FD model. However, SPH models can better simulate the phenomenon that occurs shortly after the expansion (Fig. 4). In Case 1, all models are in good agreement with the measurements regarding the wave-front propagation (Fig. 6a). Fig. 8 shows part of Channel A in axonometric view, and the water flow simulated with the coupled 2-D/3-D and the fully 3-D model with larger and smaller particle size at t = 45 s. Pink/lighter particles show the areas with higher water depth. Fig. 9 shows cross-sections at x = 32.22 m, x = 33.11 m and x = 34.00 m with the thickness of the slice d obtained from the coupled 2-D/3-D (d = 0.039 m and d = 0.02 m) and the fully 3-D (d = 0.039 m and d = 0.02 m) model at t = 45 s. These two figures show the ability of the 3-D Tis Isat model to simulate an area of supercritical flow and the oblique hydraulic jump that occurs after the expansion (Fig. 4). Both Figs. 8 and 9 show much better results using smaller particles. 3.2 Case 2 In Figs. 10 and 6b, the results of the simulations of the expanding flow over an initially wet channel bed are presented. Fig. 10 shows the water surface profile in the channel axis at t = 16.5 s and 35 s after the dam's collapse. SPH models showed good agreement with the measurements during the expanding flow immediately after the expansion (Figs. 10 and 6b). 3.3 Case 3 Figs. 11 and 6c show the significantly more complex flow in Channel B over an initially dry channel bed. Fig. 11 shows the water surface profile 20, 30 and 58 s after the dam failure. Good results were obtained only before the contraction with all models. N0------- i 0 10 20 30 x[m] 40 50 60 0 10 20 30 x[ml 40 50 60 —channel bottom - • -dam - - expansion from 0.4 m to 1.2 m • measurements —1-D finite difference model —3-D model Tis Isat —coupled 2-D/3-D model Tis Isat Fig. 10. Water surface profile in the channel axis at a) t = 16.5 s and b) t = 35 s (bottom) after the dam break for Case 2 20 a) 0 10 20 30 X[ml 40 50 60 —channel bottom - • -dam - - contraction from 0.4 m to 0.2 m ......expansion from 0.2 m to 1.2 m • measurements —1-D finite difference model —3-D model Tis Isat —coupled 2-D/3-D model Tis Isat Fig. 11. Water surface profile in the channel axis at a) t = 20 s, b) t = 30 s and c) t = 58 s (bottom) after the dam break for Case 3 The 1-D FD model shows good results after the expansion, while the water level calculated by the coupled 2-D/3-D SPH model is too low and the level calculated by the fully 3-D SPH is too high. The wave 0 10 20 30 40 50 60 0 10 20 30 x[ml 40 50 60 —channel bottom - • -dam - - contraction from 0.4 m to 0.2 m ......expansion from 0.2 m to 1.2 m • measurements —1-D finite difference model —3-D model Tis Isat —coupled 2-D/3-D model Tis Isat Fig. 12. Water surface profile in the channel axis at a) t = 19.5 s, b) t = 26.5 s and c) t = 40 s after the dam break for Case 4 front propagation is simulated well with all the models (Fig. 6c). 3.4 Case 4 The last simulation was dedicated to the flow in Channel B over an initially wet channel bed. Water surface 19.5 s after the dam failure is shown in Fig. 12a. All models show good results. Fig. 12b shows the water surface 26.5 s after the dam's collapse. In comparison to measurements, satisfactory results were achieved using SPH models. A comparison of the measured data and simulated water depth 40 s after the dam break is shown in Fig. 12c. Both SPH models show the phenomena above the contraction and below the expansion very well. The wave front propagation is simulated well with all the models (Fig. 6d). 4 DISCUSSION The results obtained by the SPH models can be divided into two groups: in the first three cases the results of the SPH simulations are aproximately as good as the results of the 1-D FD model. SPH models are able to accurately reproduce the hydraulic phenomena which occur before the contraction and after the expansion. Case 4, on the other hand, shows very good results. The agreement of the water levels with the measurements was much better than in other cases. Therefore, all four cases show that SPH simulations are suitable for this type of simulations. We also demonstrated that complex local hydrodynamic phenomena, such as the reflected wave propagation upstream from the contraction (Fig. 3) and the area of supercritical flow and non-uniform velocities ending in an oblique hydraulic jump after the expansion (Figs. 4 and 8) can easily be simulated with the SPH method. Furthermore, as the SPH method is not based on a computational grid, it can be used to effectively describe complex geometries and can therefore be used in highly complex computational domains. The main disadvantage of the SPH method is its high cost with regard to computational time. SPH models give more accurate results using a larger number of particles, which strongly increases the computational cost of simulations, as a search for all neighboring particles must be performed at each time step. Therefore, the development of new techniques for reducing computational time is very important in SPH. All simulations were performed on a quad-core processor. Details regarding the processor used and the required computational times are shown in Table 3. Case 1 was simulated with larger and smaller particles. By using an aproximately 2-fold smaller Table 3. Used processor, required computational time and reduction in computational time using the coupled 2-D/3-D SPH model in comparison to the fully 3-D SPH model Quad-core processor: Intel® Core™ i7 3.33 GHz Installed memory (RAM): 6.00 GB (2.99 GB usable) Computational time [h] Reduction [%] Case Coupled 2-D/3-D SPH model 3-D SPH model Case 1 (d = 0.039 m) 0.50 5.60 91 Case 1 (d = 0.02 m) 11.5 106 89 Case 2 2.25 4.00 44 Case 3 7.00 11.00 36 Case 4 8.25 12.00 31 particle size, computational time was increased 20fold. By using the coupled SPH model in comparison to the fully 3-D SPH model, computational time was significantly reduced in all cases (Table 3). A further decrease in computational time could be achieved in Channels A and B since it is possible to use the 2-D SPH model again in the downstream expansion part of both channels, when the flow again becomes regular and width-averaged. 5 CONCLUSIONS The simulations presented in this work were performed with a coupled 2-D/3-D and a fully 3-D SPH model. It was shown that the coupled 2-D/3-D SPH model is as good as the fully 3-D SPH model. Both SPH models can simulate the hydraulic phenomena occurring before the contraction and after the expansion. The most important advantage of using the 2-D/3-D coupling technique is the significant decrease in the computational time with no decrease in the accuracy of the results. 6 ACKNOWLEDGEMENTS The research was funded by the Ministry of Education, Science, Culture and Sport of the Republic of Slovenia, DEM Maribor, d.o.o., SENG Nova Gorica, d.o.o. and CGS plus d.o.o. Special thanks to Professor Rudi Rajar for his help with the manuscript. 7 REFERENCES [1] Gingold, R.A., Monaghan, J.J. (1977). Smoothed particle hydrodynamics. Theory and application to non-spherical stars. Monthly Notices of the Royal Astronomical Society, vol. 181, p. 375-389. [2] Lucy, L.B. (1977). A numerical approach to the testing of the fission hypothesis. Astronomical Journal, vol. 82, no. 12, p. 1013-1024, D0I:10.1086/112164. [3] Libersky, L.D., Petsheck, A.G., Carney,T.C., Hipp, J.R., Allahdadi, F.A. (1993). High strain Lagrangian hydrodynamics - a three-dimensional SPH code for dynamic material response. Journal of Computational Physics, vol. 109, no. 1, p. 67-75, D0I:10.1006/ jcph.1993.1199. [4] Liu, M.B., Liu, G.R., Lam, K.Y. (2006). Adaptive smoothed particle hydrodynamics for high strain hydrodynamics with material strength. Shock Waves, vol. 15, no. 1, p. 21-29, D0I:10.1007/s00193-005-0002-1. [5] Monaghan, J.J. (1992). Smooth particle hydrodynamics. Annual Review of Astronomy and Astrophysics, vol. 30, p. 573-574, D0I:10.1146/ annurev.aa.30.090192.002551. [6] Vesenjak, M., Ren, Z., Mullerschon, H., Matthaei, S. (2006). Computational modelling of fuel motion and its interaction with the reservoir structure. Strojniški vestnik - Journal of Mechanical Engineering, vol. 52, no. 2, p. 85-100. [7] Liu, M.B., Liu, G.R., Zong, Z., Lam, K.Y. (2003). Computer simulation of the high explosive explosion using smoothed particle hydrodynamics methodology. Computers & Fluids, vol. 32, no. 3, p. 305-322, D0I:10.1016/S0045-7930(01)00105-0. [8] Cummins, J.S., Rudman, M. (1999). An SPH projection method. Journal of Computational Physics, vol. 152, no. 2, p. 587-607, D0I:10.1006/jcph.1999.6246. [9] Petkovšek, G., Džebo, E., Četina, M., Žagar, D. (2010). Application of non-discrete boundaries with friction to smoothed particle hydrodynamics. Strojniški vestnik -Journal of Mechanical Engineering, vol. 56, no. 5, p. 307-315. [10] Martin, J.C., Moyce, W.J. (1952). An experimental study of the collapse of liquid columns on a rigid horizontal plane. Philosophical Transactions of the Royal Society of London, no. 244, p. 312-324. [11] Rajar, R. (1972). Recherche théorique et expérimentale sur la propagation des ondes de rupture de barrage dans une vallée naturelle - Theoretical and experimental research of dam-break waves in natural river valleys. PhD. thesis, University of Toulouse, Toulouse. (In French) [12] Žagar, D., Džebo, E., Četina, M., Petkovšek, G. (2008). Effects of boundary friction in SPH flow simulations. SPHERIC newsletter, no. 7, from http:// wiki.manchester.ac.uk/spheric/index.php/Newsletters, accessed on 2013-07-08. [13] Petkovšek, G., Rajar, R., Četina, M., Žagar, D. (2010). Numerical simulation of flow expansion with SPH. 1st European Division Congress IAHR. [14] Žagar, D., Džebo, E., Četina, M., Petkovšek, G. (2009). Simulations of dam-break and flow through a steep valley using SPH. 33rd IAHR Congress, p. 171-178. [15] Krzyk, M., Klasinc, R., Četina, M. (2012). Two-dimensional mathematical modelling of a dam-break wave in a narrow step stream. Strojniški vestnik -Journal of Mechanical Engineering, vol. 58, no. 4, p. 255-262, D0I:10.5545/sv-jme.2010.216. [16] Narayanaswamy, M., Crespo, A.J.C., Gômez-Gesteira, M., Dalrymple, R.A. (2010). SPHysics - FUNWAVE hybrid model for coastal wave propagation. Journal of Hydraulic Research, vol. 48, no. 1, p. 85-93, D01:10.10 80/00221686.2010.9641249. [17] Berve, S., Omang, M., Trulsen, J. (2001). Regularized smoothed particle hydrodynamics: a new approach to simulating magnetohydrodynamic shocks. The Astrophysical Journal, vol. 561, no. 1, p. 82-93, D0I:10.1086/323228. [18] Berve, S., Omang, M., Trulsen, J. (2005) Regularized smoothed particle hydrodynamics with improved multi-resolution handling. Journal of Computational Physics, vol. 208, no. 1. p. 345-367, D0I:10.1016/j. jcp.2005.02.018. [19] Lastiwka, M., Quinland, N., Basa, M. (2005). Adaptive particle distribution for smoothed particle hydrodynamics. International Journal for Numerical Methods in Fluids, vol. 47, no. 10-11, p. 1403-1409, D0I:10.1002/fld.891. [20] Feldman, J., Bonet, J. (2007). Dynamic refinement and boundary contact forces in SPH with applications in fluid flow problems. International Journal for Numerical Methods in Engineering, vol. 72, no. 3, p. 295-324, D0I:10.1002/nme.2010. [21] Vacondio, R., Rogers, B.D., Stansby, P.K., Mignosa, P., Feldman, J. (2013). Variable resolution for SPH: A dynamic particle coalescing and splitting scheme. Computer Methods in Applied Mechanics and Engineering, vol. 256, p. 132-148, D0I:10.1016/j. cma.2012.12.014. [22] Džebo, E., Žagar, D., Četina, M., Petkovšek, G. (2012). Simulation of dam-break flow in channel expansion with coupled 2-D/3-D SPH model. 7th International SPHERIC Workshop Proceedings, p. 403-408. [23] Monaghan, J.J., Lattanzio, J.C. (1985). A refined particle method for astrophysical problems. Astronomy andAstrophyspcs, vol. 149, no. 1, p. 135-143. [24] Monaghan, J.J. (1994). Simulating free surface flows with SPH. Journal of Computational Physics, vol. 110, no. 2, p. 399-406, D0I:10.1006/jcph.1994.1034. [25] Rajar, R. (1978). Mathematical simulation of dam-break flow. Journal of the Hydraulics Division, vol. 104, no. 7, p. 1011-1026. [26] Violeau, D. (2012). Fluid Mechanics and the SPH Method - Theory and Applications. 0xford University, 0xford, D0I:10.1093/acprof:o so/9780199655526.001.0001. [27] Gesteira, M.G., Rogers, B.D., Dalrymple, R.A., Crespo, A.J.C., Narayanaswamy, M. (2009). User's guide for the SPHysics code. User's Guide.