Induced flow in a T-junction

Application Challenge 3-02 © copyright ERCOFTAC 2004

Overview of CFD Simulations

Simulations of this test case with the code ESTET Version 3.2 (supplied by EDF) have been carried out for a range of the test data. These include a study of the influence of the bulk velocity in the auxiliary pipe as indicated in Table 1, a comparison of two turbulence models and grid sensitivity studies.

The turbulence models considered here are:

- the standard high Reynolds number k-epsilon model by Launder and Spalding (1974)

- a high Reynolds number Reynolds Stress Model (or SMC for “Second Moment Closure”) with isotropic dissipation, turbulent self-transport modelled by the usual gradient diffusion following Daly and Harlow (1970) and the pressure-strain correlation term recommended by Launder (1989), i.e. consisting of the sum of Rotta's return to isotropy and isotropization of production terms.

The grid dependence study has been carried out with meshes of 100000, 400000, and 1500000 nodes for each turbulence model and each auxiliary mass flow rate considered. Aspect ratios (ratio of the longest to the smallest edge in each cell) have been kept approximately between 1 and 2.

Simulation CFD1 is detailed hereafter. It corresponds to EXP1, with k-epsilon model and a coarse mesh. The 17 remaining simulations differ from CFD1 by the turbulence model, the grid and the auxiliary pipe velocity only: Table 3 CFD B provides a summary.

Remark: with the Reynolds Stress Model all conditions are identical to those used for the k-epsilon model, except for the two differences reported hereafter:

- the inlet conditions for the Reynolds stresses are derived from the k-epsilon inlet conditions assuming isotropic turbulence: the diagonal components are taken equal to 2/3 k whereas the extra-diagonal components are set to zero. The inlet conditions for the other variables remain unchanged.

- the dissipation and Reynolds stresses are cell-centered whereas with the k-epsilon model, k and epsilon are located at the vertices.

Simulation Case CFD1

Solution strategy

The flow investigated being unsteady, the convergence had to be determined from time averages. Hence, once the vortical motion established, time average values of the variables were estimated in the dead leg. The vortex penetration was then determined by visualizing secondary velocities. The calculations were stopped once the value for the mean vortex penetration was stabilized.

Computational Domain

The computational domain represents the geometry described in Figure 1 with the same sharp-edged junction but with different pipe lengths to minimize the size of the numerical model. The inlet of the main pipe has been placed relatively near from the junction to reduce the size of the numerical model. Hence, special attention has been devoted to determination of inlet boundary conditions. The outlet is sufficiently far away to prevent any recirculating structure from crossing it. The height of the auxiliary pipe has been chosen large enough to prevent the corkscrew pattern from reaching its top.

The mesh is single block structured in cartesian orthogonal coordinates. Curved boundaries are represented by special features of the code which enables the treatment of slanted boundary cells (prisms, tetrahedrons...).

Figure 3 : Computational domain.

Boundary Conditions

Default initial conditions implemented in the code have been used for each simulation.

As for the boundary conditions at the inlets (main pipe and auxiliary pipe), Dirichlet conditions are used for all transported variables. In the main pipe, the incoming flow is supposed to be fully developed. Hence, the mean velocity and turbulent quantities were obtained from preliminary periodic pipe computations with the same bulk velocity and the same cross-section. For the auxiliary pipe, due to the very small flow-rates, it is not necessary to describe so accurately the incoming flow. The mean velocity was taken constant and turbulent quantities were estimated from usual experimental laws for the friction velocity. In fully developed pipe flows with zero-roughness, and for Reynolds numbers Re (based on the bulk velocity V_A and the hydraulic diameter D) ranging between 5,000 and 30,000, the friction velocity u_t can be determined from $u_{t}=V_{A}((0.3164/Re^{0.25})/8)^{0.5}$ The turbulent kinetic energy k can then be estimated from u_t as k = u_t²/0.3 and the dissipation as ε = u_t³/(0.42D^0.1).

At the outlet (main pipe), zero gradient conditions were applied for all transported variables.

Wall functions were used at solid walls.

The velocity in the main pipe is 9.2m/s, given in Table 1. The auxiliary pipe velocity can be found in Table 3 CFD B.

Application of Physical Models

The standard high Reynolds number k-epsilon model proposed by Launder and Spalding (1974) has been used for CFD1. See Table 3 CFD B for the other simulations.

Numerical Accuracy

The numerical techniques are based on finite difference and finite volume discretizations. The equations are discretized on a 3D semi-staggered grid. Velocity is located at the vertices and pressure is cell-centered. When the k-epsilon model is used, the turbulent variables are located at the vertices (whereas they are cell-centered with the Reynolds stress model).

The algorithm for the solution of transient Navier-Stokes equations relies on a segregated velocity-pressure formulation.The advective terms are treated by a method of characteristics. The trajectory is approximated by a second order Runge Kutta scheme with a third order 3D Hermitian polynomial for interpolation. Diffusion with explicit and implicit source terms for dynamic variables is solved implicitly. For the computation of the velocity components, a third step is required in order to prescribe the mass conservation, leading to a Poisson equation for the pressure increment. In order to avoid non physical oscillations of the pressure field and the associated difficulties in obtaining a converged solution, a variant of the Rhie and Chow interpolation is used as in Méchitoua et al., (1994). The overall precision of the discretization is first order in time and second order in space.

CFD Results

Figure 4 compares the length of vortex penetration obtained by numerical simulation with experimental measurements.

Name	GNDPs	PDPs (problem definition parameters)			SPs (simulated parameters)	*Sensitivity Analysis*
	Inlet Re	Main Pipe Bulk Velocity $V_{M}$ (M/S)	Auxiliary Pipe Bulk Velocity $V_{A}$ (m/s)	Velocity Ratio $V_{A}/V_{M}$	DOAPs	Grid (number of nodes)	Turbulence model
CFD 1	895,000	9.2	0.092	1%	vortex penetration	100,000	k-epsilon
CFD 2	895,000	9.2	0.092	1%	vortex penetration	400,000	k-epsilon
CFD 3	895,000	9.2	0.092	1%	vortex penetration	1,500,000	k-epsilon
CFD 4	895,000	9.2	0.046	0.5%	vortex penetration	100,000	k-epsilon
CFD 5	895,000	9.2	0.046	0.5%	vortex penetration	400,000	k-epsilon
CFD 6	895,000	9.2	0.046	0.5%	vortex penetration	1,500,000	k-epsilon
CFD 7	895,000	9.2	0.023	0.25%	vortex penetration	100,000	k-epsilon
CFD 8	895,000	9.2	0.023	0.25%	vortex penetration	400,000	k-epsilon
CFD 9	895,000	9.2	0.023	0.25%	vortex penetration	1,500,000	k-epsilon
CFD 10	895,000	9.2	0.092	1%	vortex penetration	100,000	SMC
CFD 11	895,000	9.2	0.092	1%	vortex penetration	400,000	SMC
CFD 12	895,000	9.2	0.092	1%	vortex penetration	1,500,000	SMC
CFD 13	895,000	9.2	0.046	0.5%	vortex penetration	100,000	SMC
CFD 14	895,000	9.2	0.046	0.5%	vortex penetration	400,000	SMC
CFD 15	895,000	9.2	0.046	0.5%	vortex penetration	1,500,000	SMC
CFD 16	895,000	9.2	0.023	0.25%	vortex penetration	100,000	SMC
CFD 17	895,000	9.2	0.023	0.25%	vortex penetration	400,000	SMC
CFD 18	895,000	9.2	0.023	0.25%	vortex penetration	1,500,000	SMC