Repositorio Institucional de Documentos

2014 62 Carmelo Juez Jiménez Development of robust, physicallybased numerical models for transport processes and geomorphodynamics changes Departamento Director/es Ciencia y Tecnología de Materiales y Fluidos Murillo Castarlenas, Javier Director/es Tesis Doctoral Autor Repositorio de la Universidad de Zaragoza – Zaguan http://zaguan.unizar.es UNIVERSIDAD DE ZARAGOZA Departamento Director/es Carmelo Juez Jiménez DEVELOPMENT OF ROBUST, PHYSICALLY-BASED NUMERICAL MODELS FOR TRANSPORT PROCESSES AND GEOMORPHODYNAMICS CHANGES Director/es Ciencia y Tecnología de Materiales y Fluidos Murillo Castarlenas, Javier Tesis Doctoral Autor 2014 Repositorio de la Universidad de Zaragoza – Zaguan http://zaguan.unizar.es UNIVERSIDAD DE ZARAGOZA Departamento Director/es Director/es Tesis Doctoral Autor Repositorio de la Universidad de Zaragoza – Zaguan http://zaguan.unizar.es UNIVERSIDAD DE ZARAGOZA Development of robust, physically-based numerical models for transport processes and geomorphodynamic changes Carmelo Juez Jim´enez A thesis for the title of Doctor of Philosophy Fluid Mechanics Doctoral Program Zaragoza, March 2014 Supervisor: Dr. Javier Murillo Castarlenas Escuela de Ingenier´ ıa y Arquitectura Universidad de Zaragoza Abstract iii Development of robust, physically-based numerical models for transport processes and geomorphodynamic changes Abstract Bed changes in rivers may occur under several morphodynamics and hydrodynamics conditions. The modeling of this type of phenomena can be performed coupling the Shallow Water Equations (SWE) for the hydrodynamic part and the Exner equation for the morphodynamic part. The Exner equation states that the time variation of the sediment layer is due to the sediment transport discharge through the boundaries of the volume. Considering that sediment transport discharge are computed by means of sediment capacity formulae based on 1D experimental steady flows, the assessment of these empirical relations under unsteady 1D and 2D situations must be studied. In order to ensure the reliability of the numerical experimentation, the numerical scheme must handle correctly the coupling between the 2D SWE and the Exner equation under any condition. If possible, it is convenient to express the formulation of different empirical laws under a general framework. In consequence, a finite-volume numerical scheme that includes these two main features has been chosen as a benchmark for comparing the 1D and 2D results obtained when using several well known sediment transport formulae: Meyer-Peter and M¨uller, Ashida and Michiue, Engelund and Fredsoe, Fernandez Luque and Van Beek, Parker, Smart, Nielsen, Wong and Camenen and Larson. In addition, a new interpretation of the Smart empirical law is presented in order to cope with bed load transport over irregular beds of changing slope. Detailed results for this new modified empirical law together with the ones obtained with Meyer-Peter and M¨uller (which is the sediment capacity formula more used in hydraulic engineering) are provided for every test case analyzed. Furthermore, the Root Mean Square Error (RMSE) associated to every formula at each experimental condition is calculated with the purpose of evaluating quantitatively the overall behavior of each one. The results point out that the new interpretation of the Smart formula reaches the most accurate results in all cases, but in a genuinely 2D flow, that is, a situation involving more than one flow direction, the differences among sediment transport formulae are not as noticeable as in the 1D studied situations. Once the forecasting capacity of each sediment transport formula has been studied, another concern is the computational cost. The coupling between the SWE and the Exner equation by means of an augmented Jacobian matrix involves a high number of algebraic operations for computing the eigenvalues and the eigenvectors. Therefore, the computational cost is increased significantly, limiting the applicability of the numerical scheme to realistic situations where large domains are involved. In order to improve the computational efficiency, the coupling technique is modified, not decreasing the number of waves involved in the Riemann Problem but simplifying their definitions. The approach proposed in this thesis is a new strategy to combine concepts from hyperbolic conservation laws and conservative finite volume schemes. With the aim to control numerical stability in the most efficient form possible, a numerical eigenvalue iv Abstract is defined to control the discrete Exner equation in the explicit scheme. This bed wave celerity helps mainly to ensure conservation and to control automatically the numerical stability of the explicit scheme. The effects of the numerical coupling strategy proposed in this thesis are tested against exact solutions and 1D and 2D experimental data. The results emerging from this analysis show that efficiency and accuracy can be obtained when choosing an adequate sediment transport law and the stability condition is augmented by including a new celerity associated to the bed changes. On the other hand, in environmental and civil engineering applications, geomorphological changes are not only present in rivers but also in steep areas where massive mobilizations of poorly sorted material can occur. This sliding material is usually composed by a mixture of sand and water. For simplifying the phenomenon, dry granular flows have been considered as a starting point for the understanding of the physics involved within the landslides. The hypothesis of Saint-Venant equations are considered valid for modeling these land movements. Taking advantage of this approach, in this thesis approximate augmented Riemann solvers are formulated providing appropriate numerical schemes for mathematical models of granular flow on irregular steep slopes. Fluxes and source terms are discretized to ensure steady state configurations including correct modeling of start/stop flow conditions, both in a global and a local system of coordinates. The weak solutions presented involve the effect of bed slope in pressure distribution and frictional effects by means of the adequate gravity acceleration components. The numerical solvers proposed are first tested against 1D cases with exact solution and then are compared with 2D experimental data in order to check the suitability of the mathematical models described in this thesis. Comparisons between results provided when using global and local system of coordinates are presented. Both the global and the local system of coordinates can be used to predict faithfully the overall behavior of the landslides. The performance of the numerical scheme has been studied using novel experimental situations. These laboratory works include bidimensional configurations, the inclusion of obstacles in the flow path and a variable slope in the domain. Hence, a further step in mimicking realistic situations is obtained, since the behavior of the granular flow is affected by the presence of natural elements such as boulders or trees. Three situations have been considered. The first experiment is based on a single obstacle, the second one is performed against multiple obstacles and the third one study the influence of a dike when an overtopping situation takes place. Due to the impact of the flow against the obstacles, fast moving shocks appear, and a variety of secondary waves emerge. Comparisons between computed and experimental data are presented for the three cases. The computed results show that the numerical tool previously developed is able to predict faithfully the overall behavior of this type of complex dense granular flow. Resumen v Desarrollo de esquemas num´ericos robustos y basados en modelos f´ısicos para procesos de transporte y cambios geomorfodin´amicos Resumen Los cambios en la topograf´ıa de los r´ıos pueden ocurrir bajo diferentes condiciones hidrodin´amicas y morfodin´amicas diferentes. La modelizaci´on de este tipo de fen´omenos se puede desarrollar mediante un acoplamiento entre un modelo de aguas poco profundas (SWE) para la parte hidrodin´amica y la ecuaci´on de Exner para la parte morfodin´amica. La ecuaci´on de Exner relaciona la variaci´on temporal del nivel de fondo con los flujos de transporte de sedimento que atraviesan el volumen de control. Considerando que las f´ormulas de transporte de sedimento est´an basadas en situaciones experimentales 1D con flujo estacionario, la validaci´on de estas relaciones emp´ıricas para situaciones transitorias 1D y 2D es imprescindible. Para garantizar la confianza en los resultados computacionales obtenidos, el esquema num´erico empleado debe manejar correctamente el acoplamiento entre las ecuaciones 2D SWE y la ecuaci´on de Exner bajo cualquier situaci´on. Adem´as, es conveniente expresar la formulaci´on de las diferentes formulaciones de transporte general de forma general para que se puedan incorporar con facilidad al esquema num´erico. En consecuencia, un esquema en vol´umenes finitos que incluye ambas caracter´ısticas ha sido utilizado para comparar los resultados 1D y 2D obtenidos con diversas f´ormulas de transporte de sedimentos ampliamente conocidas: Meyer-Peter and M¨uller, Ashida and Michiue, Engelund and Fredsoe, Fern´andez Luque and Van Beek, Parker, Smart, Nielsen, Wong and Camenen and Larson. Adem´as, una nueva interpretaci´on de la f´ormula de Smart es presentada para tener en cuenta el efecto del transporte de fondo sobre topograf´ıas irregulares con pendientes cambiante. Resultados detallados para esta nueva interpretaci´on de la f´ormula junto con los obtenidos con Meyer-Peter y M¨uller (que es la f´ormula de sedimento m´as utilizada en ingenier´ıa hidra´ulica) son mostrados para cada caso analizado. Adem´as, el error cuadr´atico medio asociado a cada f´ormula para cada condici´on experimental es calculado con el prop´osito de evaluar cuantitativamente el comportamiento general de cada relaci´on emp´ırica. Los resultados demuestran que la nueva interpretaci´on de la f´ormula de Smart obtiene los resultados m´as precisos en todos los casos, aunque, en un caso genuinamente 2D, las diferencias entre las leyes de transporte de sedimento no son tan notables como en los casos 1D estudiados. Una vez analizada la precisi´on de los resultados obtenidos con cada formulaci´on de transporte de sedimento, se ha estudiado otro hecho importante como es el coste computacional del esquema num´erico empleado. El acoplamiento entre las SWE y la ecuaci´on de Exner a trav´es de una matriz Jacobiana ampliada requiere un elevado n´umero de operaciones algebraicas para calcular los valores y vectores propios. De esta manera, el coste computacional se incrementa notablemente, limitando la aplicabilidad del esquema num´erico ante situaciones realistas. Para mejorar la eficiencia computacional, la t´ecnica de acoplamiento es simplificada, pero sin reducir el n´umero de ondas involucradas en el problema a Riemann. La aproximaci´on considerada en esta tesis vi Resumen combina conceptos de las ecuaciones conservativas hiperb´olicas y de los esquemas conservativos en vol´umenes finitos. Con el prop´osito de controlar la estabilidad num´erica de la forma m´as eficaz posible, un valor propio num´erico es definido para controlar la ecuaci´on discreta de Exner en el esquema expl´ıcito. Esta celeridad asociada al fondo ayuda principalmente a garantizar la conservaci´on y a controlar autom´aticamente la estabilidad num´erica del esquema expl´ıcito. Los efectos del acoplamiento num´erico propuesto en este trabajo son verificados frente a soluciones exactas y casos experimentales 1D y 2D. Los resultados obtenidos muestran que eficiencia y precisi´on pueden obtenerse si se escoge una formulaci´on de transporte de sedimento adecuada y adem´as, se amplia la condici´on de estabilidad del esquema num´erico para considerar la nueva celeridad asociada a los cambios de fondo. Por otra parte, en la ingenier´ıa medio ambiental y civil, los cambios geomorfol´ogicos no est´an solo presentes en los r´ıos, sino tambi´en en ´areas con fuertes pendientes donde masivas movilizaciones de terreno con escasa cohesi´on pueden producirse. Este material deslizante suele estar compuesto por una mezcla de arenas y agua. Para simplificar el fen´omeno, los flujos granulares secos han sido considerados como un punto de partida para comprender la f´ısica involucrada en estos deslizamientos de terreno. Adem´as, las hip´otesis de las ecuaciones de Saint-Venant son v´alidas para modelar este tipo de movimientos t´erreos. Por ello, esquemas aproximados aumentados de tipo Riemann han sido formulados incorporando las caracter´ısticas propias de flujos que evolucionan sobre pendientes elevadas. Los flujos y t´erminos fuente son discretizados para garantizar la correcta modelizaci´on de las condiciones de parada y comienzo de movimiento tanto en coordenadas locales como en globales. Las soluciones d´ebiles presentadas tienen en cuenta los efectos de las proyecciones de la gravedad en la distribuci´on de presiones y en los t´erminos de fricci´on. Los esquemas num´ericos propuestos son primeramente testados frente a casos 1D con soluci´on exacta y luego, son comparados con casos experimentales 2D para verificar la idoneidad de los modelos matem´aticos propuestos. Los resultados obtenidos con las coordenadas locales y globales son presentados, concluyendo que ambos sistemas de coordenadas pueden ser usados para predecir adecuadamente el comportamiento de los deslizamientos de terreno. Gracias a la herramienta num´erica desarrollada para el c´alculo de deslizamientos de material granular seco, una serie de situaciones experimentales nuevas han sido estudiadas. El denominador com´un de estos ensayos de laboratorio se basa en una configuraci´on bidimensional, la incorporaci´on de obst´aculos al paso del flujo y una pendiente variable en el dominio de estudio. De esta manera, se intenta conseguir un mayor acercamiento a la realidad donde el comportamiento de los flujos granulares est´a influenciado por la presencia de elementos naturales como grandes bloques de piedra o ´arboles. Tres situaciones experimentales han sido consideradas. El primer experimento est´a basado en un ´unico obst´aculo, el segundo es realizado con varios obst´aculos y el ´ultimo, estudia el efecto que la presencia de un dique tiene sobre el flujo. Los resultados muestran una sucesi´on de r´apidos choques, que evolucionan desplegando una variedad de ondas secundarias alrededor de los obst´aculos. La comparativa con los datos experimentales es presentada. Los resultados computacionales muestran que el esquema num´erico es capaz de predecir la evoluci´on del flujo ante este tipo de situaciones complejas. CONTENTS xiii 8 Weakly-coupled numerical scheme 75 8.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75 8.2 Finite Volume Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77 8.3 Approximate Riemann Solution for the Hydrodynamic model . . . . . . 79 8.4 Approximate Riemann Solution for the Morphodynamic model . . . . . 84 8.5 Stability region . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 88 8.6 Geomorphological collapse . . . . . . . . . . . . . . . . . . . . . . . . . 88 9 Weakly-coupled scheme: results 91 9.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91 9.2 Problems with exact solutions . . . . . . . . . . . . . . . . . . . . . . . 91 9.3 One dimensional cases . . . . . . . . . . . . . . . . . . . . . . . . . . . 95 9.3.1 Dam break test cases . . . . . . . . . . . . . . . . . . . . . . . . 95 9.3.2 1D Knickpoint test case . . . . . . . . . . . . . . . . . . . . . . 101 9.4 Two dimensional cases . . . . . . . . . . . . . . . . . . . . . . . . . . . 103 9.4.1 2D Numerical modeling of dam failure . . . . . . . . . . . . . . 103 9.4.2 2D Dam break with a sudden enlargement . . . . . . . . . . . . 107 10 Weakly-coupled scheme: conclusions 113 10.1 Further research . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 114 II Mass motion over steep areas 115 11 Introduction 117 11.1 State of the art of the numerical techniques . . . . . . . . . . . . . . . 118 11.2 State of the art for the experimental works . . . . . . . . . . . . . . . . 120 11.3 Outline . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 122 12 Mathematical model and numerical scheme following local coordinates 123 12.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 123 xiv CONTENTS 12.2 Mathematical model . . . . . . . . . . . . . . . . . . . . . . . . . . . . 123 12.3 Finite Volume Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . 125 12.3.1 Definition of the Riemann problem . . . . . . . . . . . . . . . . 127 12.3.2 Integration of the bed slope source term . . . . . . . . . . . . . 130 12.3.3 Integration of the friction stress source term . . . . . . . . . . . 131 12.3.4 Consistency Condition . . . . . . . . . . . . . . . . . . . . . . . 132 12.3.5 2D first order finite volume model . . . . . . . . . . . . . . . . . 136 12.3.6 Stability region . . . . . . . . . . . . . . . . . . . . . . . . . . . 136 13 Mathematical model and numerical scheme following global coordinates 139 13.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 139 13.2 Mathematical model . . . . . . . . . . . . . . . . . . . . . . . . . . . . 139 13.3 Finite Volume Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . 141 13.3.1 Integration of the bed slope source term . . . . . . . . . . . . . 142 13.3.2 Integration of the friction stress source term . . . . . . . . . . . 143 13.3.3 Approximate solution . . . . . . . . . . . . . . . . . . . . . . . . 144 14 Results following local and global coordinates 145 14.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 145 14.2 Quiescent equilibrium and start/stop flow conditions . . . . . . . . . . 145 14.3 Dam break test cases with exact solution . . . . . . . . . . . . . . . . . 152 14.4 Experimental 1D dam break . . . . . . . . . . . . . . . . . . . . . . . . 154 14.5 Experimental spreading of cylindrical granular mass . . . . . . . . . . . 157 14.6 Experimental spreading of granular mass over a fixed rough inclined plane165 14.7 Experimental spreading of granular mass over a initially static layer . . 167 14.8 Spreading of granular mass over a rough parabolic inclined plane . . . . 174 15 Conclusions for the local and global coordinates 181 15.1 Further research . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 182 CONTENTS xv 16 Small-scale environmental problems 183 16.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 183 16.2 Experimental setup . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 183 16.3 Extra considerations about the friction law . . . . . . . . . . . . . . . . 184 17 Results for the small-scale environmental problems 187 17.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 187 17.1.1 Gravity driven flow facing up a single obstacle . . . . . . . . . . 188 17.2 Gravity driven flow facing up three obstacles . . . . . . . . . . . . . . . 202 17.2.1 Gravity driven flow facing up a dike . . . . . . . . . . . . . . . . 212 18 Conclusions for the small-scale environmental problems 221 18.1 Further research . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 222 Conclusiones generales 225 Bibliography 238 A Calculus of eigenvalues and eigenvectors for the coupled-Jacobian numerical scheme 241 B Conservation of the coupled-Jacobian numerical scheme 243 List of Figures 3.1 Interfaces in the domain . . . . . . . . . . . . . . . . . . . . . . . . . . 14 3.2 Depth averaged quantities within the layers . . . . . . . . . . . . . . . 14 3.3 Reynolds theorem applied in an arbitrary volume . . . . . . . . . . . . 15 3.4 Mass conservation in layer 1 . . . . . . . . . . . . . . . . . . . . . . . 15 3.5 Mass conservation in layer 2 . . . . . . . . . . . . . . . . . . . . . . . 17 3.6 Momentum balance in layer 1 . . . . . . . . . . . . . . . . . . . . . . . 18 3.7 Interfaces in the domain . . . . . . . . . . . . . . . . . . . . . . . . . . 20 4.1 Types of sediment transport . . . . . . . . . . . . . . . . . . . . . . . 25 6.1 Results for test A . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37 6.2 Comparison between MPM and Smart CFBS for test A . . . . . . . . 38 6.3 Errors for test A . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38 6.4 RMSE for test A . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39 6.5 Results for test B . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40 6.6 Comparison between MPM and Smart CFBS for test case B . . . . . . 41 6.7 Errors for test B . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41 6.8 RMSE for test B . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41 6.9 Comparison between MPM and Smart CFBS for test D . . . . . . . . 42 6.10 Results for test D . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43 6.11 Errors for test D . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44 6.12 RMSE for test D . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44 6.13 Comparison between MPM and Smart CFBS for test case F . . . . . . 45 xviii LIST OF FIGURES 6.14 Results for test F . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46 6.15 Errors for test F . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47 6.16 RMSE for test F . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47 6.17 Sketch for the 1D dam failure experimental setup . . . . . . . . . . . . 48 6.18 Bed level results for 1D dam failure test . . . . . . . . . . . . . . . . . 49 6.19 Water depth results for 1D dam failure test . . . . . . . . . . . . . . . 50 6.20 Overtopping results for 1D dam failure test . . . . . . . . . . . . . . . 50 6.21 Errors for 1D dam failure test . . . . . . . . . . . . . . . . . . . . . . . 51 6.22 Sand cube sketch . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52 6.23 Results for the sand cube test . . . . . . . . . . . . . . . . . . . . . . . 53 6.24 Errors for the sand cube test . . . . . . . . . . . . . . . . . . . . . . . 54 6.25 Detail of the triangular mesh for the 2D dam failure test . . . . . . . . 55 6.26 Numerical results for the 2D dam failure test . . . . . . . . . . . . . . 55 6.27 Water depth and bed evolution for 2D dam failure test . . . . . . . . . 56 6.28 Overtopping for the 2D dam failure test . . . . . . . . . . . . . . . . . 56 6.29 RMSE for bed level for the 2D dam failure test . . . . . . . . . . . . . 57 6.30 Sketch for the 2D symmetric dam break test . . . . . . . . . . . . . . 58 6.31 Location of the probes for the 2D symmetric dam break test . . . . . . 58 6.32 Results for the 2D dam break test case . . . . . . . . . . . . . . . . . . 60 6.33 Experimental results for the 2D dam break test case . . . . . . . . . . 60 6.34 Comparison between MPM and Smart CFBS for the 2D symmetric dam break test at section S1 . . . . . . . . . . . . . . . . . . . . . . . 61 6.35 Comparison between MPM and Smart CFBS for the 2D symmetric dam break test at section S2 . . . . . . . . . . . . . . . . . . . . . . . 61 6.36 RMSE for the sections for the 2D symmetric dam break test . . . . . . 61 6.37 Probes comparison (U1-U4) for the 2D symmetric dam break test . . . 62 6.38 Probes comparison (U5-U8)for the 2D symmetric dam break test . . . 63 6.39 RMSE for the probes (U1-U4) for the 2D symmetric dam break test . 63 6.40 RMSE for the probes (U5-U8) for the 2D symmetric dam break test . 64 LIST OF FIGURES xix 6.41 Sketch for the 2D dam break test with a sudden enlargement . . . . . 65 6.42 Triangular mesh for the 2D dam break test with a sudden enlargement 66 6.43 Results for 2D dam break test with a sudden enlargement . . . . . . . 68 6.44 Probes comparison for 2D dam break test with a sudden enlargement . 68 6.45 Sections comparison (S1-S7) for 2D dam break test with a sudden enlargement . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69 6.46 Sections comparison (S9) for 2D dam break test with a sudden enlargement . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70 8.1 Cell parameters . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78 8.2 Riemann problem in 2D along the normal direction to a cell side . . . 79 8.3 Integration control volume for the hydrodynamic model . . . . . . . . 81 8.4 Integration control volume for the morphodynamic model . . . . . . . 85 9.1 Exact and computed solution for Test A . . . . . . . . . . . . . . . . . 93 9.2 Exact and computed solution for Test B . . . . . . . . . . . . . . . . . 93 9.3 Exact and computed solution for Test C . . . . . . . . . . . . . . . . . 94 9.4 Results for test A . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 96 9.5 Time step evolution for test A . . . . . . . . . . . . . . . . . . . . . . 97 9.6 Results for test A without imposing CFL limitation . . . . . . . . . . 97 9.7 Results for test F . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 99 9.8 Time step evolution for test F . . . . . . . . . . . . . . . . . . . . . . 100 9.9 Results for test F without imposing CFL limitation . . . . . . . . . . . 100 9.10 Sketch for the knickpoint test . . . . . . . . . . . . . . . . . . . . . . . 101 9.11 Results for knickpoint test . . . . . . . . . . . . . . . . . . . . . . . . 102 9.12 Time step evolution for knickpoint test . . . . . . . . . . . . . . . . . 103 9.13 Detail of the triangular mesh for the 2D dam failure test . . . . . . . . 104 9.14 Temporal evolution for the 2D dam failure test . . . . . . . . . . . . . 104 9.15 Results for the 2D dam failure test . . . . . . . . . . . . . . . . . . . . 105 9.16 Time step evolution for the waves celerities for the 2D dam failure test 106 xx LIST OF FIGURES 9.17 Time step comparison with the CJ and the WC scheme for the 2D dam failure test . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 106 9.18 Location of probes and sections with the 2D dam break test with a sudden enlargement . . . . . . . . . . . . . . . . . . . . . . . . . . . . 107 9.19 Computed results for the 2D dam break test with a sudden enlargement 109 9.20 Results for the 2D dam break test with a sudden enlargement . . . . . 110 9.21 Sections comparison for the 2D dam break test with a sudden enlargement111 12.1 Sketch for local and global coordinates . . . . . . . . . . . . . . . . . . 124 12.2 Cell parameters . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 125 12.3 2D Riemann problem . . . . . . . . . . . . . . . . . . . . . . . . . . . 126 12.4 Integration control volume . . . . . . . . . . . . . . . . . . . . . . . . 129 12.5 Frictionless quiescent equilibrium in local coordinates . . . . . . . . . 130 12.6 Quiescent equilibrium involving Coulomb stress in local coordinates . . 131 12.7 Integration control volume with source terms . . . . . . . . . . . . . . 134 13.1 Relation among local and global coordinates . . . . . . . . . . . . . . 140 13.2 Frictionless quiescent equilibrium in global coordinates . . . . . . . . . 142 13.3 Quiescent equilibrium involving Coulomb stress in global coordinates . 143 14.1 Initial condition for the numerical test . . . . . . . . . . . . . . . . . . 146 14.2 Results for numerical test when using GC . . . . . . . . . . . . . . . . 147 14.3 Results for numerical test when using LC . . . . . . . . . . . . . . . . 148 14.4 Module of flow velocity for numerical test . . . . . . . . . . . . . . . . 148 14.5 Results with different meshes for numerical test when using GC . . . . 149 14.6 Results with different meshes for numerical test when using LC . . . . 150 14.7 Module velocity for numerical test . . . . . . . . . . . . . . . . . . . . 150 14.8 Mesh refinement for numerical test . . . . . . . . . . . . . . . . . . . . 151 14.9 1D Results for the exact solution test . . . . . . . . . . . . . . . . . . 153 14.10 2D Results for the exact solution test . . . . . . . . . . . . . . . . . . 154 14.11 Sketch for the 1D dam break test . . . . . . . . . . . . . . . . . . . . . 154 LIST OF FIGURES xxi 14.12 3D contour plot of the free surface level when using GC . . . . . . . . 155 14.13 3D contour plot of the free surface level when using LC . . . . . . . . 156 14.14 Results for different meshes for test A . . . . . . . . . . . . . . . . . . 158 14.15 Results for test A using M3. . . . . . . . . . . . . . . . . . . . . . . . 160 14.16 Results for test A using M1and M2A. . . . . . . . . . . . . . . . . . . 161 14.17 Results for test A using M2Band M3. . . . . . . . . . . . . . . . . . 162 14.18 Temporal evolution of maximum modulus of flow velocity for test A . 162 14.19 Results for different meshes for test B . . . . . . . . . . . . . . . . . . 163 14.20 Results for test B using M3. . . . . . . . . . . . . . . . . . . . . . . . 163 14.21 Temporal evolution of maximum modulus of flow velocity for test B . 164 14.22 Results for test C using M3. . . . . . . . . . . . . . . . . . . . . . . . 164 14.23 Temporal evolution of maximum modulus of flow velocity for test C . 165 14.24 3D Initial condition for the inclined plane test . . . . . . . . . . . . . . 165 14.25 Longitudinal initial condition for the inclined plane test . . . . . . . . 166 14.26 Results with different meshes for the inclined plane test . . . . . . . . 168 14.27 Results for the inclined plane test when using LC and GC . . . . . . . 169 14.28 Results 2 for the inclined plane test when using LC and GC . . . . . . 170 14.29 Temporal evolution of thickness contours using GC . . . . . . . . . . . 170 14.30 Results with different friction angles when using LC and GC . . . . . 171 14.31 Results for the inclined plane with an initial layer test when using LC and GC . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 172 14.32 3D results for the inclined plane with and without an initial layer . . . 173 14.33 Sketch for the parabolic chute test . . . . . . . . . . . . . . . . . . . . 174 14.34 3D contour views for the parabolic chute test when using GC . . . . . 175 14.35 Results for the parabolic chute test when using GC . . . . . . . . . . . 177 14.36 Results for the parabolic chute when using LC . . . . . . . . . . . . . 178 14.37 Results comparison for the parabolic chute test . . . . . . . . . . . . . 179 14.38 Results comparison 2 for the parabolic chute test . . . . . . . . . . . . 179 xxii LIST OF FIGURES 16.1 Sketch for the experimental setup . . . . . . . . . . . . . . . . . . . . 184 17.1 Probes location for the experimental work . . . . . . . . . . . . . . . . 188 17.2 Initial configuration for Experiment 1 . . . . . . . . . . . . . . . . . . 189 17.5 Final stage for Experiment 1 with and without considering gravity projections . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 191 17.6 Final stage for Experiment 1 when including Manning’s law . . . . . . 191 17.3 Flow structures for Experiment 1 . . . . . . . . . . . . . . . . . . . . . 194 17.4 Final stage for Experiment 1 with different friction angles . . . . . . . 194 17.7 3D results for Experiment 1 . . . . . . . . . . . . . . . . . . . . . . . . 195 17.8 Velocity field for Experiment 1 . . . . . . . . . . . . . . . . . . . . . . 196 17.9 Comparison for the sand depth for Experiment 1 (t= 540-700 ms) . . 197 17.10 Comparison for the sand depth for Experiment 1 (t= 1000-2000 ms) . 198 17.11 Computational error for Experiment 1 . . . . . . . . . . . . . . . . . . 199 17.12 Probes comparison for Experiment 1 . . . . . . . . . . . . . . . . . . . 200 17.13 Longitudinal profile for Experiment 1 . . . . . . . . . . . . . . . . . . 201 17.14 Initial configuration for Experiment 2 . . . . . . . . . . . . . . . . . . 202 17.15 Flow structures in Experiment 2 . . . . . . . . . . . . . . . . . . . . . 204 17.16 3D results for Experiment 2 . . . . . . . . . . . . . . . . . . . . . . . . 205 17.17 Velocity field for Experiment 2 . . . . . . . . . . . . . . . . . . . . . . 206 17.18 Comparison for the sand depth for Experiment 2 (t= 460-640 ms) . . 207 17.19 Comparison for the sand depth for Experiment 2 (t= 740-1500 ms) . 208 17.20 Computational error for Experiment 2 . . . . . . . . . . . . . . . . . . 209 17.21 Probes comparison for Experiment 2 . . . . . . . . . . . . . . . . . . . 210 17.22 Longitudinal profile for Experiment 2 . . . . . . . . . . . . . . . . . . 211 17.23 Initial condition for Experiment 3 . . . . . . . . . . . . . . . . . . . . 212 17.24 Flow structures for Experiment 3 . . . . . . . . . . . . . . . . . . . . . 213 17.25 3D results for Experiment 3 . . . . . . . . . . . . . . . . . . . . . . . . 215 17.26 Comparison for the sand depth for Experiment 3 (t= 490-710 ms) . . 216 1.1 Goal 3 1.1 Goal The physic of the problems described above needs to be described in order to define suitable mathematical models. A theoretical framework is established for the two phenomena studied in this work: the sediment transport in alluvial channels and the geomorphodynamic flow over steep areas. The relevant formulation of these 2D phenomena derives from the depth-averaged equation of bulk mass conservation, mixture momentum conservation and conservation of the mass of the different sediments. The main goal of this thesis is the development of numerical models able to handle with the mathematical model studied. The milestones overcame in this work for achieving this goal are detailed below •Numerical assessment of several closure equations for bed-load transport Usual formulations for morphodynamic bed load transport are given by empirical sediment transport laws. The different sediment transport capacity formulae, used worldwide to control the erosion and deposition rates that deform the bed in transient cases, are based on equilibrium closure equations obtained from experimental observation in 1D steady cases. Their general applicability to 2D unsteady problems requires a careful analysis to assess whether they are able to predict sediment transport in complex transient flows. This point is of paramount importance and requires a well tested and robust numerical method. For this reason, the numerical scheme proposed in Murillo and Garc´ıa-Navarro (2010a) and based on a Jacobian-coupled model has been chosen as a numerical benchmark for analyzing the relative performance of several well-known formulations under different morphodynamic conditions in 1D and 2D configurations. •Development of a efficient and robust weakly-coupled numerical technique for the shallow water equations and the Exner equation Recent advances in free surface flows over mobile bed have shown that accurate and stable results in realistic problems can be provided if an appropriate coupling between the shallow water equations (SWE) and the Exner equation is performed. This coupling can be done if using a suitable Jacobian matrix, as the one employed in the previous milestone. However, when considering this coupling option despite that the SWE are enhanced by only considering one extra conservation law, i.e, the sediment mass conservation, the computational cost may become unaffordable in situations where the initial SWE for rigid bed can be used involving large time and space scales without giving up to the adequate level of mesh refinement. In order to restore the computational efficiency, the coupling technique has been studied and simplified, not decreasing the number of waves involved in the Riemann Problem but simplifying their definitions. The effects of the approximations made have been tested against experimental data which include transient problems over erodible bed. The simplified model has been formulated under a general framework able to insert any desirable discharge solid load formula. 4Introduction •Development of a 2D dry granular flow solver Landslides, rockfalls and debris avalanches take place when a mixture of mud, sand and rocks slide down a slope together. As suggested by Denlinger and Iverson (2004) the study of granular flows constitutes an starting point for the understanding of the more complex mass movement phenomena mentioned before. For this reason a numerical scheme following Murillo and Garc´ıa-Navarro (2010b) is developed, taking into account the particularities which arises in this type of flows and ensuring quiescent equilibrium stages. Taking advantage on the reliability of the numerical scheme developed, a series of novel experimental cases which represent small-scale up-to-date environmental problems have been studied for delving into the physics of this type of phenomena. 1.2 Outline The outline of the present document is structured in two parts. Part I is devoted to the sediment motion involving the presence of water. Part II addresses the mass motion over steep areas, i.e. landslides. Within Part I the milestones 1.1,1.1 are developed. In Chapter 2a brief introduction of the science of sediments dealing with water bodies is provided. In Chapter 3different mathematical models are described in order to clarify the assumptions made in these type of flow, leading to as suitable description of the problem by means of a reduced set of partial differential equations. In Chapter 4the bed load formulations employed in this work are described, and are written using a differentiable expression. The study of the numerical assessment of several closure equations for bed-load transport, milestone 1.1, is developed within Chapter 5and Chapter 7. The development of a efficient and robust weakly-coupled numerical technique for the shallow water equations and the Exner equation, milestone 1.1, is performed from Chapter 8to Chapter 10. For the achievement of each milestone, several 1D and 2D experimental test cases presenting transient states, complex geometries and wet/dry boundaries have been been compared with the computational results. The particular conclusions and future research line are also included. Part II brings together all the features regarding the mass motion over steep areas. Therefore the milestone 1.1 is addressed here. Chapter 11 provides an introduction of the nature of landslides and the previously laboratory work developed for its study. Chapters 12 and 13 describe the numerical schemes particularized for this phenomena when using both local and global coordinate systems. Chapter 14 presents the verification of the numerical models and Chapter 15 is devoted to the conclusions and further research. Once the reliability of the numerical scheme has been proven, a series of small-scale environmental problem are studied in Chapter 16. Chapter 17 shows the results obtained. Finally, Chapter 18 summarizes the conclusions of these small-scale environmental problem and future research lines. Part I Sediment motion in alluvial channels Chapter 2 Sediment transport The science of sediment transport deals with the interrelationship between flowing water and sediment particles. Despite having been studied since the 1950s and being widely employed in real-life engineering (Nielsen,1992;Julien,1998), the develop in the improvement of sediment management remains at present one of the most active topics in the field of hydraulic research The hydraulic and sediment systems are not static even under nature conditions, the cyclic flooding events cause an imbalance on the sediment processes leading to changes in river and coastal morphology. These differences on the sediment transport behavior can be largely augmented by human activities such as river regulation, agriculture, forestry, dredging, coastal and port construction and soil degradation. For this reason there is a a general agreement about the significance of sediment management in rivers, estuaries and coastal areas. Additionally, sediment not only affects to the morphodynamic changes. Sediment is also a key part of the ecosystem and directly concerns the biodiversity: it is the responsible of the habitat formation and of the adequate ecological and chemical quality of the water volumes. In order to develop a sustainable use of river, coastal and marine environments several practical solutions have been proposed. However, some of this sediment management actions have been only focused on the initial domains of concern, leading to local positive effects, but causing unforeseen negative consequences in other places. For all these reasons, the inclusion of sediment management into river, coastal management programs has became a reality and allows to get closer to the holistic idea which represents the understanding of the water bodies, where the same level of importance should be paid to the whole, the local river, coastal and marine environment, and to the interdependence of each part, i.e. the hydrology, the hydraulic, the sediment, the biology and the pollutant issues. In response to the necessity of the integration of sediment in the water bodies manage- 8Sediment transport ment a computational tool is required for the analysis and prediction of these complex systems. The forecasting capacities of the numerical technology allow to obtain proactive solutions and provides not only a local but also, a global feed-back on how a man-made action or a natural event may alter a particular domain. In this fashion this thesis pushes the development of numerical morphodynamic models for increasing accuracy and efficiency. Since the sediment environmentals (river, coastal and marine) are wide an each one has its own particularities, this document is focused on the river ones. 2.1 State of the art It is generally accepted that two of the fundamental concerns in modern sediment hydraulic engineering practice is the need for accurate and, in the same level of importance, efficient schemes for computing the shallow water equations together with the movement of sediment particles. The numerical strategy proposed must mimic the principal phenomenae observed in the flow field and in the movable bed. In the search for capturing this physically significant processes Hudson and Sweby (2002,2005) studied the influence of steady and unsteady approaches in the mathematical model when computing free surface flows considering a bed-load transport. It was commanded to consider the unsteady system contrary to what was assumed in earlier works (De Vriend et al.,1993;Abderrezzak and Paquier,2011). Ignoring unsteady hydrodynamical effects means that the time scales of the morphodynamics changes are smaller in comparison with the morphodynamic ones and only nearly steady process where the bed changes are generated in a slow way could be computed. Focusing on the numerical techniques employed for obtaining the solution, a classification between asynchronous and synchronous strategies can be established (Aric`o and Tucciarelli,2008). Asynchronous procedures imply that the changes in the bed level are not of enough importance for affecting the hydrodynamic equations during a computational time step. This way, the continuity and momentum equations for the fluid phase are decoupled of the sediment continuity equation. They are also known as uncoupled models. On the other hand, numerical methods which solve at the same time step the hydrodynamic and morphodynamic equations are called synchronous and also, coupled. De Vriend et al. (1993) justified that asynchronous/uncoupled techniques were only valid for a limited range of hydrodynamic regimes governed by low Froude numbers and weak interactions between the flow and bed dynamics. For this reason, other authors, Holly and Rahuel (1990); Cao et al. (2002); Wu and Wang (2004); Xia et al. (2010); Cordier et al. (2011), have studied synchronous/coupled procedures, able to handle a wider range of hydrodynamic and morphodynamic situations. In some of those previous works, despite considering an extra equation for computing the sediment dynamics no additional conditions to the classical Courant-Friedrichs-Lewy (CFL) were provided for controlling the numerical stability. In particular, the lack of knowledge of an automatic numerical stability condition in Wu et al. (2012) has driven to calibrate, 2.1 State of the art 9 by trial and error a CFL condition for obtaining a stable solution to each particular case. In order to overcome the challenge when building a self-stable numerical scheme, several strategies have been proposed: ones are based on the development of the exact form of the eigenvalues through the mathematical model (Kassem and Chaudry,1998;Tassi et al.,2008;Lyn and Altinakar,2002;Cao et al.,2006;Gouti`ere et al.,2008) and other in the numerical treatment of the whole set of equations (Hudson and Sweby,2002, 2005). This work is focused on this last idea. In Hudson and Sweby (2002,2005) thanks to the Riemann theory and using a Roe’s approximate Jacobian matrix of the whole system of equations was developed. Hence, the hydrodynamic and morphodynamic equations were not only solved at the same time step but also the wave celerities, which participate in the stability condition, incorporated information from both phases: water and sediment. The term coupled-Jacobian will be used for that model from now on. The main drawback of this Jacobian matrix was a strong dependence on the bed reference level. Additionally, this Jacobian matrix included the definition of the sediment transport formula through the Grass law, Grass (1981). This formula is based in a power law of the velocity, which is nicely differentiable, and in a global calibration parameter, which is unique for all the computational domain and must be tuned in each particular problem. Following with the Jacobian-coupled strategy, other schemes have been proposed and extended to 2D triangular meshes more recently. In Castro Diaz et al. (2009) the identification of the approximate Jacobian matrix was achieved by means of the distribution theory (Dal Maso et al.,1995). However, this numerical technique needs to select families of paths that cannot be generalized. In Soares-Frazao and Zech (2010) a first order HLLC scheme was proposed and a novel wave-speed estimator was provided for the Exner equation. The results were affected by numerical diffusion and a fine mesh was required by obtaining accurate results. The work in Rosatti et al. (2008a) described a Roe solver for a two-phase problem where the attention was devoted to the non-linear relations between primitive and conserved variables. Only the 1D approach of the problem was studied. In Canestrelli et al. (2010); Siviglia et al. (2013) high order numerical techniques were explained over fixed and mobile beds. However, no clear evidence of the behavior of the numerical scheme under a real and experimental case is provided, since only a laboratory test case is studied in the second word. Additionally, the high computational cost of such schemes is not addressed. In Murillo and Garc´ıa-Navarro (2010a) a novel coupled-Jacobian model was proposed and the Jacobian matrix was built with independence of the bed level reference. Regarding the calibration coefficient of Grass law, the uniqueness of this parameter in all the problem was avoided, (Murillo and Garc´ıa-Navarro,2010a), by writing the law in terms of several bed-load sediment transport formulae. Numerical solutions obtained probed to be robust and accurate. Nevertheless, the applicability of this numerical scheme to a real situation, where the domain contains kilometers of river and several types of sediment, is in somehow limited by the computational cost, which is prohibitively expensive. The computational time is highly penalized by the number of 10 Sediment transport algebraic operations need for computing the eigenvectors and eigenvalues of the augmented Jacobian matrix. In order to overcome this huge numerical effort in Serrano et al. (2012) a partially coupled model was proposed, although the quality of the results were compromised by the poor sediment transport law employed. Furthermore, no clear evidence of the effect of the bed wave speed in the time step restriction was provided. 2.2 Outline Following the previous effort made by the authors mentioned above, the concern of this part of the work is twofold: accuracy and efficiency. The first one is related with the sediment transport law employed. Several well known capacity formulae based on 1D experimental steady flows have been analyzed under unsteady 1D and 2D situations. Moreover, a new interpretation of the Smart (1984) empirical law is presented in order to cope with bed load transport over irregular beds of changing slope. Detailed results for this new modified empirical law together with the ones obtained with Meyer-Peter and M¨uller (1948) (which is the sediment capacity formula more used in hydraulic engineering) are provided for every test case analyzed. Furthermore the Root Mean Square Error (RMSE) associated to every formula at each experimental condition is calculated with the purpose of evaluating quantitatively the overall behavior of each one. In order to ensure the reliability of the numerical experimentation the coupledJacobian model previously developed and tested in Murillo and Garc´ıa-Navarro (2010a) has been used. The second objective, the efficiency, has been addressed studying a novel weaklycoupled numerical strategy for coupling the hydrodynamic and the morphodynamic models. The sediment transport law employed in the simulations have been the most accurate one chosen from the previous analysis. Some of the experimental test cases employed when studying the accuracy among the sediment discharge formulae are also employed in this part of the work. Results of the computational cost between the coupled-Jacobian scheme and the weakly-coupled scheme are provided. Furthermore, a clear evidence of how a non-carefully treatment of the stability condition can ruined the numerical results is provided. The outline of this part is as follows: in Chapter 3different mathematical models are described in order to clarify the assumptions made in these type of flow, leading to a suitable description of the problem by means of a reduced set of partial differential equations. In Chapter 4the bed load formulations employed in this work are described, and are written using a differentiable expression. Chapter 5presents the coupledJacobian numerical scheme used to study the differences among the bed-load discharge formulae and a novel numerical discretization of the Smart one is provided. Chapter 6displays the numerical results obtained with the coupled-Jacobian model in 1D and 2D test cases when using several sediment discharge formulae. Chapter 7is devoted to the conclusions about the bed-load formulae studied and further research. The 2.2 Outline 11 novel weakly-coupled numerical scheme proposed in this thesis is showed in Chapter 8. Results obtained with this numerical strategy are depicted in Chapter 9and conclusions and future research lines are written in Chapter 10. 3.2 Two layer model 19 Sediment mass, layer 1 ∂(ρ1h1φ1) ∂t +∂(ρ1h1u1φ1) ∂x −ρ1Ψnet s2,1= 0 (3.17) Sediment mass, layer 2 ∂(ρ2h2φ2) ∂t +∂(ρ2h2u2φ2) ∂x −ρbΨnet s3,2+ρ1Ψnet s2,1= 0 (3.18) Momentum of the mixture, layer 1 ∂(h1ρ1u1) ∂t +∂(h1ρ1u2 1) ∂x +g∂ ∂x 1 2ρ1h2 1+g∂(zρ1h1) ∂x = = (ρ2u2Ψ2,1−ρ1u1Ψ2,1)−τ2,1(3.19) Momentum of the mixture, layer 2 ∂(h2ρ2u2) ∂t +∂(h2ρ2u2 2) ∂x +g∂ ∂x 1 2ρ2h2 2+ρ1h1h2+g∂(zρ2h2) ∂x = =−(ρ2u2Ψ2,1−ρ1u1Ψ2,1+ρ2u2Ψ3,2) + (τ2,1−τ2,3) (3.20) In case that the granular material is not homogeneous in size or density, the subscript pshould be employed to distinguish among cases with non uniform size and specific weight distributions, inside each liquid-granular layer. Thus the formulation becomes: Total mass, fraction p, layer 1 ∂(ρ1ph1) ∂t +∂(ρ1ph1pu1) ∂x +ρ1pΨnet 2,1= 0 (3.21) Total mass, fraction p, layer 2 ∂ρ2ph2p ∂t +∂(ρ2ph2pu2) ∂x −ρ2pΨnet 3,2+ρ2pΨnet 2,1= 0 (3.22) Sediment mass, bed ∂(ρbzφb) ∂t +ρbΨnet s3,2= 0 (3.23) Sediment mass, fraction p, layer 1 ∂(ρ1ph1pφ1p) ∂t +∂(ρ1ph1pu1φ1p) ∂x −ρ1pΨnet s2,1= 0 (3.24) Sediment mass, fraction p, layer 2 ∂(ρ2ph2pφ2p) ∂t +∂(ρ2ph2pusφ2p) ∂x −ρbΨnet s3,2+ρ1pΨnet s2,1= 0 (3.25) 20 Mathematical model Momentum of the mixture, layer 1 ∂(h1pρ1pu1) ∂t +∂(h1pρ1pu2 1) ∂x +g∂ ∂x 1 2ρ1ph2 1p+g∂(zρ1ph1p) ∂x = = (ρ2pu2Ψ2,1−ρ1pu1Ψ2,1)−τ2,1(3.26) Momentum of the mixture, layer 2 ∂(h2pρ2pu2) ∂t +∂(h2pρ2pu2 2) ∂x +g∂ ∂x 1 2ρ2ph2 2p+ρ1ph1ph2p+g∂(zρ2ph2p) ∂x = =−(ρ2pu2Ψ2,1−ρ1pu1Ψ2,1+ρ2pu2Ψ3,2) + (τ2,1−τ2,3) (3.27) There are (5 + 2)Npequations, being Npthe number of size fractions pwhich had the bed material, and (5 + 2)Npindependent variables: the flow depth for each layer, h1p, h2p; the depth averaged velocity in layer 1, u1and in layer 2, u2; the bottom elevation, zand finally the sediment concentration in layer 1 of fraction p,φ1pand in layer 2, φ2p. Several closure equations are required to express the shear stress between layers, τ2,1 and τ2,3and the sediment fluxes, Ψ2,1and Ψ3,2in terms of the independent variables. This represents such a complex task that it justifies further simplification of the model. Next section is devoted to discuss this. 3.3 One layer model The one layer model, Figure (3.7), is built upon a set of assumptions in relation to the two layer model: (i) a unique layer of depth his considered, which includes previous layer 1 and 2, (ii) continuity approach, assuming the same velocity for the liquid and for the solid phase, u, which leads to continuity of momentum and consequently to a continuity of shear stresses, avoiding the necessity of calculating τij between interfaces. Q z, φb h1 - Fluid layer 2 - Solid layer Figure 3.7: Interfaces in the domain 3.3 One layer model 21 3.3.1 Conservation equations The relevant formulation of the model derives from the depth-averaged equation of bulk mass conservation, mixture momentum conservation and conservation of the mass of the different constituents. The term φprepresents the scalar depth-averaged volumetric concentration of component p, with p= 1, ..., Npand Npthe number of different components transported. The mixture density is given by ρm=ρwrwhere ρwis the density of the water and rmeans the relative density of the bulk mixture with respect the clean water r= 1 + Np X p=1 ∆pφp(3.28) where ∆p= (ρp−ρw)/ρwis the relative density of the solid phase p. It is assumed that dissolved species with low concentration do not change bulk density ∆p= 0. In case of having an unique specie the relative density of the bulk mixture becomes r= 1 + ∆φ. Mass conservation Considering a generic control volume for a horizontal flow over a mobile bed where the velocity is depth averaged, defined in Figure 3.7, the Reynolds transport for mass conservation at the liquid layer leads to: Zx2 x1 ∂ ∂t(ρmh)dx +Zx2 x1 ∂ ∂x(ρmhu)dx +Zx2 x1 ρmΨnet sdx = 0 (3.29) and to the sediment balance mass Zx2 x1 ∂ ∂t(h Np X p=1 ρpφp)dx +Zx2 x1 ∂ ∂x(hu Np X p=1 ρpφp)dx +Zx2 x1 Np X p=1 ρpΨnet pdx = 0 (3.30) Following the same procedure for mass sediment balance at the bottom and considering, φbp= (1 −pp), being ppthe porosity of each sediment, drives to Zx2 x1 ∂ ∂t(z Np X p=1 ρp(1 −pp))dx −Zx2 x1 Np X p=1 ρpΨnet pdx = 0 (3.31) 22 Mathematical model The term Ψnet s, which appears in the above set of equations, includes the vertical sediment flux, suspension transport, and the horizontal sediment flux, bed load transport. Ψnet s= Ψload + Ψsusp (3.32) Momentum equation For the xdirection, and considering Figure 3.7, the momentum conservation equation for the mixing layer, which is the unique zone where there exists velocity, and with the xcomponent of the gravity mass force equal to 0, (fv)x= 0, leads to Zx2 x1 ∂ ∂t(ρmhu)dx +Zx2 x1 ∂ ∂x(ρmhu2)dx =Zx2 x1 pbdx −Zx2 x1 τbdx (3.33) The term of superficial forces, (fs)x, has been split in its two components, the hydrostatic pressure, pb, and the friction term exerted over the bed, τb. fs=pb−τb(3.34) Summary of conservation equations The set of developed differential equations is newly written below in terms of the relative density r. Total mass ∂(hr) ∂t +∂(hur) ∂x = Np X p=1 ∆Ψnet p(3.35) Sediment mass, bed ∂(z) ∂t + Np X p=1 Ψnet p (1 −pp)= 0 (3.36) Sediment mass of the mixing layer for specie p ∂(hφp) ∂t +∂(huφp) ∂x = Ψnet p(3.37) Momentum of the mixing layer ∂(hur) ∂t +∂[hu2r+ (1/2)gh2r] ∂x =pb ρw−τb ρw (3.38) 3.3 One layer model 23 There are 3+Npequations and 3+Npvariables, being Npthe number of species: the flow depth, h, the mean flow velocity, u, the bed level, z, and the depth averaged sediment concentration, φp. Furthermore two closure equations are still needed, one for the bed shear stress, τb, and another one for the formulation of the sediment flux between flow and bed, Ψnet p. In the search for the simplest model involving the minimum number of closure relations, the above formulation can be transformed into Exner equation, next presented. 3.3.2 Exner equation The above set of equations may be manipulated leading to a simpler model. Inserting Ψpfrom (3.36) in (3.37) leads to the following sediment mass conservation, ∂z ∂t +1 (1 −pp) ∂(huφp) ∂x =−1 (1 −pp) ∂(hφp) ∂t (3.39) The second term on the left hand side of (3.39) is the derivative of transported sediment flow qs,x =huφpalong the xcoordinate, whereas the term on the right side contains information about the temporal evolution of the bed level due to vertical fluxes of material in cases of suspended material. They become the Exner equation (Kalinske, 1947), expressed as follows ∂z ∂t +ξ∂qs,x ∂x =ξωs(Es−cb) (3.40) with ξ=1 1−pp,ωsthe settling velocity of the sediment particles, Esa dimensionless factor accounting for the sediment material entering the volume by suspension and φb is the suspended material concentration. Both terms of qs,x and ξωs(Es−φb) can be estimated if using empirical closure formulae, that depend on the flow conditions. Regarding the bulk density, it can be evaluated assuming that the volumetric concentration is given by the closure formulae themselves (Rosatti et al.,2007). In many environmental problems, the bulk density remains almost constant and furthermore low concentrations of transported material are present. This means that further simplications over liquid phase mass and momentum conservation equations are admissible, allowing the elimination of the dependence with the relative density of the mixture, r. Alternatively, assuming that the sediment material presents low concentration and does not change the bulk density, the relative density of the mixture, rcan be made constant and equal to 1. Gathering the depth averaged set of equations which governed the 1D flow and the sediment dynamics leads to 24 Mathematical model Mass ∂(h) ∂t +∂(hu) ∂x = 0 (3.41) Momentum ∂(hu) ∂t +∂[hu2+ (1/2)gh2] ∂x =pb ρw−τb ρw (3.42) Sediment mass, bed ∂z ∂t +ξ∂qs,x ∂x =ξωs(Es−cb) (3.43) The extension of the formulation of the shallow water equations to unsteady 2D flow over mobile bed using the Exner equation approach is: Mass ∂(h) ∂t +∂(hu) ∂x +∂(hv) ∂y = 0 (3.44) Momentum in xdirection ∂(hu) ∂t +∂[hu2+ (1/2)gh2] ∂x +∂(huv) ∂y =pbx ρw−τbx ρw (3.45) Momentum in ydirection ∂(hu) ∂t +∂(huv) ∂x +∂[hv2+ (1/2)gh2] ∂y =pby ρw−τby ρw (3.46) Sediment mass, bed ∂z ∂t +ξ∂qs,x ∂x +ξ∂qs,y ∂y =ξωs(Es−cb) (3.47) with (u, v) the depth averaged components of the velocity vector along the (x, y) coordinates. Considering that the present work is focused on the numerical simulation of bed load transport since the influence of the suspended load is assumed negligible, the Exner equation turns into a reduced form: ∂z ∂t +ξ∂qs,x ∂x +ξ∂qs,y ∂y = 0 (3.48) Chapter 4 Bed load transport 4.1 Introduction Sediment transport includes suspended and bed-load sediment transport. Suspended sediment is present when the flux is intense enough for allowing the sediment grains to move away from the bed. Bed-load transport is the kind of sediment motion where the grains roll, slide or even jump over the bed, Figure 4.1. In this work it is faced the study of bed-load sediment transport and the suspended transport is neglected. Rolling Sliding Jumping Figure 4.1: Types of sediment transport As it has been depicted in the previous chapter, when using the Exner equation, horizontal solid fluxes can be evaluated using capacity formulae. In this chapter, different formulations empirically proposed for the modeling of non-cohesive granular material flows are presented and written following a unified expression. In this work mass exchange fluxes associated to suspended load will be considered negligible in comparison with bed load transport, and therefore will not be included in the mathematical model. 4.2 Description of bed load formulation Considering a bidimensional flow where the solid transport is focused on the bed load, the Exner equation can be written as 26 Bed load transport ∂z ∂t +ξ∂qs,x ∂x +ξ∂qs,y ∂y = 0 (4.1) The formulation of the bed load discharge qscan be based on deterministic laws (MeyerPeter and M¨uller,1948), (Camenen and Larson,2005), (Smart,1984) or in probabilistic methods (Kalinske,1947), (Einstein,1950), always supported by experimentation. Grass (Grass,1981) discussed one of the most basic sediment transport laws that in 2D can be written as (Hudson,2001) qs,x =Aguu2+v2qs,y =Agvu2+v2(4.2) This deterministic formulation is well suited for the modeling of non-cohesive granular material and, as a basic feature, this model does not involve any sediment movement threshold but assumes that the flow is always able to mobilize the bed. The model requires a dimensional calibration constant Ag, accounting for the effects associated to the grain size and the kinematic viscosity. Ranging typically from 0 to 1, it represents a stronger interaction between flow and sediment as it approaches 1. Following the idea presented in (Murillo and Garc´ıa-Navarro,2010a), Agcan be determined by using the empirical deterministic formulae avoiding the necessity of expressing this quantity as a calibration constant in each particular problem. To do this, several empirical formulations for sediment transport will be analyzed assuming that it is possible to write them all as Ag=Ag(h, qs,x, qs,y) (4.3) The bed load transport is often represented by the following dimensionless parameter, Φ = |qs| pg(s−1)d3 m (4.4) where s=ρp/ρwis the ratio between solid material (ρp) and water densities, and dm is the median diameter. The dimensionless bottom shear stress or Shields parameter, can be expressed as: θ=|Tb| g(ρs−ρw)dm (4.5) where Tb= (τb,x, τb,y) is the shear stress at the bottom due to the steady flow, that written in terms of the Manning-Strickler’s coefficient (16.2) can be expressed as 4.2 Description of bed load formulation 27 Formula Φ Meyer-Peter and M¨uller (1948) 8 (θ−θc)3/2 Ashida and Michiue (1972) 17 (θ−θc)(√θ−√θc) Engelund and Fredsoe (1976) 18.74 (θ−θc)(√θ−0.7√θc) Luque and van Beek (1976) 5.7 (θ−θc)3/2 Parker (1979) fit to Einstein (1950) 11.2θ3/2(1 −θ/θc)9/2 Smart (1984) 4 (d90/d30)0.2S0.6 oCθ1/2(θ−θS c) Nielsen (1992) 12 θ1/2(θ−θc) Wong (2003) 4.93 (θ−θc)1.6 Wong (2003) 3.97 (θ−θc)3/2 Camenen and Larson (2005) 12 θ3/2exp (−θ/θc) Table 4.1: Summary of sediment formulae τb,x ρw=ghSf,x Sf,x =n2u√u2+v2 h4/3 τb,y ρw=ghSf,y Sf,y =n2v√u2+v2 h4/3 (4.6) This allows to express |Tb|as |Tb|=qτ2 b,x +τ2 b,y =q(ρwghSf,x)2+ (ρwghSf,y)2(4.7) leading to the following expression for the Shields parameter: θ=n2 (s−1)dmh1/3(u2+v2) = n2 (s−1)dmh1/3|u|2(4.8) Different commonly applied empirical deterministic formulae are written in terms of Φ and θ. The formulae tested in this work are gathered in Table 4.1, where d90 and d30 are the grain diameter for which 90% and 30% of the weight of a nonuniform sample is finer respectively, Cis the flow resistance factor C=u/(ghSf)0.5,Sois the bed slope, θcis the critical Shields parameter, Table 4.2, expressing the sediment movement threshold, and θS c=θccos φ1−tan φ tan ψ(4.9) with φthe angle of the bed slope and ψthe angle of repose of saturated bed material. Using (4.8) and (4.4) the transport formulae in (4.1) can be expressed as |qs|=K0K1(u2+v2)3/2=Ag|u|3(4.10) 28 Bed load transport Formula θc Meyer-Peter and M¨uller (1948) 0.0470 Ashida and Michiue (1972) 0.0500 Engelund and Fredsoe (1976) 0.0500 Fern´andez Luque and Van Beek (1976) 0.037–0.0455 Parker (1979) fit to Einstein (1950) 0.030 Smart(1984) 0.0470 Nielsen (1992) 0.0470 Wong (2003) 0.0470 Wong (2003) 0.0495 Camenen and Larson(2005) 0.0400 Table 4.2: Summary of threshold of non dimensional shear stress with Ag=K0K1,K0=g1/2n3 (s−1)h1/2and K1varying in each case as displayed in Table 4.3. Formula K1 Meyer-Peter and M¨uller (1948) 8 (1 −θc/θ)3/2 Ashida and Michiue (1972) 17 (1 −θc/θ)(1 −pθc/θ) Engelund and Fredsoe (1976) 18.74 (1 −θc/θ)(1 −0.7pθc/θ) Fern´andez Luque and Van Beek (1976) 5.7 (1 −θc/θ)3/2 Parker (1979) fit to Einstein (1950) 11.2 (1 −θ/θc)9/2 Smart(1984) 4 (d90/d30)0.2S0.6 oC(1 −θc/θ) Nielsen (1992) 12 (1 −θc/θ) Wong (2003) 4.93 (1 −θc/θ)3/2(θ−θc)0.1 Wong (2003) 3.97 (1 −θc/θ)3/2 Camenen and Larson(2005) 12 exp (−θ/θc) Table 4.3: Summary of Grass coefficients written for sediment formulae These more complex definitions provided for Ag, (4.10), allows to standardize sediment transport formulae and perform a study about their relative behavior under different hydrodynamic and morphodynamic conditions. Chapter 6 CJ scheme: numerical results 6.1 Introduction This Chapter gathers 1D and 2D cases with experimental data in order to study the relative behavior of the numerical results predicted when using different sediment transport formulae. These closure laws were derived from 1D experimental steady flows and are going to be tested in order to verify their capacity of prediction in unsteady situations. In addition, a novel numerical discretization of the classical Smart formula is also tested. Firstly 1D results are presented. A series of sudden dam break test cases are presented, with a combination of morphodynamic and hydrodynamic situations. In the next test case, dam erosion in time due to flow overtopping is considered. In all these numerical experiments the flow finds different regions under subcritical or supercritical regime. The last experiment considers a case of fully subcritical flow, with an important discontinuity at the bottom. The second section of this chapter is devoted to 2D hydro-morphodynamic changes. The dam failure is the first test case studied. Then, two test cases of dam break over a channel with a a/symmetric enlargement are analyzed. 6.2 One dimensional cases 6.2.1 Dam break test cases These experiments were performed in a flume designed at the UCL Civil Engineering Department (Spinewine and Zech,2007). The flume had a length of 6 m, 3 m on both sides of a central gate simulating an idealized dam. The channel width was set constant and equal to 25 cm. The bed material was uniform coarse sand with the 36 CJ scheme: numerical results Test hLhRzLzR A 0.35 0.00 0.00 0.00 B 0.40 0.00 -0.05 0.00 D 0.25 0.00 0.10 0.00 F 0.25 0.10 0.10 0.00 Table 6.1: Summary of dam break test cases following properties: particle sizes ranging from 1.2 to 2.4 mm, with d50 = 1.82 mm, density ρs= 2683 kg m−3, a friction angle ϕ= 30o, negligible cohesion, porosity p= 0.47 and was characterized by a Manning roughness factor n= 0.0165 sm−1/3. Table 6.1 summarizes the set of experiments selected in this work. The regions upstream and downstream the gate were filled with sediments and different water depths. The three first test cases, A, B, and D, have been chosen to guarantee the correct performance of the numerical scheme in combination with a discharge formulation, in cases where morphological changes are produced in presence of dry bed and null, adverse or in favorable slope. Case F allows checking if the numerical scheme in combination with a discharge formulation is able to handle with the different type of waves that may arise in a dam break case over wet bed. Numerical simulations have been performed using ∆x= 0.01 m and CFL = 1.0. In all the simulations the bed domain is considered deformable and no boundary condition is imposed at the downstream section. 6.2.2 Test A Test A is a dam break over dry bed with an initially plain bed level. The flow evolves in time leading to a left moving rarefaction wave upstream the gate ending in a flooding front dominated by friction. The experimental results are close to those ones obtained for dam break cases over dry and fixed bed (Dressler,1954). In Figure 6.1 numerical results and experimental data have been plotted for test case A, for times ranging from 0 to 1.5 seconds. The front wave is numerically well reproduced in space and time when using Smart CFBS. Figure 6.2 shows the numerical results and experimental data for the dam break test case A using MPM (left) and Smart CFBS (right). In this case little scour is produced and both formulations provide indistinguishable results. The Smart CFBS formulation provides a correct tracking of the advance velocity, bed level and water level surface in time, as shown in Figure 6.1. Considering that the numerical scheme is conservative, differences among measured and computational data are expected to be produced by the lack of an infiltration parameter in the numerical model. 6.2 One dimensional cases 37 -0.1 0 0.1 0.2 0.3 0.4 -1 0 1 2 3 h+z (m) x (m) -0.1 0 0.1 0.2 0.3 0.4 -1 0 1 2 3 h+z (m) x (m) -0.1 0 0.1 0.2 0.3 0.4 -1 0 1 2 3 h+z (m) x (m) -0.1 0 0.1 0.2 0.3 0.4 -1 0 1 2 3 h+z (m) x (m) -0.1 0 0.1 0.2 0.3 0.4 -1 0 1 2 3 h+z (m) x (m) -0.1 0 0.1 0.2 0.3 0.4 -1 0 1 2 3 h+z (m) x (m) Figure 6.1: Numerical results and experimental data for the dam break test case A at times t= 0.025, 0.050, 0.075, 0.100, 0.125 and 1.5 s, using a variable value of Ag computed using Smart CFBS: measured water level surface (− • −), measured bed level surface (−◦−), computed water level surface (−4−), computed bed level surface (−N−) 38 CJ scheme: numerical results -0.1 0 0.1 0.2 0.3 0.4 -1 0 1 2 3 h+z (m) x (m) -0.1 0 0.1 0.2 0.3 0.4 -1 0 1 2 3 h+z (m) x (m) Figure 6.2: Numerical results and experimental data for the dam break test case A at t= 1.5 s, using a variable value of Agcomputed using MPM (left) and Smart CFBS (right): measured water level surface (− • −), measured bed level surface (− ◦ −), computed water level surface (−4−), computed bed level surface (−N−) The similarity among computational results for the different discharge formulations is clear when observing Figure 6.3, that displays the modulus of the water level surface error (left) and bed level (right) error in xfor the different formulations at t= 1.5 s. The RMSE (Root median square error) for the different formulations plotted at Figure 6.4, confirms that in plain bed, accurate results are given by all formulas. 0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 1.2 Error (m) x (m) MPM ASHIDA-MICHIUE ENGELUND-FREDSOE FERNANDEZ LUQUE-VAN BEEK PARKER SMART NIELSEN WONG (4.93) WONG (3.97) CAMENEN-LARSON SMART CFBS 0 0.005 0.01 0.015 0.02 0.025 0.03 0.035 0.04 0.045 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 Error (m) x (m) MPM ASHIDA-MICHIUE ENGELUND-FREDSOE FERNANDEZ LUQUE-VAN BEEK PARKER SMART NIELSEN WONG (4.93) WONG (3.97) CAMENEN-LARSON SMART CFBS Figure 6.3: Modulus of the water level surface error (left) and bed level error (right) in xfor the different formulations at t= 1.5 s in test A 6.2 One dimensional cases 39 0 0.005 0.01 0.015 0.02 0.025 0.03 0.035 0.04 0.045 MPM ASHIDA-MICHIUE ENGELUND-FREDSOE FERNANDEZ LUQUE-VAN BEEK PARKER SMART NIELSEN WONG (4.93) WONG (3.97) CAMENEN-LARSON SMART CFBS RMSE (m) 0 0.005 0.01 0.015 0.02 MPM ASHIDA-MICHIUE ENGELUND-FREDSOE FERNANDEZ LUQUE-VAN BEEK PARKER SMART NIELSEN WONG (4.93) WONG (3.97) CAMENEN-LARSON SMART CFBS RMSE (m) Figure 6.4: RMSE for water level surface (left) and bed level surface (right) with different formulas at t= 1.5 s in test A 6.2.3 Test B Test B is a case of advance front over dry bed and adverse discontinuity. The flow evolves leading to a left moving rarefaction wave ending in front wave dominated by friction. Figure 6.5 shows how front wave celerity is well reproduced in time when using the Smart CFBS formula. Figure 6.6 shows the numerical results and experimental data for the dam break test case B when using MPM (left) and Smart CFBS (right). In both cases, the most relevant difference with measured data is observed over the step, due to the lack of erosion with respect to experimental data. Upstream and downstream the step both numerical simulations provide identical results, being able to reproduce accurately the free surface level in space. The lack of precision over the upward step is observed for all discharge formulations if observing Figure 6.7, that provides level errors in space. The rest of the domain presents an acceptable error. The RMSE for water level surface (left) and bed level surface (right) at t= 1.5 s plotted at Figure 6.8 shows that in this test case there is not clearly a more advantageous formula. 40 CJ scheme: numerical results -0.1 0 0.1 0.2 0.3 0.4 -1 0 1 2 3 h+z (m) x (m) -0.1 0 0.1 0.2 0.3 0.4 -1 0 1 2 3 h+z (m) x (m) -0.1 0 0.1 0.2 0.3 0.4 -1 0 1 2 3 h+z (m) x (m) -0.1 0 0.1 0.2 0.3 0.4 -1 0 1 2 3 h+z (m) x (m) -0.1 0 0.1 0.2 0.3 0.4 -1 0 1 2 3 h+z (m) x (m) -0.1 0 0.1 0.2 0.3 0.4 -1 0 1 2 3 h+z (m) x (m) Figure 6.5: Numerical results and experimental data for the dam break test case B at times t= 0.025, 0.050, 0.075, 0.100, 0.125 and 1.5 s, using a variable value of Ag computed using Smart CFBS: measured water level surface (− • −), measured bed level surface (−◦−), computed water level surface (−4−), computed bed level surface (−N−) 6.2 One dimensional cases 41 -0.1 0 0.1 0.2 0.3 0.4 -1 0 1 2 3 h+z (m) x (m) -0.1 0 0.1 0.2 0.3 0.4 -1 0 1 2 3 h+z (m) x (m) Figure 6.6: Numerical results and experimental data for the dam break test B at t= 1.5 s, using a variable value of Agcomputed using MPM (left) and Smart CFBS (right): measured water level surface (− • −), measured bed level surface (− ◦−), computed water level surface (−4−), computed bed level surface (−N−) 0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 1.2 Error (m) x (m) MPM ASHIDA-MICHIUE ENGELUND-FREDSOE FERNANDEZ LUQUE-VAN BEEK PARKER SMART NIELSEN WONG (4.93) WONG (3.97) CAMENEN-LARSON SMART CFBS 0 0.005 0.01 0.015 0.02 0.025 0.03 0.035 0.04 0.045 -0.4 -0.2 0 0.2 0.4 Error (m) x (m) MPM ASHIDA-MICHIUE ENGELUND-FREDSOE FERNANDEZ LUQUE-VAN BEEK PARKER SMART NIELSEN WONG (4.93) WONG (3.97) CAMENEN-LARSON SMART CFBS Figure 6.7: Modulus of the water level surface error (left) and bed level error (right) in xfor the different formulations at t= 1.5 s in test B 0 0.005 0.01 0.015 0.02 0.025 0.03 0.035 0.04 0.045 MPM ASHIDA-MICHIUE ENGELUND-FREDSOE FERNANDEZ LUQUE-VAN BEEK PARKER SMART NIELSEN WONG (4.93) WONG (3.97) CAMENEN-LARSON SMART CFBS RMSE (m) 0 0.005 0.01 0.015 0.02 MPM ASHIDA-MICHIUE ENGELUND-FREDSOE FERNANDEZ LUQUE-VAN BEEK PARKER SMART NIELSEN WONG (4.93) WONG (3.97) CAMENEN-LARSON SMART CFBS RMSE (m) Figure 6.8: RMSE for water level surface (left) and bed level surface (right) with different formulas at t= 1.5 s in test B 42 CJ scheme: numerical results 6.2.4 Test D Test D represents a reservoir partially filled with sediments and includes a downward step. In this case, once flow passes through the gate location accelerates and decelerates in the friction dominated front. Figure 6.9 shows numerical results and experimental data for the dam break using MPM (left) and Smart CFBS (right). Smart CFBS formulation is able to handle perfectly with this kind of bed discontinuity, tracking the water level surface and redrawing correctly the bed level. Different time instants captured in Figure 6.10 allow appreciating the accuracy and the grade of detail of the computational results in time. Free surface and bed levels are correctly captured for both rarefaction wave and advance front wave, as well as, the bed level at the discontinuity point. Figure 6.11 shows how Smart CFBS formulation provides the lowest level for bed level (right) and free surface (left) error in space at t= 1.5 s if compared with the rest of formulations. Also, the RMSE for water level surface (left) and bed level surface (right) displayed in Figure 6.12 confirms that Smart CFBS formulation gives the better results. Compared with test cases A and B, error is drastically reduced with the proposed formulation. -0.1 0 0.1 0.2 0.3 0.4 -1 0 1 2 3 h+z (m) x (m) -0.1 0 0.1 0.2 0.3 0.4 -1 0 1 2 3 h+z (m) x (m) Figure 6.9: Numerical results and experimental data for the dam break test case D at t= 1.5 s, using a variable value of Agcomputed using MPM (left) and Smart CFBS (right): measured water level surface (− • −), measured bed level surface (− ◦ −), computed water level surface (−4−), computed bed level surface (−N−) 6.2 One dimensional cases 43 -0.1 0 0.1 0.2 0.3 0.4 -1 0 1 2 3 h+z (m) x (m) -0.1 0 0.1 0.2 0.3 0.4 -1 0 1 2 3 h+z (m) x (m) -0.1 0 0.1 0.2 0.3 0.4 -1 0 1 2 3 h+z (m) x (m) -0.1 0 0.1 0.2 0.3 0.4 -1 0 1 2 3 h+z (m) x (m) -0.1 0 0.1 0.2 0.3 0.4 -1 0 1 2 3 h+z (m) x (m) -0.1 0 0.1 0.2 0.3 0.4 -1 0 1 2 3 h+z (m) x (m) Figure 6.10: Numerical results and experimental data for the dam break test case D at times t= 0.025, 0.050, 0.075, 0.100, 0.125 and 1.5 s, using a variable value of Ag computed using Smart CFBS: measured water level surface (− • −), measured bed level surface (−◦−), computed water level surface (−4−), computed bed level surface (−N−) 44 CJ scheme: numerical results 0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 1.2 Error (m) x (m) MPM ASHIDA-MICHIUE ENGELUND-FREDSOE FERNANDEZ LUQUE-VAN BEEK PARKER SMART NIELSEN WONG (4.93) WONG (3.97) CAMENEN-LARSON SMART CFBS 0 0.005 0.01 0.015 0.02 0.025 0.03 0.035 0.04 0.045 -0.4 -0.2 0 0.2 0.4 Error (m) x (m) MPM ASHIDA-MICHIUE ENGELUND-FREDSOE FERNANDEZ LUQUE-VAN BEEK PARKER SMART NIELSEN WONG (4.93) WONG (3.97) CAMENEN-LARSON SMART CFBS Figure 6.11: Modulus of the water level surface error (left) and bed level error (right) in xfor the different formulations at t= 1.5 s in test D 0 0.005 0.01 0.015 0.02 0.025 0.03 0.035 0.04 0.045 MPM ASHIDA-MICHIUE ENGELUND-FREDSOE FERNANDEZ LUQUE-VAN BEEK PARKER SMART NIELSEN WONG (4.93) WONG (3.97) CAMENEN-LARSON SMART CFBS RMSE (m) 0 0.005 0.01 0.015 0.02 MPM ASHIDA-MICHIUE ENGELUND-FREDSOE FERNANDEZ LUQUE-VAN BEEK PARKER SMART NIELSEN WONG (4.93) WONG (3.97) CAMENEN-LARSON SMART CFBS RMSE (m) Figure 6.12: RMSE for water level surface (left) and bed level surface (right) with different formulas at t= 1.5 s in test D 6.2 One dimensional cases 51 (a) 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0 20 40 60 80 100 120 Error (m) Time (s) MPM ASHIDA-MICHIUE ENGELUND-FREDSOE FERNANDEZ LUQUE-VAN BEEK PARKER SMART NIELSEN WONG (4.93) WONG (3.97) CAMENEN-LARSON SMART CFBS (b) 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0 20 40 60 80 100 120 Error (m) Time (s) MPM ASHIDA-MICHIUE ENGELUND-FREDSOE FERNANDEZ LUQUE-VAN BEEK PARKER SMART NIELSEN WONG (4.93) WONG (3.97) CAMENEN-LARSON SMART CFBS (c) 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0 20 40 60 80 100 120 Error (m) Time (s) MPM ASHIDA-MICHIUE ENGELUND-FREDSOE FERNANDEZ LUQUE-VAN BEEK PARKER SMART NIELSEN WONG (4.93) WONG (3.97) CAMENEN-LARSON SMART CFBS (d) 0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 MPM ASHIDA-MICHIUE ENGELUND-FREDSOE FERNANDEZ LUQUE-VAN BEEK PARKER SMART NIELSEN WONG (4.93) WONG (3.97) CAMENEN-LARSON SMART CFBS RMSE (m) Probe 2 Probe 3 Probe 4 Figure 6.21: Modulus of bed level error in time at (a) station SA, (b) station SB and (c) station SC with different formulas. RMSE for bed level zwith different formulas in time (d) 52 CJ scheme: numerical results 6.2.7 Sand cube The last experiment studied in this paper is a test where the flow has a subcritical regime in opposition to previous tests. The experiment was made in a 15 m long channel, with a cross section of 0.5 x 0.5 m2, at the Hydraulics Laboratory of the Civil Engineering School of the University of A Coru˜na (Spain) (Pe˜na et al.,2008). The bottom of the flume was characterized by uniform slope, 0.00052, and a sediment layer 4.5 cm height, was placed in the central part, between 4.5 and 9 m from its upstream end. A sketch can be appreciated in Figure 6.22. The sand employed had the following properties: ρs= 2680kgm−3,d50 = 1mm (uniform size), ϕ= 30o, negligible cohesion, porosity p= 0.5 and was characterized by a Manning roughness factor n= 0.015 sm−1/3. Initial conditions used were a water surface level downstream set to 0.115 m and a flow value enforced to be 21.8 l/s. Numerical simulations were performed using cells ∆x= 0.05 m and CFL = 1. 0.115 m 4.5 m 4.5 m 0.045 m 15 m S0= 0.00052 Figure 6.22: Sand cube sketch Figure 6.23 shows experimental data and numerical results calculated using MPM (left) and Smart CFBS (right) at different times. The bed evolution in time is well described with Smart CFBS. In the first part of the simulation there is an important mobilization of material up to time t= 40 min, when the sediment bed tends to stabilize. Most relevant differences between numerical and experimental data appear downstream the cube. This difference is more noticeable at time t= 120 min and is attributable to the fact that in the sediment transport model suspended load is not considered. A careful data analysis of the measured bed level reveals that at this time, the initial mass associated to the cube is not conserved, may be due to suspension effects. On the other hand, the numerical scheme used in this work is exactly mass conservative, so differences between numerical and experimental data downstream the cube are expectable. The results provided by MPM formula are unable to gather information correctly, leading to a poor bed level prediction as time increases. Correct performance of Smart CFBS in comparison with the rest of sediment discharge formulae is well appreciated in Figure 6.24, where the modulus of bed level error in x 6.2 One dimensional cases 53 (left) and RMSE for bed level surface at time t= 120 min (right) are plotted. Smart CFBS presents the more accurate results. (a) 0 0.01 0.02 0.03 0.04 0.05 0.06 0 2 4 6 8 10 12 14 16 Bed level (m) x (m) (b) 0 0.01 0.02 0.03 0.04 0.05 0.06 0 2 4 6 8 10 12 14 16 Bed level (m) x (m) (c) 0 0.01 0.02 0.03 0.04 0.05 0.06 0 2 4 6 8 10 12 14 16 Bed level (m) x (m) (d) 0 0.01 0.02 0.03 0.04 0.05 0.06 0 2 4 6 8 10 12 14 16 Bed level (m) x (m) (e) 0 0.01 0.02 0.03 0.04 0.05 0.06 0 2 4 6 8 10 12 14 16 Bed level (m) x (m) (f) 0 0.01 0.02 0.03 0.04 0.05 0.06 0 2 4 6 8 10 12 14 16 Bed level (m) x (m) Figure 6.23: Results for the sand cube test case. Initial bed level (···), measured bed and water level (−• −) and computed (−4−) using MPM at times (a) t= 10 min, (c) t= 40 min, (e) t= 120 min, and using Smart CFBS at times (b) t= 10 min,(d) t = 40 min, (f) t= 120 min 54 CJ scheme: numerical results 0 0.005 0.01 0.015 0.02 0.025 0.03 0.035 0.04 4 5 6 7 8 9 10 11 Error (m) x (m) MPM ASHIDA-MICHIUE ENGELUND-FREDSOE FERNANDEZ LUQUE-VAN BEEK PARKER SMART NIELSEN WONG (4.93) WONG (3.97) CAMENEN-LARSON SMART CFBS 0 0.002 0.004 0.006 0.008 0.01 0.012 0.014 MPM ASHIDA-MICHIUE ENGELUND-FREDSOE FERNANDEZ LUQUE-VAN BEEK PARKER SMART NIELSEN WONG (4.93) WONG (3.97) CAMENEN-LARSON SMART CFBS RMSE (m) Figure 6.24: Results for the sand cube test case. Modulus of the bed level error in x for the different formulations after 120 min (left) and RMSE for bed level surface with different formulae at 120 min (right) 6.3 Two dimensional cases 6.3.1 2D Numerical modeling of dam failure The test case studied above, in section 6.2.6, is reproduced in a 2D mesh. Being the flow mostly onedimensional in this case, it is important to check the performance of the numerical discretization of the empirical formulations in a 2D mesh to ensure that numerical results are not influenced by the grid definition. This case is of great interest, as it allows a direct comparison between 1D and 2D simulations in a wide variety of flow conditions. 2D numerical simulations have been performed using a coarse unstructured triangular mesh, with a maximum cell size of 0.01m2, Figure 6.25. The CFL is retained equal to 0.5. Figure 6.26 displays the numerical results obtained using Smart CFBS formulation for both the water level and the bed level. During the first seconds the erosion rate reduces drastically the height of the crest and downstream the dam a hydraulic jump appears. At the final stage of the simulation, a large wedge has been developed. Also, the presence of incipient antidunes is observed. Figures 6.27 (a) and (b) show the water and bed level surface computed after 120s using MPM and Smart CFBS formulations respectively. The bed level evolution recorded in time at the three stations SA, SB and SC, located downstream from the edge of the original dam crest, are plotted in Figures 6.27 (c) and (d). The evolution of the measured and computed water reservoir level is depicted in Figures 6.27 (e) and (f). In all cases the Smart CFBS formulation presents accurate results, while the MPM formulation shows noticeable discrepancies with respect to the experimental data. Figures 6.28 (a) and (b) display the measured and computed overtopping discharge just upstream the breach using MPM and Smart CFBS formulations respectively. It 6.3 Two dimensional cases 55 Figure 6.25: Detail of the triangular mesh (a) 0.052 0.156 0.259 0.363 0.467 0.571 0.674 0.778 0.05 0.15 0.25 0.35 0.45 0.55 0.65 0.75 (b) 0.052 0.156 0.259 0.363 0.467 0.571 0.674 0.778 0.05 0.15 0.25 0.35 0.45 0.55 0.65 0.75 Figure 6.26: Numerical results of water level (top image) and bed level (bottom level) in the dike at 0s (a) and 120s (b) using Smart CFBS formulation. is observed that the experimental overtopping discharge is better tracked with Smart CFBS while MPM predictions are quite far from experimental data. The relative performance of the different formulations in terms of RMSE is plotted in Figure 6.29 at the three stations SA, SB and SC, showing important differences among numerical results depending of the experimental law selected. The Engelund and Fredsoe sediment transport relation was derived for a wide range of slopes, and Figure 6.29 shows how this formulation leads to low values of RMSE. The Smart formula was derived for a set of experimental cases with steep slopes, therefore it can be expected that in this case any numerical discretization would provide accurate predictions. Contrarily, numerical simulation shows that the Smart FS discretization leads to less accurate results if compared with those given by the Smart CFBS discretization. The rest of formulations, derived from experiments ranging from low to medium slopes provide higher RMSE. When comparing the numerical results of the 2D simulation with those obtained of a 56 CJ scheme: numerical results (a) 0 0.2 0.4 0.6 0.8 1 0 2 4 6 8 10 12 14 Elevation (m) x (m) (b) 0 0.2 0.4 0.6 0.8 1 0 2 4 6 8 10 12 14 Elevation (m) x (m) (c) 0.4 0.5 0.6 0.7 0.8 0.9 1 0 20 40 60 80 100 120 140 z (m) Time (seconds) (d) 0.4 0.5 0.6 0.7 0.8 0.9 1 0 20 40 60 80 100 120 140 z (m) Time (seconds) (e) 0 0.2 0.4 0.6 0.8 1 0 20 40 60 80 100 120 140 Water level (m) Time (seconds) (f) 0 0.2 0.4 0.6 0.8 1 0 20 40 60 80 100 120 140 Water level (m) Time (seconds) Figure 6.27: Initial bed level (- - -), computed water level surface (−4−) and bed level surface (−N−) at t= 120 s using (a) MPM and (b) Smart CFBS. Bed level surface evolution in time measured at stations SA (−◦−) (−−),SB (−•−), and SC (−4−) and computed at stations SA (−?−),SB (−−), and SC (−−) using (c) MPM and (d) Smart CFBS. Evolution in time of the measured water reservoir level (−◦−) and computed water reservoir level (−•−) using (e) MPM and (f) Smart CFBS. (a) 0 10 20 30 40 50 0 20 40 60 80 100 120 140 Q (l/s) Time (seconds) (b) 0 10 20 30 40 50 0 20 40 60 80 100 120 140 Q (l/s) Time (seconds) Figure 6.28: Evolution in time of the measured (−◦−) and computed (•−) overtopping discharge using (a) MPM and (b) Smart CFBS. 6.3 Two dimensional cases 57 0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 MPM ASHIDA-MICHIUE ENGELUND-FREDSOE FERNANDEZ LUQUE-VAN BEEK PARKER SMART NIELSEN WONG (4.93) WONG (3.97) CAMENEN-LARSON SMART CFBS RMSE (m) Probe 2 Probe 3 Probe 4 Figure 6.29: RMSE for bed level zat stations SA, SB and SC (−4−) with different formulas in time. 1D discretization, section 6.2.6, it can be observed that the RMSE is slightly bigger in the 2D cases and that 2D results follow closely the tendencies given by the 1D formulation. 6.3.2 Symmetric configuration for 2D dam break flow over erodible bed This experiment was designed at the laboratory of UCL (Soares-Frazao et al.,2012) consisting of a dam break over a 3.6 m wide and about 36 m long flume. The gate was connected to an upstream reservoir and was 1 m wide. The sand was extended over 9 m downstream the gate and 1 m upstream the gate, having a thickness of 0.085 m. A complete sketch of the set up of the experiment is shown in Figure 6.30. The properties of the sand were ρs= 2630 kg m−3,d50 = 1.61 mm, ϕ= 30o, negligible cohesion, porosity p= 0.40 and was characterized by a Manning roughness factor n= 0.019 sm−1/3. Initial conditions used were: upstream, the water level was imposed to 0.047 m, and downstream, a control section at the end of the flume with the same height as the sand layer, 0.085 m. The measurements carried out during the experiments consisted of recording the water level evolution for the first 20 s at different probes, Figure 6.31, and the longitudinal bed profiles measured from x= 0.5 m to x= 8 m at two ycoordinates, Table 6.2, and at t= 100 s. The domain was discretized on a non-uniform triangular mesh, with a higher density downstream the widening, being the total number of cells equal to 12500. The CFL used was imposed to 0.5. Figures in 6.32 show a sequence of plant views of the computed bed evolution in time predicted by the Smart CFBS discretization, characterized by fast morphodynamic changes. Figure 6.32 (a) at t= 10 s shows how the flow generates a wavefront which 58 CJ scheme: numerical results 0.47 m 0.085 m Sediment bed Clean water z x Gate y x 1.8 m 1 m 1.8 m 10 m 1 m 9 m 1 m A A C C B B S1 S2 S3 3.6 m 0.34 m 0.155 m 1 m Figure 6.30: Experimental set up: transversal sketch, 2D sketch and cross sections (AA-CC and BB) e e e e e e e e Gate 0.64 m 1.94 m 0.66 m 0.33 m 0.33 m 0.335 m 0.33 m 0.335 m U1 U2 U3 U4 U5 U6 U7 U8 Figure 6.31: Position of probes in the experiment 6.3 Two dimensional cases 59 Section Y coordinate (m) S1 0.20 S2 0.70 Table 6.2: Position of the sections causes an important erosion process in the enlargement zone of the channel. While the flooding wave advances the sand particles grabbed in this process are carried out to the wavefront and to the wall, where they tend to sediment, as shown in Figure 6.32 (b) at t= 20s respectively. Symmetric elongated sedimentary bodies appear on the right and left banks of the channel, that grow in time to merge generating a diamondshaped erosion region at t= 40s, shown in Figure 6.32 (c). At t= 60 s most of the morphodynamic changes have taken place, and the drainage of the water contained in the upstream reservoir smooths the bed surface, attenuating the bed forms previously generated. For longer times, no more important morphodynamic changes happen. At t= 100 s, Figure 6.32 (f) shows how only the diamond-shaped erosion region in the enlargement zone, generated by the sudden change in flow direction after the opening of the gate, remains in time. The rest of the bed surface becomes almost planar. Figure 6.33 displays the final bed surface at t= 100 s obtained with the experimental data (left) and with the numerical results using Smart CFBS formula (right). Numerical results follow correctly the tendency of the final bed morphology although they tend to underestimate the length of the diamond-shaped body and the thickness of the eroded layer, resulting in smaller heights for the deposition forms. In the experimental data the length of the bed-form zones is bigger than the one provided by the numerical simulation. This may be explained, if considering that, due to the underestimation of erosion rates along the numerical simulation, the magnitude of the bed forms is smaller, and consequently, they are more easily eroded. Also, differences between numerical and experimental bed surfaces can be justified by two important points: i) the 2D SW model neglects the vertical accelerations and decreases the erosion/deposition rate and ii) errors associated to the reconstruction of the experimental bed surface, which was generated through the interpolation of measured bed profiles. The results shown in, Figures 6.34,6.35, display the experimental bed level against the computed one using the MPM and the Smart CFBS formulae at the two control sections. The first one, section S1, which is placed to study the effect of the flow over the bottom in the enlargement zone presents differences between both load discharge formulae. The Smart CFBS formula obtains a better tracking of the sedimentary process, getting more accurate results for the maximum erosion position, x= 1.4 m, and in the maximum deposition position, x= 2.6 m. At the second control section, section S2, differences are also noticeable between both sediment transport formulae, being the Smart CFBS the formula which achieves a better averaged bed level. The computed results obtained with MPM show a zone at 60 CJ scheme: numerical results (a) 0 5 X Y Z (b) 0 5 X Y Z (c) 0 5 X Y Z (d) 0 5 X Y Z (e) 0 5 X Y Z (f) 0 5 X Y Z 0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 Bottom (m) Figure 6.32: Numerical results of bed level in the enlargement zone at 10s (a), 20s (b), 40s (c), 60s (d), 80s (e) and 100s (f) using Smart CFBS formula 0 5 X Y Z 0 5 X Y Z 0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 Bottom (m) Figure 6.33: Experimental results (left) and numerical results using Smart CFBS formula (right) of bed level in the enlargement zone at 100s 6.3 Two dimensional cases 67 Computed results at probe U1, which is placed within the channel, where the flow is mostly one dimensional, provided accurate results with respect the experimental data. Numerical simulations at probe U3, located closer to the enlargement zone, where erosion is of maximum importance, reproduce less accurately the measured water surface level if compared with the rest of probes. Numerical results for probes U6 and U7 located downstream the widening zone lead to accurate predictions of water level surface. Numerical predictions using MPM, Smart, Wong (3.97), Fernandez Luque and Van Beek and Smart CFBS obtain closer results to the experimental data. Smart formula provides less accurate results than the Smart CFBS one, although it achieves in tracking the general trend of temporal evolution. The measured bed level after the dam break event and the numerical predictions at cross sections S1, S3, S5 and S7, S9 are plotted in Figures 6.45 (left) and 6.46 (left). The RMSE obtained with every sediment transport discharge formula at cross sections S1, S3, S5 and S7, S9 are plotted in Figures 6.45 (right) and 6.46 (right) respectively. All sediment transport formulations are able to describe the deposition of material on the left bank and the erosion on the right bank. More noticeable differences appear among them for the predicted bed level at the right bank, where deposition processes take place. On the left bank (y=0, looking upstream) of section S1, located close to the widening zone, all sediment transport formulations predict a bed profile that follows closely the pattern given by the experimental data. Smart CFBS, Wong (4.93) and Wong (3.97) formulae provide the most accurate bed elevations levels. Regarding the right bank (y=0.5), all formulations generate a less sharp slope than the one given by the experiments and Wong (4.93) and Wong (3.97) also obtained the most accurate results. Sections S3 and S5 show that the numerical results track correctly the bed level surface for both left and right banks, giving similar results and RMSE values. The Smart CFBS formula provides the least error. Section S7 shows that on the left bank the level of erosion is well captured with independence of the formulae. Noticeable differences among sediment discharge formulae appear in the stagnation flow region, located at the right wall, where Smart CFBS obtained a better prediction for the bed slope shape. At Section S9, Figure 6.46, which is the cross section placed farthest from the enlargement location, numerical results present the lowest values of RMSE. The bed level on the left bank is newly well tracked by all the formulations but the key zone close to the right bank is only well predicted by Smart CFBS and Wong (4.93) and Wong (3.97). 68 CJ scheme: numerical results (a) 3456 X Y Z (b) 3456 X Y Z (c) 3456 X Y Z (d) 3456 X Y Z 0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 Bottom (m) Figure 6.43: Numerical results of water level (top image) and bed level (bottom level) in the enlargement zone at 2s (a), 3s (b), 5s (c) and 20s (d) (a) 0.08 0.1 0.12 0.14 0.16 0.18 0.2 0.22 0 2 4 6 8 10 Water level (m) t (s) EXPERIMENTAL MPM ASHIDA-MICHIUE ENGELUND-FREDSOE FERNANDEZ LUQUE-VAN BEEK PARKER SMART NIELSEN WONG (4.93) WONG (3.97) CAMENEN-LARSON SMART CFBS (b) 0.08 0.1 0.12 0.14 0.16 0.18 0.2 0.22 0 2 4 6 8 10 Water level (m) t (s) EXPERIMENTAL MPM ASHIDA-MICHIUE ENGELUND-FREDSOE FERNANDEZ LUQUE-VAN BEEK PARKER SMART NIELSEN WONG (4.93) WONG (3.97) CAMENEN-LARSON SMART CFBS (c) 0.08 0.1 0.12 0.14 0.16 0.18 0.2 0.22 0 2 4 6 8 10 Water level (m) t (s) EXPERIMENTAL MPM ASHIDA-MICHIUE ENGELUND-FREDSOE FERNANDEZ LUQUE-VAN BEEK PARKER SMART NIELSEN WONG (4.93) WONG (3.97) CAMENEN-LARSON SMART CFBS (d) 0.08 0.1 0.12 0.14 0.16 0.18 0.2 0.22 0 2 4 6 8 10 Water level (m) t (s) EXPERIMENTAL MPM ASHIDA-MICHIUE ENGELUND-FREDSOE FERNANDEZ LUQUE-VAN BEEK PARKER SMART NIELSEN WONG (4.93) WONG (3.97) CAMENEN-LARSON SMART CFBS Figure 6.44: Numerical results and experimental data of water level for probes U1 (a), U3 (b), U6 (c), U7 (d) 6.3 Two dimensional cases 69 (a) 0.06 0.08 0.1 0.12 0.14 0.16 0.18 0.2 0 0.1 0.2 0.3 0.4 0.5 Bed level (m) y (m) EXPERIMENTAL MPM ASHIDA-MICHIUE ENGELUND-FREDSOE FERNANDEZ LUQUE-VAN BEEK PARKER SMART NIELSEN WONG (4.93) WONG (3.97) CAMENEN-LARSON SMART CFBS 0 0.005 0.01 0.015 0.02 0.025 0.03 0.035 0.04 MPM ASHIDA-MICHIUE ENGELUND-FREDSOE FERNANDEZ LUQUE-VAN BEEK PARKER SMART NIELSEN WONG (4.93) WONG (3.97) CAMENEN-LARSON SMART CFBS RMSE (m) (b) 0.06 0.08 0.1 0.12 0.14 0.16 0.18 0.2 0 0.1 0.2 0.3 0.4 0.5 Bed level (m) y (m) EXPERIMENTAL MPM ASHIDA-MICHIUE ENGELUND-FREDSOE FERNANDEZ LUQUE-VAN BEEK PARKER SMART NIELSEN WONG (4.93) WONG (3.97) CAMENEN-LARSON SMART CFBS 0 0.005 0.01 0.015 0.02 0.025 0.03 0.035 0.04 MPM ASHIDA-MICHIUE ENGELUND-FREDSOE FERNANDEZ LUQUE-VAN BEEK PARKER SMART NIELSEN WONG (4.93) WONG (3.97) CAMENEN-LARSON SMART CFBS RMSE (m) (c) 0.06 0.08 0.1 0.12 0.14 0.16 0.18 0.2 0 0.1 0.2 0.3 0.4 0.5 Bed level (m) y (m) EXPERIMENTAL MPM ASHIDA-MICHIUE ENGELUND-FREDSOE FERNANDEZ LUQUE-VAN BEEK PARKER SMART NIELSEN WONG (4.93) WONG (3.97) CAMENEN-LARSON SMART CFBS 0 0.005 0.01 0.015 0.02 0.025 0.03 0.035 0.04 MPM ASHIDA-MICHIUE ENGELUND-FREDSOE FERNANDEZ LUQUE-VAN BEEK PARKER SMART NIELSEN WONG (4.93) WONG (3.97) CAMENEN-LARSON SMART CFBS RMSE (m) (d) 0.06 0.08 0.1 0.12 0.14 0.16 0.18 0.2 0 0.1 0.2 0.3 0.4 0.5 Bed level (m) y (m) EXPERIMENTAL MPM ASHIDA-MICHIUE ENGELUND-FREDSOE FERNANDEZ LUQUE-VAN BEEK PARKER SMART NIELSEN WONG (4.93) WONG (3.97) CAMENEN-LARSON SMART CFBS 0 0.005 0.01 0.015 0.02 0.025 0.03 0.035 0.04 MPM ASHIDA-MICHUE ENGELUND-FREDSOE FERNANDEZ LUQUE-VAN BEEK PARKER SMART NIELSEN WONG (4.93) WONG (3.97) CAMENEN-LARSON SMART CFBS RMSE (m) Figure 6.45: Numerical results and experimental data of bed level for sections S1 (a), S3 (b), S5 (c), S7 (d), and its corresponding RMSE obtained with every sediment transport formula 70 CJ scheme: numerical results 0.06 0.08 0.1 0.12 0.14 0.16 0.18 0.2 0 0.1 0.2 0.3 0.4 0.5 Bed level (m) y (m) EXPERIMENTAL MPM ASHIDA-MICHIUE ENGELUND-FREDSOE FERNANDEZ LUQUE-VAN BEEK PARKER SMART NIELSEN WONG (4.93) WONG (3.97) CAMENEN-LARSON SMART CFBS 0 0.005 0.01 0.015 0.02 0.025 0.03 0.035 0.04 MPM ASHIDA-MICHIUE ENGELUND-FREDSOE FERNANDEZ LUQUE-VAN BEEK PARKER SMART NIELSEN WONG (4.93) WONG (3.97) CAMENEN-LARSON SMART CFBS RMSE (m) Figure 6.46: Numerical results and experimental data of bed level for section S9, and its corresponding RMSE obtained with every sediment transport formula . Chapter 7 CJ scheme: conclusions Several well-known sediment discharge formulae have been studied and included in a general form in a coupled-Jacobian model for the shallow water equations and the Exner equation. Additionally, the Smart formula, that includes the bed slope, has been formulated in both in 1D, 2D configurations and extended to distinguish among different situations. In special, the possibility to evaluate situations where the flow encounters an adverse slope, and consequently has a less erosion capacity as well as other situations where the flow reaches a favorable slope and has a bigger erosion rate has been included. This discretization has been called Smart CFBS (Combined Friction and Bed Slope). 1D test cases In the first set of test cases, the bed load formulae have been applied to solve dam break flows over dry/wet initial conditions. Whilst advance front celerity has been well captured in the dam break cases over dry bed with independence of the type of capacity formula used, noticeable differences appear in the bed level predictions, except in the dam break test case A with initially flat bed level and in the case B with adverse slope. In experiment D, over favorable slope and dry bed, erosion produces a meaningful variation of the initial bed step, leading to a rate of erosion and deposition only well captured in time and space if using Smart CFBS formula. In test case F, where both sides are initially filled with water, and a favorable slope is present, only Smart CFBS formula leads to a correct erosion evolution in time, that becomes negligible in the final stage of the experiment. From this set of test cases it can be concluded that Smart CFBS formula can be recommended for dam break test cases with null, adverse and favorable slopes, and wet/dry problems. When numerically modeling dam erosion and failure it has been found that Smart CFBS formula is applicable in all cases analyzed in this study. Also, Engelund and Fredsoe capacity formula provides correct results in bed level predictions, although Smart CFBS formula estimates much better the maximum discharge values reached in all experiments. It is also worth mentioning that, for downstream steep slopes, the 72 CJ scheme: conclusions computational time associated to the peak discharge value is calculated earlier. The computed sand cube test showed that the best agreement between experimental and numerical data are obtained with Smart (computing slope as friction slope) and Smart CFBS. This can be explained considering that, in this case, unsteady hydrodynamic effects are a quasi-steady process of slowly varying bed-load, and friction slope is adapted to bed slope. 2D test cases For the first experiment, the dike failure by overtopping, characterized by a onedimensional flow, it was clearly stated that the Smart CFBS formulation provides the most accurate results in time and in space. In the second experiment, a symmetric dam break over a mobile bed in a channel with an enlargement zone was numerically reproduced. In this case a two-dimensional flow is generated and differences among different sediment formulations are less noticeable. Numerical results follow the tendency of the final bed morphology, underestimating the length of the diamond-shaped body and the thickness of the eroded layer. In the third experiment, where an erodible channel with a sudden enlargement produces a two-dimensional flow, the computational results provided good agreement with experimental values for the different sediment formulae. Comparing with previous results from other authors (Spinewine and Zech,2004;Abderrezzak and Paquier,2011;Wu and Wang,2007,2008) it can be stated that the numerical scheme used in this work allows to clarify the differences among different formulations which were derived by 1D stationary laboratory experiments. The Smart CFBS discretization reaches the more accurate results in all cases, although in a genuinely 2D flow, that is, a situation involving more than one flow direction, the differences between sediment transport formulae are not as noticeable as in the 1D situations. The results and conclusions of this part of the document have been published in Juez et al. (2013b). 7.1 Further research Numerical experimentation is necessary to include non-equilibrium state formulation in the mathematical model. Hence, it is necessary to develop the one layer model derived of mass conservation equations, section 3, including a non uniform density along the longitudinal profile. Several authors have suggested to include the difference between the actual transported material and the equilibrium sediment transport capacity by means of the definition of an adaptation length. 7.1 Further research 73 The inclusion of the suspended transport coupled with the bed load model developed here is other natural follow-up of this work. For this phenomenon, both mathematical and numerical model are still to be studied. Another important feature that should be addressed is the numerical assessment of the the bed load sediment transport discharge formulae under a bed composed of a non-uniform sand particles. The interaction between grains of different diameter requires an additional effort. The study of the mathematical and numerical properties of more complex friction laws for the definition of the shear stress at the bottom is also necessary. This feature is oriented to the definition of a hyperconcentrated model, where the rheology of the flow presents a pseudo plastic behavior due to high values of depth averaged concentrations. Chapter 8 Weakly-coupled numerical scheme 8.1 Introduction The development of the novel numerical strategy proposed for coupling the hydrodynamic and the morphodynamic models is described in this Chapter. The weaklycoupled numerical scheme will be noted as WC from now on. The numerical scheme departs from the previous system of equations presented in the one layer model for the shallow water (3.41) and for the Exner model (8.8) which are written separately as follows: Hydrodynamic model Neglecting diffusion of momentum due to viscosity and turbulence, wind effects and the Coriolis term, the two-dimensional shallow water equations are formulated as, ∂U ∂t +∂F(U) ∂x +∂G(U) ∂y =S(U, x, y) (8.1) where U= (h, qx, qy)T(8.2) are the conserved variables with hrepresenting the water depth, qx=hu and qy=hv, with (u, v) the depth averaged components of the velocity vector ualong the (x, y) coordinates respectively. The fluxes of these variables are given by: F=qx,q2 y h+1 2gh2,qxqy hT ,G=qy,qxqy h,q2 y h+1 2gh2T (8.3) where gis the acceleration of the gravity. The source terms of the system are 76 Weakly-coupled numerical scheme S=0,pb,x ρw−τb,x ρw ,pb,y ρw−τb,y ρwT (8.4) which express the x-component and y-component of: i) the pressure force along the bottom line, pb,x and pb,y, being ρwthe water density, that in differential form are expressed as a function of the bed slope, So pbx ρw=ghSo,x So,x =−∂z ∂x pby ρw=ghSo,y So,y =−∂z ∂y (8.5) and ii) the bed shear-stress, τb,x and τb,y. System (8.1) is time dependent, non-linear, and is non-homogeneous due to the presence of source-terms. The pure shallow water model is hyperbolic since the eigenvalues of its Jacobian matrices are always real. The presence of the source-terms leads to a non-strictly hyperbolic system. However, it is assumed that under the hypothesis of dominant advection it can be classified and numerically dealt with as belonging to the family of hyperbolic systems. Hence, the mathematical properties of (8.1) include the existence of a Jacobian matrix, Jn, of the flux normal to a direction given by the unit vector, n,En=Fnx+Gny, defined as Jn=∂En ∂U=∂F ∂Unx+∂G ∂Uny(8.6) whose components are Jn=  0nxny (gzh−u2)nx−uvnyvny+ 2unxuny (gzh−v2)ny−uvnxvnxunx+ 2vny (8.7) The eigenvalues of this Jacobian matrix (λ1=un −c,λ2=un and λ3=un +c, with c=√gh) constitute the wave speeds in the linearized problem and provide information about directions in which the information travels. Morphodynamic model Sediment dynamics are assumed to be well modeled through the bed-load Exner equation, ∂z ∂t +ξ∂qs,x ∂x +ξ∂qs,y ∂y = 0 (8.8) where zis the bed elevation, ξ=1 1−p,pis the material porosity, qs,x and qs,y denote the solid transport discharge along the (x, y) coordinates respectively, influenced by the water depth hand the depth averaged velocities uand v. 8.3 Approximate Riemann Solution for the Hydrodynamic model 83 e Λk=  e λ10 0 0e λ20 0 0 e λ3  k (8.38) Difference in vector Uacross the grid edge is projected onto the matrix eigenvectors basis δUk=e PkAk(8.39) where Ak= ( α1α2α3)T kcontains the set of wave strengths. Furthermore, for linking the source terms to the set of eigenvalues they are also projected onto the matrix eigenvectors basis (¯ Sn)k= (e PB)k(8.40) with Bk= (β1, β2, β3)T k. Using (12.28), (12.48) and (12.56) matrix δMkcan be expressed as δMk=e Jn,kδUk−e Pk(B)k= 3 X m=1 e λ θαeem k(8.41) with θm k=1−β e λαm k (8.42) or in matrix form δMk= (e Pe ΛΘe P−1)kδUk(8.43) Therefore the value for the desired matrix Ln,k in (12.47) is Ln,k = (e Pe ΛΘe P−1)k(8.44) where Θis a diagonal matrix Θk=  θ10 0 0θ20 0 0 θ3 k (8.45) 84 Weakly-coupled numerical scheme that relates fluxes and source terms and that becomes equal to the identity matrix in absence of source terms. 2D first order finite volume for the Hydrodynamic Model Definition of matrix Ln,k allows to define directly right-going and left-going wave propagations. This way the flux δMkis written in splitting way as follows δMk=δM− i,k +δM+ j,k (8.46) with δM− i,k = (e Pe Λ−Θe P−1)kδUkδM+ j,k = (e Pe Λ+Θe P−1)kδUk(8.47) the flux splitting version of the Godunov first order method is Un+1 i=Un i− NE X k=1 δM− i,k ∆t lk Ai (8.48) Superindex −becomes necessary to distinguish from outcoming fluxes to cell iat edge k, that will be referred to as δM+ j,k, as they update the adjacent jcell sharing the k edge. 8.4 Approximate Riemann Solution for the Morphodynamic model Coming back to the equation which governs the sediment dynamics (8.8), the same steps followed with system (8.1) can be applied. A local 1D RP is obtained projecting the sediment fluxes onto the normal direction nkof each kedge of each cell ∂z ∂t +ξ∂(qsn) ∂x0= 0 (8.49) Using the integral form of (8.49) the weak solutions of the RP can be found. For this purpose a suitable control volume, Figure 8.4, is integrated over the following time interval [0,∆t] and the space interval [−∆x0,∆x0], with x0sufficiently large, Z+∆x0 −∆x0 z(x0, t = ∆t)dx0= ∆x0(zi+zj)−ξδqsn∆t(8.50) 8.4 Approximate Riemann Solution for the Morphodynamic model 85 t x0 ∆t(qsn)i(qsn)j ∆x0∆x0 zizj z(t > 0) x0=0 x0 zi zj λb Figure 8.4: Integration control volume defined by a time interval [0,∆t] and a space interval [−∆x0,∆x0] Again, the piecewise representation of the variables is hypothesized and the first order Godunov method is used for updating the averaged quantities. Consistency condition for the Morphodynamic Model Following the philosophy employed for the hydrodynamic model a Roe approach is going to be used, i.e., the solution of each RP is obtained from the exact solution of a locally linearized problem defined by an approximate solution ˆz(x, t). This constant linear problem is based on the definition of an approximate wave speed of the non-linear sediment flux, qsn. This following equivalent equation is written ∂ˆz ∂t +e λbn,k ∂ˆz ∂x0= 0 (8.51) with the following initial conditions ˆz(x0,0) = ziif x0<0 zjif x0>0(8.52) The approximate solution must fulfill the Consistency Condition (Leveque,2002), forcing the integral of the exact solution (8.49) and the integral of the locally linearized solution, (8.51) to be the same. Thanks to this constraint it is possible to obtain the following expression for the wave speed which updates the bed level, 86 Weakly-coupled numerical scheme e λbn,k =δ(ξqsn,k) δz (8.53) with δz =zj−ziand δqsn,k =qsn,j −qsn,i. Regarding equation (4.2) it is necessary to compute the Grass coefficient for defining the bed load discharge in each cell. Following Murillo and Garc´ıa-Navarro (2010a) as the coefficient Agis not a constant but varies from cell to cell, at every edge ka local Ag,k value is defined as an arithmetic mean between neighboring cells. Consequently, the term δqsn,k is written as δqsn,k =Ag,kδunk. Additionally, when applying numerical modeling techniques under a flat bottom situation, the bed level difference is null and consequently the bed wave speed is not defined. In order to overcome this difficulty the computation of the friction slope, Sf,k Murillo and Garc´ıa-Navarro (2010b), is proposed. The friction slope is commonly used in a high number of sediment transport empirical laws, Meyer-Peter and M¨uller (1948); Smart (1984); Ashida and Michiue (1972); Camenen and Larson (2005), as these formulae were derived from 1D steady solid transport experiments. Additionally, its employment is coherent with the fact that transport process implies a loss of energy through the interaction between the sediment and the flow (Smart,1984;Whittaker and Davies,1982). Also it is worth noting that the linearization of e λbn,k in cases of almost flat bottom can lead to unphysical huge values of the bed wave speed. This is avoided by imposing a lower threshold for the bed level difference between cells: up to grain size, ds, the approximation of the friction slope will be considered. This limitation ensures coherent values in the estimation of the bed wave speed, and wave celerity in (8.53) is approximated by e λbn,k =ξδqsn,k δz0(8.54) with δz0=δz if δz0> ds −Sf,kdnif δz0< ds (8.55) being dnthe normal distance between cell centers (Murillo and Garc´ıa-Navarro,2010b). 2D first order finite volume for the Morphodynamic Model The evaluation of the wave speed, e λbn,k as in (8.54), brings the opportunity of splitting the sediment flux difference δqsn,k in right-going and left-going wave propagations. Consequently the Godunov first order method can be written as 8.4 Approximate Riemann Solution for the Morphodynamic model 87 δqsn,k =δqsn+ i,k +δqsn− j,k (8.56) with δqsn+ i,k =e λ+ bn,kδzkδqsn− j,k =e λ− bn,kδzk(8.57) and e λ± bn,k =1 2(e λbn,k ±|e λbn,k|). Therefore, zn+1 i=zn i− NE X k=1 δqsn− i,k ∆t lk Ai− NE X k=1 δqsnIi,k ∆t lk Ai (8.58) where the second term of the right side in (8.58) evaluates the flux in the cell edge and the third term completes the updating formula to consider the spatial variation of Ag, as it was justified in Murillo and Garc´ıa-Navarro (2010a). Another possibility for defining the Godunov first order method is through a flux scheme, considering outcoming and incoming fluxes through the edges of the cell. Hence the bed level is updated as zn+1 i=zn i− NE X k=1 ξq∗ sn,k ∆t lk Ai (8.59) where q∗ sn,k =(qsn,i if e λbn,k >0 qsn,j if e λbn,k <0(8.60) being qsn,i and qsn,j the bed load discharge computed in the cell iand in the cell j. Although both numerical schemes (8.58) and (8.59) are completely equivalent, it must be stressed that the flux version is computationally more efficient, as minor algebraic operations are need. Additionally, with the flux form of the numerical scheme in (8.59), ghost cells must be considered in the boundary cells in order to complete the information over the entire cell, Leveque (2002). It is worth noting that the application of ghost cells almost does not penalize the computational effort. In this fashion, since the computational cost when using the flux scheme in (8.59) is less, this alternative has been chosen for obtaining the results displayed in the next sections. 88 Weakly-coupled numerical scheme 8.5 Stability region Updated values of Un+1 iand zn+1 iare defined cell averaging the contributions of the local RPs, and in consequence the time step ∆thas to be taken small enough so that there is no interaction of waves from the kneighboring RPs. In the 2D framework, considering unstructured meshes, the relevant distance, that will be referred to as χi in each cell imust consider the volume of the cell and the length of the shared kedges (Murillo and Garc´ıa-Navarro,2010b) χi=Ai maxk=1,NE lk (8.61) Considering that each kRP is used to deliver information to a pair of neighboring cells of different size, the distance min(Ai, Aj)/lkis relevant, so in case that the water depth is greater than zero in all the regions of the RP solution the time step is limited by ∆t≤CFL ∆t e λ∆t e λ=min(χi, χj) max |e λm|(8.62) with CFL=1 in case of 1D meshes, CFL=1/2 in case of 2D structured or unstructured meshes (Toro,1997) and being e λmthe wave speeds. When the advection structure of the problem is all contained in the system matrices, i.e. coupled-Jacobian approach (Murillo and Garc´ıa-Navarro,2010a;Castro Diaz et al., 2009;Soares-Frazao and Zech,2010;Siviglia et al.,2013), the linearised wave speeds provided by the eigenvalues allow to define a suitable CFL condition, retaining the sediment transport part of the system. However, when using uncoupled/asynchronous (De Vriend et al.,1993) or coupled/synchronous models (Holly and Rahuel,1990;Cao et al.,2002;Wu and Wang,2004;Xia et al.,2010), it has been considered traditionally that since the wave speeds associated to water surface and bed level present different magnitudes, not straightforward limitation has to be considered in the stability condition. Nevertheless, this is no longer admissible when the celerities are in the same order of magnitude. Therefore, an extra limitation linked to the bed wave speed is required ∆t≤CFL ∆t e λ∆t e λ=min(χi, χj) |e λm,e λb|(8.63) 8.6 Geomorphological collapse The same strategy proposed in 5.4 for modeling the geotechnical equilibrium bank characteristics is employed here: a simple mass conservative mechanism of slope sliding 8.6 Geomorphological collapse 89 failure, assuming that the angle of repose of submerged material of the bed can be approximated by the friction angle. Chapter 9 Weakly-coupled scheme: results 9.1 Introduction This Chapter gathers the validation tests that allow to show the assessment of the numerical schemes described in Chapter 8. Numerical results have been compared with experimental data and exact solution considering 1D and 2D situations. The bed-load discharge law employed for computing the bed evolution is the Smart CFBS, which was introduced in Chapter 5. Furthermore, in all the simulations a conservative mechanism of slope sliding failure has been considered, 5.4. The experimental tests employed for comparing with the numerical results are the same as the ones proposed in Chapter 6. Therefore the detailed description is omitted and the reference to each particular section from Chapter 6is given. 9.2 Problems with exact solutions In this section the numerical solutions for three bidimensional test cases that will be named A, B and C, as summarized in Table 9.1 are presented. The tests are Riemann problems for the movable bed equations, in which the friction shear-stress has been neglected in the momentum equation. The two first test cases, A and B, have been chosen to assess the performance of the numerical scheme assuming that morphological changes can be characterized using constant values of Ag. The third one, Test C, assumes that the value of Agdepends on the water depth. These exact solutions were firstly reported in Murillo and Garc´ıa-Navarro (2010a). The exact solutions were built by nesting several waves, departing from a left state until reaching to define the right state. The CFL condition is equal to 1.0, the mesh size is x= 0.1mand the simulation is computed up to t= 2s. The value of the parameter Agfor the Grass law is considered as 92 Weakly-coupled scheme: results Ag=Ag,o hr(9.1) being Ag,o = 0.01 in all cases, r= 0 in test cases A and B, and r= 1 in test case C. Test hLhRuLuRvLvRzLzR A 2.0 2.0 0.25 2.3247449 0.05 0.04 3.0 2.846848 B 2.25 1.18868612 0.20 2.4321238 0.045 0.02 5.0 5.124685 C 6.0 5.2 0.3 15.167196 0.015 0.04 3.0 4.631165 Table 9.1: Summary of dam break test cases with exact solution In order to compare the accuracy of the weakly-coupled technique (WC) proposed in this work, the results obtained with the coupled-Jacobian technique used in Murillo and Garc´ıa-Navarro (2010a) (CJ) are also plotted. TestA: the solution proposed in this test case is based on two outcoming rarefactionwaves and a central shock together with a contact wave evolving downstream, Figure 9.1. The shock and the contact wave move slowly, compared with the other two waves. The shock absorbs most of the initial step in bed profile. The numerical solution is able to capture the general trend of the flow behavior, without arising numerical problems at the step area. The unit sediment discharge in both directions is also displayed. TestB: the second solution analyzed is built through two-rarefaction waves, a contact wave and a shock, Figure 9.2. he central wave is a quite slowly-moving rarefaction, bearing most of the initial step of the bed profile, whereas the shock is very weak, almost ineffective for the bed. The computed results are able to depict the moving waves in all the wet domain with an adequate level of accuracy. 9.3 One dimensional cases 99 -0.1 0 0.1 0.2 0.3 0.4 -1 0 1 2 3 h+z (m) x(m) z Experimental h+z Experimental z Numerical h+z Numerical -0.1 0 0.1 0.2 0.3 0.4 -1 0 1 2 3 h+z (m) x(m) z Experimental h+z Experimental z Numerical h+z Numerical -0.1 0 0.1 0.2 0.3 0.4 -1 0 1 2 3 h+z (m) x(m) z Experimental h+z Experimental z Numerical h+z Numerical -0.1 0 0.1 0.2 0.3 0.4 -1 0 1 2 3 h+z (m) x(m) z Experimental h+z Experimental z Numerical h+z Numerical -0.1 0 0.1 0.2 0.3 0.4 -1 0 1 2 3 h+z (m) x(m) z Experimental h+z Experimental z Numerical h+z Numerical -0.1 0 0.1 0.2 0.3 0.4 -1 0 1 2 3 h+z (m) x(m) z Experimental h+z Experimental z Numerical h+z Numerical Figure 9.7: Numerical results and experimental data for the dam break test case F at times t= 0.25, 0.50, 0.75, 1.0, 1.25 and 1.5 s, using a variable value of Agcomputed using Smart CFBS: measured water level surface (−•−), measured bed level surface (−◦−), computed water level surface (−4−), measured bed level surface (−N−) 100 Weakly-coupled scheme: results 0.0000 0.0005 0.0010 0.0015 0.0020 0.0025 0.0030 0.0035 0 0.2 0.4 0.6 0.8 1 1.2 1.4 ∆t (s) t (s) ∆t Exner ∆t SWE Figure 9.8: Time step evolution in test case F for the water waves speed as in (8.62), (−•−), and for the bed wave speed as in (8.63), (−◦−) during time simulation -0.1 0 0.1 0.2 0.3 0.4 -1 0 1 2 3 h+z (m) x(m) z Experimental h+z Experimental z Numerical h+z Numerical -0.1 0 0.1 0.2 0.3 0.4 -1 0 1 2 3 h+z (m) x(m) z Experimental h+z Experimental z Numerical h+z Numerical Figure 9.9: Numerical results and experimental data for the dam break test case F at times t= 1.0 and 1.5 s, when CFL limitation related to the bed speed is removed 9.3 One dimensional cases 101 9.3.2 1D Knickpoint test case Morphological changes due to the transition between two planes with different slope (knickpoint) were measured in (Bellal et al.,2004). This test case is useful to compare the capacity of the numerical schemes to handle with a sudden flow transition from subcritical regime over a mild slope to supercritical regime over a steep slope. A sketch of the experiment, with the initial conditions of bed slope, is shown in Figure 9.10. The knickpoint is defined as the point of abrupt change in the longitudinal bottom profile of the channel. This experiment was carried out using a coarse and uniform size sand with the following properties ρs= 2680kgm−3,d50 = 1.65mm,ϕ= 30o, negligible cohesion, porosity p= 0.42 and was characterized by a Manning roughness factor n= 0.0165 sm−1/3. Initial conditions employed are: upstream, water level surface (0.028 m) and discharge (9.8 l/s); downstream, a known water surface level at the end of the flume (0.11 m). The domain, 7.4 meters long, is divided using ∆x= 0.05 m. In all simulations CFL = 1. S01 = 0.0057 S02 = 0.024 0.0274 m0.063 m 0.001 m 0.109 m0.109 m 6.3 m 1.1 m Figure 9.10: Knickpoint sketch Bed level variation in the longitudinal profile was recorded in time and is compared with the predictions supplied by the numerical schemes in Figure 9.11. The computed solution describes a good trend when comparing with the experimental solution. The erosion located in the knickpoint is predicted at the same rate as the experiment and the final bottom is also well achieved. Since in this experimental case an important change in the bottom morphology takes place, Figure 9.12 shows the more restrictive time step associated to the wave speeds of water and bed in time simulation. Bed time step imposes a harder restriction than the fluid flow and for this reason has to be considered in (8.62) for preserving the numerical stability of the numerical scheme. 102 Weakly-coupled scheme: results 0.1 0.12 0.14 0.16 0.18 0.2 0.22 0 1 2 3 4 5 6 7 Levels (m) x (m) Initial h+z Numerical z Numerical h+z Experimental z Experimental 0.1 0.12 0.14 0.16 0.18 0.2 0.22 0 1 2 3 4 5 6 7 Levels (m) x (m) Initial h+z Numerical z Numerical h+z Experimental z Experimental 0.1 0.12 0.14 0.16 0.18 0.2 0.22 0 1 2 3 4 5 6 7 Levels (m) x (m) Initial h+z Numerical z Numerical h+z Experimental z Experimental 0.1 0.12 0.14 0.16 0.18 0.2 0.22 0 1 2 3 4 5 6 7 Levels (m) x (m) Initial h+z Numerical z Numerical h+z Experimental z Experimental 0.1 0.12 0.14 0.16 0.18 0.2 0.22 0 1 2 3 4 5 6 7 Levels (m) x (m) Initial h+z Numerical z Numerical h+z Experimental z Experimental Figure 9.11: Results for the knickpoint test case. Initial bed level (···), measured bed and water level (−•−) and computed (−4−) at times t= 165, 223, 345, 589 and 851 s with variable value of Agcomputed using Smart CFBS 9.4 Two dimensional cases 103 0.0000 0.0050 0.0100 0.0150 0.0200 0.0250 0.0300 0.0350 0 100 200 300 400 500 600 700 800 ∆t (s) t (s) ∆t Exner ∆t SWE Figure 9.12: Time step evolution for the water waves speed, as in (8.62), (−•−), and for the bed wave speed as in (8.63), (−◦−) during time simulation 9.4 Two dimensional cases 9.4.1 2D Numerical modeling of dam failure This experiment has been previously defined in 6.2.6 and it was studied by Tingsanchali and Chinnarasri (2001). Following prior work developed in Juez et al. (2013b) the 2D numerical simulation has been performed using a coarse unstructured triangular mesh, with a maximum cell size of 0.01m2. The mesh together with the initial water depth is displayed in Figure 9.13. CFL is imposed equal to 0.5. Free boundary condition is considered at the outflow section. Figure 9.14 displays the bed level evolution when using Smart CFBS formulation. At the crest of the dike strong erosion occurred because of the strong initial discontinuity of water depth and the severe slope downwards the gate. The granular material of the dike is completely mobilized within a short period of time and it is grabbed downstream the dam by the flow. Figure 9.15(a) shows the water and bed level surface computed after 120 s when using Smart CFBS formulation. As the bed level was recorded in time at three stations SA, SB and SC, located downstream from the edge of the original dam crest, the comparison between experimental data and computed results are displayed in Figure 9.15(b). Numerical results are able to handle the strong morphodynamics changes which take place without displaying numerical oscillations and additionally, well tracking the experimental data. On the other hand, the evolution of the measured and computed water reservoir level is depicted in Figure 9.15(c). Figure 9.15(d) displays the measured and computed overtopping discharge just upstream the breach. Both measurements provide high quality and useful information about this type of phenomena. Numerical schemes allows to obtain a good detail of forecasting capacity for the bed and water level evolution together with an efficient computational cost. 104 Weakly-coupled scheme: results Figure 9.13: Detail of the triangular mesh and initial condition for the water depth Figure 9.14: Computed results of the bed level evolution when using a variable value of Agbuilt with Smart CFBS and at times t= 0, 30, 80 and 120 s 9.4 Two dimensional cases 105 (a) 0 0.2 0.4 0.6 0.8 1 0 2 4 6 8 10 12 14 h+z (m) x (m) Initial z Numerical h+z Numerical (b) 0.4 0.5 0.6 0.7 0.8 0.9 1 0 20 40 60 80 100 120 140 z (m) Time (seconds) (c) 0 0.2 0.4 0.6 0.8 1 0 20 40 60 80 100 120 140 Water level (m) Time (seconds) (d) 0 10 20 30 40 50 0 20 40 60 80 100 120 140 Q (l/s) Time (seconds) Figure 9.15: (a) Initial bed level (- - -), computed water level surface (−4−) and bed level surface (−N−) at t= 120 s. (b) Bed level surface evolution in time measured at stations SA (− ◦ −) (−−), SB (− • −), and SC (−4−) and computed at stations SA (−?−),SB (−−), and SC (−−). (c) Evolution in time of the measured water reservoir level (− ◦ −) and computed water reservoir level (−• −). (d) Evolution in time of the measured (−◦−) and computed (−•−) overtopping discharge For this test case, the time step evolution associated to each wave speed is also studied, Figure 9.16. Initially, heavier restrictions are required by the water flow, as the overtopping event has not provoked yet the dike failure. However, as time advances and the geomorphic changes become more severe, time step restrictions come from the bed celerity. At the end of time simulation, where most of the sediment particle movement has occurred, the time step is newly governed by flow characteristics. In view of these results, it is proved the efficiency of the solver, as only when important bed changes exist the classical time step of water flow is decreased. Additionally to the study of the time step evolution this test case has been chosen also for comparing the computational time cost with respect to the coupled-Jacobian technique used in Murillo and Garc´ıa-Navarro (2010a) (CJ) and the weakly-coupled scheme (WC) proposed in this work. For this purpose three meshes with increasing number of elements are considered. In Table 9.3 are displayed the ratio between the computational cost when employing Murillo and Garc´ıa-Navarro (2010a) and when considering the procedure explained in this work. Results plotted above belongs to the second mesh. Noticeable computational efficiency is achieved, being more important as the level of mesh refinement is increased. The computational cost time with the 106 Weakly-coupled scheme: results 0.0005 0.0010 0.0015 0.0020 0.0025 0.0030 0.0035 0 20 40 60 80 100 120 ∆t (S) t (s) ∆t Exner ∆t SWE Figure 9.16: Time step evolution for the water waves speed, as in (8.62), (−•−), and for the bed wave speed, as in (8.63), (−◦−) during time simulation CJ is penalized by the high number of algebraic operations need for computing the eigenvalues and eigenvectors. In order to support this fact and employing the second mesh, the time step evolution, associated to the CJ and to WC is displayed in Figure 9.17. Despite of presenting a bigger time step on average when using the CJ, the computational cost is higher. 0.0005 0.0010 0.0015 0.0020 0.0025 0.0030 0.0035 0 20 40 60 80 100 120 ∆t (s) t (s) ∆t CJ ∆t WC Figure 9.17: Time step evolution following the CJ technique in Murillo and Garc´ıaNavarro (2010a), (−◦−), and the WC technique explained in this work, (−•−) during time simulation Together with the computational cost time, the RMSE (Root median square error) for the three stations SA, SB and SC obtained when using the CJ and the WC technique, is displayed in Table 9.4. The weakly-coupled technique provides computational results close to the experimental ones whilst the computational time is decreased. 9.4 Two dimensional cases 107 N. of elements Ratio of computational cost time =CJ/WC 2000 8.46 4100 10.15 8300 13.72 Table 9.3: Summary of ratios of computational cost time when using the JCM technique and the WC technique N. of elements RMSE(m)SA RMSE(m)SB RMSE(m)SC CJ WC CJ W C CJ W C 2000 0.065 0.039 0.042 0.034 0.058 0.037 4100 0.043 0.021 0.028 0.019 0.038 0.023 8300 0.028 0.014 0.019 0.012 0.025 0.015 Table 9.4: Summary of the RMSE associated to each station when using the JC technique and the WC technique 9.4.2 2D Dam break with a sudden enlargement This experiment was performed at the laboratory of the Civil and Environmental Engineering Department of the UCL (Palumbo et al.,2008;Gouti`ere et al.,2011) and has been previously numerically reproduced in 6. During the development of the experiment the water level evolution was recorded at different points as well as the final bed surface at several cross sections, Figure 9.18. An unstructured mesh is considered and CFL condition is imposed equal to 0.5. Figure 9.18: Plan view of the experimental flume. Locations of the probes (left) and the cross sections (right) This experimental case represents a complete challenge as it gathers several highlighted situations which can occur in the real engineering life: an area where the flow is genuinely one-dimensional, an abrupt expansion which provokes the change to a twodimensional flow, important velocity gradients which create a recirculating area, moving shocks close to the wall zone and moreover a severe local erosion together with a noticeable sediment deposition area. It constitutes the perfect benchmark for checking the assessment of the numerical schemes against sudden and strong changes in the flow and the bed. Due to these characteristics other authors have also studied recently this test case (Soares-Frazao and Zech,2010;Xia et al.,2010;Siviglia et al.,2013). 108 Weakly-coupled scheme: results A series of computed bed surface evolutions are shown in Figure 9.19. The bed deformation is very sensitive since the flow evolves over a initially dry bed: sediment particles start to bounce as soon as the water reaches their position creating a kind of ripples or dunes, at time t= 2 s. As the water overtakes the corner of the channel the flow expands, causing the water depth to decrease and the bed level suffers a dramatic local erosion, at time t= 4 s. Close to the wall area the flow tends to slow down and the material grabbed upstream is settled. In this zone of the channel the loss of energy is so strong that a bed sharp surface emerges, at times t= 4 and 6 s. Downstream, the sediment grains are pushed outward the domain and eventually intersects driving to settling zones, times t= 6 and 10 s. At the last time, t= 20 s, the drainage of water leads to soften the bed surface although the minimum and maximum sediment peak areas are clearly identified. It is worth noting that the bed ripples plotted have a twofold nature. Firstly they have a numerical origin, since they are generated by the mesh topology as the numerical technique employed in this work make use of the Riemann theory, which is built considering the edges of each cell. Additionally, the bed ripples have also a physical nature, as the flow evolves over a dry bottom. Once the experiment has been qualitatively described, computed and experimental data are faced. Comparison between the water level measured and the numerical solution is showed in Figure 9.20. The majority of the probes achieve a good trend in relation with the experimental data. Probes U3 and U4 are the ones which provide less accurate results. This is justified by the fact that they are located close to the expansion (probe U3) and close to the wall (probe U4), where three dimensional flow structures are generated due to the sudden expansion and the shock against the lateral side. With the present mathematical model, where the set of equations is depth-averaged, the vertical accelerations are neglected and consequently, this flow behavior cannot be properly treated (Xia et al.,2010). Figure 9.21 gathers the measured bed level after the dam break event and the numerical predictions at control sections S1, S2, S3, S4 and S5. In all the sections the computed bed surface is able to follow the measured evolution. Section S1 which is the closest to the expansion does not obtain neither the maximum nor the minimum of sediment peaks, although the prediction follows the sediment movement pattern: particles are grabbed from left and settled to the right bank. In control sections, S2, S3 and S4, the computed bed surface follows correctly the tendency of the final bed morphology although the final bed slopes are less sharp than the ones recorded after the experiment. As it has been noted before, since the mathematical model is depth-averaged the vertical accelerations are not considered. Consequently, the erosion/deposition rates are decreased and differences in the granular material lying close to the right wall are expected. Section S5, positioned far away from the area of stronger influence, obtains a good tendency when comparing with the experimental data. Part II Mass motion over steep areas Chapter 11 Introduction Several catastrophic events during the past decades have showed that important flood induced by a dam-break is in relation with a strong erosion in the bed and in the banks. These rapid and variably geomorphodynamic processes affect significantly to the flow behavior compounding the harmful effects of the flooding waves. For this reason, the numerical modeling of severe transient geomorphic flows is an active topic in the research field. In addition, the study of these geomorphic flows/landslides and their movement constitutes an important environmental issue as they play a key role in landscape evolution. Currently, the triggering mechanisms, mechanical properties and assessment of likelihood and consequences as well as the development of measures to limit their impact, is an active topic in the field of the geophysical flows research. As this phenomena involved a mixture of mud, sand and water sliding down a slope together, the study of granular flows constitutes an starting point for the understanding of the more complex mass movement phenomena mentioned before (Denlinger and Iverson,2004). Therefore, several experiments on granular dry flows have been carried out in the past (Savage and Hutter,1989;Iverson and Denlinger,2001;Pouliquen and Forterre,2002;Lajeunesse et al.,2004;Mangeney et al.,2010) as the initial target to be overcame by the numerical modeling tools. Since the computational models open a wide range of possibilities for handling with this type of phenomena a numerical scheme has been developed for the analysis of dry granular flows. Moreover, taking advantage of the numerical experimentation an extra work has been carried out in order to push the knowledge about the dry granular behavior. 118 Introduction 11.1 State of the art of the numerical techniques Granular dry flows show fluid-like behavior where the front of the avalanche moves as a thin layer along high distances. Well known approaches for describing geophysical flows consider the Saint-Venant equations as an starting point. Depth averaged equations were first employed for solving geomorphologic flows by Savage and Hutter (1989), where granular mass sliding was modeled including Coulomb-like basal frictions, and assuming a cohesionless Mohr-Coulomb type material. Since then, new mathematical models have appear in the literature. In Denlinger and Iverson (2004) an extensive and complete review of predictive models for geomorphologic flows was provided and special attention was devoted to the computation of the Coulomb stresses conjugated to the deformation in solid-like behavior avalanches. Contrary to the Saint-Venant equations, defined in a Cartesian coordinate, the Savage-Hutter model uses a curvilinear coordinate along the topography. Denlinger and Iverson (2004) formulated the depth-averaged governing equations referenced to a rectangular Cartesian coordinate system (with Zvertical) and in their new approach the estimated frictional stresses were defined with independence of the orientation of the coordinate system. The model was tested against analytical solutions and experimental data. Bouchut et al. (2003) introduced an extra term in the original Savage-Hutter mathematical model, related to the curvature of the bottom, which is usually neglected when compared in terms of magnitude, in order to ensure the equilibrium at rest of the mass whatever the flow conditions (topography, friction coefficients, etc). Some phenomenas, such as landslides, where the curvature terms play an important role can be found in Favreau et al. (2010); Moretti et al. (2012). In Bouchut and Westdickenberg (2004) this model was extended to consider an arbitrary coordinate system for shallow flow over a 2D topography, retaining the curvature terms. In recent works concerning geomorphologic flows over 2D irregular bed topographies (Pirulli et al.,2007;Pirulli and Mangeney,2008) this term was omitted and promising computational results were obtained. Following Pirulli et al. (2007); Pirulli and Mangeney (2008) curvature terms related with the geometry are not considered here and the rheology of the material will be described using a Coulomb-type friction law. Once a mathematical model is selected, another separate issue is the numerical scheme used. Considering the hyperbolic nature of the depth averaged equations, Godunov type schemes are commonly used in literature (Denlinger and Iverson,2004;MangeneyCastelnau et al.,2003). Godunov type schemes can be constructed departing from the definition of approximate solvers of the Riemann problem (RP). Approximate solvers provide a comprehensible definition of the conserved variables in the inner states of the same RP. Among the most successful and disseminated approximate solvers, Roe’s method (Roe,1986) and the HLL method (Harten et al.,1983), were defined to approximate solutions for hyperbolic system of equations without source terms. When including the presence of source terms in the system of equations, it is possible to extend the numerical schemes defined for the homogeneous case using point-wise explicit or implicit discretizations of the source terms. In Mangeney-Castelnau et al. (2003), internal stresses represented by the Coulomb friction law were discretized us- 11.1 State of the art of the numerical techniques 119 ing a point-wise implicit discretization. Numerical experimentation has shown that point-wise discretizations lead to undesirable results, as non uniform discharge values in steady solutions (Burguete et al.,2008b;Murillo et al.,2008). A proper discretization of the frictional source terms must ensure a correct balance among fluxes and source terms (well-balanced property). Following previous work in well balanced numerical schemes for the shallow flow equations in presence of bed variations (Hubbard, M. E. and Garc´ıa-Navarro, P.,2000), in Denlinger and Iverson (2004), the Roe scheme was applied in combination with an upwind technique applied to the internal stresses. The apparent topography method to deal with generic source terms in Bouchut and Westdickenberg (2004), based on the well-balanced property, was applied successfully to the simulation of the spreading of a granular column over a rough horizontal plane in Mangeney-Castelnau et al. (2005) and over a over an inclined plane in MangeneyCastelnau et al. (2007). Contrary to shallow water flows, where quiescent flow is given in cases of horizontal water level surface, in granular flows, steady state configurations include correct modeling of starting and stopping flow conditions (Bouchut and Westdickenberg,2004;Mangeney-Castelnau et al.,2007;Pirulli et al.,2007;Pirulli and Mangeney,2008). The presence of source terms leads to non-strictly hyperbolic systems of equations and, as a consequence, they have an impact in the solution of the RP. In the shallow water equations with variable topography, different approximations to the Riemann problem have been presented in the literature (Alcrudo and Benkhaldoun,2001;Chinnayya et al.,2004;LeFloch and Thanh,2007;Bernetti et al.,2008;Rosatti and Begnudelli, 2010;LeFloch and Thanh,2011). The properties of these RP solutions not only must guarantee the well-balanced property, but also, the associated numerical scheme must ensure convergence to the solution. Convergence to the solution is not an easy task, as in problems with source terms the total number of waves can be larger than the number of characteristic fields (LeFloch and Thanh,2011). Two augmented solvers which consider intrinsically the presence of source terms, named ARoe (Augmented Roe) and HLLCS (HLL with Contact wave and Source terms), were presented in Murillo and Garc´ıa-Navarro (2010b) and Murillo and Garc´ıa-Navarro (2012b) respectively. Both schemes include an extra static wave for considering the effect of the source terms in the stability region. The ARoe solver was exploited in Murillo and Garc´ıa-Navarro (2012a) allowing correct approximate solutions of wave Riemann problems involving complex rheology when using depth average equations. An accurate and robust first order finite volume scheme, able to handle correctly transient problems including modeling of starting and stopping flow conditions was presented in Murillo and Garc´ıa-Navarro (2012a). Then, in contrast with prior works, Bouchut and Westdickenberg (2004); Pirulli et al. (2007), where numerical fluxes were constructed to ensure well-balanced arguments, in Murillo and Garc´ıa-Navarro (2012a), the definition of the complete approximate solution ensured correct integral estimations of the source terms under all type unsteady of flow conditions. It is worth mentioning, that, in general, only in cases of quiescent equilibrium, the source terms can be integrated exactly. In any other case, the approximate solver provides the rules to avoid unphysical results, allowing the correction of the estimations made for the source terms if necessary. This result is 120 Introduction of utmost importance when modeling of starting and stopping flow conditions and can be applied with independence of the type of rheological model selected. In presence of steep slopes the usual hypothesis of hydrostatic pressure in the vertical direction Zin the shallow water equations is not longer admissible. This fact has consequences when deriving the mathematical model. In global coordinates (X, Y, Z) (Figure 12.1), the gravity vector has a simple form g=−(gX, gY, gZ)T= (0,0,−g)T(11.1) and the bed vector Socan be written as So= (tan θ, tan γ) = −∂Zb ∂X ,−∂Zb ∂Y (11.2) where Zbis the bed level surface in global coordinates. If a coordinate system linked to the topography (x, y, z) is preferred, the gravity vector must be projected following the new system of coordinates. By means of two rotations around the angles γand θ, as in Pirulli (2005), g=−(gx, gy, gz) = −gsin θ, −gsin γcos θ, −gcos θcos γT(11.3) The correct discretization of numerical fluxes and source terms is an important issue in both coordinate systems, when the bed slopes may change from cell to cell. Therefore, in order to preserve steady state configurations, and with independence of the numerical solver selected, it is necessary to ensure exact balance among fluxes and source terms, including the correct modeling of starting and stopping flow conditions. 11.2 State of the art for the experimental works Due to the fact that avalanches are initiated on steep slopes, pioneer experimental studies concerning granular flows were focused on the grain movement over constant inclined planes with slopes larger than the material repose angle Wieland et al. (1999); Pouliquen (1999); Pouliquen and Forterre (2002); Mangeney et al. (2010). This type of movement is governed by the gravity component along the slope direction. Experiments developed in Pouliquen (1999); Mangeney et al. (2010) were performed over a genuine 1D configuration whilst Wieland et al. (1999); Pouliquen and Forterre (2002) were devoted to 2D events. All of them brought the opportunity of studying unstable granular masses, focusing on the maximum spreading or the avalanche front and tail speeds. dditionally to these prior laboratory works, in Lajeunesse et al. (2004); Boutreux and deGennes (1997) other type of configuration was experimentally addressed: the sudden 11.2 State of the art for the experimental works 121 release of a surface over a quasi-horizontal surface. In such case, the weight of the sand grains was the responsible for the onset of the movement, while the frictional forces were in charge of the stopping condition. These experiments, being free from the influence of the topography, were of utmost importance, since they provided results concerning the quantity of mass mobilized by the flow, the final shape and the maximum spreading of the granular mass. Another important configuration which has been recently mimicked in the laboratory consists of granular flows traveling over erodible topography, Mangeney et al. (2010); Roche et al. (2011). This phenomena is easily found in nature, as under certain circumstances landslides can move over deposits built up by earlier events. The strong effects of erosion processes can significantly increase the mobility of avalanches, changing drastically the final distribution of the granular mass, Mangeney-Castelnau et al. (2005); Bouchut et al. (2008); Mangeney et al. (2010). The study of granular flows in combination with obstacles has also acquired prominence during the last years. The impact of the obstacle in the flow behavior needs to be understood for a better design of civil engineering elements such as mast of electrical power lines, buildings, ski lifts, dams and other man-made structures. Several works have dug on this active research field, some of them analyzing the flow overtopping on dike elements Hakonardottir et al. (2003); Faug et al. (2008) and other ones focusing on the shock waves generated by the impact between the flow and the single obstacle Gray et al. (2003); Hakonardottir and Hogg (2005); Hauksson et al. (2007). Following the previous effort made by the authors mentioned above Gray et al. (2003); Hakonardottir et al. (2003); Hakonardottir and Hogg (2005); Hauksson et al. (2007), one of the main concern of this work is in relation with the study of the variable nature of the moving shocks and their complex birth and propagation. Since we want to get closer to the phenomenology which takes place in nature, a series of laboratory experiments have been carried out for studying novel an real-life configuration: 2D spread of the granular mass over variable topography with a changing slope and multiple shock waves derived from the presence of multiple obstacles. The experimental avalanche is triggered by a simple mechanism: a granular mass which is suddenly released from a semi-spherical container. The full description of the experimental facility, methods and cases is found in Caviedes et al. (2014). To provide a physical insight into these phenomena, the spatial and temporal spreading dynamics and the morphology of the resulting shape are investigated and discussed in this work through the wave theory (Roe,1986). For this purpose, the computational scheme developed for dry granular flow in this work as previous goal is the basis for performing the numerical experimentation. Since the computational scheme has been previously tested and validated the numerical results obtained are free of distorting numerical effects allowing to study the physical features involve in the granular flows. 122 Introduction 11.3 Outline In this work, the results presented in Murillo and Garc´ıa-Navarro (2010b,2012a) are extended to provide appropriate numerical schemes for mathematical models of 2D granular flow written in global and local system of coordinates. Taking advantage of this fact a numerical experimentation has been performed through the basis of a novel experimental work. In Chapter 11 weak solutions for Riemann problems over steep and variable slopes are presented focusing in the analysis and definition of the correct numerical discretization of both fluxes and source terms in local coordinates. Also, a detailed definition of the numerical fluxes based on the analysis of the approximate solution is provided. The results are extended to global coordinates in Chapter 12. The augmented approximate solver in Murillo and Garc´ıa-Navarro (2010b,2012a) is modified in both systems of coordinates to consider the effect of bed slope in pressure distribution and frictional effects. In Chapter 13, the numerical solvers are tested against 1D and 2D experimental data in order to check the suitability of the mathematical models described in this work following global and local system of coordinates. Conclusions and further research are written in Chapter 14. Regarding the experimental work, in Chapter 15 it is brought together a summary of the laboratory setup where the experiments have been carried out and in addition, it is also depicted extra consideration about the friction law employed. Chapter 16 is devoted to the comparison between the experimental data and the computed results, focusing on the discussion of the physics involved in the granular flow behavior. Finally, the conclusions of the laboratory work and future research line are written in Chapter 17. Chapter 12 Mathematical model and numerical scheme following local coordinates 12.1 Introduction In this Chapter, a mathematical model for describing granular flow in local coordinates, LC from now on, is presented. It assumes that the flow is oriented in a predominantly longitudinal direction and is confined to a layer which is thin compared to the scale of interest. Hydrostatic pressure distribution in the normal direction to the bed is assumed and a Coulomb type bed friction is used to model the basal stress. The set of depth averaged equations developed in Chapter 3for the one layer model are also valid in this phenomena. Additionally, a numerical scheme is developed regarding the particularities of the problem. 12.2 Mathematical model Bearing in mind the expression for the gravity vector in (11.3), the depth averaged equations expressing volume and momentum conservation are written as follows ∂U ∂t +∂F(U) ∂x +∂G(U) ∂y =Sτ+Sb(12.1) where U=h, hu, hv T(12.2) are the conserved variables, with hrepresenting granular material depth in the zcoordinate and (u, v) the depth averaged components of the velocity vector along x, y 124 Mathematical model and numerical scheme following local coordinates coordinates. The fluxes are given by F=hu, hu2+1 2gzh2, huvT G=hv, huv, hv2+1 2gzh2T (12.3) and the source terms of the system are split in two kind of terms. The term Sτ represents the frictional effects in the bed, and is defined as Sτ=0,−τb,x ρ,−τb,y ρT (12.4) with τb,x, τb,y the bed shear stress in the xand ydirection respectively and ρthe density of the fluid. These tangential forces are evaluated in this work through a Coulomb law Sτ= (0,−gzhtan θb,−gzhtan θb)T(12.5) being θbthe dynamic friction angle between the bed and the flowing mass. The term Sbis defined as Sb= (0,−gxh, −gyh)T(12.6) and expresses the variation of the pressure force in the xand ydirection respectively. Note that no bed derivatives are included. Instead the relevant quantity is the gravity acceleration projection. Figure 12.1 shows a 1D sketch of the relative position of local and global coordinates. g X Z x Zb(X) h(x) z θ Figure 12.1: 1D sketch of global and local coordinates System (12.1) is time dependent, non linear, and contains source terms. Under the hypothesis of dominant advection it can be classified and numerically dealt with as C Bibliography Abderrezzak, K. K. and Paquier, A. 2011. Applicability of Sediment Transport Capacity Formulas to Dam-Break Flows over Movable Beds. Journal of Hydraulic Engineering 137, 209–221. Alcrudo, F. and Benkhaldoun, F. 2001. Exact solutions to the Riemann problem of the shallow water equations with a bottom step. Computers and Fluids 30, 643– 671. Aric` o, C. and Tucciarelli, T. 2008. Diffusive Modeling of Aggradation and Degradation in Artificial Channels. Journal of Hydraulic Engineering 134(8), 1079– 1088. Ashida, K. and Michiue, M. 1972. Study on hydraulic resistance and bedload transport rate in alluvial streams. Transactions, Japan Soc. Civil Eng. 206, 569– 589. Bellal, M.,Iervolino, M.,and Zech, Y. 2004. Knickpoint migration process: experimental and numerical approaches. Proc., 12th Conf. on ”Sediment and Sedimentation Particles”. Prague, Czech Republic , –. Berm´ udez, A. and V´ azquez-Cend´ on, M. E. 1994. Upwind methods for hyperbolic conservation laws with source terms. Computers and Fluids 23, 1049–1071. Bernetti, R.,Titarev, V.,and Toro, E. 2008. Exact solution of the Riemann problem for the shallow water equations with discontinuous bottom geometry. Journal of Computational Physics 227, 3212–3243. Bilanceri, M.,Beux, F.,Elmahi, L.,Guillard, H.,and Salvetti, M. 2012. Linearized implicit time advancing and defect correction applied to sediment transport simulations. Computers and Fluids 63, 82–104. Bouchut, F.,Fern´ andez-Nieto, E. D.,Mangeney, A.,and Lagr´ ee, P. Y. 2008. On new erosion models of Savage-Hutter type for avalanches. Acta Mechanica 199, 1-4, 181–208. Bouchut, F.,Mangeney-Castelnau, A.,Perthame, B.,and Vilotte, J. 2003. A new model of Saint-Venant and Savage-Hutter type for gravity driven shallow water flows. C. R. Acad. Sci. Paris Ser. I, 336, 531–536. 230 BIBLIOGRAPHY Bouchut, F. and Westdickenberg, M. 2004. Gravity driven shallow water models for arbitrary topography. Community of Math Sciences 2, 3, 359–389. Boutreux, T. and deGennes, P. 1997. Evolution of a step in a granular material: the Sinai problem. Comptes rendus de l academie des sciences. Serie II Fascicule B-Mechanique Physique Chimie Astronomie 325, 2, 85–89. Burguete, J.,Garc´ ıa-Navarro, P.,and Murillo, J. 2008a. Friction term discretization and limitation to preserve stability and conservation in the 1D shallowwater model: Application to unsteady irrigation and river flow. International Journal of Numerical Methods in Fluids 54, 403–425. Burguete, J.,Garc´ ıa-Navarro, P.,and Murillo, J. 2008b. Preserving bounded and conservative solutions of transport in one-dimensional shallow-water flow with upwind numerical schemes: Application to fertigation and solute transport in rivers. International Journal of Numerical Methods in Fluids 56, 1731–1764. Camenen, B. and Larson, M. 2005. A general formula for non-cohesive bed load sediment transport. Estuarine, Coastal and Shelf Science 63, 249–260. Canestrelli, A.,Dumbser, M.,Siviglia, A.,and Toro, E. 2010. Well-balanced high-order centred schemes on unstructured meshes for shallow water equations with fixed and mobile bed. Advances in Water Resources 33, 291–303. Cao, Z.,Day, R.,and Egashira, S. 2002. Coupled and decoupled numerical modeling of flow and morphological evolution in alluvial rivers. Journal of Hydraulic Engineering 128, 306–321. Cao, Z.,Pender, G.,and Carling, P. 2006. Shallow water hydrodynamic models for hyperconcentrated sediment-laden flows over erodible bed. Advances in Water Resources 29(4), 546–557. Castro Diaz, M.,Fernandez Nieto, E.,Ferreiro, A.,and Pares, C. 2009. Two-dimensional sediment transport models in shallow water equations. A second order finite volume approach on unstructured meshes. Computer Methods in Applied Mechanics and Engineering 198, 2520–2538. Caviedes, D.,Juez, C.,Murillo, J.,and Garc´ ıa-Navarro, P. 2014. Experimental study of 2D granular flow over a complex topography with obstacles using a consumer-grade RGB-D sensor. Under review in Journal of Fluid Mechanics –, –. Chinnayya, A.,LeRoux, A.,and Seguin, N. 2004. A well-balanced numerical scheme for the approximation of the shallow water equations with topography: the resonance phenomenon. International Journal On Finite Volumes 1, 1–33. Cordier, S.,Le, M.,and Morales de Luna, T. 2011. Bedload transport in shallow water models: Why splitting (may) fail, how hyperbolicity (can) help. Advances in Water Resources 34, 980–989. BIBLIOGRAPHY 231 Cui, X. and Gray, J. 2013. Gravity-driven granular free-surface flow around a circular cylinder. Journal of Fluid Mechanics 720, 314–337. Da Cruz, F.,Emam, S.,Prochnow, M.,Roux, J.,and Chevoir, F. 2005. Rheophysics of dense granular materials: Discrete simulation of plane shear flows. Physical Review 72, 021309. Dal Maso, G.,LeFloch, P.,and Murat, F. 1995. Definition and weak stability of nonconservative products. . Math. Pures Appl. 74, 483–548. De Vriend, H.,Zyserman, J.,Nicholson, J.,Roelvink, J.,Pechon, P.,and Southgate, H. 1993. Medium-term 2DH coastal area modelling. Journal of Coastal Engineering 21, 193–224. Denlinger, R. and Iverson, R. 2004. Granular avalanches across irregular three dimensional terrain: 1. Theory and computation. Journal of Geophysical Research 109, F1, F01014. Douady, S.,Andreotti, B.,and Daerr, A. 1999. On granular surface flow equations. The European Physical Journal B 11, 131–142. Dressler, R. F. 1954. Comparison of theories and experiments for the hydraulic dam-break wave. Int. Assoc. Sci. Hydrology 3, 319–328. Einstein, H. 1950. The bed-load function for sediment transportation in open channel flows. Tech. Rep. Engelund, F. and Fredsoe, J. 1976. Sediment transport model for straight alluvial channels. Nordic Hydrology 7:5, 293–306. Faccanoni, G. and Mangeney, A. 2012. Exact solution for granular flows. International Journal for Numerical and Analytical Methods in Geomechanics 26, –. Faug, T.,Gauer, K.,Lied, K.,and Naaim, M. 2008. Overrun length of avalanches overtopping catching dams: cross-comparison of small-scale laboratory experiments and observations from full-scale avalanches. Journal of Geophysical Research 113, F03009. Favreau, P.,Mangeney, A.,Lucas, A.,Crosta, G.,and Bouchut, F. 2010. Numerical modeling of landquakes. Geophysical Research Letters 37, L1530. Forterre, Y. and Pouliquen, O. 2003. Long-surface-wave instability in dense granular flows. Journal of Fluid Mechanics 486, 21–50. Fraccarollo, L. and Capart, H. 2002. Riemann Wave description of erosional dam-break flows. Journal of Fluid Mechanics 461, 115–133. Garegnani, G.,Rosatti, G.,and Bonaventura, L. 2013. On the range of validity of the Exner-based models for mobile-bed river flow simulations. Journal of Hydraulic Research 51(4), 380–391. 232 BIBLIOGRAPHY Gouti` ere, L.,Soares-Frazao, S.,Savary, C.,Laraichi, T.,and Zech, Y. 2008. One-dimensional model for transient flows involving bed-load sediment transport and changes in flow regimes. Journal of Hydraulic Engineering 134(6), 726–735. Gouti` ere, L.,Soares-Frazao, S.,and Zech, Y. 2011. Dam-break flow on mobile bed in abruptly widening channel: experimental data. Journal of Hydraulic Research 49(3), 367–371. Grass, A. 1981. Sediments transport by waves and currents. SERC London Cent. Mar. Technol, Report No. FL. Gray, J.,Tai, Y.,and Noelle, S. 2003. Shock waves, dead zones and particle-free regions in rapid granular free-surface flows. Journal of Fluid Mechanics 491, 161–81. Gray, J.,Wieland, K.,and Hutter, K. 1999. Gravity-driven free surface flow of granular avalanches over complex basal topography. Proc. Royal Soc. London Ser. A, 455, 1841. Hakonardottir, K. and Hogg, A. 2005. Oblique shocks in rapid granular flows. Physics of Fluids 17, 077101. Hakonardottir, K.,Hogg, A.,Batey, J.,and Woods, A. 2003. Flying avalanches. Geophysical Research Letters 30, 2191. Harten, A.,Lax, P.,and Leer, B. V. 1983. On upstream differencing and Godunov-type schemes for hyperbolic conservation laws. SIAM Review 25:1, 35– 61. Hauksson, S.,Pagliardi, M.,Barbolini, M.,and Johannesson, T. 2007. Laboratory measurements of impact forces of supercritical granular flow against mast-like obstacles. Cold Regions, Science and Technology 49, 54–63. Holly, F. M. and Rahuel, J. L. 1990. New numerical/physical framework for mobile-bed modelling. I: Numerical and physical principles. Journal of Hydraulic Research 28(4), 401–416. Hubbard, M. E. and Garc´ ıa-Navarro, P. 2000. Flux difference splitting and the balancing of source terms and flux gradients. Journal of Computational Physics 165, 89–125. Hudson, J. 2001. Numerical techniques for morphodynamic modelling. Ph.D. thesis, Department of Mathematics, The University of Reading, Whiteknigths, Reading. Hudson, J. and Sweby, P. K. 2002. Formulations for Numerically Approximating Hyperbolic Systems Governing Sediment Transport. Journal of Scientific Computing 19, 225–251. Hudson, J. and Sweby, P. K. 2005. A high-resolution scheme for the equations governing 2D bed-load sediment transport. International Journal of Numerical Methods in Fluids 47, 1085–1091. BIBLIOGRAPHY 233 Iverson, R. and Denlinger, R. 2001. Flow of variably fluidized granular masses across three-dimensional terrain. A Coulomb mixture theory. Journal of Geophysical Research 106, B1, 537–552. Juez, C.,Murillo, J.,and Garc´ ıa-Navarro, P. 2013a. 2D simulation of granular flow over irregular steep slopes using global and local coordinates. Journal of Computational Physics 255, 166–204. Juez, C.,Murillo, J.,and Garc´ ıa-Navarro, P. 2013b. Numerical assesment of bed load discharge formulations for transient flow in 1D and 2D situations. Journal of Hydroinformatics -, In press. Julien, P. 1998. Erosion and Sedimentation. Cambridge University Press. Kalinske, A. 1947. Movement of sediment as bed load in rivers. Trans. AGU 28, 615–620. Kassem, A. A. and Chaudry, M. H. 1998. Comparison of coupled and semicoupled numerical models for alluvial channels. Journal of Hydraulic Engineering 124(8), 794–802. Kerswell, R. 2005. Dam break with Coulomb friction: A model for granular slumping? Physics of Fluids 17. Lajeunesse, E.,Mangeney-Castelnau, A.,and Villote, J. P. 2004. Spreading of a granular mass on a horizontal plane. Physics of Fluids 16, 2371–2381. LeFloch, P. and Thanh, M. 2007. The Riemann problem for shallow water equations with discontinuous topography. Community of Math Sciences 5, 865–885. LeFloch, P. and Thanh, M. 2011. A Godunov-type method for the shallow water equations with discontinuous topography in the resonant regime. Journal of Computational Physics 230, 7631–7660. Leveque, R. 2002. Finite Volume Methods for Hyperbolic Problems. Cambridge University Press, New York. Luque, R. F. and van Beek, R. 1976. Erosion and transport of bedload sediment. Journal of Hydraulic Research 14, 127–144. Lyn, D. and Altinakar, M. 2002. St. Venant-Exner equations for near-critical and transcritical flows. Journal of Hydraulic Engineering 128(6), 579–587. Mangeney, A.,Roche, O.,Hungr, O.,Mangold, N.,Faccanoni, G.,and Lucas, A. 2010. Erosion and mobility in granular collapse over sloping beds. Journal of Geophysical Research 115, F03040. Mangeney, A. and Heinrich, Ph. and Roche, R. 2000. Analytical and numerical solution of the dam-break problem for application to water floods, debris and dense snow avalanches. Pure and Applied Geophysics 157, 1081–1096. 234 BIBLIOGRAPHY Mangeney-Castelnau, A.,Bouchut, F.,Thomas, N.,Vilotte, J. P.,and Bristeau, M. 2007. Numerical modeling of self-channeling granular flows and of their levee-channel deposits . Journal of Geophysical Research 112, F02017. Mangeney-Castelnau, A.,Bouchut, F.,Vilotte, J. P.,Lajeunesse, E., Aubertin, A.,and Pirulli, M. 2005. On the use of Saint-Venant equations for simulating the spreading of a granular mass . Journal of Geophysical Research 110, B09103. Mangeney-Castelnau, A.,Vilotte, J. P.,Bristeau, M. O.,Perthame, B., Bouchut, F.,Simeoni, C.,and Yernini, S. 2003. Numerical modeling of avalanches based on Saint-Venant equations using a kinetic scheme. J. Geophys. Res. 108, 2527. Manning, R. 1895. On the flow of water in open channels and pipes. Transactions of the Institution of Civil Engineers of Ireland 20, 161–207. Meyer-Peter, E. and M¨ uller, R. 1948. In: Report on the 2nd Meeting International Association Hydraulic Structure Research. Stockholm, Sweden. Moretti, L.,Mangeney, A.,Capdeville, Y.,Stutzmann, E.,Huggel, C., Schneider, D.,and Bouchut, F. 2012. Numerical modeling of the Mount Steller landslide flow history and of the generated long period seismic waves. Geophysical Research Letters 39, L16402. Murillo, J. and Garc´ ıa-Navarro, P. 2010a. An Exner-based coupled model for two-dimensional transient flow over erodible bed. Journal of Computational Physics 229, 8704–8732. Murillo, J. and Garc´ ıa-Navarro, P. 2010b. Weak solutions for partial differential equations with source terms: Application to the shallow water equations. Journal of Computational Physics 229, 4327–4368. Murillo, J. and Garc´ ıa-Navarro, P. 2011. Improved Riemann solvers for complex transport in two-dimensional unsteady shallow flow. Journal of Computational Physics 230, 7202–7239. Murillo, J. and Garc´ ıa-Navarro, P. 2012a. A Riemann solver for unsteady computation of 2D shallow flows with variable density. Journal of Computational Physics 231:4, 1963–2001. Murillo, J. and Garc´ ıa-Navarro, P. 2012b. Augmented versions of the HLL and HLLC Riemann solvers including source terms in one and two dimensions for shallow flow applications. Journal of Computational Physics 231, 6861–6906. Murillo, J.,Garc´ ıa-Navarro, P.,and Burguete, J. 2008. Time Step Restrictions For Well Balanced Shallow Water Solutions In Non-Zero Velocity Steady States. International Journal of Numerical Methods in Fluids 56, 661–686. BIBLIOGRAPHY 235 Murillo, J.,Garc´ ıa-Navarro, P.,and Burguete, J. 2009. Conservative Numerical Simulation of Multicomponent Transport in Two-Dimensional Unsteady Shallow Water Flow. Journal of Computational Physics 228, 5539–5573. Nielsen, P. 1992. Coastal Bottom Boundary Layers and Sediment Transport. Advanced Series on Ocean Engineering. World Scientific Publishing. Palumbo, A.,Soares-Frazao, S.,Goutiere, L.,Pianese, D.,and Zech, Y. 2008. Proc., River Flow 2008 International Conference on Fluvial hydraulics, Cesme. Parker, G. 1979. Hydraulic geometry of active gravel rivers. Journal of Hydraulic Engineering 105:9, –. Pe˜ na, E.,Fe, J.,S´ anchez-Tembleque, F.,Puertas, J.,and Cea, L. 2008. Experimental validation of a sediment transport two-dimensional depth-averaged numerical model using PIV and 3D Scanning technologies. Journal of Hydraulic Research Vol. 46, 489–503. Pirulli, M. 2005. Numerical modelling of landslide runout. Ph.D. Degree in Geotechnical Engineering , . Pirulli, M.,Bristeau, M.,Mangeney-Castelnau, A.,and Scavia, C. 2007. The effect of the earth pressure coefficients on the runout of granular material. Environmental Modelling and Software 22, 1437–1454. Pirulli, M. and Mangeney, A. 2008. Results of Back-Analysis of the Propagation of Rock Avalanches as a Function of the Assumed Rheology. Rock Mechanics and Rock Engineering 41, 59–84. Pouliquen, O. 1999. Scaling laws in granular flows down rough inclined planes. Physics of Fluids 11, 542–548. Pouliquen, O. and Forterre, Y. 2002. Friction law for dense granular flows: application to the motion of a mass down a rough inclined plane. Journal of Fluid Mechanics 453, 133–151. Pouliquen, O. and Forterre, Y. 2008. Flows of Dense Granular Media. Annual Review of Fluid Mechanics 40, 1–24. Ritter, A. 1892. Die Fortpflanzung der Wasserwelle. Vereine Deutscher Ingenieure Zeitswchrift 36, 947–954. Roche, O.,Attali, M.,Mangeney, A.,and Lucas, A. 2011. On the run-out distance of geophysical gravitational flows: Insight from fluidized granular collapse experiments. Earth and Planetary Science Letters 311, 3, 375–385. Roe, P. 1986. Numerical Methods in Fluid Dynamics. Vol II. Oxford University Press, Oxford. 236 BIBLIOGRAPHY Rosatti, G. and Begnudelli, L. 2010. The Riemann Problem for the onedimensional, free-surface Shallow Water Equations with a bed step: theoretical analysis and numerical simulations. Journal of Computational Physics 229, 760–787. Rosatti, G.,Murillo, J.,and Fraccarollo, L. 2007. Generalized Roe schemes for 1D two-phase, free-surface flows over a mobile bed. Journal of Computational Physics 54, 543–590. Rosatti, G.,Murillo, J.,and Fraccarollo, L. 2008a. Generalized Roe schemes for 1D, two-phase, free-surface flows over a mobile bed. Journal of Computational Physics 227(4), 10058–10077. Rosatti, G.,Murillo, J.,and Fraccarollo, L. 2008b. Generalized Roe schemes for 1D two-phase, free-surface flows over a mobile bed. Journal of Computational Physics 227, 10058–10077. Savage, S. and Hutter, K. 1989. The motion of a finite mass of granular material down a rough incline. Journal of Fluid Mechanics 199, 177–215. Serrano, A.,Murillo, J.,and Garc´ ıa-Navarro, P. 2012. Finite volumes for 2D shallow-water flow with bed-load transport on unstructured grids. Journal of Hydraulic Research 50(2), 154–163. Siviglia, A.,Stecca, G.,Vanzo, D.,Zolezzi, G.,Toro, E.,and Tubino, M. 2013. Numerical modelling of two-dimensional morphodynamics with applications to river bars and bifurcations. Advances in Water Resources 52, 243–260. Smart, G. 1984. Sediment transport formula for steep channels. Journal of Hydraulic Engineering 3, 267–276. Soares-Frazao, S.,Canelas, R.,Cao, Z.,Cea, L.,Chaudhry, H.,Moran, A.,Kadi, K.,Ferreira, R.,Fraga-Cadorniga, I.,Gonzalez-Ramirez, N., Greco, M.,Huang, W.,Imran, J.,Coz, J. L.,Marssoli, R.,Paquier, A., Pender, G.,Pontillo, M.,Puertas, J.,Spinewine, B.,Swartenbroekx, C.,Tsubaki, R.,Villaret, C.,Wu, W.,Yue, Z.,and Zech, Y. 2012. Dambreak flows over mobile beds: experiments and benchmark tests for numerical models. Journal of Hydraulic Research 50:4, 364–375. Soares-Frazao, S. and Zech, Y. 2010. HLLC scheme with novel wave-speed estimators appropiate for two-dimensional shallow-water flow on erodible bed. International Journal of Numerical Methods in Fluids 66(8), 1019–1036. Spinewine, B. and Zech, Y. 2004. Proc., 4th Workshop of IMPACT Project. Alkema, Rotterdam, The Netherlands. Spinewine, B. and Zech, Y. 2007. Small-scale laboratory dam-break waves on movable beds. Journal of Hydraulic Research 45, 73–86. Appendix B Conservation of the coupled-Jacobian numerical scheme For the sake of clarity the necessity of the last term in the numerical scheme, B.1, is going to be explained. The lack of this term in the computed method may lead to a non conservative solution, which has been misunderstood by some authors as a diffusivity problem, Un+1 i=Un i− NE X k=1 4 X m=1 (e λ−α−β−)m keem JI,klk ∆t Ai− NE X k=1 δEIi,knklk ∆t Ai (B.1) In case of having a set of equations ∂U ∂t +∂F(U) ∂x +∂G(U) ∂y =S(U,x,y) (B.2) which can be manipulated as follows ∂U ∂t +Mn∂U ∂x +∂U ∂y −Hn∂U ∂x +∂U ∂y =Ss(U,x,y) (B.3) where Mnis the flux normal to a direction given by the unit vector n,En =Fnx+Gny, defined as Mn=∂(En) ∂U(B.4) and Hnis the flux associated to the bed slope, projected onto the unit vector n, Tbn=Sbnx+Sbny 244 Conservation of the coupled-Jacobian numerical scheme Hn=∂(Tbn) ∂U(B.5) It is possible to define the following Jacobian matrix, Jn, through the definitions in (B.4) and (B.5) Jn=Mn−Hn(B.6) which will allow us to define system (B.2) as belonging to the family of hyperbolic systems. ∂U ∂t +Jn∂U ∂x +∂U ∂y =Ss(U,x,y) (B.7) Due to the non linearity of the flow En, the Jacobian matrix has to be approximated in order to generate a local linearization. Roe (1986), proposed an approximated Jacobian matrix, e Jn,k, for clean water, imposing that e Jn(Ui,Ui) = Jn(Ui) (B.8) In Murillo and Garc´ıa-Navarro (2010a) an augmented Roe solver is proposed to allow the inclusion of the sediment transport terms. Hence, it was developed the building of an approximate Jacobian matrix e Jn,k at each kedge of each cell combining the normal flux En=Fnx+Gnywith the bed slope source term Tn,b at each cell edge, (δE−Tb)knk= (f Mn−e Hn)k(Uj−Ui) (B.9) (δE−Tb)knk=e Jn,kδUk(B.10) with δ(En)k= (Ej−Ei)nk,δUk=Un j−Un i, and Un iand Un jthe initial values at cells iand jsharing edge k. Figure B.1 collects the previous information in order to solve the Riemann problem in a 2D situation. On the other hand, the main difficulties in the definition of the approximate Jacobian matrix (B.9) for sediment transport that allows the construction of approximate solutions arises from the presence of a local Agvalue, which is variable within each cell, see Figure B.2. For this reason the net exchange of flow between internal walls must include the difference between the value of the cell in its centroid, Ag,i, and the value close to the wall within the cell Ag,I. This fact leads to new definitions of Eat each pair of cells iand j, connected through edge k, that will be referred to as 245 Un iUn j U Un i Un j nk x0 x0 x0= 0 Figure B.1: Riemann problem in 2D along the normal direction to a cell side EI=E(Ui, Ag,k)EJ=E(Uj, Ag,k) (B.11) An g,I,k An g,J,k An g,i An g,j Figure B.2: Linear representation by cells. Analyzing the flux term δEn drives to δEn =Ej−Ei=EJ−EI+ (Ej−EJ)−(Ei−EI) = (B.12) =δEJI +δEjJ −δEiI = = (δEJI +δEIi)−δEJj = = (δE− JI +δEIi) + δE+ JI −δEJj = = (δE− JI +δEIi)−(δE− IJ +δEJj) 246 Conservation of the coupled-Jacobian numerical scheme As it can be appreciated not only a flux crossing the wall, δE− JI, is necessary to ensure a conservative numerical scheme. A flux with information of the variation of AgI,i is also necessary and this is the additional term which appears in the numerical scheme.