Numerical Study of the Effects of an Obstacle and its Position in a Curved Channel in a Lock-Exchange Flow

Researchers have recently shown interest in complex flow patterns, especially the effects of secondary flows, in a curved channel. This paper studied the obstacle effects in a gravity current in a channel with a symmetrical 120° bend using the OpenFoam toolbox and the Realizable k – ε turbulence model for simulation. The models studied included no-obstacle curved channel, curved channel with obstacle in 30° position, curved channel with obstacle in 60° position and curved channel with obstacle in 60° position with increased radius. Results showed that the obstacle directed the concentration towards the banks with its maximum value tending from the outer to the inner bank, especially in the tail. Although the post-obstacle head height did not change, that of the tail did (fell); the tail longitudinal velocity was maximized near the channel bed in areas far from the obstacle, and in the outer bank in areas near it. The secondary flow was so reduced that its lowest and most different pattern was observed around the obstacle. In displacing the latter, if the front was at a certain distance from it, the secondary flow did not change much, but if it was at the channel end, the post-obstacle secondary flow increased as the obstacle neared the lock.

In displacing the latter, if the front was at a certain distance from it, the secondary flow did not change much, but if it was at the channel end, the post-obstacle secondary flow increased as the obstacle neared the lock.

INTRODUCTION
Gravity currents are two-phase flows formed due to the density difference between two fluids because of engineering reasons, natural phenomena (avalanche, sea breeze, sandstorm, etc.) or human interference (discharging pollutants into rivers and dams' reservoirs, tunnel fires, etc.) (Simpson, 1999;Meiburg et al., 2015). Such factors as temperature variations, chemical compounds and solid particles can cause density differences, and since solid particles are difficult to simulate, use can be made of high-density salt water solutions for related studies (Ooi et al., 2007;Paik et al., 2009). Ungarish (2009) describes formation conditions of the gravity current (c > a) under Lock-exchange topics two dimensionally at position zero ( 0 ) ( Fig. 1). According to Fig.1, the interface under consideration in the present configuration has two components. The first is a vertical front (also referred to as nose) and the second, and major, part of the interface is the inclined, sometimes piecewise inclined-horizontal, surface (Ungarish, 2009).
With the gate removal at 0 , denser current begins to have a sharp front somewhere near the floor wall, forms along the channel base and becomes more along the flow (Hacker et al., 1996). released from behind a lock (Ungarish, 2009) Hacker et al. (1996) was the first to do gravity-current lab studies in high-Reynolds Lockexchange topics where the two fluids' density difference was 12%, the denser fluid's volume was 0.0328 m3, and Reynolds number was 19700. Paik et al. (2009), Ooi et al. (2007) and Mahdinia et al. (2012) are among many researchers who based their works on this experiment; the present study has used the configuration of the Hacker test. Mahdinia et al. (2012), who simulated, for the first time, the gravity-current in bend under Lock-exchange, believed that such currents could occur in rivers, at the entrance of reservoirs or in any other natural environment. Following the latter work, the present study has used an obstacle in a curved channel to study gravity-currents' hydraulic properties.
Examples of comprehensive studies on gravity currents in direct channels under Lockexchange conditions include those of Firoozabadi et al. (2009) who simulated the velocity/density of channels' inflows, Lam et al. (2018) and  who simulated lock dimensions variations and density difference of two fluids, Zhao et al. (2018), Lombardi et al. (2015) and  who simulated the channel gradient flow hydrodynamics, and so on; Longo's (2015) too is among limited studies on the effects of the cross-section type on the flow. The modified RANS model introduced by Paik et al. (2009) observed near-wall flows well. Ooi's (2007) study results showed that it is better, in numerical studies, to use the no-slip boundary condition for walls so that density lines could have better validations. Among those who studied obstacle and roughness effects on the gravity current in direct channels, Tokyay et al. (2012) investigated the roughness effects on the channel bed wall shear stress in one study, and the flow regime due to the obstacle height variations in another (Tokyay and Constantinescu, 2015).
Lock-exchange experiments were used to examine the temporal entrainment properties and mixing processes of sediment-laden turbidity currents interacting with a rectangular obstacle.
During the slumping process, currents of varying densities were studied over smooth and rough substrates (Wilson et al. 2017).
The entrainment and mixing dynamics of a gravity current produced by a lock-exchange mechanism were investigated experimentally by simultaneously measuring velocity and density fields. Near the front of the gravity present, the values of entrainment coefficient were found to be the largest (Balasubramanian and Zhong, 2018).
To find out how three-dimensional canopy geometry affects the front propagation of an incoming gravity current under a given initial forcing, large-eddy simulations were conducted by Zhou and Venayagamoorthy, 2020. It was shown that when applied to buoyancy-driven flows, the traditional geometrical parameter of submerged canopies in constant-density flows is misleading.
In a general geometry-based grouping, Brice et al. (1974) divides meanders into simple symmetric, simple asymmetries, symmetric combination and asymmetric combination. Based on his studies and morphological analyses, Magdaleno et al. (2011) believes that rivers naturally tend to have a symmetrical form.
Bend-geometry open-channel fellow studies are many; Abad and Garcia (2009a) measurements in an asymmetric lab channel show that the maximum velocity occurs at the flow surface in the inner bank of the bend, and at the channel bed in the outer bank of the bend; in the channel width, the highest velocity difference is observed at the bend vertex. Next in his research, Abad added the erodible bed to the test conditions and studied movements of sediment particles (Abad and Garcia, 2009b). In a numerical study, Farshi et al (2018) presented relationships for secondary flows in 120° symmetric bends in a flood channel and showed, by the flow pattern, that the secondary current increased after the bend vertex and reached its maximum near the outer bank; the pattern was, then, repeated from this vertex to the next.
Flow patterns in many studies (Kassem et al., 2004;Islam and Imran., 2008;Ezz et al., 2013 andSumner et al., 2014) on submarine density currents, formed like a bend under Coriolis force effects of sinusoidal currents, have similarities and differences compared to those in open channel bends.
Since Mahdinia et al. (2012) is the only one who studied gravity currents in curved channels under Lock-exchange conditions and showed that mid-flow (somewhere between head and tail) has features similar to submarine density currents, his study is of particular importance to the present work. For instance, the flow, in both cases, increases in the outer bank and secondary currents occur at the bend vertex near the canal bed. In his study, while there is almost no secondary flow at the head, the one near the bed has a direction toward the outer bank. Near the tail, the pattern is different due, maybe, to the reduced density/velocity that has no definite direction across the channel width; the highest longitudinal velocity too occurs in the bend's outer bank. He, then, studied the effects of the bend radius and showed that its decrease increased the secondary flow because the centrifugal force and bend radius are inversely related.
According to the literature, as no studies have been done so far on the obstacle effects in curved channels under gravity currents, this paper has used an obstacle in a 120° curved channel to study the hydraulic features (e.g., gravity and velocity) of the flow and location effects of the obstacle by some numerical modellings.

NUMERICAL METHODS
For modeling, this research has used OpenFOAM, which is an open-source and uses the C++ language to solve CFD problems, and the twoliqiudMixingFoam solver, which is used in Lock exchange cases and is defined for 2-phase currents where both phases are incompressible and have different densities.
For velocity and coupling it with pressure, use has been made of the noSlip boundary condition [Ooi et al., 2007] and the PIMPLE algorithm, respectively, and the boundary condition for pressure at all boundaries is the fixedFluxPressure. Wall functions are used for , k, and ε.

GOVERNING EQUATIONS AND TURBULENCE MODELING
The twoLiquidMixingFoam is an open-source solver form miscible and incompressible fluid. The solver finds the volume fraction αo and α1= 1-αo, and then the velocity U by alpha diffusion equation, the continuity, and the momentum equations (Xie and Chu, 2019): Where Dab is the molecular diffusivity and Sc is the turbulent Schmidth number. More details of this solver and its equations in numerical simulation of related problems can be found in (Zhang et al, 2016 andXie andChu, 2019). where is the earth gravitational acceleration, is the fluid density inside the lock, is that in the ambient fluids and time 0 equals ℎ/ .
From among RANS models, this study has selected the realizable k-ε turbulence model wherein the kinetic energy (k) and energy dissipation (ε) are as follows (Shaheed et al., 2019): Where (Eq. 6) increases the kinetic energy due to mean velocity gradients and is found as follows (Shaheed et al., 2019): Where is the strain-rate tensor, and are diffusivity effects for K and ε, 2 is taken to be 1.9 (constant), 1 is not constant and is a turbulent Prandtl number assumed to be 1.2. An important feature of the realizable k -ε model is that it does not assume to be constant; it calculates it for the turbulence viscosity (Shaheed et al., 2019).
Wall functions are boundary conditions and in OpenFOAM all of them inherit from the abstract class FvPatchField (Fangqing, 2016). In this paper epsilonWallFunctions, kqRWallFunctions, and nutWallFunctions are used.

PROBLEM SETUP
The B test Hacker et al., (1996) was modeled in a direct channel and called SC to validate the present study. The length (L), width (b) and height (H) of the model geometry were 3480, 200 and 500 mm, respectively, and those of the Lock at its beginning (x 0 , b and h - Fig. 2a) were 400, 200, and 400 mm, respectively. In this paper, the fluids' density difference in the SC and research models ( Fig. 2b) was 12%. The latter contained a symmetrical 120° bend in the middle of the obstacle channel. The bend radius from O to the geometry center (where z = 0) was considered as 1 m to ensure the centrifugal force efficacy from center to the channel. In   (Hacker, 1996) used for validation, b) Curved channel used in present study's models In all models, obstacle dimensions are fixed, its height is 0.3h (Tokyay and Constantinescu, 2015) and a hypothetical line in the middle of its width makes a 180° angle with the bend radius; Table (1) summarizes all the geometric parameters. In this paper, BOSD60 is the symbol for bending channel with obstacle in section 60° and BC means a no-obstacle curved channel. BOSD30 (obstacle in section 30°) has been modeled to study the effects of the obstacle distance from the lock because the bend is symmetrical and its radius is constant along the channel length. Finally, the bending channel with obstacle in section 60° with increased radius has been modeled with an obstacle at 60° position (BOSD60IR) and the results have been compared with those of BOSD60. Since the lab modeling of a meandering river with varying-radius curves along the channel length is quite costly, this study has assumed an asymmetric bend made with different-radius circular pieces to study the effects of the obstacle positions on an asymmetric bend by comparing BOSD60IR and BOSD60 models.
The t/t 0 in Table 1 is the run time in each model wherein the flow front reaches the channel end. During the solution, Δt is a variable (in present models, it varies in a 0.0005-0.009" range) and the Courant maximum is used in the solver's "controlDict" file to ensure its value will not exceed.

SOLVER VALIDATION
The lab channel Hacker et al. (1996) is straight (Fig. 2a) and its dimensions are shown in same in both curved and straight channel models. Organized meshes generated in the validation model by the "blockMesh" program, get finer towards the channel bed by a factor of 4. To ensure the mesh quality, the model was run with 1, 2 and 3 million cells in all of which the results, including the data of the flow front position versus time and the general flow shape were equal, but the 3-million cell computational network (Fig. 3) performed better in presenting the mixing of vortices.  (Simpson, 1972) because of the numerical models' boundary conditions. Although results of the LES turbulence model simulation Ooi et al. (2006) were similar, he attributed it to lab errors when removing the gate in the lab. For more details see Paik et al. (2009).
Solving micro-current vortices by Mahdinia's (2012) LES model cannot be a defect and selecting a RANS turbulence model can be a good option here because of the research objectives and the models' high computational costs. As shown in Fig. (3d), the present model has nicely modeled the flow front collapse near the channel bed. In Mahdinia's (2012) simulation model, the tail of the dense current is at a lower height than its head. It is worth noting that in the lab model, the flow tail vortex is slightly further back along the channel length due, maybe, to the lab gate removal error. In this figure, the bases to compare the tail heights with those of the lab model are points where these heights are the maximum.

GRID SENSITIVITY STUDY
Since the flow pattern varies in three directions in curved channels causing complicated computational network and increased related costs due to the obstacle presence in the channel, high Reynolds number and 2-phase problem, this research has performed grid sensitivity study for BC, BOSD30, BOSD60 and BOSD60IR models.
Grids have been generated by the blockMesh program that produces organized background grids. Here, the latter are three degrees finer towards the channel bed wall and are completed by the snappyHexMesh program that creates the final computational grid by generating different shape/size cells.
During processing, when the flow front was around the obstacle, the reduced value of Δt caused the solution to diverge due, maybe, to the intricate-shape grids around the obstacle. In other words, the snappyHexMesh program generates, for the solver, intricate-and (sometimes) incomprehensible-shape cells around the obstacle due to the channel curvature making convergence difficult; this problem will be solved with some changes in the background mesh.
All the research models of the present study have been solved by keeping a maximum value of 0.3 for the Courant number. Table (2) lists the specifications of the grid independence runs. The front position versus time was the first parameter studied and the results were almost equal in all different-volume models, i.e. closer results even in smaller grid volumes. Next, such parameters as pressure, concentration and velocity in different directions were investigated where the highest difference in results was observed in concentration contours, especially in the front position. As shown in Fig. 5, results are, in general, quite close with no major difference among the three computational grids. However, since the highest number of cells in each model had closer results than the two models with the lowest number of cells, the highest number of cells in each model was selected as the final grid. It is worth mentioning that  In the present study, the sensitivity data of the velocity along the channel width, especially near the bed, are very close, but since its main objective is the analysis of the secondary flow, its results are shown in Fig. 6; results in the longitudinal velocity and pressure parameters are almost the same. normalized potential energy, kinetic energy, and dissipation rate of gravity currents before, during, and after the impact. They identified various stages of propagation, including the impact, transient, and quasi-steady stages.

DISCUSSION AND RESULTS
This section is to study the results of the no-obstacle curved channel model (BC) and that with obstacle at 60° position (BOSD60), and compare the results to check the obstacle effects.
Comparing BOSD60 and BOSD60IR will enable the studying of the increased radius effects (with obstacle in channel), and comparing BOSD30 and BOSD60 (when obstacle nears the lock) will enable the studying of the obstacle displacement effects on the flow.

Flow Front Position
In curved channels, = where is the front position, r is bend radius and is front position in the curved channel. As shown in Fig. 7, the front in the BC model moves at a constant speed and the line gradient is constant. In BOSD60, BOSD30 and BOSD60IR models, the gradient is equal to that in the BC model in the pre-and post-obstacle modes (where the front is not affected by the obstacle). Under the obstacle effect, the gradient is reduced and reaches zero (m = 0) where the nose falls ahead in the outer bank. This means sometimes, while the flow maintains its dynamicity, the front does not fall ahead; rather, it is the nose that does so.

Fig. 7. Front position versus time in research models
Discussed in detail by Mahdinia et al. (2012), in the BC model, the nose falls ahead in the inner bank in the first half of the channel and then tends to do so in the outer bank after the obstacle if there is one. Since this is the case in BOSD30 model, it can be concluded that the mere obstacle causes the nose to fall ahead in the outer bank in some part of its length and the obstacle position at the bend is irrelevant; time intervals when the gradient is zero are different in each model. These results are also in a good agreement with Wu and Ouyang (2020) numerical research that investigated flow morphology in bottom-propagating gravity currents over immersed obstacles in a direct channel. Leong et al. (2006) also showed the similar results in a direct channel and mentioned that the velocity of the advancing front fluctuated as the current moved along.

Obstacle in the Channel Bend
The lock-released gravity current is a time-dependent deformable "obstacle," whose shape interacts with the waves it produces in the ambient and analysis and investigation of the field In a direct channel, entrainment was shown to increase downstream of the obstacle when the obstacle was present as shown by Wilson et al (2017). In the same way, our present results show this flow feature in a bend. Fig. 8 shows contours of the concentration in the channel cross section; in 8a (5/6 position), the head height has found its shape, is almost equal in both models, is not affected by the obstacle, concentration is highly reduced due to the obstacle and is almost maximum near the channel bed on both inner and outer banks. In 8b (4/6 position), obstacles reduce the tail height and maximum concentration is in the inner bank. In the BC model, concentration increases in the tail in the outer bank and reaches its maximum at the channel bed by moving away from the flow front; this maximum is usually in the inner bank if there is an obstacle in different positions. The post-obstacle tail height is almost balanced across the channel due, maybe, to the low flow velocity (8c). Being almost equal across the channel, the pre-obstacle tail height is maximized in the outer bank as it approaches the lock (8d).
A similar conclusion was reached by Wilson et al. (2018) who showed that in a turbidity currents in a direct channel, interacting with an obstacle, the obstacle reduce the maximum velocity and turbulence intensity downstream of it. Fig. 9 shows contours of the longitudinal velocity across the channel. The flow front is at channel position 11/12 and cross-section positions are as in Fig. 8. The longitudinal velocity direction is negative at the flow head, but changes by moving away from the head to the tail.
An obstacle causes this to happen sooner; in 9a, the maximum longitudinal velocity at the flow head occurs in the outer bank nearly around y/h = 0.3. In 9b, the velocity decrease is quite clear in the tail reaching its maximum near the channel bed. In 9c, the post-obstacle longitudinal velocity distribution is irregular due, maybe, to the vortices and flow turbulence. In 9d, the preobstacle longitudinal velocity direction is negative in the channel bed due, clearly, to a direction change. As shown, the longitudinal velocity is maximized in the inner bank with a quite opposite flow pattern in the channel upper half. In the BC model, the longitudinal velocity is maximized before the vertex in the outer bank and is balanced in the channel width after the vertex. especially before, the obstacle (Fig. 10d). According to 10a and 10b, both in post-obstacle positions, the secondary flow direction is similar to that in the BC model with the difference that it decreases in the outer bank and Increases in the inner one. In 10c and 10d, the flow has a different pattern near the obstacle due to vortices and flow turbulence; in the channel bed, the secondary flow direction is from the outer bank to the inner, changes above it and again from the outer to the inner bank in the channel upper half. Near the lock, the direction is totally reverse; from outer to inner near the channel bed, and from inner to outer near the interface.

Obstacle effects when channel radius increases
As mentioned before, with an obstacle, an increase in the radius will not change the gradient of the flow front motion graph. Mahdinia et al. (2012) had already shown that in a no-obstacle channel, a change in the radius would not change the front position versus time graph gradient and concluded that the secondary flow would decrease with an increase in the radius. Despite the presence of an obstacle, the present study too has had similar flow pattern results to be described in the following parts. Figure 11 shows the BOSD60IR model contours in sections 5/6 and 4/6 when the flow front is at position 11/12 to be compared with those of the BOSD60 model. As shown, an increase in the bend radius in the former reduces the centrifugal force and, hence, the velocity compared to the latter. With this increase, the channel-across flow components tend to increase in the banks as well as in other positions. Fig. 11a shows that the heights of flow tail and head have decreased.

Pre-vertex obstacle effects
According to Mahdinia et al. (2012), the initial-moment maximum longitudinal flow velocity is observed in the inner bank where the flow nose is further ahead in the first half of the channel. Since these results have also been observed in the present study, the obstacle effects, in the BOSD30 model, have also been studied in this half the results of which are shown in Fig. 12.  increased. Observations (Fig. 12c) show that the secondary flow amount/pattern does not change much compared to the BOSD60 model.

CONCLUSIONS
This paper has investigated the effects of the obstacle presence/position in a gravity-flow in a symmetrical curved channel. Settings presented for the numerical model are for good validation results. The mesh sensitivity analyses results show that the transverse velocity parameter is less sensitive to concentration than to the computational network. At the channel beginning, the front moves at a constant speed, and the nose is farther in the inner bank. After hitting the obstacle, its headway speed is reduced, but in post-obstacle conditions, while maintaining its dynamism, the flow front does not move along the channel for a while when its nose gets ahead in the outer bank after which it continues moving at a constant speed.
With an obstacle, the flow-head density is reduced, but the head height does not change.
The density is maximized in most channel positions in the inner bank, but in the no-obstacle case, the density increases in the outer bank, especially in the tail. As the radius increases, the head and tail heights of the flow are reduced and the density tends towards the banks. With an obstacle, the longitudinal velocity is highly reduced in the tail and has a maximum value in the channel bed; in the outer bank, it is maximized around the obstacle and in the inner bank it does so far from it.
In general, obstacles reduce the transverse velocity along the channel, but cause different secondary flow patterns at different positions. The highest reduction is observed around, especially before the obstacle. Results show that displacing obstacles enables the secondary flow and its min/max values to be controlled and managed along the bend. According to the results of the present study, when the flow front is at the channel end and the obstacle is near the lock, there will be more secondary flow in the lower half of the channel (where the flow direction is toward the outer bank) and its maximum value occurs in the outer bank, but when it is at a certain distance from the obstacle, the secondary flow pattern/amount does not change much.

DATA AVAILABILITY STATEMENTS
The data that supports the findings of this study are available within the article.