Numerical Simulation of Asymmetric Merging Flow in a Rectangular Channel ()
Received 8 November 2015; accepted 19 April 2016; published 22 April 2016
1. Introduction
The mixing of two different fluids flows and concentration at a junction of two channels constitutes the field of concentration. It is also called merging flow. It has great interest in physical application point of view. For example [1] in his consideration of problem of the mixing of two rivers assumed that the concentration field was completely unmixed at the point of confluence. It has a lot of use in physiological flows, internal machinery dynamics, lakes, estuaries, and rivers, see for example: [2] [3] . This type of geometry is used in quadrupole magnetic cell sorter [4] and is also dominant in respiratory flow [5] [6] .
In this paper, we consider the steady, two-dimensional flow of incompressible viscous fluid that constitutes the merging flow in which two oncoming channel flows are forced to join into one. The geometry of problem consists of an infinite long straight channel |y¢|£h, −¥ < x¢ < ¥ of width 2h containing a disposed straight splitter plate placed at y¢ = ah, −1 < α < 1 within the channel as shown in Figure 1.
Here (x¢, y¢) represent dimensional rectangular coordinates. The position of splitter is controlled by a splitter position parameter a. If a = 0 then the splitter plate will coincide with the x-axis while the plate will lie above and below the x-axis if a > 0 and a < 0 respectively. The inlet and outlet boundary flow assumed to be same as the plane Poiseuille flow appropriate to the single width channel for downstream while double channels for upstream of splitter. Although for symmetric case, the various attempts had been made by several authors for example; [2] [7] - [13] yet owing to the complexity for asymmetric flow which resembles actually to the real life problems, unfortunately due attention was not given.
Obviously a question arises how does our study differ with previous attempts? Following are few facts which cover the answer of this question.
We consider non-uniform rectangular meshing while they used mostly square grid.
We use higher order scheme while previously the authors used mostly 2nd order scheme.
We consider asymmetric flow while previously the symmetric flow cases were examined.
We study the low, moderate, and high Reynolds numbers flow while previous attempts are mostly based either on low, moderate or high Reynolds number see for example [7] .
The main aim of this attempt is to explore all the aspects of the flow behaviour and examine the understanding of the complicated features such as boundary layer separation, growth and decay of vortices, reattachment, trailing-edge and leading-edge properties and variation of skin friction on the solid walls and splitter plate for asymmetric case. In addition to this, the flow behaviour is also simulated for various values of Reynolds number R and splitter position parameter a. We also observe how asymmetric study differs from the previously attempted symmetric cases. We examine the variation of eddies when the splitter plate is moved upward or downward within the channel. The extrapolated results obtained are simulated graphically and are compared with those of previously attempted symmetric cases, like [2] . These are given in Section 5. The formulation and basic flow analysis are presented in Section 2 while Sections 3 and 4 contain the detail of numerical scheme so adopted and computational procedure used respectively.
2. Basic Analysis
The continuity and the Navier-Stokes equations for steady flow of incompressible fluid in dimensional form, in the absence of body force are: [2]
(1)
(2)
where 
is the fluid velocity vector, r is the constant density, m is the viscosity coefficient, and 
is the pressure.
The mechanics of the problem under consideration can briefly be stated as; flow is through an infinite length uniform channel with finite width parallel walls. The upstream region is divided into two channels with the help of semi-infinite splitter plate while downstream is flow between two plates without any splitter plate. The walls and splitter are parallel to each other and x-axis. The upstream of splitter plate is separated by 2ah while downstream is separated by 2h, as shown in Figure 1. Moreover, the flow is being considered as steady, two-dimen- sional, laminar, and asymmetric about x-axis in general. The flow is generated by the uniform pressure gradient analogous to that of plane Poiseuille flow on inlet and outlet flow to preserve the continuity.
On deforming Equations (1)-(2) into vorticity transport equation in dimensionless form by normalizing velocity, space coordinates and vorticity vector by U (main stream velocity), h, and h/U respectively, we will obtain:
(3)
(4)
where 
is the vorticity vector, R the Reynolds number which is equal to the ratio Uh/n and n being the kinematic viscosity.
If we introduce 
as velocity components in terms of stream function y relation into Equations (3)-(4), we get,
(5)
where
(6)
For the geometry of Figure 1 the boundary conditions for these equations are given as follows:
(a) No slip at the walls;
i.e.
, on all the solid walls and the splitter plate. (7a)
(b)
(7b)
(c)
(7c)
3. Numerical Scheme
To overcome the problem of rate of convergence and stability of the numerical scheme, we adopt the numerical scheme which consists of the following two methods.
The grid size along x and y directions, will be taken by H, K1 respectively, and around a typical internal grid point 
we adopt the convention that quantities at![]()
, ![]()
, ![]()
, ![]()
, ![]()
, ![]()
, ![]()
, ![]()
, and ![]()
are represented by the subscripts 0, 1, 2, 3, 4, 5, 6, 7 and 8 respectively.
3.1. Fourth order Finite-Difference Method
If we employ the method similar to [14] , then governing Equations (5)-(6) are deformed as:
![]()
(8)
![]()
(9)
where
![]()
and
![]()
.
3.2. Special Finite-Difference Method
In order to approximate linear partial differential equation (6), we employ second order standard central difference formulation at point “0”, we obtain: [15]
![]()
(10)
To deal with Equation (5), we split it equation into following two equations as:
![]()
(11)
![]()
(12)
where A is an arbitrary function of x and y while ![]()
Equation (11) is then approximated along the grid line y = y0, over the range![]()
, by applying a local transformation for E, given as follows:
![]()
(13)
where![]()
.
On the introduction of Equation (13), into Equation (11) and approximation of the derivatives involved by 2nd order central differences, we get
![]()
(14)
Similarly, Equation (12) is approximated along the grid line x = x0, over the range![]()
, by applying a local transformation for E, given as follows:
![]()
(15)
where![]()
.
On using Equation (15) into Equation (12), after the similar approximations as above, we shall obtain
![]()
(16)
On eliminating A0 between Equations (14) and (16), we have
![]()
(17)
In order to express l and m back in terms of E, it can be done from the definitions. It is found that
![]()
(18)
where i = 1, 3 and j = 2, 4.
Now, we expand above exponential in powers of their arguments and keep the truncation error of order H4 and![]()
. After some simplifications under above arguments and using the Taylor’s theorem during calculations, we obtain
![]()
(19)
![]()
(20)
On using Equations (19) and (20) into Equation (17), we get
![]()
(21)
These approximations give rise to an associated matrix that is always diagonally dominant. The present formulation was found to work well for all values of Reynolds number R, whereas results could only be obtained for relatively small values of R by using standard central-difference approximations to Equations (5) and (6). Finally the condition for E is required at grid points on the solid walls. Here we use the same approximation as that is originally due to [16] and is given by
![]()
(22)
where the subscripts w1, w2 and w3 represent for the value at the approximate boundary points on upper, lower, and splitter plate respectively while the subscripts a, b and c signify to the internal grid points most immediate to w1, w2 and w3 respectively.
4. Computational Procedure
This numerical experiment consists of two steps. In the first step, the Equations (10) and (21) are solved iteratively with respect to the boundary conditions (7) and (22) by SOR-method [17] until the convergence is achieved according to the criterion:
![]()
(23)
In the second step, Equations (8) and (9) are solved numerically with respect to boundary conditions (7) and (22) by SOR-method on using the data obtained from step 1 as an initial estimation. This procedure is repeated until the convergence is attained according to the above criterion (23). The solutions obtained are further extrapolated by Richardson’s extrapolation technique [18] .
5. Results and Discussions
Two parameters R and α have a vital role in the general fluid motion. The ranges for R and a are 0 £ R £ 105 and −1 < a < 1 but, for the purpose of presentation and discussion, we consider R Î {1, 10, 100, 500, 1000, 100,000} while a Î {−0.5, 0, 0.2, 0.5} throughout the study. The numerical calculations are carried out on grid sizes H Î {1/20, 1/40, 1/80, 1/160, 1/320} and K1 Î {1/60, 1/120, 1/240, 1/480} and the presented results are also extrapolated using the Richardson Extrapolation method. In all the cases the infinite length of channel is numerically approximated and is taken as four.
The streamlines for aforementioned R and a are presented in Figures 2-6. These figures depict that the vortices are produced in the vicinity of the trailing edge of the splitter plate. The size of eddies increases as a increases or if the splitter plate moves upward towards upper plate as shown in Figures 2-6. An interesting situation exists for R = 100,000 and a = 0.5 where whole fluid below the splitter plate is rotating clockwise as shown in Figure 6(d).
The trailing edge of the splitter plate needs special attention due to the infinite or singular behaviour of vorticity at this point. To deal with this singular point two special treatments are being adopted during calculations for all values of R and a: (i) the method, which is very similar to [7] [8] , is used; (ii) the refine meshing is taken in the vicinity of the splitter plate and specially its trailing edge, as shown in Figure 7. More clustering near the singular point is needed as R increases.
![]()
Figure 2. Streamlines for L' = 4, R = 1, (a) a = −0.5; (b) a = 0; (c) a = 0.2; (d) a = 0.5.
![]()
Figure 3. Streamlines for L' = 4, R = 100, (a) a = −0.5; (b) a = 0; (c) a = 0.2; (d) a = 0.5.
![]()
Figure 4. Streamlines for L' = 4, R = 500, (a) a = −0.5; (b) a = 0; (c) a = 0.2; (d) a = 0.5.
Here it is necessary to comment on Figure 2 of [2] that our results for ua (velocity at centre line) for symmetric flow match well for low Reynolds number and slightly differ for high R near downstream end of the splitter (in our case it is x = 2) where the matching boundary condition
![]()
Figure 5. Streamlines for L' = 4, R = 1000, (a) a = −0.5; (b) a = 0; (c) a = 0.2; (d) a = 0.5.
![]()
Figure 6. Streamlines for L' = 4, R = 105, (a) a = −0.5; (b) a = 0; (c) a = 0.2; (d) a = 0.5.
![]()
is applicable and on the velocity will be![]()
. For symmetric case, the velocity on the central line (ucl in
![]()
Figure 7. Clustering in the vicinity of splitter plate for different values of a considered.
[2] and ua: a = 0 in our case) will be obviously 1.5 for all R. This fact is missing in [2] and this deviation is corrected here and is displayed in Figure 8(b). In Figure 8, the distribution for velocity ua along the splitter plate is given for various values of Reynolds number and splitter position parameter. When a = 0 it is observed that the velocity along the splitter increases in downstream of trailing edge of the splitter as we move away from it and curvature of velocity profile increases as R increases. This situation will exist till R = 1000 but it reverses dramatically when R > 1000 as shown in Figure 6(b). It is noted that the optimum value of velocity lies as x increases far from the trailing edge of the splitter, for symmetric case but it is not so for asymmetric case as shown in Figures 6(a)-(d). Moreover, the velocity increases at once after passing the trailing edge for asymmetric case which is clear in Figure 6.
The vortex effects (or in other words the variation of the skin friction) on the upper and lower plates of the channel are also examined here and are displayed in Figure 9 and Figure 10 for various values of aforementioned values of R and a. One effect can be noted that the vorticity increases as R increases and it decreases as a increases. These effects occur on both the walls of the channel as shown in Figure 9 and Figure 10. Figures 11(a)-(d) depict the variation of dimensionless vorticity on the splitter for various values of R and a, so considered. Figure 11(b) belongs to symmetric flow case and it coincides and compares well with Figure 1 of [2] . This verifies our results. Moreover, the notable thing is that the decay of vorticity to its downstream value takes place over an increasingly longer scale of x as R increases for symmetric case but it is not so for asymmetric one, as shown in Figure 11.
The coefficient of skin friction Cf has vital role in the motion of the fluids and loading-effects on the solid walls. It generally depends upon vorticity, Reynolds number, viscosity, and other geometric aspects. This effect has also been studied here for various values of R and a so considered and its variation with respect to x is also displayed graphically in Figure 12. This figure indicates that the shear stress on the splitter increases as x in- creases and as R decreases. The maximum shear stress will occur at the singular point, which resembles with the physical reality. Figure 12 represents the comparison of our numerical solutions with analytical estimation of [2] . The curves for analytical results are displayed as dotted line. The comparison is seemed to be satisfactory. Moreover, it is observed that; as the splitter plate is moved upward then the shear stress on the splitter plate
![]()
Figure 8. The variation of ucl for various values of R and (a) a = −0.5; (b) a = 0; (c) a = 0.2; (d) a = 0.5.
![]()
Figure 9. The variation of wall vorticity on upper plate of channel for various values of R and (a) a = −0.5; (b) a = 0; (c) a = 0.2; (d) a = 0.5.
increases as shown in Figure 12.
For comparison purpose, Figure 13 indicates the variation of vorticity on the splitter for symmetric case
![]()
Figure 10. The variation of wall vorticity on lower plate of channel for various values of R and (a) a = −0.5; (b) a = 0, (c) a = 0.2; (d) a = 0.5.
![]()
Figure 11. The variation of wall vorticity on splitter of channel for various values of R and (a) a = −0.5; (b) a = 0; (c) a = 0.2; (d) a = 0.5.
![]()
Figure 12. Comparison of present numerical solution for shear stress on the splitter for symmetric case with analytic solution of [1] (shown here dotted) when (a) R = 1; (b) R = 100; (c) R = 500; (d) R = 1000.
![]()
Figure 13. Comparison of present results with analytical results of Badr [1] (here shown as *) for the variation of wall vorticity on splitter of channel for various values of R and (a) a = −0.5; (b) a = 0; (c) a = 0.2; (d) a = 0.5.
for the various values of R so considered. The curve with “*” represents the asymptotic solution of while other curves with solid lines belong to our numerical solutions.