Strojniški vestnik - Journal of Mechanical Engineering 64(2018)9, 536-542 © 2018 Journal of Mechanical Engineering. All rights reserved. D0l:10.5545/sv-jme.2017.5183 Original Scientific Paper Received for review: 2017-12-28 Received revised form: 2018-05-16 Accepted for publication: 2018-06-06 Analysis of Flow in a Curved Channel Using the Curvilinear Orthogonal Numerical Mesh Mario Krzyk* - Matjaž Četina University of Ljubljana, Faculty of Civil and geodetic Engineering, Slovenia Along the flow in a curved channel, energy losses occur due to external influences on the considered water body and due to the internal friction caused by the turbulent flow of water. In practice, free surface flow is often calculated by using two-dimensional depth-averaged mathematical models. The mathematical model PCFLOW2D-ORTHOCURVE was built to analyse the flow in steep curved streams. It is based on the orthogonal curvilinear numerical mesh. Before the application of the model on examples of complex geometry natural flows, its accuracy was verified in the case of flow in a semi-circular curved channel. A 20 m wide channel with a curvature radius of 30 m in the axis and the horizontal bottom of the channel was applied. There were straight channel sections of 100 m before and behind the curve, respectively. The roughness coefficient of the channel was negligibly small, thus eliminating the impact of bottom and banks on the flow and energy losses. The flow depth distribution was calculated for different average flow velocities in the curve of approximately 1 m/s to about 9 m/s. The critical depth, the depth of water at the upstream and downstream boundaries of the straight section of the channel, and the mean velocity of the flow in the curved section were determined, by means of which the theoretical value of the transverse slope of free surface was calculated. This was compared with the numerically calculated value. Keywords: flow in a curved channel, two-dimensional mathematical model, orthogonal curvilinear co-ordinates, depth averaged flow, PCFLOW2D-ORTHOCURVE Highlights • Two-dimensional depth averaged mathematical model was used. • Mathematical model based on Cartesian coordinate system was transformed to the mathematical model based on curvilinear orthogonal numerical mesh. • Flow in semi-circular channel with rectangular cross-section was analyzed. • Transversal slope of the free surface in the curved section of the channel calculated by mathematical model was compared with theoretical solution. • Very good agreement was reached. 0 INTRODUCTION Flows in steep natural channels are usually very dynamic, which is a consequence of the channels' distinctly irregular shape and nonuniform course. The mathematical modelling of such flows, using the computational mesh based on the Cartesian coordinate system, is rather irrational and can cause inaccurate results and instability of the calculation procedure [1]. The use of curvilinear coordinates allows for better adjustments to the channel shape and profile. The PCFLOW2D mathematical model [2], based on the Cartesian numerical mesh, was upgraded by using the hydrodynamic equations derived for the orthogonal curvilinear coordinate system. This expanded the model's applicability to solving practical problems that require very fine discretisation due to the morphology of the channel bed and banks. There are two more advantages to using such an approach: • it allows the impact of the centrifugal force in curved sections to be taken into account, • because of the improved alignment of the numerical mesh with the flow direction it is possible to avoid computational errors due to excess numerical diffusion. This paper shows flow equations and the applicable turbulence model used in the PCFLOW2D-ORTHOCURVE model. Prior to applying the model on examples of complex geometry natural flows, its accuracy was verified on the case of a flow in a semi-circular curved channel with a rectangular cross-section. The results were compared against the theoretical equation for calculating the transverse water surface slope in a curve. In recent literature, examples of flow analysis in the bend using an analytic depth averaged two-dimensional (2D) model [3] or using various numerical approaches, including a gene expression programming model [4], can be found. These are examples of flow in the bend of 90°. Several analysis were also performed on flow in a confluent meander bend channel using large eddy simulations [5]. Due to the length of flow in such curved channel, under 536 *Corr. Author's Address: University of Ljubljana, Faculty of Civil and Geodetic Engineering, Jamova 2, Ljubljana, Slovenia, mkrzyk@fgg.uni-lj.si Strojniski vestnik - Journal of Mechanical Engineering 64(2018)9, 536-542 certain hydrodynamic conditions it is not possible to establish a parallel depth along the inner and outer edge of the channel. Therefore, a channel with a curve of 180° was used. Also, in such channel geometry it is easier to avoid the influence of the inlet and outlet hydrodynamic characteristics on the flow in the bend. 1 MATHEMATICAL MODEL EQUATIONS As analysis was focused on 2D depth-averaged flow, the definition of two coordinate axes in the horizontal plane (%, n) will be sufficient. The meaning of Lame's coefficients m% and mn [6] used in the equations is evident from Fig. 1. The basic hydrodynamic equations used in the model based on the orthogonal curvilinear coordinate system are given below. The continuity, Eq. (1) and momentum equations, (Eqs. (2) and (3)) describing a two-dimensional unsteady depth averaged flow are written in the conservative form. d(hu^) d(hvn ) hu4 dmn dt m^d^ mndn m^mv 3Ç hv„ dm, = 0, m^mv dn (1) dhuP d(hu- ) d(hu-vn) hu-vn dm-dt m-d-/2 2 ) h dmv mdn n m m dn mn m-d- - gh m-d- m-d- dzb - ghn2 1 d d- . m- d- j dn h4'3 ^ mP du, hü r ----- V mn dn j (2) dhvn , d(hu4Vn) d(H) hu^n dmn dt m£dâ, mndn m4mn dâ, (2 2\ h dmâ yn — U )- m4 mndn 4 n -- gh- mndn mndn r 2 2 V JU, +V b - ghn2 W 4 n 1 d m4 d4 J dn 4/3 hv,^ mn dn J (3) Fig. 1. The meaning of Lame's coefficients in the orthogonal curvilinear (natural) coordinate system The 2D mathematical model PCFLOW2D for simulating velocity and water depth fields, based on the Cartesian coordinate system, contains the basic differential equations governing the time-averaged flow, where the impact of Reynolds (turbulent) stresses is captured using the k- s turbulence model. Such an approach was used in developing the PCFLOW2D-ORTHOCURVE mathematical model, where the use of the k- s model was kept in the curvilinear coordinates. Therefore, the transport equations for turbulent kinetic energy per unit mass and the degree of its dissipation had to be converted to the natural coordinate system. The equation for turbulent kinetic energy per unit mass k is: dhk d(hu^k) d(hvnk) hu^k dmn hvnk dm^ dt mÇdÇ m^dn mm dÇ mm dn 1 mem„ d dÇ hUeLm±dk_ mÇ dÇ _d_ dn h ujLm^dk_ mn dn +hG - clhs + hPv, and the equation for its dissipation s: (4) dhs d(hu— d(hv—) hu— dmn hv— dmÇ __ + + + + m.dÇ 1 d m^dn mm dÇ mm dn dÇ hUeLm±d± °s m _d_ dn Oe mn dn +c1 — hG - c2 —) 1 k 2 k hP . (5) The expression of G (the production of k due to horizontal velocity gradients of the basic flow) in + + Analysis of Flow In a Curved Channel Using the Curvilinear Orthogonal Numerical Mesh 537 Strojniski vestnik - Journal of Mechanical Engineering 64(2018)9, 536-542 orthogonal curvilinear coordinates is calculated using the following equation: G = u ef du^ m.dÇ duÇ dvn V mndn j dvn V mndn m^dÇ (6) Effective viscosity value vef is derived using Eq. (7), [7] as the sum of coefficient of laminar and coefficient of turbulent viscosity: (7) The values for Pkv and Pev are calculated using the following Eqs. (8) [7], (9) and (11) [7]: 3 4 P = P = Pkv Ck h ' P¡¡v C¡¡ h2 ' =Vcf(+v ), (8) (9) = Cs = CsChKAH. (10) ylcf cf The value for the constant cer is 3.6 for the laboratory channels. In the Table 1 used values for the other empirical constants proposed by Launder and Spalding [8] are presented. Table 1. Values of the constants in the k- e model Cl C2 Gk G 1.44 1.92 0.09 1.00 1.21 2 NUMERICAL METHOD Finite volume numerical scheme proposed by [9] is used for solving partial differential equations presented in the previous section. Characteristics of the method are staggered control volumes and "hybrid" scheme. The "hybrid" scheme is a combination of the upwind and central difference scheme. The first order upwind scheme assures simplicity and robustness [9] and it remains stable even with very complex geometry, relatively coarse numerical grid and complicated boundary conditions. The use of a curvilinear numerical mesh also contributes to the stability of the method. For the depth corrections an iterative procedure is used. A fully implicit scheme is used for time integration providing stable and accurate solutions. Numerical method could be used to simulate subcritical and supercritical flows [9]. 3 FLOW IN THE CURVED SECTION Free-surface flow in a curved section is much more complex than that in a straight section. One of its basic characteristics is the secondary flow phenomenon, which is generated as a result of the imbalance between the surface gradient in the transversal flow direction and the centrifugal force. The pressure gradient at the bottom in the direction towards the inner side of the bend dominates over the centrifugal force, generating the flow at the bottom towards the centre of the bend and the flow at the surface towards the outer bank. The transverse velocity component is a size class smaller than the longitudinal or primary velocity component. This results in additional energy losses. The energy loss depends on the relationship between the water depth and the radius of the curve. With smaller wall roughness, due to the impact of transverse circulation, energy losses can even exceed half the value of the total loss [10] and [11]. For the cases of curved channels, where acting centrifugal forces cause the surface gradient to occur in the direction transversal to the flow, we can theoretically calculate the surface height difference between the outer and inner banks, Dh, using the following theoretically derived formula: v R Ah = — ■ ki-ot. g R_. (11) The equation refers to flows in curves where there are no energy losses due to friction along channel walls or changes in channel cross-section geometry in the transversal and longitudinal direction (e.g. narrowing, sills, or bed dredging). For testing the 2D mathematical model based on the orthogonal curvilinear numerical mesh, a semi-circular smooth channel with a rectangular cross-section was used, which on both ends extended into long straight sections. This channel shape was chosen to remove the impact of boundary conditions on the flow in the curve; we obtained results that are conditioned only by the channel geometry and the appropriate discharge or flow velocity. A rectangular channel of a width of 20 m with a semi-circular curvature and a channel curvature radius of 30 m in the axis was used. The straight inflow section and the 2 2 + 2 c 538 Krzyk, M. - Cetina, M. Strojniski vestnik - Journal of Mechanical Engineering 64(2018)9, 536-542 final section of the model are L = 100 m long and the channel bottom is horizontal. The channel geometry is shown in Fig. 2. Manning's roughness coefficient n = 1x10-6 s/m1/3 was assumed for the channel bottom and banks. The channel geometry in the mathematical model is described by 7104 points; 24 points in the direction transversal to the flow and 296 points in the direction of the flow (n). In the straight section of the channel the mesh was uniform, with a distance of computation points 1 m x 1 m and in the curved section adjusted to the requirements of the orthogonal curvilinear mesh. 4 RESULTS AND DISCUSION The hydraulic conditions were calculated for six discharges. The discharge rate varied from Table 2. Comparison of mathematical model result and theoretical solution of transversal slope of water surface Q hcrit hups ^downs ^aver AhCurv (model) Ahcurv (theory) [m3/s] [m] [m] [m] [m/s] [m] [m] 100 1.366 5.002 4.999 1.000 0.07 0.07 300 2.841 5.017 4.984 3.025 0.67 0.66 500 3.994 5.321 5.160 4.980 1.81 1.75 700 4.998 6.228 4.519 6.870 3.42 3.34 900 5.910 4.600 5.125 8.020 4.64 4.55 1100 6.756 5.200 5.805 8.675 5.52 5.32 Discharge Fig. 3. Height differences between the inner and outer edges of the curve from the mathematical model and theory Analysis of Flow In a Curved Channel Using the Curvilinear Orthogonal Numerical Mesh 539 Strojniski vestnik - Journal of Mechanical Engineering 64(2018)9, 536-542 Q = 100 m3/s to 1100 m3/s, with a step of 200 m3/s. For these discharges, with the average flow velocity in the curve ranging from approx. 1 m/s to approx. 9 m/s, the water level formation in the channel was calculated. In the cases investigated, flow with a horizontal surface and uniform velocity distribution was created at the upstream and downstream model boundaries. In the curved section, the water surface appropriately leaned towards the inner bank of the channel. Table 2 summarises the results of the calculated values: the flow rate (Q) for each case, calculated critical depth hcrit, water depth at the upstream section of the straight channel hups, water depth at the downstream section of the channel hdowns, average flow velocity in the curve vaver, between channel stationing 100 and 194 m used to calculate the theoretical value of the transverse slope of the surface. The next two columns give the height differences between the inner and outer edges of the curve Dhcurv, derived from the mathematical modelling results and calculated theoretically using Eq. (9). Those values are presented graphically in Fig. 3. Figs. 4 to 6 show the water surface formation at outer and inner edges of the channel for discharges 100 m3/s, 700 m3/s, and 1100 m3/s. At Q = 100 m3/s and the average water depth in the channel of around 5 m (Fig. 4), flow velocity of approx. 1 m/s is created. The flow along the entire channel is in the subcritical regime and, because of the low flow velocity, the transverse slope of the free surface at the curve is low as well (Dhcurv = 7 cm). At both edges the surface has a falling trend in the flow direction. Under such conditions the average velocity along the channel does not change considerably; therefore, to calculate the water surface slope in the bend the velocity at any cross-section can be taken without causing any considerable differences between the theoretically calculated and the modelled free-surface slope in the bend. The shape of the surface line at both edges is similar for the discharges up to 700 m3/s and follows Stationing | Fig. 4. Water surface profile along channel walls at Q = 100 m3/s 75 100 125 150 175 200 225 250 Stationing [m] Fig. 5. Water surface profile along channel walls at Q = 700 m3/s 540 Krzyk, M. - Cetina, M. Strojniski vestnik - Journal of Mechanical Engineering 64(2018)9, 536-542 \ \ -Outer edee \ 'inner edge • * • 50 75 100 125 150 175 200 225 250 Stationing [m] Fig. 6. Water surface profile along channel walls at Q = 1100 m3/s the shape shown in Fig. 4 at discharge Q = 100 m3/s. By comparing the offset of the water level at the channel edges and the mean depth it could be find that the lowering of the water level at the inner edge of the curve is much more pronounced than the water level rise at the outer edge. The flow rate Q = 700 m3/s is characteristic because in the downstream flow section, the flow regime changes into supercritical flow. The surface shape along the edges is of a somewhat different character than that with an unchanged regime (Figure 5). Because of inertia at increased velocity under the supercritical regime the surface at the end of the curve slightly fluctuates in the transversal direction. In the initial part of the curve the water surface resembles that in the subcritical regime, i.e. the water level rise at the outer edge is smaller than the water surface lowering at the inner edge. By the end of the curve this difference almost approaches zero, while the water surface rise is almost equal to the lowering at the mean water surface level in the straight part of the downstream section. By increasing the discharge and velocity, the flow regime at the upstream boundary changes as well. At Q = 1100 m3/s, supercritical flow conditions occur throughout the channel length. The calculation results are shown in Fig. 6. At the start of the curve the impact of flow inertia is evident, thus the surface slightly fluctuates. In the remaining part of the curve the surfaces at the edges are almost horizontal, with a considerable water surface superelevation compared to the mean depth at the outer edge and a water surface lowering by almost a half at the inner edge of the curve. Under these conditions, surface fluctuations at the outflow from the curve occur, which is even more pronounced at high velocities. Because of energy losses, which are higher on examples with the higher velocities, and because of the supercritical flow conditions in the channel, the surface in the part of the channel following the curve is higher than that in the straight upstream section. Fig. 7 shows an axonometric illustration of the calculated water surface on the model of the curve at Q = 1100 m3/s. Analysis of Flow In a Curved Channel Using the Curvilinear Orthogonal Numerical Mesh 541 Strojniski vestnik - Journal of Mechanical Engineering 64(2018)9, 536-542 5 CONCLUSIONS The analysis of the flow in a semi-circular curve using the 2D PCFLOW2D-ORTHOCURVE model, where the influence of channel walls was considered negligible, revealed the following: • the form of the free-surface flow in the curved section depends on the flow regime (subcritical, transitional or supercritical); • in subcritical flow conditions, the surface level rise at the outer edge is smaller by one half than the water level lowering at the inside edge; • in supercritical flow conditions, the rise in the surface level at the outer edge is almost two times larger than the lowering of the surface level at the inner edge; • at flow through a bend, where the regime transitions from subcritical to supercritical, the water surfaces at both curve boundaries are almost parallel and follow the linear change in the average depth; • by taking into account the average water flow velocity across the curve, which is the result of the mathematical model, we obtained a very good fit between the theoretically calculated transverse slope of the free surface and the surface slope in the mathematical model. 6 NOMENCLATURE £, n curvilinear system coordinates, I, J succession numbers of numeric cel in the £ and n directions, m(, Lame's coefficients in the £ and n directions t time, [s] h water depth, [m] u^ velocity component in the £ direction, [m/s] velocity component in the n direction, [m/s] vn zb bottom level, [m] n Manning's friction coefficient, [s/m1/3] g acceleration due to gravity, [m/s2] v laminar coefficient of viscosity, [m2/s] vt turbulent coefficient of viscosity, [m2/s] vef effective coefficient of viscosity, [m2/s] k turbulent kinetic energy per unit mass, [m2/s2] s degree of k dissipation, [m2/s2] Pfo>, Psv source terms due to the bottom friction for the k and e respectively, [m2/s3] Cf Chezy friction coefficient, [m1/2/s] v average flow velocity in the channel, [m/s] Rout radii of the outer edge of the curved channel, [m] Rins radii of the inner edge of the curved channel, [m] 7 REFERENCES [1] Roache, P.J. (1972). Computational Fluid Dynamics, Hermosa Publishers, Socorro. [2] Cetina, M. (1988). Mathematical Modelling of Two-Dimensional Turbulent Flows, MSc Thesis, University of Ljubljana, Faculty of Civil and Geodetic Engineering, Ljubljana. (in Slovene) [3] Tang, X., Knight, D.W. (2015). The lateral distribution of depth-averaged velocity in a channel flow bend. Journal of Hydro-environment Research, vol. 9, no. 4, p. 532-541, D0l:10.1016/j.jher.2014.11.004. [4] Gholam, A., Bonakdari, H., Zaji, A.F., Akhtari, A.A., Khodashenas, S.R. (2015). Predicting the velocity field in a 90° Open channel bend using a gene expression programming model. Flow Measurement and Instrumentation, vol. 46, part A, p. 189-192, D0I:10.1016/j.flowmeasinst.2015.10.006. [5] Sui, B., Huang, S. (2017). Numerical analysis of flow separation zone in a confluent meander bend channel, Journal of Hydrodynamics, vol. 29, p. 716-723, DOI: 10.1016/ S1001-6058(16)60783-7. [6] Ivanovic, D.M. (1963). Vector analysis, Scientific Book, Beograd. (in Serbian) [7] Rodi, W. (1993). Turbulence Models and Their Application in Hydraulics, 3rd edition, A.A. Bakema, Rotterdam, Brookfield. [8] Launder, B.E., Spalding, D.B. (1972). Lectures in Mathematical Models of Turbulence, Academic Press, London/New York. [9] Patankar, S.V. (1980). Numerical Heat Transfer and Fluid Flow, McGraw-Hill Book Company, New York, D0I:10.1201/9781482234213. [10] Bridge, J.S., Jarvis, J. (1976). Flow and sedimentary processes in the meandring river South Esk, Glen Clova, Scotland. Earth Surface Processes, vol. 1, no. 4, p. 303-336, D0I:10.1002/ esp.3290010402. [11] Bhowmik, N.G. (1979). Hydraulics of Flow in the Kaskaskia River, Illionis, Report of Investigation 91, Illionis State Water Survey, Urbana, Illionis. m 542 Krzyk, M. - Cetina, M.