INTERNATIONAL JOURNAL FOR NUMERICAL METHODS IN ENGINEERING Int. J. Numer. Meth. Engng. 2000; 47:2019}2037 Development of physically based meshes for two-dimensional models of meandering channel #ow M. S. Horritt*,s School of Geographical Sciences, ;niversity of Bristol, Bristol BS8 1SS, ;.K. SUMMARY The choice of mesh generation and numerical solution strategies for two-dimensional "nite element models of #uvial #ow have previously been based chie#y on experience and rule of thumb. This paper develops a rationale for the "nite element modelling of #ow in river channels, based on a study of #ow around an annular reach. Analytical solutions of the two-dimensional Shallow Water (St. Venant) equations are developed in plane polar co-ordinates, and a comparison with results obtained from the TELEMAC-2-D "nite element model indicates that of the two numerical schemes for the advection terms tested, a #ux conservative transport scheme gives better results than a streamline upwind Petrov}Galerkin technique. In terms of mesh discretization, the element angular deviation is found to be the most signi"cant control on the accuracy of the "nite element solutions. A structured channel mesh generator is therefore developed which takes local channel curvature into account in the meshing process. Results indicate that simulations using curvature-dependent meshes o!er similar levels of accuracy to "ner meshes made up of elements of uniform length, with the added advantage of improved model mass conservation in regions of high channel curvature. Copyright ( 2000 John Wiley & Sons, Ltd. KEY WORDS: mesh generation; shallow water #ow; model error INTRODUCTION Finite element models solving the two-dimensional shallow water equations are currently the most sophisticated tools capable of predicting free-surface #ows in compound channels at the reach scale (typically up to 30 km long). While these models are capable of making distributed predictions of free-surface height and depth-averaged velocities, the modelling process is fraught with uncertainties. The governing equations are incapable of capturing three-dimensional processes, and therefore do not describe the #ow physics completely. The high Reynolds numbers associated with such #ows mean that some practical technique of dealing with turbulence is also * Correspondence to: M. S. Horritt, School of Geographical Sciences, University of Bristol, University Road, Bristol B58 1SS, U.K. s E-mail: Matt.Horritt@bris.ac.uk Contract/grant sponsor: U.K. Natural Environment Research Council; contract/grant number GR3 CO/030 Copyright ( 2000 John Wiley & Sons, Ltd. Received 5 February 1999 Revised 18 June 1999 2020 M. S. HORRITT needed. The "nite element technique only delivers an approximate solution to these equations, and the accuracy of this solution may be a function of "nite element mesh discretization and the exact formulation of the numerical solver used. There is also uncertainty present in the data input to these models [1] (boundary friction parameterization, boundary conditions, bed topography [2]) which will a!ect model output. Fluvial process studies [3] can shed light on the "rst of these issues, by identifying the more signi"cant #ow mechanisms and structures and assessing their importance. Recent advances in the provision of distributed input data (especially through remote sensing techniques [4]) can also reduce uncertainty in the modelling process, and may encourage a move away from calibration studies, where these parameters are adjusted to give a best "t between model predictions and real-world observations, to more physically based parameterizations. Fluvial #ow models have been calibrated and validated through comparison with environmental data [5], and this will tend to result in "nal model uncertainties which stem from all three sources listed above. It is thus di$cult to pinpoint exactly where a model is going wrong, and a more appropriate approach would be to separate these three sources of uncertainty and hence quantify each of them. The aim of this paper is thus an assessment of the accuracy of "nite element solutions, i.e. the ability of the "nite element model to give accurate solutions of the shallow water equations. Fundamental to this problem is that of the "nite element mesh discretization. While much literature concerned with the generation of unstructured triangular meshes [6}8] exists, and techniques for adapting meshes to the spatial properties of the "nite element solutions (for example [9]) are available, mesh generation for #uvial applications has been governed by practical considerations. The need for at least three elements across a channel to capture, even in a crude fashion, variations in velocity and the channel bathymetry, generally leads to meshes of a prohibitively high number of elements if such triangular meshing strategies are adopted for reaches which can be up &10 km in length, but with channels only &10 m wide. A typical method of reducing the computational burden is to discretize the channel with elements elongated in the #ow direction, thus allowing a reasonable representation of channel cross-section and cross-channel variations in predicted variables without an excessive number of elements. A region of mesh of this type is shown in Figure 1, where a channel made up of elements with an aspect ratio of 5 : 1 is seen to be embedded in a #oodplain of more equilateral elements [5]. The choice of this aspect ratio is essentially arbitrary, and its e!ect on the quality of model solutions for #ow in the channel has yet to be assessed. Figure 1. Part of a "nite element mesh for #ood #ow modelling. The channel is described by elongated elements (with a length-to-breadth ratio of 5) embedded in the approximately equilateral elements of the #oodplain. Copyright ( 2000 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng. 2000; 47:2019}2037 DEVELOPMENT OF PHYSICALLY BASED MESHES 2021 The mesh size a!ects the model's ability to represent second-order process gradients, which can a!ect bulk #ow properties, as seen in [10]. It is unclear, however, whether the use of a "ner mesh merely a!ects the detail at which the solutions are represented or can make a more fundamental systematic di!erence to model predictions. The choice of mesh resolution is often based on computational considerations, opting for the "nest mesh feasible within the constraints of the hardware available and reasonable simulation times. The choice of numerical solver used has also been seen to a!ect bulk #ow properties [10], but previous studies have tended to concentrate upon the solver's computational e$ciency (again a pragmatic consideration) and intercomparisons between the results from di!erent solvers or "nite element models [11], rather than using objective tests of model accuracy. This has mainly been due to a lack of detailed validation data against which model results can be compared. To overcome this lack of validation data, we here develop analytical solutions of the model governing equations to aid an assessment of the e!ects of mesh discretization and solver formulation when compared with the "nite element solution. The development of analytical solutions of the shallow water equations raises a dichotomy with respect to the choice of test strategy. A modelling situation is required which is on the one hand complex enough to be non-trivial (and interesting) but on the other simple enough to be analytically soluble for the set of simultaneous non-linear partial di!erential equations. The analytical approach does, however, have the advantage of being able to invoke a number of symmetry arguments which will not be apparent in the model formulation. For example, a transformation of the governing equations into plane polar co-ordinates can reveal rotational symmetry in the #ow, a symmetry that will not be apparent to the "nite element model working in an essentially cartesian system. Thus, the problem of #ow around a channel which is rectangular in cross-section but annular in planform may succumb to analytical methods and yet provide enough variations in velocity magnitude and direction and free surface height to make an interesting "nite element problem. The result will be a rationale for the "nite element discretization of the channel region for #uvial #ow models, based on the accuracy of model solutions rather than rule of thumb. The technique also isolates and quanti"es one of the sources of uncertainty in the modelling process, an advantage over comparisons with "eld or laboratory data, where numerical e!ects may be masked by poor process representation and model parameterization. MODELLING AND ANALYSIS OF AN ANNULAR REACH The TELEMAC-2-D model [10,12] is used throughout this study, having undergone extensive development in its application to #uvial #ow problems and validation/calibration with respect to point hydrometric data [5]. The technique will, however, be applicable to any model solving the two-dimensional shallow water equations. The model aims to solve the two-dimensional shallow water equations, which can be written as a vector equation for the conservation of momentum and a scalar equation for the conservation of mass: Lv l gDvDv #(v ) +)v#g+(z #h)! t +(h+ ) v)# "0 0 Lt h C2h (1) Lh #+(vh)"0 Lt (2) Copyright ( 2000 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng. 2000; 47:2019}2037 2022 M. S. HORRITT where v describes the two-dimensional depth-averaged velocity, h is the water depth, z the bed 0 elevation, g the acceleration due to gravity, l the turbulent kinematic eddy viscosity, and C is 5 Chezy's coe$cient of bed friction. In this study, in order to obtain a comparison with an analytical solution, the equations are simpli"ed by neglecting the turbulent momentum di!usion term. This may seem to be an oversimpli"cation, but the e!ect of l on reach length studies has 5 been found to be small [10], and the value of l poorly de"ned [5]. A further simpli"cation is that 5 this study is limited to steady-state #ows. In this case, the time dependence of Equation (1) is used as a surrogate for an iterative procedure to "nd steady-state solutions, and Equation (2) reduces to a statement of the non-divergence of the mass #ux. To achieve an analytic solution, these equations (simpli"ed by assuming steady state and l "0) can be written in polar form in terms of the depth h and v , v , the radial and tangential 5 r h velocity components. The assumption of radial symmetry (for depth, velocity and bed elevation) means all derivatives with respect to the polar angle, h, are eliminated, and the continuity equation becomes 1 d (rv h)"0 r dr r (3) With all variables now simply functions of the radius, r, the partial derivatives are replaced with normal derivatives. The derivation of the polar versions of the spatial derivatives found in Equations (1) and (2) can be found in [13]. If the side walls of the channel are impermeable, v will r be zero at the walls, and from Equation (3) zero across the whole cross-section (assuming non-zero depth). Thus, all radial velocities can be eliminated, and the equations for radial and tangential momentum become dh v2 ! h "0 dr r (4) g dz gv2 0# h "0 C2h r dh (5) g These two equations can be interpreted physically: the "rst equates the free surface slope in a radial direction with the centripetal acceleration required to maintain the circular motion of the #uid, the second balances down reach slope against bed friction. If we assume constant friction and down reach slope, v2 is easily eliminated and the equations integrated to give the analytical solutions: h C C2 dz 0 h"D exp ! gr dh C D dz DC2 C2 dz 0 exp ! v2" 0 h dh r gr dh (6) D (7) where D is a constant of integration. It will be useful to derive the volumetric #ow rate, Q, down the reach from these expressions r1 Pr hvh dr Q" (8) 0 Copyright ( 2000 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng. 2000; 47:2019}2037 DEVELOPMENT OF PHYSICALLY BASED MESHES 2023 where the limits r and r are the inner and outer radii of the channel. The expression de"es 0 1 analytic integration, but can be integrated numerically (using the trapezium rule) or can be written as a Taylor expansion in powers of the dimensionless channel width, u"(r !r )/r . The 1 0 0 power series converges to within 1 per cent with terms up to and including O(u3). This enables the constant of integration D to be derived from a given #ow rate Q. The experimental methodology runs as follows. Eleven "nite element meshes were generated of channels 20 m in width and of inner radii 50 and 80 m (Figure 2). These parameters were taken as being representative of a small U.K. meandering river, such as the upper Thames. Meshes both of unstructured equilateral elements and structured elongated elements are included. The down reach slope was set at 0.001, giving values of dz /dh of !0.06 m for the smaller meshes and 0 !0.09 m for the larger. The Chezy friction coe$cient was set at 40 m1@2 s~1 and the #ow rate at 40 m3 s~1, a typical "gure for bankful #ow in such a channel. The boundary conditions for the TELEMAC-2-D model were the imposed #ow rate at the upstream end (#ow is in an anticlockwise direction) and periodic boundary conditions were applied at the downstream end, the water depths at this boundary being "xed as equal to those at the equivalent nodes on the upstream boundary. These boundary conditions leave the model free to "x both velocities and water depths, but still gives a well-posed numerical problem. The TELEMAC-2-D model incorporates a number of options for the numerical solution of the shallow water equations [12], and two of these are investigated here. TELEMAC-2-D generally uses a fractional step technique: advection terms are solved initially, with propagation, di!usion and source terms solved in a second step, with these steps being repeated to achieve iterative convergence. The "rst scheme uses a conservative streamline-upwind Petrov}Galerkin (SUPG) technique for the advection of depth in the continuity equation, and the method of characteristics for the advection of velocity. This scheme solves the continuity equation exactly (no fractional step method is used), and should, in theory, give exact mass conservation. The second scheme uses a #ux conservative transport (FCT) method for the advection both of depth and velocity, and although not guaranteed to be totally mass conservative, mass error can be reduced by using a number of sub-iterations of the advection solving step. The mass conservative nature of the SUPG scheme without the need for sub-iterations makes it a popular choice for reach length simulations [5]. The combination of the 11 computational meshes and the 2 advection schemes gives a total of 22 simulations. Figures 3 and 4 show the modelling results using the SUPG and FCT schemes in comparison to the analytical solutions for mesh 3, the "nest mesh used in the study. The FCT scheme gives markedly better results, both in terms of the absolute errors in depth and velocity and with respect to the cross channel velocity gradient. The poorest results were those from mesh 6 (Figures 5 and 6), where #ow depth has been over-predicted by '50 per cent for the SUPG scheme. In terms of mesh discretization, it was found that the ratio of mesh size to channel planform radius (or the elemental angular deviation) was the most signi"cant control on model prediction errors, having a far greater in#uence than the number of elements across the channel, the choice of structured or unstructured mesh or the channel radius. Figures 7 and 8 summarize this result, demonstrating that both schemes converge with decreasing angular deviation, but with FCT giving the lower overall errors. The dependence of model accuracy on element angular deviation could be a result of two e!ects: the angular deviation of the streamline over an element may a!ect the "nite element model, and the representation of the channel itself is di!erent. Each mesh represents the boundary of the circular channel as two concentric polygons, and for the "ner meshes, these polygons coincide more closely with the true circular boundaries. Copyright ( 2000 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng. 2000; 47:2019}2037 2024 M. S. HORRITT Figure 2. Finite element meshes for modelling #ow around an annular reach. Copyright ( 2000 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng. 2000; 47:2019}2037 DEVELOPMENT OF PHYSICALLY BASED MESHES 2025 Figure 2. Continued. Two obvious conclusions can be drawn from these results. Firstly, the FCT scheme gives more accurate solutions than the SUPG scheme, despite the better mass conservation properties theoretically achievable with the SUPG scheme. This may in part be a property of this relatively simple test case, and when dynamic e!ects and turbulence come to be taken into account, the signi"cance of the advection terms may be reduced, and the two schemes may give more similar results. Given the lack of validation data for more complex #ows, however, this is mere Copyright ( 2000 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng. 2000; 47:2019}2037 Copyright ( 2000 John Wiley & Sons, Ltd. Figure 4. (a) Comparison between model predicted depth and true depth for the FCT scheme operating on mesh 3. (b) Comparison between model predicted velocity and true velocity for the FCT scheme operating on mesh 3. Figure 3. (a) Comparison between model predicted depth and true depth for the SUPG scheme operating on mesh 3. (b) Comparison between model predicted velocity and true velocity for the SUPG scheme operating on mesh 3. 2026 M. S. HORRITT Int. J. Numer. Meth. Engng. 2000; 47:2019}2037 Copyright ( 2000 John Wiley & Sons, Ltd. Figure 6. (a) Comparison between model predicted depth and true depth for the FCT scheme operating on mesh 6. (b) Comparison between model predicted velocity and true velocity for the FCT scheme operating on mesh 6. Figure 5. (a) Comparison between model predicted depth and true depth for the SUPG scheme operating on mesh 6. (b) Comparison between model predicted velocity and true velocity for the SUPG scheme operating on mesh 6. DEVELOPMENT OF PHYSICALLY BASED MESHES 2027 Int. J. Numer. Meth. Engng. 2000; 47:2019}2037 2028 M. S. HORRITT Figure 7. Height prediction error as a function of element angular deviation for the SUPG and FCT schemes. Figure 8. Velocity prediction error as a function of element angular deviation for the SUPG and FCT schemes. speculation, and the poor treatment of the advection term by the SUPG scheme would indicate that it should be abandoned in favour of the FCT scheme. Secondly, the strong dependence of solution accuracy upon the elemental angular deviation would imply that local channel curvature should be taken into account in mesh generation strategy, rather than simply using meshes with uniform element lengths. This is explored further in the next section. Copyright ( 2000 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng. 2000; 47:2019}2037 DEVELOPMENT OF PHYSICALLY BASED MESHES 2029 A CURVATURE-DEPENDENT MESH GENERATION STRATEGY The results given above can be used to formulate a channel meshing strategy given the level of accuracy required of the "nite element solution. Mesh generation with the constraint of uniform angular deviation would seem to be a better use of computational resources than elements of uniform length, concentrating computational resources in bend regions where the solutions are likely to be less accurate due to the high local streamline curvature. To this end, a structured "nite element channel mesh generator is developed, which calculates element lengths according to local channel curvature. The local curvature, i of an arc can be de"ned as dh i" ds (9) where h is the angle between the tangent to the curve and some arbitrary direction, and s is the arc length parameter along the curve. i can thus be thought of as factor for converting between lengths and angular deviations, and a mesh generator should aim to produce elements not of equal length, l, but with equal DiDl. Since information on the channel planform is likely to be in the form of a series of nodes de"ning the channel banks, a more useful "nite di!erence approximation of curvature at a node is given by (v1]v2) 2 i" ' Dv1 D Dv2 D (Dv1 D#Dv2 D) (10) where v and v are vectors linking the previous and successive nodes to the one in question. The 1 2 "rst term gives the angular deviation dh (using the small angle approximation sin h&h), the second approximates the change in arc length ds. The channel form is likely to be subject to a certain degree of noise (uncertainty in the position of the individual nodes) due to the digitization process and errors in the surveying of the channel, and this may generate high values of DiD, obscuring the underlying channel curvature. Coherent smoothing of i (observing the sign of i) along the channel can eliminate this noise, but this can have the drawback of producing zero curvature in cross-over regions, where the contributions from convex and concave sections before and after cancel. Further incoherent smoothing (DiD is smoothed rather than i) is used to avoid this problem, allowing DiD to di!use into these cross-over regions and promoting the use of shorter elements therein. This smoothing procedure is illustrated in Figure 9. Having established an appropriate measure of curvature along each of the banks of the channel, the meshing algorithm proceeds as follows. A cross-section moves along the reach, in the direction of a vector determined from the mean of the direction of the two banks. The distance from the start of the reach, l, is measured, and if the quantity DiDl exceeds a threshold, the end of the "rst set of elements is inserted. The process is then repeated, inserting elements whenever the length from the last set, l is such that DiDl exceeds the threshold. A second threshold value for l itself is also set, to prevent overly long elements forming in straight sections of the reach. Copyright ( 2000 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng. 2000; 47:2019}2037 2030 M. S. HORRITT Figure 9. Schematic showing action of coherent and incoherent smoothing on the curvature measure for a meandering channel. This mesh generation strategy is tested on a &1 km reach of the upper Thames, in southern England, the results being shown in Figure 10, along with the results of applying two other mesh generation strategies commonly adopted in #uvial #ow studies [5]. Meshes 1 and 2 were formed using the older strategy of equal size elements (e!ectively "xing i"0 along the reach) with lengths 10 and 30 m. The elemental angular deviations are up to 903 for mesh 2. Meshes 3 and 4 use the new curvature-dependent strategy with angular deviations between elements of 0.4 and 0.1 rad. The technique ensures that the meander regions of high curvature are represented by smaller elements, and that the angular deviation between sets of elements remains approximately constant (Figure 11). Finally, two meshes (5 and 6) were also generated using unstructured near-equilateral elements of size 2 and 10 m. These parameters, along with the number of nodes and elements for each mesh, are summarised in Table I. The six meshes thus represent coarse and "ne examples of the results of each of the three generation strategies. Copyright ( 2000 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng. 2000; 47:2019}2037 DEVELOPMENT OF PHYSICALLY BASED MESHES 2031 Figure 10. Finite element meshes of the Thames reach. Mesh 1 is shown in its entirety, for the others a detail of a meander region is shown. PERFORMANCE OF CURVATURE-DEPENDENT MESHES Steady state simulations were performed using the six meshes of Figure 10 at approximately bankful discharge (40 m3 s~1) with both the SUPG and FCT schemes. Tests were also performed Copyright ( 2000 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng. 2000; 47:2019}2037 2032 M. S. HORRITT Figure 11. Detail of mesh 4, demonstrating the constant element angular deviation property. Sets of edges between elements from regions of the channel with di!erent curvatures subtend similar angles at points a and b. Table I. Parameters for the six Thames meshes. Mesh 1 2 3 4 5 6 Element size (m) Angular deviation (rad) 10 30 0.4 0.1 2 10 Nodes Elements 378 130 146 470 7685 522 562 190 214 700 14370 770 on meshes of intermediate resolution, but the results were predictable from those from the coarse and "ne meshes, and have therefore been omitted for the sake of clarity. The results, in the form of water depths as functions of down reach distance, are shown in Figures 12 and 13. Cross-reach variations in water level were eliminated by taking the mean depth across the channel. The predicted water levels are more consistent between meshes for the FCT scheme, where some level of convergence can be seen (the similarity between prediction from meshes 4 and 5 for example). The "gures demonstrate the advantages of the new curvature dependent meshing scheme: the results from mesh 4 are similar to those from mesh 5 but with fewer elements, and the advantage of shorter computation times. This advantage will become greater when such channel meshes are embedded in larger meshes for the simulation of #oodplain #ows. Because process gradients on the #oodplain are lower than those in the channel, generally larger elements can be used, but there is often an upper limit imposed on their size by the need to couple with the mesh in the channel. The curvature-dependent scheme will impose small #oodplain elements on the inside of meander bends, but should make for larger elements surrounding straighter sections of the reach. Copyright ( 2000 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng. 2000; 47:2019}2037 DEVELOPMENT OF PHYSICALLY BASED MESHES 2033 Figure 12. Water depths (averaged across the channel) predicted using the SUPG scheme as a function of along reach distance, for the six "nite element meshes. Figure 13. Water depths (averaged across the channel) predicted using the FCT scheme as a function of along reach distance, for the six "nite element meshes (note the change in scale of the vertical axis). A measurement of mass conservation error is often used as a measure of con"dence in model predictions in situations such as these where no validation data is available. Mass conservation may be de"ned on a global basis, in terms of mass #uxes at the boundaries and the rate of change of volume of #uid within the domain (which will be zero at steady-state), but it may be more instructive to look at mass balance errors on a local basis. For a steady-state simulation, the divergence of the mass #ux vector (vh) should be zero, and its actual value will be a measure of local mass balance error. This can be calculated for each element by adding the #uxes across its Copyright ( 2000 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng. 2000; 47:2019}2037 2034 M. S. HORRITT Table II. Mass balance errors for the Thames reach simulations. Mesh/advection scheme 1 1 2 2 3 3 4 4 5 5 6 6 Final d</dt (m3 s~1) : +(vh) dA (m3 s~1) [1 :(+(vh))2 dA]1@2 A (cm s~1) 0.16 0.90 !0.65 0.011 !0.65 0.45 0.22 !0.15 (0.001 (0.001 (0.001 (0.001 1.36]10~3 !5.68]10~1 !2.76]10~2 4.71]10~1 !5.82]10~2 !7.72]10~1 !3.13]10~3 !3.50]10~1 5.52]10~4 !8.85]10~1 !8.27]10~3 !1.23]10~1 2.38 2.70 4.24 4.02 4.11 3.80 2.52 3.04 0.96 1.03 1.27 1.28 SUPG FCT SUPG FCT SUPG FCT SUPG FCT SUPG FCT SUPG FCT three boundaries and dividing by the element area (according to the divergence theorem). This is essentially a "nite volume measure of mass balance. The integral of the divergence over the area of the model domain is presented in Table II. Some of this error may be assigned to the fact that the simulations are not true steady state, but are really dynamic simulations with steady-state boundary conditions, run for long enough for transients to die down. Rates of change of volume in the domain are also given in Table II, measured over the last two time steps of each simulation. For a simulation that conserved mass, d</dt and the integral of the divergence would be equal in magnitude and opposite in sign (a result obtained by integrating the continuity equation over the domain), but there appears to be no such correlation in Table II. This means that some of the mass balance error is genuine, and not due to transients still present in the simulation. A second measure of the mass balance error can be found by calculating the root-mean-square value of the divergence over the reach, and is also shown in Table II. The "gure can be interpreted as the volume error per unit area over the reach, or the depth of water lost or created averaged over the mesh. All the errors are of the same order, but the errors for the triangular meshes are lower than those for the regular meshes. As a further demonstration of potential problems with mass conservation, Figures 14 and 15 show down reach #ow (which should be uniformly 40 m3 s~1) for meshes 1 (10 m elements) and 4 (angular deviation of 0.1 rad). The #ow for mesh 1 varies by up to 5 per cent (2 m3 s~1), mainly in regions of high curvature. Mesh 4 performs better in the main section of the reach, but there are large #uctuations in #ow towards the upstream boundary, perhaps associated with the large elements in this part of the reach. Elsewhere, the discrepancy is generally less than 0.5 m3 s~1. Both numerical advection schemes (SUPG and FCT) give similar mass balance errors. The better mass conservation properties of the curvature dependent mesh are a further advantage of its use, but some special treatment at the upstream end of the reach may be required if adequate mass conservation is required in this region. The #uctuations in along reach #ow are mainly due to velocity variations, the depths predicted by TELEMAC-2-D being relatively much smoother and free from node-to-node oscillations (the SUPG and FCT schemes are very good at damping oscillations in depth). Both numerical schemes seem to encounter problems dealing with the non-linear velocity advection term, Copyright ( 2000 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng. 2000; 47:2019}2037 DEVELOPMENT OF PHYSICALLY BASED MESHES 2035 Figure 14. Flow rates as functions of along reach distance for the two numerical schemes operating on mesh 1. Figure 15. Flow rates as functions of along reach distance for the two numerical schemes operating on mesh 4. a problem not encountered with the study of the annular reach. The problem only arises in the Thames reach probably because irregularities in channel planform excite these modes of oscillation, especially in regions of high curvature and at the upstream boundary of the domain. The modes of node to node oscillations will tend to be aligned with the channel direction for meshes 1}4 and at approximately 303 to the channel for the triangular meshes (Figure 16). This may explain the poorer mass conservation properties (at least according to the RMS measure) of the regular meshes, as the oscillations are able to propagate freely along the channel. The poor performance of even the conservative SUPG scheme is partly due to the "nite volume method used to assess mass balance in this study. The conservative scheme ensures mass Copyright ( 2000 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng. 2000; 47:2019}2037 2036 M. S. HORRITT Figure 16. (a) a possible mode of node to node oscillation for a regular mesh (e.g. meshes 1}4) of elongated elements. (b) the mode of oscillation for approximately equilateral elements is no longer aligned with the channel direction. conservation integrated over the "nite element test functions, allowing velocity variations to cancel and conserve mass on a global basis (which is the usual measure used). The results do, however, demonstrate that even though the conservative SUPG scheme gives much better results in terms of global mass balance, both schemes perform similarly in terms of local mass balance. This in some way negates one of the advantages of the conservative SUPG scheme, in that its globally conservative property is not re#ected in the local mass balance, and even simulations that appear to conserve mass at the boundaries of the domain may have local mass balance problems, caused by unphysical spatial #uctuations in the predicted velocity "eld. CONCLUSIONS Uniform #ow around an annular reach has proved to be a good test of the TELEMAC-2-D model when compared with analytical solutions of the shallow water equations, discriminating between di!erent numerical schemes and "nite element meshing strategies. The #ux conservative transport scheme for the advection of velocity and depth gives more accurate solutions than the streamline upwind Petrov}Galerkin technique when compared to the analytical solutions. The solution accuracy is a strong function of the elemental angular deviation, and this measure should be used to drive any mesh generation strategy, rather than simply prescribing element length. The analysis of this simple scenario has thus given a rationale for the choice of numerical scheme and "nite element mesh. With this rationale in place, a mesh generator for meandering #uvial channels has been developed, with element sizes dependent on a measure of local channel curvature. Meshes generated using this technique appear to give better results than those of uniform element size, when compared to irregular meshes of small triangular elements. Despite the lack of validation data for this reach, the FCT scheme still seems to o!er an advantage over the SUPG scheme, in terms of the apparent convergence of the solutions as progressively "ner meshes are used. Both numerical schemes su!er from mass conservation errors due to unphysical node-to-node oscillations in the predicted velocity "eld, but these are smaller in magnitude for the curvature-dependent meshes, o!ering a further incentive for their use. It must be noted that this modelling study makes certain simpli"cations to facilitate the analytical solution of the shallow water equations. Turbulence has been neglected, but its e!ect on #oodplain #ow studies is as yet only poorly understood, and this o!ers some justi"cation for its omission. The assumption of a constant eddy viscosity is in itself a gross oversimpli"cation, and the use of a more sophisticated model (such as the k}e model) may have a bearing on the mesh generation strategy. Model accuracy may change when #ow through more complex channel Copyright ( 2000 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng. 2000; 47:2019}2037 DEVELOPMENT OF PHYSICALLY BASED MESHES 2037 geometries (changes in curvature, variable widths, complex bathymetry) are simulated, and the modelling of overbank #ows (with the implication of complex interactions between channel and #oodplain #uid) may introduce further inaccuracies. The simple relationship between model accuracy and streamline curvature found here may then no longer hold. It is also unclear from this study whether model accuracy is solely dependent on streamline curvature, or whether the representation of the model boundary is as important, but either way, the curvature-dependent mesh generation strategy is still valid. Dynamic simulations may pose more problems, as the "nite element mesh may have an e!ect on the subtle changes that will need to be represented to model the passing of a #oodwave that may last hours or days. These results will have some practical implications for future #uvial modelling studies. Given the error as a function of angular deviation (as in Figures 7 and 8), and the accuracy required of model predictions, an appropriate mesh resolution can be selected. Conversely, using a given elemental angular deviation, the uncertainty in the "nite element solutions due to the approximate way in which the governing equations are solved can be estimated. In studies comparing modelled results to laboratory or "eld data, this numerical error can be discounted, and the remaining discrepancies assigned to poor process representation and model parameterization. This will shed light on how well the governing equations describe the physics of the problem, and aid in the calibration and validation process. ACKNOWLEDGEMENTS This research has been supported by the UK Natural Environment Research Council grant GR3 CO/030. Thanks go to Paul Bates (University of Bristol) and Jean-Michel Hervouet (Laboratoire National d'Hydraulique, ElectriciteH de France) for comments on the text. REFERENCES 1. Anderson MG, Bates PD. Evaluating data constraints on two dimensional "nite element models of #oodplain #ow. Catena 1994; 22:1}15. 2. Bates PD, Horritt MS, Hervouet J-M. Investigating two-dimensional, "nite element predictions of #oodplain inundation using fractal generated topography. Hydrological Processes 1998; 12:1257}1277. 3. Bathurst JC. Flow processes and data provision for channel #ow models. In Modelling Geomorphological Systems, Anderson MG (ed.). Wiley: Chichester, 1988; 127}150. 4. Bates PD, Horritt MS, Smith CN, Mason DC. Integrating remote sensing observations of #ood hydrology and hydraulic modelling. Hydrological Processes 1997; 11:1777}1795. 5. Bates PD, Stewart MD, Siggers GB, Smith CN, Hervouet J-M, Sellin RJH. Internal and external validation of a two-dimensional "nite element code for river #ood simulations. Proceedings of the Institution of Civil Engineers, =ater, Maritime and Energy 1998; 130:127}141. 6. Frey WH. Selective re"nement: a new strategy for automatic node placement in graded triangular meshes. International Journal for Numerical Methods in Engineering 1987; 24:2183}2200. 7. Rebay S. E$cient unstructured mesh generation by means of Delaunay triangulation and Bowyer}Watson algorithm. Journal of Computational Physics 1994; 110:23}38. 8. Frey WH, Field DA. Mesh relaxation: a new technique for improving triangulations. International Journal for Numerical Methods in Engineering 1991; 31:1121}1133. 9. Knupp P. Mesh generation using vector "elds. Journal of Computational Physics 1995; 119:142}148. 10. Bates PD, Anderson MG, Price DA, Hardy RJ, Smith CN. Analysis and development of hydraulic models for #oodplain #ow. In Floodplain Processes, Anderson MG, Walling DE, Bates PD (eds). Wiley: Chichester, 1996; 215}254. 11. Bates PD, Anderson MG, Hervouet J-M. An initial comparison of two 2-dimensional "nite element codes for river #ood simulation. Proceedings of the Institution of Civil Engineers, =ater, Maritime and Energy, 1995; 112:238}248. 12. Hervouet J-M, Van Haren L. Recent advances in numerical methods for #uid #ow. In Floodplain Processes, Anderson MG, Walling DE, Bates PD (eds). Wiley: Chichester, 1996; 183}214. 13. Arfken G. Mathematical Methods for Physicists (3rd edn). Academic Press: London, 1985. Copyright ( 2000 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng. 2000; 47:2019}2037

1/--страниц