Modeling Asymmetric Flow of Viscoelastic Fluid in Symmetric Planar Sudden Expansion Geometry Based on User-Defined Function in FLUENT CFD Package

Abstract

Through embedding an in-house subroutine into FLUENT code by utilizing the functionalization of user-defined function provided by the software, a new numerical simulation methodology on viscoelastic fluid flows has been established. In order to benchmark this methodology, numerical simulations under different viscoelastic fluid solution concentrations (with solvent viscosity ratio varied from 0.2 to 0.9), extensibility parameters ( $100 \leq L^{2} \leq 500$ ), Reynolds numbers (0.1 ≤ Re ≤ 100), and Weissenberg numbers (0 ≤ Wi ≤ 20) are conducted on unsteady laminar flows through a symmetric planar sudden expansion with expansion ratio of 1: 3 for viscoelastic fluid flows. The constitutive model used to describe the viscoelastic effect of viscoelastic fluid flow is FENE-P (finitely extensive nonlinear elastic-Peterlin) model. The numerical simulation results show that the influences of elasticity, inertia, and concentration on the flow bifurcation characteristics are more significant than those of extensibility. The present simulation results including the critical Reynolds number for which the flow becomes asymmetric, vortex size, bifurcation diagram, velocity distribution, streamline, and pressure loss show good agreements with some published results. That means the newly established method based on FLUENT software platform for simulating peculiar flow behaviors of viscoelastic fluid is credible and suitable for the study of viscoelastic fluid flows.

1. Introduction

Adding polymer or some certain surfactant into water or other Newtonian fluid, the solution has the behaviors of both viscous Newtonian fluid and elastic solid simultaneously and is called viscoelastic fluid. The presence of elasticity brings viscoelastic fluid special rheological characteristics which are different from those of Newtonian fluid, such as shear thinning in viscosity, nonzero normal stress difference, Weissenberg effect (rod-climbing), extrusion swell, tubeless siphon, elastic recoil [1], and turbulent drag reduction effect [2 –7]. Among these characteristics, turbulent drag reduction is an exciting phenomenon and has great potential in energy saving for industrial application systems. Therefore, since Toms [8] reported the drag-reducing effect in viscoelastic fluid turbulent flow for the first time, this has motivated researchers to conduct studies on it by theoretical, experimental, and numerical methods. However, the up-to-date experimental techniques cannot obtain the information of the molecular or micelle conformation tensor and elastic stress field in viscoelastic fluid flow, which is a key parameter for studying the mechanism of turbulent drag reduction. Besides, the development of computer technology prompts researchers to pay more attention to numerical simulation studies on viscoelastic fluid turbulent drag-reducing flows. Direct numerical simulation (DNS), as one of numerical simulation methods solving the governing equations directly, can obtain the most detailed information of the turbulent flow field and microstructure conformation field, and thus it is widely adopted to study the turbulent drag reduction mechanisms of viscoelastic fluid flow. However, the extremely large temporal and spatial resolutions induced by the multiscale characteristics of turbulence bring about extremely high demands for computer's memory and time consumption. According to the exponential relationship between the number of grid and Reynolds number, the computable Reynolds number for DNS is confined to be small or moderate, which cannot satisfy the need of engineering practice. Compared with DNS, the current commercial computational fluid dynamics (CFD) software using RANS (Reynolds averaged Navier-Stokes) method such as FLUENT (used in the present study as the major tool to realize our method for numerical simulation on viscoelastic fluid flow) can provide mature, reliable, and robust numerical simulation algorithms for turbulent flows, which is suitable for engineering applications. However, there is no constitutive model for viscoelastic fluid in those commercial CFD codes for turbulent flow simulation. The functionalization of user-defined function (UDF) provided by the FLUENT package makes it achievable to embed user-defined scalar equations and the source terms of governing equations into the calculations. With considerations of the above-mentioned issues, we have been then motivated to develop a new numerical simulation methodology for viscoelastic fluid flow based on FLUENT CFD package. Our final goal is, mainly for the engineering purpose, to realize numerical simulation on turbulent drag-reducing flows of viscoelastic fluid at Reynolds number as high as in a real application system. We have developed an in-house subroutine modeling viscoelastic fluid flow and embedded this subroutine into the FLUENT code through utilizing the functionalization of UDF in the software package. Before the development of subroutine modeling turbulent flow for the next step, the established viscoelastic fluid subroutine for laminar flow without consideration of turbulence modeling needs to be benchmarked, which is the purpose of the present paper.

Flow through abrupt expansion geometry is a classical problem and involves fundamental characteristics including separation, recirculation, and reattachment. It attracts many numerical and experimental studies because of its practicality in many areas such as the design of water channels, ducts, and the air flow around the buildings. For viscoelastic fluid, sudden expansion flow is a common phenomenon in chemical and food industries, such as extrusion process, mold filling, and food processing. Exploring and understanding the different flow behaviors of viscoelastic fluid from that of Newtonian fluid in the expansion geometry is of great importance in practical application. Therefore, the flows of Newtonian and non-Newtonian fluids through sudden expansion geometry have attracted considerable attentions, especially for laminar flow through the expansion of planar pattern, where the asymmetric flows are known to occur above a critical Reynolds number even though the expansion is symmetric. However, the published numerical researches for sudden expansion flow of viscoelastic fluid are scarce. One of the earlier numerical researches on laminar flow of viscoelastic fluid through symmetric planar sudden expansion geometry was performed by Townsend and Walters [11] with adopting the linearised form of Phan-Thien and Tanner (PTT) viscoelastic constitutive model and using finite element method. Flows through two-dimensional expansion with expansion ratio of 3: 40 at low Reynolds number were simulated for comparison with some cases in their experiments, and satisfactory agreement was obtained. Baloch et al. [12] used the same method and viscoelastic constitutive model as used in Townsend and Walters [11] to simulate the two-dimensional expansion flows with expansion ratios of 3: 40 and 1: 80 at low Reynolds numbers (Re ≤ 4). The results for the flows with expansion ratio of 3: 40 showed good qualitative agreement with some experimental cases conducted by Townsend and Walters [11]. Darwish et al. [13] and Missirlis et al. [14] applied upper-convected Maxwell (UCM) model and developed finite volume method to simulate the creeping flows of viscoelastic fluid in a 1: 4 sudden expansion at Re = 0.1. Poole et al. [15, 16] numerically investigated the creeping flow of viscoelastic fluid with three different models (UCM, Oldroyd-B, and the linear form of PTT model) through a 1: 3 planar sudden expansion and studied the effect of expansion ratio on the creeping flow of viscoelastic fluid obeying UCM model. The results of numerical simulations with different constitutive models showed that the reduction of both the length and intensity of the vortex, which was caused by the effect of elasticity, is much lower than that reported in previous studies, and a significant region of recirculation still exists for all the models at high Deborah number [15]. The study of expansion ratio effect on the flow characteristics showed that both the length and intensity of the recirculation zone are increased with increasing expansion ratio at the same Deborah number, and for large expansion ratios the recirculation size and strength are decreased with increase of Deborah number, while for small expansion ratios the recirculation length is initially decreased at low Deborah numbers, followed by an increase as Deborah number increases [16]. Since the Reynolds number in above numerical researches is lower than a certain critical Reynolds number, the flow in sudden expansion geometry is symmetry. Oliveira [10] employed modified finitely extensive nonlinear elastic-Chilcott-Rallison (FENE-CR) model to simulate the behaviors of viscoelastic fluid flows in a symmetric planar expansion geometry with expansion ratio of 1: 3 under laminar state and is likely to be the first one who reported the asymmetric flows at high Reynolds numbers. Some detailed studies were conducted on the effects of Weissenberg numbers, Reynolds numbers, concentration parameters, and extensibility parameters on the flow characteristics. The results showed that extensibility had a rather small effect upon the flow pattern, while concentration had a very strong effect that the asymmetry was completely removed with the increase of solvent viscosity ratio. The effect of elasticity mainly occurred for small Weissenberg numbers with tendency to reduce the vortex asymmetry. The increase of inertia prompted the formation of asymmetry with larger onset Reynolds number than that of Newtonian fluid. Afterwards, Rocha et al. [17] employed the modified FENE-CR model to study the bifurcation phenomena in viscoelastic fluid flow through a symmetric planar sudden expansion with expansion ratio of 1: 4 by numerical methods. Similar results to those in Oliveira's research were obtained for the effects of Weissenberg number, Reynolds number, concentration parameter, and extensibility parameter on the flow characteristics. Besides, the smaller the onset Re was, larger and more intense the vortices acquired at larger expansion ratio by comparison with Oliveira's result were, and a new recirculation region was formed at the downstream of small vortex with Re higher than 64 and 73.5 for the Newtonian and viscoelastic fluids, respectively.

In this paper, considering the comprehensive discussions and analyses on the flow behaviors in the numerical study performed by Oliveira [10], the unsteady viscoelastic fluid and Newtonian fluid laminar flows through a symmetric planar sudden expansion with expansion ratio of 1: 3 with main focus on the flow behavior of asymmetric flow are selected as the test case to benchmark the subroutine constituting viscoelastic fluid embedded into FLUENT software. The constitutive model used to simulate viscoelastic fluid flow is finitely extensive nonlinear elastic-Peterlin (FENE-P) model. By comparing our simulation results with that by Oliveira [10], the established new simulation method for viscoelastic fluid flows using FLUENT software is verified. Besides, some special phenomena induced by elasticity in viscoelastic fluid are also studied. In Section 2, the description of the flow problem and numerical method, including the algorithm and the way to embed the constitutive model of viscoelastic fluid into FLUENT software, is presented. The results, for Newtonian fluid flow at different Reynolds number and for viscoelastic fluid flow with one changing parameter and three fixed parameters among the Reynolds number, Weissenberg number, concentration parameter, and extensibility parameter, are given in Section 3. The analyses of flow characteristics and the effect of those four parameters on the flow behavior, comparison between the results of Newtonian fluid case, and that of viscoelastic fluid case and comparison with Oliveira's results are also elucidated. Finally, the main conclusions are given in Section 4.

2. Numerical Simulation Procedures

2.1. Problem Description

The asymmetric flow phenomenon in symmetric planar sudden expansion geometry is schematically shown in Figure 1. For both Newtonian fluid and viscoelastic fluid, there exists a critical Reynolds number at which the symmetric flow in the symmetric expansion geometry becomes asymmetric with a larger and a smaller recirculation zone behind the step change. This bifurcation effect is possibly caused by Coanda effect as explained by Oliveira in which fluid tends to flow along the protruding surface, and hence any perturbation will push the main flow to one side where larger velocity and lower pressure are brought about, and ultimately the asymmetry is naturally formed [10]. In order to compare with Oliveira's earlier work, the symmetric abrupt expansion channel with expansion ratio E = D/d = 3 is selected as the computational model. Herein, d and D are the height of the smaller channel and the larger channel, respectively, as shown in Figure 1. According to the recommendation from the study of Hawa and Rusak [18] and Oliveira [10], the lengths of the smaller channel and the larger channel are set to be L₁ = 2d and L₂ = 50d, respectively. With L₂ = 50d, the flow at the outlet will not affect the flow in the expansion geometry, and the accuracy and facticity of flow are guaranteed. Considering that too small velocity at the inlet is not conducive to the calculation for the case at small Reynolds number, small value of the height of the smaller channel d is adopted (d = 1 mm is used herein). Besides, the set of coordinate system is also shown in Figure 1. Numerical simulations are conducted on two-dimensional isothermal flows of incompressible fluid for Newtonian and viscoelastic fluid cases, respectively. The flow characteristics under different Reynolds numbers, Weissenberg numbers, concentration parameters, and extensibility parameters are obtained and analyzed, including the recirculation zone, bifurcation effect, streamline, velocity distribution, and flow resistance. The differences between Newtonian fluid and viscoelastic fluid and those between the present and previously reported results are compared.

Figure 1:

The schematic diagram of asymmetric flow in symmetric sudden expansion geometry.

2.2. Numerical Method

Nonuniform collocated mesh, which is the same as used by Oliveira [10], is adopted in the computational domain, as shown in Figure 2. The meshes behind the step change are arranged to be denser. The grid size along the x-direction is gradually increased with a constant factor for x ≥ 0, and the mesh along the y-direction is uniform and of the same size with the minimum cell size along the x-direction.

Figure 2:

Partial view of mesh (−2d ≤ x < 4d, −1.5d ≤ y ≤ 1.5d).

For a two-dimensional unsteady isothermal flow of incompressible fluid, the continuity and momentum equations are as follows:

\frac{\partial u_{i}}{\partial x_{i}} = 0,

(1)

ρ \frac{\partial u_{i}}{\partial t} + ρ u_{j} \frac{\partial u_{i}}{\partial x_{j}} = - \frac{\partial p}{\partial x_{i}} + \frac{\partial τ_{i j}^{s}}{\partial x_{j}} + \frac{\partial τ_{i j}^{p}}{\partial x_{j}},

(2)

where u_i is the velocity, ρ is the fluid density and constantly equals 1000 kg/m³, p is the pressure, τ_ij^s (i = 1, 2; j = 1, 2) is the stress caused by the solvent, and τ_ij^p is the stress induced by the elasticity. For Newtonian fluid; the third term on the right side of (2) disappears. Consider τ_ij^s = μ_s (∂ u_i/ ∂ x_j + ∂ u_j/ ∂ x_i), where μ_s is the solvent viscosity and constantly equals 0.001 Pa · s. For viscoelastic fluid, FENE-P model is adopted as the constitutive equation for elastic stress:

τ_{i j}^{p} = n (x, t) \frac{μ_{p}}{λ} [f (r) C_{i j} - δ_{i j}],

(3)

where n(x, t) is the unit concentration and equals 1 when the concentration of polymer or surfactant is homogeneous (assumed in this paper), μ_p is the viscosity of viscoelastic polymer or surfactant solution, λ is the relaxation time of viscoelastic fluid, δ_ij is the Kroneker symbol, f(r) is the Peterlin function and for FENE-P model it is described as

f (r) = \frac{L^{2} - 3}{L^{2} - r^{2}},

(4)

where r is the molecular or micellar length under equilibrium state and equals the square root of the trace of conformation tensor and L is the molecular or micellar extensibility in relation to its equilibrium size. In (3) C_ij is the conformation tensor of polymer molecule or surfactant micelle; for i ≠ j, there exists C_ij = C_ji, and its transport equation is as follows:

\frac{\partial C_{i j}}{\partial t} + u_{k} \frac{\partial C_{i j}}{\partial x_{k}} = C_{i k} \frac{\partial u_{j}}{\partial x_{k}} + C_{k j} \frac{\partial u_{i}}{\partial x_{k}} - \frac{1}{λ} [f (r) δ_{i k} C_{k j} - δ_{i j}] + κ \frac{\partial}{\partial x_{k}} (\frac{\partial C_{i j}}{\partial x_{k}}) .

(5)

The last term in (5) is the artificial viscosity term, which is added to enhance the stability and convergence of calculation. κ is a nonzero value in the calculation, while, in the present study, converged solution can be obtained with κ = 0 by adjusting the underrelaxation factors, and so the artificial viscosity term is actually not used.

The Reynolds number (Re) is defined as follows:

Re = \frac{ρ U d}{μ_{s} + μ_{p}},

(6)

where U is the average velocity at the inlet. The decisive similarity criterion for viscoelastic flow is the Weissenberg number (Wi) defined as follows:

Wi = λ \frac{U}{d} .

(7)

As mentioned above, the viscoelastic constitutive equation is embedded into the FLUENT code through UDF functionality. In UDF, the general user-defined scalar transport equation is shown as follows, in which there are four terms to be customized: transient term, convection term, diffusion term, and source term, respectively:

\frac{\partial ρ ϕ}{\partial t} + \frac{\partial}{\partial x_{i}} (F_{i} ϕ - Γ \frac{\partial ϕ}{\partial x_{i}}) = S_{ϕ},

(8)

where ϕ is the scalar, corresponding to components of the conformation tensor C_ij in (5) herein, and Γ is diffusion coefficient, corresponding to the coefficient of artificial viscosity term κ, and equals 0 in the present calculation as mentioned above. In (5), the transient term corresponds to the first term on the left side; the convection term is the second term on the left side; the remaining terms in (5) belong to the source term S_ϕ. The customizations of these terms are realized by writing corresponding functions with C language, connecting the functions with interfaces through the corresponding settings in the code. By customizing these terms, the viscoelastic constitutive model is then embedded into the FLUENT code, and the viscoelastic fluid flow problem can be solved.

The boundary condition of entrance is set to be velocity inlet with a fully developed velocity profile (as shown in Figure 1) which is defined by the following expression:

u = 1.5 U [1 - {(\frac{y}{d / 2})}^{2}] .

(9)

The exit is set as outflow condition, and all the walls have no slip boundary conditions. The algorithm provided by FLUENT software is based on finite volume method. Double precision is employed for our calculations. Pressure implicit with splitting of operators (PISO) scheme, which is more suitable for the calculation of transient flow, is adopted to solve the coupling of velocity and pressure; standard scheme is utilized to discretize the pressure equation; QUICK scheme is used to discretize the momentum equation and user-defined scalar equations. The convergence criteria (residuals below which the calculation is considered to be converged) and underrelaxation factors are shown in Tables 1 and 10, respectively. Note that different underrelaxation factors are adopted for different conditions. The time step is set to be 0.001 s for transient calculation. The velocity of whole flow field is initialized by the average velocity at the inlet U, and the conformation tensor is initialized to be C_xx = 1, C_yy = 1, and C_xy = 0.

Table 1:

Convergence criteria.

Variable	Continuity	u	v	C _xx	C _yy	C _xy

Residual	10⁻⁸	10⁻⁶	10⁻⁶	10⁻⁶	10⁻⁶	10⁻⁶

3. Results and Discussions

Numerical simulations are performed on the flow through sudden expansion geometry for Newtonian fluid at different Reynolds numbers, while, for viscoelastic fluid, three more parameters, Weissenberg number, concentration (measured by the solvent viscosity ratio defined as $β = μ_{s} / (μ_{s} + μ_{p})$ ), and extensibility parameters (L²), are taken into consideration. The numerical simulations are conducted at one changing parameter and the other three unchanged parameters to study the effects induced by these four parameters individually. Note that the existence of the dynamic viscosity of viscoelastic polymer or surfactant solution μ_p in the definition of Reynolds number makes the average velocity at the inlet U different for Newtonian fluid and viscoelastic fluid cases at the same Reynolds number.

3.1. Mesh Independence Verification

In order to guarantee the numerical accuracy, the mesh independence is studied at first. Three meshes are selected, including the coarse mesh with 1550 grids and minimum cell size of 0.1d, the medium mesh with 6200 grids and minimum cell size of 0.05d, and the fine mesh with 24800 grids and minimum cell size of 0.025d, which are the same as those in Oliveira's work. Numerical simulations are conducted at Re = 60 for Newtonian fluid flow and at Re = 60, Wi = 2, β = 0.9, and L² = 100 for viscoelastic fluid flow with three different meshes. To demonstrate the mesh independence of solutions, Richardson's extrapolation technique [10] is applied to analyze the recirculation size (X_r shown in Figure 1), as listed in Table 2.

Table 2:

Verification of mesh independence.

	Newtonian case		Viscoelasticcase
	X_r1/d	X_r2/d	X_r1/d	X_r2/d

Coarse mesh	3.9314	7.5749	4.1452	7.4470
Medium mesh	3.9498	7.5820	4.0891	7.4800
Fine mesh	3.9515	7.6044	4.0588	7.5136
Richardson's extrapolation	3.9521	7.6119	4.0487	7.5248
Discretization error (%)	0.06	0.39	1.00	0.60

The basic principle of Richardson's extrapolation is that any result which is directly evaluated from the numerical solution or any integrated result, φ, is assumed to converge with mesh refinement as φ = φ_h + Ch^p, where h is the mesh size, C is a constant, p is the order of the scheme, and φ_h is the value of φ obtained based on the mesh with size of h. If the cell size of fine mesh is denoted as h, then the cell size of medium mesh would be 2h, and there exists φ = φ_2h + C(2h)^p. Therefore, the constant C equals $Δ φ / [h^{p} (2^{p} - 1)]$ (Δφ = φ_h – φ_2h), and the corresponding discretization errors will be $ε_{h} = φ - φ_{h} = Δ φ / (2^{p} - 1)$ on fine mesh and $ε_{2 h} = φ - φ_{2 h} = 2^{p} Δ φ / (2^{p} - 1)$ on medium mesh. By assumption of second-order accuracy (p = 2), the extrapolated values are calculated from $φ = (4 φ_{h} - φ_{2 h}) / 3$ .

Note that the discretization error in Table 2 is based on the medium mesh. It can be obtained that the discretization errors are all within 1%. Considering the economy of the calculation and convenience of comparison with Oliveira's results at the same time, the medium mesh with 6200 grids is used in the present numerical simulations.

3.2. Results for Newtonian Fluid Cases

The bifurcation effect (recirculation size asymmetry DX = (X_r2 – X_r1) versus Re) obtained in this paper, together with that of Fearn et al. [9] and Oliveira [10], is shown in Figure 3. Recirculation sizes for different Reynolds numbers are given in Table 3, and the comparison with Oliveira's results is shown in Figure 4. From Table 3, we can obtain that the critical Reynolds number at which the flow becomes asymmetric is 55. It is in good agreement with those of previous numerical researches with different numerical methods, which are shown in Table 4. The bifurcation effect (Figure 3) and recirculation size (Figure 4) also show excellent agreement with previous researches.

Table 3:

Recirculation size for Newtonian fluid cases.

Re	X_r1/d	X_r2/d

0.1	0.502	0.502
1	0.548	0.548
10	1.195	1.195
20	2.096	2.096
30	3.065	3.065
40	4.057	4.057
50	5.062	5.062
52	5.265	5.265
54	5.467	5.467
55	4.841	6.197
56	4.480	6.628
58	4.132	7.175
60	3.950	7.582
64	3.770	8.230
70	3.670	9.011
75	3.650	9.572
80	3.656	10.074
85	3.677	10.531
90	3.706	10.950
95	3.741	11.334
100	3.779	11.686

Table 4:

Critical Reynolds numbers obtained by previous researches.

Sources	Methods	Re_c

Oliveira [10]	Finite volume method	54
Hawa and Rusak [18]	Vorticity/stream function finite difference	53.8
Fearn et al. [9]	Finite element method	53.9
Drikakis [19]	Fourth-order finite difference method	53.3
Mishra and Jayaraman [20]	Finite element method	54

Figure 3:

Comparisons of the bifurcation effect for Newtonian fluid case with the experimental data from Fearn et al. [9] and the numerical results from Oliveira [10].

Figure 4:

Comparison of the recirculation size for Newtonian fluid case with results in Oliveira [10].

Besides, the velocity distribution in the flow field is compared with previous work at Re = 80 in Figure 5. The profiles of dimensionless velocity in the x-direction at four different positions along the streamwise direction are plotted. It can be obtained that the results in this paper are closely consistent with experimental data of Fearn et al. [9] and numerical results of Oliveira [10].

Figure 5:

Comparison of velocity profiles with experimental data of Fearn et al. [9] and numerical results of Oliveira [10].

For further verification, we consider a case with Re = 187 at which a new recirculation zone (measured by X_r3 at detachment point and X_r4 at reattachment point) appears adjacent to the smaller recirculation zone behind the step change, and the streamlines obtained in this paper are plotted in Figure 6. The result, sizes of three recirculation zones, is compared with previous results and shows excellent agreement with them, as shown in Table 5.

Table 5:

Comparisons of the sizes of recirculation zones at Re = 187.

Sources	X_r1/d	X_r2/d	X_r3/d	X_r4/d

Present results	4.48	15.24	12.50	25.46
De Zilwa et al. [21]	5.0	15.2	12.6	26
Oliveira [10]	4.5	15.1	12.4	25.3

Figure 6:

Streamline for Newtonian fluid case at Re = 187.

3.3. Results for Viscoelastic Fluid Cases

In this section, in order to study the effect induced by elasticity, extensibility, concentration, and inertia on the sudden expansion flow, numerical simulations are systematically conducted on viscoelastic fluid flows through sudden expansion geometry at different Weissenberg numbers, extensibility parameters, solvent viscosity ratios, and Reynolds numbers, respectively.

3.3.1. Effect of Elasticity: Weissenberg Number Dependence

In order to study the effect of elasticity on the flow, numerical simulations are conducted at changing Weissenberg numbers and constant Reynolds number, concentration, and extensibility. In order for comparison, the base values of the three constant parameters are selected to be the same with those in Oliveira [10]; that is, Re = 60, β = 0.9, and L² = 100. As analyzed by Oliveira, the extensibility parameter of L² = 100 is typical for FENE models, and solvent viscosity ratio of β = 0.9 corresponds to the concentration parameter of c ≈ 0.1 (concentration parameter c is defined as $c = μ_{p} / μ_{s}$ , and hence $β = 1 / (1 + c)$ ), which is typical for dilute solutions. Reynolds number is set to be 60, which is not far away from the critical Reynolds number 55 for Newtonian fluid flow, and is chosen for the convenience of studying the effect on the flow characteristics, especially on the flow asymmetry. However, because the convergence of calculation cannot be guaranteed by only adjusting the underrelaxation factors for larger Weissenberg numbers (Wi > 20), the Weissenberg numbers in the present simulations only vary from 0 to 20, while in Oliveira's work Weissenberg number is up to 100. Study on the viscoelastic fluid flow through the sudden expansion geometry for larger Weissenberg numbers will be performed and reported in the next step of our continuous studies.

The effects of elasticity on the recirculation size and flow asymmetry are given in Table 6. And the comparison with Oliveira's results is shown in Figure 7. One can obtain that, for viscoelastic fluid, the size of the larger vortex slightly increases with the increase of Wi, while the size of the smaller vortex decreases. The increase of the elasticity does not remove the asymmetry of the flow, which is consistent with Oliveira's result qualitatively. But the variation trend of the sizes of two recirculation zones is different from that of Oliveira's result, especially for the size of the larger vortex, as shown in Figure 7. For this Wi-dependent phenomenon appeared in the asymmetry of planar sudden expansion flow, elaborated experiments are called for in order to quantify the sizes of two corner vortices as a function of viscoelasticity.

Table 6:

Sizes and asymmetry of recirculation zone for different Weissenberg numbers.

Wi	X_r1/d	X_r2/d	DX/d

0	3.950	7.582	3.632
1	4.094	7.397	3.303
2	4.089	7.480	3.391
3	4.083	7.590	3.507
4	4.087	7.679	3.592
5	4.099	7.747	3.648
10	4.149	7.951	3.802
15	4.108	8.105	3.997
20	4.041	8.182	4.141

Figure 7:

Comparison of recirculation sizes at different Weissenberg numbers with results of Oliveira [10].

3.3.2. Effect of Extensibility: L² Dependence

Numerical simulations on the effect of extensibility are conducted at Re = 60, Wi = 1, and β = 0.9. The extensibility parameter L² varies from 100 to 500 with increment step of 100. The sizes and asymmetry of vortices are given in Table 7 and are compared with those of Oliveira's results in Figure 8. It shows that the extensibility has little influence on the variation of vortex length. With the increase of extensibility, the flow asymmetry decreases, but, in Oliveira's research, the asymmetry size (DX/d) is almost the same with just a slight decrease except for the point with L² = 500. And the variation trend of the size of larger recirculation zone is opposite to that in Oliveira's work though there is only small quantitative difference between them. The sizes of smaller recirculation zone are in both qualitative and quantitative agreement with Oliveira's results, as shown in Figure 8.

Table 7:

Sizes and asymmetry of recirculation zone for different extensibility parameters.

L ²	X_r1/d	X_r2/d	DX/d

100	4.094	7.397	3.303
200	4.153	7.351	3.198
300	4.189	7.329	3.140
400	4.215	7.315	3.100
500	4.236	7.303	3.067

Figure 8:

Comparison of recirculation sizes at different extensibility parameters with results of Oliveira [10].

3.3.3. Effect of Concentration: β Dependence

Numerical simulations are conducted at fixed Re = 60, Wi = 2, and L² = 100 and changing solvent viscosity ratios from 0.9 to 0.2 (the corresponding concentration of polymer or surfactant solution increases) to study the effect of concentration. Note that we do not reduce the time step with the decrease of β as Oliveira did. For each concentration, the time steps in the numerical simulations are all set to be 0.001 s. The underrelaxation factors are adjusted to be smaller for smaller β (higher concentration parameter c) since the coupling between the momentum and constitutive equations is enhanced with the increase of c. The sizes and asymmetry of recirculation zone and comparisons with Oliveira's resulet are shown in Table 8 and Figure 9, respectively. We can obtain that concentration exhibits much stronger effect on the flow by comparing it with that of extensibility. The elongation of the smaller vortex and the reduction of the larger vortex make the flow asymmetry shrink gradually with the reduction of β from Newtonian case (β = 1), which is induced by the enhancement of elasticity and in accordance with Oliveira's result. But the asymmetry is not removed and we do not capture the change from asymmetry to symmetry as Oliveira did in the presently investigated range of β. However, the variation trend of the asymmetry of recirculation zone can still illustrate that the increase of concentration facilitates the evolvement from asymmetric state to symmetric state.

Table 8:

Sizes and asymmetry of recirculation zone for different concentrations.

β	X_r1/d	X_r2/d	DX/d

0.9	4.089	7.480	3.391
0.8	4.221	7.358	3.137
0.7	4.363	7.218	2.855
0.6	4.511	7.067	2.556
0.5	4.669	6.907	2.238
0.4	4.834	6.744	1.910
0.3	4.976	6.601	1.625
0.2	5.146	6.409	1.263

Figure 9:

Comparison of recirculation sizes at different concentrations with results of Oliveira [10].

3.3.4. Effect of Inertia: Reynolds Number Dependence

In this section, we perform the numerical simulations at fixed Wi = 2, L² = 100, and β = 0.5 and changing Reynolds number from 1 to 100 to study the effect of inertia on the flow characteristics of viscoelastic fluid. The results, including the recirculation sizes at different Re and the bifurcation behaviors, are shown in Table 9 and Figures 10 and 11, respectively, together with Oliveira's numerical results for comparison. Note that all of the flows simulated in this paper are finally developed into steady state, except for the case with Re = 1, Wi = 2, L² = 100, and β = 0.5 in this section. In this case, the flow will change from steady state to unsteady state again and ultimately keep the periodical variation, which is possibly caused by the elastic instability and will be studied in our further research. The detailed results given in this paper are only those for steady states. One can obtain that the onset Reynolds number at which the flow becomes asymmetric is delayed to Re_c = 58 compared with that of Newtonian fluid (Re_c = 55), which is also obtained in Oliveira's research. But the critical Reynolds number for viscoelastic fluid in this paper is smaller than that in Oliveira's work (Re_c = 64). The delay of onset Re for viscoelastic fluid is induced by elasticity which acts as a stabilizing factor for the occurrence of bifurcation. The recirculation sizes and bifurcation effect are in qualitative agreement with Oliveira's results: the sizes of two vortices and asymmetry extent for viscoelastic fluid are smaller than those for Newtonian fluid at the same Re for the whole range of Re (except for the size of smaller vortex at higher Reynolds numbers above its critical onset value), as shown in Figures 10 and 11. The difference between the vortex sizes of viscoelastic fluid and those of Newtonian fluid can be considered as a swelling-like phenomenon similar to extrusion swell as explained by Oliveira [10].

Table 9:

Recirculation sizes for different Reynolds numbers.

Re	X_r1/d	X_r2/d

1	0.544	0.544
10	0.850	0.851
20	1.902	1.902
30	2.931	2.931
40	3.936	3.936
50	4.933	4.933
55	5.434	5.434
56	5.535	5.535
57	5.635	5.635
58	5.322	6.110
59	4.910	6.589
60	4.669	6.907
65	4.192	7.912
70	4.035	8.638
75	3.974	9.261
80	3.956	9.814
85	3.961	10.327
90	3.979	10.799
95	4.006	11.238
100	4.039	11.643

Table 10:

Underrelaxation factors.

Variable	Pressure	Density	Body forces	Momentum	C _xx	C _yy	C _xy

Value	0.005~0.3	0.02~1	0.02~1	0.01~0.7	0.02~0.9	0.02~0.9	0.02~0.9

Figure 10:

Comparisons of recirculation sizes at different Reynolds numbers for viscoelastic fluid with those of Newtonian fluid in this paper and results of Oliveira [10] for viscoelastic fluid.

Figure 11:

Comparisons of bifurcation effect for viscoelastic fluid with that of Newtonian fluid in this paper and results of Oliveira [10] for viscoelastic fluid.

Figure 12 shows the streamlines at different Re for viscoelastic fluid and clearly depicts the changes of flow state, recirculation sizes, and bifurcation with the increase of Re. For small Reynolds numbers of Re = 1 and Re = 10, we capture the streamline's convergence towards the centerline before the expansion, followed by a divergence towards the walls of larger channel after the expansion, which was also captured by Oliveira with viscoelastic case at Re = 0.1 [10], although the phenomenon is not so obvious in the figure. Furthermore, the profiles of velocity in the x-direction along the centerline at those two Reynolds numbers are plotted in Figure 13. And an overshoot of the profile with an undershoot at the downstream is observed at the expansion for viscoelastic fluid flow, while there is no such kind of phenomena for Newtonian fluid flow. The cause of this difference is that inertia is not necessary for diverging flow behavior to occur, which is proposed and explained by Alves and Poole [22] for the contraction flow, and the sudden expansion flow of viscoelastic fluid is somewhat similar to the extrusion swell, as mentioned above, with enhancing swell effect as Re decreases. Besides, the overshoot reduces with the increase of Re, as shown in Figure 13, and finally disappears.

Figure 12:

Streamlines at different Reynolds numbers for viscoelastic fluid.

Figure 13:

Profile of velocity in the x-direction along the central plane (y = 0) for Newtonian fluid and viscoelastic fluid at small Reynolds numbers.

We consider a localized loss coefficient C_I which is independent of the channel length to measure the energy losses, and its definition is described in detail by Oliveira [10]. Note that for the part of fully developed pressure drop coefficient (C_F) in C_I, the fully developed regions $Δ X (Δ X = Δ x / d$ ) in the definition of Fanning friction factor $f (f = Δ P / Δ X$ and $Δ P = Δ p / 0.5 ρ U^{2}$ ), whose choice is important for the accurate calculation of C_I, are selected to be – 2 < x/d < – 1 for smaller upstream channel and 30 < x/d < 50 for larger downstream channel, respectively. The length of – 2 < x/d < – 1 in the upstream channel is not enough for fully developed flow; thus C_I is an estimation value. The relationship between C_I and Re is plotted in Figure 14 for viscoelastic fluid and Newtonian fluid flows and compared with the numerical data of Oliveira [10]. It is seen that C_I decreases with the increase of Re for viscoelastic fluid flow, while, for Newtonian fluid case, C_I decreases at first and then increases slightly and attains the minimum at the critical Reynolds number. C_I of viscoelastic fluid flow is larger than that of Newtonian fluid flow in the whole range of Re. These results are all in qualitative agreement with those obtained by Oliveira. For Newtonian fluid, the results are also quantitatively consistent with Oliveira's. That means the choice of fully developed region within the upstream channel in this paper does not affect the accuracy of C_I so much, for the length of the smaller upstream channel in Oliveira's numerical simulations on Newtonian cases is L₁ = 20d, which can guarantee the accurate evaluation of C_I.

Figure 14:

Change of localized loss coefficient with Reynolds number for Newtonian fluid and viscoelastic fluid and comparison with results of Oliveira [10].

Apart from localized loss coefficient, the flow resistance coefficient fRe, another measure of pressure loss, is also studied in this paper. The definition of Fanning friction factor is given by Oliveira [10]. Figure 15 depicts the changes of fRe with the increase of Re for Newtonian fluid and viscoelastic fluid flows and results of Oliveira. From this figure, we can see that our results show qualitative agreement with those of Oliveira. The flow resistance coefficient decreases with increasing Re. fRe for viscoelastic fluid flow is greater than that for Newtonian fluid flow in the whole range of Re. The change of slope is also captured at critical Reynolds number for both Newtonian fluid and viscoelastic fluid flows. The quantitative difference from Oliveira's results for Newtonian cases is possibly induced by the difference of the length of the upstream channel.

Figure 15:

Variation of flow resistance coefficient with Reynolds number for Newtonian fluid and viscoelastic fluid and comparison with results of Oliveira [10].

3.3.5. Summary of the Effects of Four Parameters and Discussions

According to the results of numerical simulations on viscoelastic fluid flows, the concentration and inertia effects on the flow characteristics are greatly dependent on the results of recirculation sizes and bifurcation effect. While the effect of extensibility is small in the range of extensibility parameter discussed in this paper.

Through comparisons, the numerical simulation results for Newtonian fluid are in both qualitative and quantitative agreements with previous experimental data and numerical results. While for viscoelastic fluid, most of the results are qualitatively consistent with Oliveira's numerical results. The cause of quantitative discrepancy and other qualitative difference is possibly due to the difference of viscoelastic fluid constitutive model except for the reason mentioned above. The constitutive model of FENE-CR adopted by Oliveira in the numerical simulations on viscoelastic fluid flows is an empirical extension of FENE-P model (adopted in this paper), which has a constant shear viscosity for a simple shear flow.

4. Conclusion

Numerical simulations are conducted on the flows through a symmetric planar sudden expansion geometry with an expansion ratio of 1: 3 at different Reynolds numbers for Newtonian fluid and at different Weissenberg numbers, extensibility parameters, concentrations, and Reynolds numbers for viscoelastic fluid modeled by FENE-P model. The results are compared with previous published ones, especially with numerical data of Oliveira [10]. The main conclusions drawn from this study are given below.

The critical Reynolds number at which the flow becomes asymmetric for Newtonian fluid is Re_c = 55, in agreement with the previous experimental and numerical results. All of the results for Newtonian fluid are in excellent agreement with previous published results except for slight quantitative difference of flow resistance coefficient fRe from Oliveira's numerical data, which is possibly caused by the difference of the upstream channel length.

The critical Reynolds number at which bifurcation occurs for viscoelastic fluid with fixed Wi = 2, L² = 100, and β = 0.5 is delayed to Re_c = 58, which is smaller than that in Oliveira's research. Some of the results for viscoelastic fluid in this paper are in qualitative agreement with numerical data of Oliveira with some quantitative discrepancy. The differences are possibly induced by the different constitutive models for viscoelastic fluid.

The comparisons with the previous experimental and numerical results with present numerical results prove the correctness and feasibility of the established numerical simulation methodology for viscoelastic fluid flow by embedding user-defined function of viscoelastic fluid constitutive model into FLUENT software in the present paper.

Footnotes

Nomenclature

Acknowledgments

This work is supported by Foundation for Innovative Research Groups of the National Natural Science Foundation of China (51121004), National Natural Science Foundation of China (51076036, 51276046), and Specialized Research Fund for the Doctoral Program of Higher Education of China (20112302110020). The authors would like to appreciate the valuable discussions with the members of Complex Flow and Heat Transfer Lab and thank Dr. J. L. Yin of Shanghai Jiaotong University for his help and guidance on writing user-defined function.

References

Bird

R. B.

Armstrong

R. C.

Hassager

, and Curtiss

C. F.

, Dynamics of Polymer Liquids—Vol. 1: Fluid Mechanics, John Wiley & Sons, 1987.

F.-C.

Kawaguchi

Wei

J.-J.

, and Hishida

, “Experimental study of drag-reduction mechanism for a dilute surfactant solution flow,” International Journal of Heat and Mass Transfer, vol. 51, no. 3–4, pp. 835–843, 2008.

Cai

W.-H.

F.-C.

Zhang

H.-N.

, “Study on the characteristics of turbulent drag-reducing channel flow by particle image velocimetry combining with proper orthogonal decomposition analysis,” Physics of Fluids, vol. 21, no. 11, Article ID 115103, 2009.

Cai

W.-H.

F.-C.

, and Zhang

H.-N.

, “DNS study of decaying homogeneous isotropic turbulence with polymer additives,” Journal of Fluid Mechanics, vol. 665, pp. 334–356, 2010.

Wei

J. J.

Kawaguchi

F. C.

, “Drag-reducing and heat transfer characteristics of a novel zwitterionic surfactant solution,” International Journal of Heat and Mass Transfer, vol. 52, no. 15–16, pp. 3547–3554, 2009.

Tamano

Itoh

Inoue

Kato

, and Yokota

, “Turbulence statistics and structures of drag-reducing turbulent boundary layer in homogeneous aqueous surfactant solutions,” Physics of Fluids, vol. 21, no. 4, Article ID 045101, 2009.

Tamano

Itoh

Kato

, and Yokota

, “Turbulent drag reduction in nonionic surfactant solutions,” Physics of Fluids, vol. 22, no. 5, Article ID 055102, 2010.

Toms

B. A.

, “Some observations on the flow of linear polymer solutions through straight tubes at large Reynolds numbers,” in Proceedings of 1st International Congress on Rheology, Mysels

M. J.

, Ed., vol. 2, pp. 135–138, Springer, Amsterdam, The Netherlands, 1948.

Fearn

R. M.

Mullin

, and Cliffe

K. A.

, “Nonlinear flow phenomena in a symmetric sudden expansion,” Journal of Fluid Mechanics, vol. 211, pp. 595–608, 1990.

10.

Oliveira

P. J.

, “Asymmetric flows of viscoelastic fluids in symmetric planar expansion geometries,” Journal of Non-Newtonian Fluid Mechanics, vol. 114, no. 1, pp. 33–63, 2003.

11.

Townsend

and Walters

, “Expansion flows on non-newtonian liquids,” Chemical Engineering Science, vol. 49, no. 5, pp. 748–763, 1994.

12.

Baloch

Townsend

, and Webster

M. F.

, “On vortex development in viscoelastic expansion and contraction flows,” Journal of Non-Newtonian Fluid Mechanics, vol. 65, no. 2–3, pp. 133–149, 1996.

13.

Darwish

M. S.

Whiteman

J. R.

, and Bevis

M. J.

, “Numerical modelling of viscoelastic liquids using a finite-volume method,” Journal of Non-Newtonian Fluid Mechanics, vol. 45, no. 3, pp. 311–337, 1992.

14.

Missirlis

K. A.

Assimacopoulos

, and Mitsoulis

, “A finite volume approach in the simulation of viscoelastic expansion flows,” Journal of Non-Newtonian Fluid Mechanics, vol. 78, no. 2–3, pp. 91–118, 1998.

15.

Poole

R. J.

Alves

M. A.

Oliveira

P. J.

, and Pinho

F. T.

, “Plane sudden expansion flows of viscoelastic liquids,” Journal of Non-Newtonian Fluid Mechanics, vol. 146, no. 1–3, pp. 79–91, 2007.

16.

Poole

R. J.

Pinho

F. T.

Alves

M. A.

, and Oliveira

P. J.

, “The effect of expansion ratio for creeping expansion flows of UCM fluids,” Journal of Non-Newtonian Fluid Mechanics, vol. 163, no. 1–3, pp. 35–44, 2009.

17.

Rocha

G. N.

Poole

R. J.

, and Oliveira

P. J.

, “Bifurcation phenomena in viscoelastic flows through a symmetric 1:4 expansion,” Journal of Non-Newtonian Fluid Mechanics, vol. 141, no. 1, pp. 1–17, 2007.

18.

Hawa

and Rusak

, “The dynamics of a laminar flow in a symmetric channel with a sudden expansion,” Journal of Fluid Mechanics, vol. 436, pp. 283–320, 2001.

19.

Drikakis

, “Bifurcation phenomena in incompressible sudden expansion flows,” Physics of Fluids, vol. 9, no. 1, pp. 76–87, 1997.

20.

Mishra

and Jayaraman

, “Asymmetric flows in planar symmetric channels with large expansion ratio,” International Journal for Numerical Methods in Fluids, vol. 38, no. 10, pp. 945–962, 2002.

21.

De Zilwa

S. R. N.

Khezzar

, and Whitelaw

J. H.

, “Flows through plane sudden-expansions,” International Journal For Numerical Methods in Fluids, vol. 32, pp. 313–329, 2000.

22.

Alves

M. A.

and Poole

R. J.

, “Divergent flow in contractions,” Journal of Non-Newtonian Fluid Mechanics, vol. 144, no. 2–3, pp. 140–148, 2007.

Modeling Asymmetric Flow of Viscoelastic Fluid in Symmetric Planar Sudden Expansion Geometry Based on User-Defined Function in FLUENT CFD Package

Abstract

1. Introduction

2. Numerical Simulation Procedures

2.1. Problem Description

2.2. Numerical Method

3. Results and Discussions

3.1. Mesh Independence Verification

3.2. Results for Newtonian Fluid Cases

3.3. Results for Viscoelastic Fluid Cases

3.3.1. Effect of Elasticity: Weissenberg Number Dependence

3.3.2. Effect of Extensibility: L2 Dependence

3.3.3. Effect of Concentration: β Dependence

3.3.4. Effect of Inertia: Reynolds Number Dependence

3.3.5. Summary of the Effects of Four Parameters and Discussions

4. Conclusion

Footnotes

Nomenclature

Acknowledgments

References

3.3.2. Effect of Extensibility: L² Dependence