air diffusion jet flame

Computers & Fluids 39 (2010) 1381–1389

Contents lists available at ScienceDirect

Computers & Fluids j o u r n a l h o m e p a g e : w w w . e l s e v i e r . c o m / l o c a t e / c o m p fl u i d

Fully explicit implementation of direct numerical simulation for a transient near-field methane/air diffusion jet flame Zhihua Wang *, Yu Lv, Pei He, Junhu Zhou, Kefa Cen State Key Laboratory of Clean Energy Utilization, Zhejiang University, Hangzhou 310027, PR China

a r t i c l e

i n f o

Article history: Received 11 October 2009 Received in revised form 11 April 2010 Accepted 13 April 2010 Available online 24 April 2010 Keywords: DNS Methane Diffusion Jet flame

a b s t r a c t Direct numerical simulation (DNS) is a kind of ultimate numerical simulation tool for studying fundamental turbulent flows, mixing, chemical reactions and interactions among them. In the present work, a fully explicit method of implementing DNS is presented for investigating transient multi-component methane/air jet flame in the near field. The detailed methodology, enclosing non-dimensional governing equations, inlet velocity disturbance, chemical scheme and fluid property, was discussed. An explicit eighth-order finite-difference scheme was used combined with an explicit tenth-order filter. Conservative variables are temporally advanced in two segmented stages that handle Euler terms and viscous terms respectively. A modified non-reflecting boundary condition was used, which has better performance about the characteristic waves on boundary planes. The developed code was firstly tested with an air jet and evaluated in terms of accuracy and parallel efficiency. Then a methane/air combusting jet was simulated to study the characteristics of the chemical heat-release in turbulence. Ó 2010 Elsevier Ltd. All rights reserved.

1. Introduction In real industrial applications, such as gas turbine, boiler furnace and internal-combustion engine, etc., combustion is generally processed combining with the effect of turbulence. Turbulence not only affects the combustion efficiency, but also partially determines the features like temperature distribution and pollutant emissions. Therefore, how to predict those effects has been always one of the central topics in combustion and energy utilization fields. In many cases, Reynolds-averaged Navier–Stokes (RANS) method [1–3] with submodels can provide correct predictions in macroscale, which are usually good enough for specific engineering designs and optimization problems [4–7]. However, the RANS models and different submodels are usually case dependent and omnifarious modes have been developed for varied conditions. More precise and universal models are always needed for different problems. But, the model development usually greatly depends on the understanding of microsale flow mechanisms and turbulenceflame interactions. RANS can offer no help for model development, because the flame wrinkle, species mixing and chemical reaction all take place in smaller physical scales beyond the resolution of RANS. As for large-eddy simulation (LES) [8–10], it has just become realistic to simulate more complex flow and combustion phenomenon, but the subgrid model still need validation from experiment or direct numerical simulation (DNS). * Corresponding author. Tel.: +86 571 87953162; fax: +86 571 87951616. E-mail address: [email protected] (Z. Wang). 0045-7930/$ - see front matter Ó 2010 Elsevier Ltd. All rights reserved. doi:10.1016/j.compfluid.2010.04.007

In DNS, all the length and time scale are resolved directly without any arbitrary models. DNS has been proven to be one of the best tools for fundamental turbulence and combustion phenomenon research. It has been successfully applied to discover the autoignition phenomenon of non-premixed fuel/air mixture [11– 13], the mechanism of vortex-caused strain rate and curve effect on flame propagation [14–16], NOx [17] and soot formation [18] and also the flame-wall interaction [19,20], etc. Chen et al. [16], Poinsot et al. [20] and Vervish and coworker [10] have done a lot of excellent works in this area. The work reported here is propelled by the curiosity on reactive flow in turbulence. In more practical conditions, combusting mixture may include oil droplets and char particles injected into the combustion chambers at the same time, which will raise many unresolved problems on DNS methodology. However, DNS of multi-phase, multi-component, reactive flow with more detailed chemical kinetics and higher Reynolds number is always the targets for us and all DNS researchers which need continuous efforts. As a rudimentary step, a DNS code that can deal with 3D jet flow with combusting mixture is accordingly developed and tested here. DNS tends to resolve all the energy containing and dissipation eddies in flow field with very fine grid system which should be in the order of Kolmogrov length and time scale [21]. The huge computational load is still the main obstacle restricting its application. In the early stage, DNS was usually implemented based on 2D turbulence, simple chemistry and periodic flow assumptions. However, with the great development of computer technology, parallel calculation provides a good path for DNS advancement

1382

Z. Wang et al. / Computers & Fluids 39 (2010) 1381–1389

with spatial and temporal resolved 3D simulation. Additionally, DNS is quite sensitive to the boundary treatment and the proposed non-reflecting condition [22] that assumes characteristic waves always perpendicularly penetrate computational boundaries may not be well adaptable in 3D flow configuration. Therefore, a modified boundary condition should be initially employed. The paper was organized in three parts: firstly, a fully explicit method to implement 3D DNS with multi-component and combustion was given in a detailed and comprehensive description; secondly, parallelized DNS code was tested and evaluated in terms of computational efficiency and accuracy; finally, a methane/air jet flame was simulated and some perceptions about chemistry characteristics in turbulence was discussed.

Then non-dimensional governing equations are listed as:

  n n @ qn @ q uj þ ¼0 @t n @xnj     n @ qn uni unj @ qn uni @pn 1 @ sij þ ¼ nþ n n n @xj Re @xj @xj @t   !   n @ qn Y ns unj @ qn Y ns 1 @ k @Y ns lr þ þ ¼ x @xnj @t n RePr @xnj cnp Les @xnj qr Y r ur s ! n n n n n @ qn T n @ q T uj pn @uj cr 1 @ n @T þ ¼ ð c  1Þ þ k r @t n @xnj cnv @xnj RePr cnv @xnj @xnj

2. DNS implementation

þ In this section, the governing equations for combusting mixture are firstly presented, followed by the treatment of the relevant fluid properties. Afterwards, numerical methods involved here are specified, including the spatial difference, time advancement and filtering, etc. A modified non-reflecting boundary condition for 3D jet flow is specifically highlighted. 2.1. Fluid dynamics 2.1.1. Governing equation and dimensionless method The governing equations are compressible, multi-component Navier–Stokes equations, consistent with those used in previous studies [11–20]. In essence they comprise balance equations for mass, momentum, energy and species concentrations, as shown below:

@ q @ðquj Þ ¼0 þ @xj @t @ðqui Þ @ðqui uj Þ @p @ sij ¼ þ þ @t @xj @xj @xj   @ðqY s Þ @ðqY s uj Þ @ k @Y s þ xs ¼ þ @t @xj @xj cp Les @xj   @ qT @ qTuj @uj @ @T @ui X þ sij ¼ p þ k  h s xs cv þ cv @t @xj @xj @xj @xj @xj s k @T X cps @Y s þ cp @xj s Les @xj

sij

  @ui @uj 2 @uk ¼l þ  dij @xj @xi 3 @xk

ð1Þ

ð2Þ

ð4Þ

v

@xnj



ð8Þ

ð9Þ

ðcr  1Þt r 1 X h s xs cnv s pr ð10Þ

Here, superscript ‘‘n” denotes the non-dimensional variables and subscript ‘‘r” denotes the reference values. Re and Pr are Reynolds number and Prandtl number, respectively, which can be used to characterize the simulated jet flow. The production rate of the sth species xs is solved in subroutine in dimensional forms. In light of the practical consideration, energy equation is rewritten in temperature form, rather than in total energy form, so that it has a more clear relationship with other variables and would be much more easily handled. 2.1.2. Fluid property The enthalpy hs and specific heat at constant pressure cps are both temperature T dependent and calculated through temperature fitted polynomials. The coefficients for those polynomial fitting were taken from CHEMKIN thermal data-package [24]. The specific heat cp of the entire mixture was calculated by weighted-average cps and mass fraction:

X

Y s cps

ð11Þ

s

Then the conductivity of the mixture can be deduced though the following equation [25]:

 0:7 k T ¼ 2:58  105 cp 298

ð12Þ

As for the species viscosity ls, they are all determined by the n exponents equations:

ð5Þ

There are some assumptions underlying these equations. Because the research aims to evaluate turbulent effect on flame, no gravity and buoyancy are specified and gas phase radiation is also neglected here. Soret and Dufour effect [23], considered to be negligible compared with used heat and mass sources, are precluded as well. The above equations should be solved in non-dimensionalized forms avoiding computational truncation errors. According to Wang et al. [21], non-dimensional variables are defined as:

xj ui q t  ur p ; uni ¼ ; qn ¼ ; tn ¼ ; pn ¼ ; qr lr ur lr qr u2r T Ys k l n T n ¼ ; T r ¼ pr =ðqr Rr Þ; Y ns ¼ ; k ¼ ; ln ¼ ; Tr kr Yr lr cp cps cv cnp ¼ ; cnps ¼ ; cnv ¼ ; cv r ¼ cpr  Rr cpr cpr cv r q ur lr l cpr cpr Re ¼ r ; Pr ¼ r ; cr ¼ lr kr cv r

Re

cn

n Y r cr k @T n X cnps @Y ns þ RePr cnp cnv @xnj s Les @xnj

cp ¼ ð3Þ

cr  1 snij @uni

ð7Þ

xnj ¼

ð6Þ



ls T  T0 l0s

n ð13Þ

And the mixture viscosity was calculated as:

l¼

X

Y s ls

ð14Þ

s

In the present study, Lewis number Le was used to quantify the diffusivity of a certain species. Furthermore, all Lewis numbers are assumed to be unity, because oxygen, methane and nitrogen, together, nearly preserve the entire mass of the mixture and the Lewis numbers of these three species are approximate to unity. Here it should be pointed out that unity approximation do cause somewhat effect on flame predication, but the essence of the flame-turbulent interaction does not change, which is exactly the primary topic of this research [26–29]. 2.1.3. Inlet velocity disturbance The conventional way for generating turbulent inflow profile is to impose random fluctuations on the mean inflow velocity. It has

1383

Z. Wang et al. / Computers & Fluids 39 (2010) 1381–1389

been proved that this method does not work well because turbulence kinetic energy will be equally distributed over the whole wave number ranges and cannot represent the essential feature of the modeled turbulence. Here, the way to generate fluctuations on inflow velocity profile follows the procedure proposed by Klein et al. [30]. A digital filtering technique was used to generate time dependent three-dimensional pseudo turbulent inflow velocities with a specified length scale and the magnitude of the fluctuation was chosen to be 4.0% of the mean velocity.

0 u0j;k

¼@

Nx X

Ny X

1

Nz X

bl;m;n  cl;jþm;kþn A  0:04U m

ð15Þ

l¼N x m¼Ny n¼N z

cl,j+m,k+n are a series of random number between 0 and 1; bl,m,n is the array of the coefficients of a digital linear non-recursive filter, which is obtained by the convolution of three dimension components bl,m,n = bl  bn  bm. And,

~= bl  b l

N X

!1=2 ~2 b 1

2

;

1¼N

~ ¼ exp  pl where b l 2n2

! ð16Þ

In the formula n refers to the turbulence scale in local mesh and is deduced from:

n¼

d 5 Dx

ð17Þ

where d is the jet diameter and Dx is the size of local mesh. Detailed explanation of this method can be found in the relevant literature. 2.1.4. Chemical kinetic mechanism The kinetics used here is a six steps methane combustion mechanism developed by Chen and coworker [31] which was reduced based on GRI_Mech 1.2. The simplification is implemented mainly by imposing quasi-steady-state assumption on the properly identified species so that the reduced mechanism only includes 10 species and six global steps. The set of investigated species is comprised of O, O2, H, OH, H2, CH4, CO, CO2, H2O and N2. The performance of the mechanism has been strictly tested in several different flame configurations including perfectly stirred reactor, steady laminar 1D premixed flame and counterflow diffusion flame in Tsuji Geometry, etc. [31]. 2.2. Numerical method 2.2.1. Spatial discretization Explicit eight-order center-difference scheme was adopted to evaluate the local derivatives. Thus, for the interior nodes, the first-order derivatives are computed according to the following equation:

4ðfiþ1  fi1 Þ ðfiþ2  fi2 Þ 4ðfiþ3  fi3 Þ ðfiþ4  fi4 Þ fi0 ¼  þ  5 Dx 5Dx 105Dx 280Dx

ð18Þ

It should be noted that the spatial difference for boundary point should be particularly treated:

1 ð11f 1 þ 18f 2  9f 3 þ 2f 4 Þ 6 Dx 1 f20 ¼ ð2f 1  3f 2 þ 6f 3  f4 Þ 6 Dx 1 f30 ¼ ðf1  8f 2  8f 4  f5 Þ 12Dx 1 f40 ¼ ðf1 þ 9f 2  45f 3 þ 45f 5  9f 6 þ f7 Þ 60Dx

f10 ¼

ð19Þ ð20Þ ð21Þ

Table 1 Coefficients for five-steps fourth-order Runge–Kutta marching scheme. A1 0

A2 73838075589  153778244506

A3  6194124222391 4410992767914

A4  875753879689 434298951106

A5  375951423649 355864889808

B1 69922298306 679754813293

B2

B3

B4

B5

151102391305 203957027379

168159696081 226431017777

189406283338 403426599621

184809359603 982122979183

2.2.2. Time integration The solution of N–S equations is obtained using a segmented time-marching strategy. In this strategy, Euler term and viscous term are handled with different stepping schemes, respectively, due to their obvious difference in quantity. For Euler term including convection and pressure terms, an explicit five-step fourth-order Runge–Kutta scheme was used, followed by an explicit one-order Euler scheme to temporally advance viscous term which involves viscous dissipation, thermal conduction, species diffusion and chemical reaction source item, etc. The Runge–Kutta scheme is symbolized as:

dfj ¼ Aj dfj1 þ Dt  RHSðfj1 Þ fj ¼ fj1 þ Bj dfj

where RHS(f) represents the right-hand sides of differential equation concerning variable f; Dt is the time step determined by Courant Friedrichs Lewy (CFL) criterion, 0.1 was used here. The coefficients, Aj and Bj, are set with the concerns of computational stability and truncation errors and given in Table 1. It has been demonstrated by Carpenter [32] that this time-advancing scheme is superior to the traditional fourth-order Runge–Kutta scheme in terms of saving computational storage and ensuring stability. 2.2.3. Filtering Spatial filtering is quite necessary to suppress numerical fluctuation especially at high wave numbers and maintain stableness for DNS. The numerical fluctuation usually results from truncation errors by exerting finite-difference and will be amplified to pollute the computational solution during time marching. Spatial filtering creates an artificial dissipation to hinder this kind of instability growth. Considering the computational efficiency, here an explicit tenth-order filter introduced by Kennedy and Carpenter [33] was employed which was more computationally economical than Lele’s implicit filter [34]. For the interior points:

^f ¼ f  220  ½252f  210ðf þ f Þ þ 120ðf þ f Þ i i iþ1 i1 iþ2 i2 i  45ðfiþ3 þ fi3 Þ þ 10ðfiþ4 þ fi4 Þ  ðfiþ5 þ fi5 Þ

ð25Þ

^f 2 ¼ f2  220  ð5f þ 26f  55f þ 60f  35f þ 10f  f7 Þ 1 2 3 4 5 6 ð26Þ ^f 3 ¼ f3  220  ð10f  55f þ 126f  155f þ 110f 1 2 3 4 5  45f 6 þ 10f 7  f8 Þ

ð27Þ

^f 4 ¼ f4  220  ð10f þ 60f  155f þ 226f  205f 1 2 3 4 5 þ 120f 6  45f 7 þ 10f 8  f9 Þ ^f 5 ¼ f5  2

As for the second-order derivatives, they are all calculated based on the first-order derivatives in the same way.

ð24Þ

for the points near boundary:

^f 1 ¼ f1  220  ðf1  5f þ 10f  10f þ 5f  f6 Þ 2 3 4 5

20

ð22Þ

ð23Þ

ð28Þ

 ð5f 1  35f 2 þ 110f 3  205f 4 þ 251f 5

 210f 6 þ 120f 7  45f 8 þ 10f 9  f10 Þ

ð29Þ

The numbering is according to the listing sequence by referring from the boundary.

1384

Z. Wang et al. / Computers & Fluids 39 (2010) 1381–1389

2.2.4. Boundary conditions Proper boundary conditions are of great importance for DNS, especially when jet flow with comparatively higher Reynolds number is modeled. Stronger turbulence tends to produce smaller vortex and results in characteristic wave with higher wave numbers. If boundary conditions are not properly imposed, this type of waves will be reflected when passing through the computational boundaries. Consequently, the interior domain will be contaminated and sequentially results will be distorted. To ensure the non-reflecting effect, Poinsot and Lele [22] firstly proposed Navier–Stokes characteristics boundary condition (NSCBC) developed based on the characteristic wave analysis to Euler equations. Later, Baum et al. [35] deduced it to be adaptable for multi-component reactive flows. However, NSCBC has only proven to perform well when flows vertically pass through boundaries, because the contribution of transverse terms [36] to characteristic wave is given no consideration. In the present study, the proposed NSCBC is further modified to enclose this contribution. By truncating all the viscous terms, the governing equations at the boundary normal to x1 were expressed though characteristic wave analysis as follows:

 @q 1 1 þ 2 n2 þ ðn5 þ n1 Þ  w1 ¼ 0 2 @t c @u1 1 þ ðn  n1 Þ  w2 ¼ 0 2qc 5 @t @u2 þ n3  w3 ¼ 0 @t @u3 þ n4  w4 ¼ 0 @t  @T T c1 þ 2 n2 þ ðn5 þ n1 Þ  w5 ¼ 0 @t qc 2 @Z þ n6  w6 ¼ 0 @t

ð30Þ ð31Þ ð32Þ ð33Þ ð34Þ ð35Þ

where ni denotes the magnitude of characteristic waves that go through the mentioned boundary, while wi is the transverse contribution in the boundary plane. Each characteristic wave has one characteristic velocity symbolized as ki. Moreover, in terms of quantification for x1 direction, k1 and k5 are the velocity of sound wave traveling in negative and positive, respectively; k2 is the convection velocity of the entropy wave; k3 and k4 are the advection velocities of u2 and u3 in x1 axis; finally k6 indicates the advection velocity of the other scalar, such as species mass. In summary,

k1 ¼ u 1  c k2 ¼ k3 ¼ k4 ¼ k6 ¼ u1

ð36Þ

k5 ¼ u 1 þ c 2

where c is the frozen speed of sound wave calculated by c = cRT. The expressions of ni and wi are given by

1 1  qc @u @x1 B 2  C C B @q qc2 @T C B k2 c ðcc1Þ @x  cT @x1 C B 1 C B @u2 C B k 3 @x C 1 n¼B C B @u3 C B k4 @x1 C B   C B 2 2 @q q c @u c @T 1 C B k5 þ þ q c c @x1 cT @x1 @x1 A @ @Z k6 @x1 1 0 k  @m @xk C B k C B uk @u C B @xk C B B uk @u2  c2 @T  c2 @ q C B @xk cT @x2 cq @x2 C w¼B Cðk ¼ 2; 3Þ B uk @u3  c2 @T  c2 @ q C B @xk cT @x3 cq @x3 C C B B uk @T  Tðc  1Þ @uk C @xk @xk A @ @Z uk @x 0

k1



c2 @ q c @x1

2

þ qccT

@T @x1

k

ð37Þ

For the flow inlet, a subsonic inflow boundary condition is imposed and the specified variables are temperature, velocity and species mass fraction. According to Thompson [37]’s analysis for Euler equation, n1 is the amplitude of the wave traveling out of the computational domain and was evaluated from interior points based on one-sided difference scheme. In addition to the calculation of w in inlet boundary, n5 and n2 were obtained from Eq. (37) and Eq. (40) in sequence. So in the boundary the density that only needs to be advanced in time was calculated from Eq. (36). Except for the flow inlet, other boundaries are all set with subsonic non-reflecting outflow boundaries. No variables are fixed and only the pressure is set close to the infinity pressure p1. All variables at boundaries were advanced in this type of boundary treatment. In this condition, the characteristic waves, including sound wave exiting the computational domain, entropy convection wave, velocity and species mass advection waves, are moving out through the boundary and were determined by the information in the interior points. For the only incoming wave n1, the wellposed method is used according to Poinsot and Lele [22], which is evaluated based on the pressure information, combined with the special treatment of transverse terms proposed by Yoo [38]:

n1 ¼ r

cð1  M2 Þ ðp  p1 Þ þ ð1  bÞðw5  qcw2 Þ L

ð38Þ

where M is the maximum Mach number for the flow, L is a characteristic length of the domain, and r and b are two constants. This ingoing wave can assure that the pressure adequately approaches the specified value and the transverse terms are properly damped in the outflow boundaries. It has been demonstrated by Poinsot and Lele that r = 0.25 is a reasonable value for ducted shear flows, which was also used in the present study. The damping parameter b of transverse terms is chosen to be the typical Mach number to avoid exerted distortion on the flow field, as suggested by Lodato et al. [36]. After all the ni and wi are calculated, Eqs. (36)–(41) are solved in the outflow boundary to prevent acoustic wave reflecting from the outflow boundary back into the computational domain. In addition, with the above inviscid conditions, the following viscous conditions are complementally adopted to implement the application of Navier–Stokes boundary conditions:

  @ s12 @ s13 @q1 @ @Z ¼0 ¼ ¼ ¼ qD @x1 @x1 @x1 @x1 @x1

ð39Þ

Table 2 The computation procedure of our DNS code. Initialization: Obtain one data block which will be handled in local CPU Set initial profiles for density q, velocity ui, temperature T, species mass fraction Ys for each time-step do for each node in local block do Compute and store species reaction rate xs by calling chemical subroutine end for Exchange necessary data with every neighbor blocks for each node in local block do Compute Euler terms (ET) for q, qui, qT, qYs for each Runge–Kutta substep n do Temporal stepping: D(qn) = An  D(qn1) + Dt  ET(qn), qn = qn1 + Bn  D(qn) Update qui, qT, qYs in the same way end for end for Exchange necessary data with every neighbor blocks for each node in local block do Compute viscous terms (VT) for q, qui, qT, qYs Temporal stepping: qm = qm1 + Dt  VT(qm1) Update qui, qT, qYs in the same way Compute new ui, T, Ys from updated q, qui, qT, qYs Calculate new P based on new ui, T, Ys end for end for

1385

Z. Wang et al. / Computers & Fluids 39 (2010) 1381–1389

2.3. Pseudo-code

4.0

The code is necessarily developed in parallel way due to the considerable computational expense of DNS. To parallelize the code, the computational domain was divided into many blocks which contain different segments of the total data and each block was designated to an independent CPU. The necessary data at boundary of each block were transferred to all spatially contiguous blocks through message passing interface (MPI) functions. A detailed pseudo-code is given in Table 2 to outline the solution procedure of our DNS code. Based on this scheme it is easy to found that our code is suitable to simulation at low Mach number because shock wave is not considered. However, to avoid modeling error and offer authentic DNS, compressible effect is still preserved.

3.5

2×4×8 3×4×5

3.0

Speed-up

2.0

2×4×6

1.5

2×3×10

1.0 0.5 0.0 48

52

56

60

64

Number of CPUs used

3. Test for parallelized code

Fig. 1. Parallel computing efficiency.

1.2

Exp. (Panchapakesan and Lumley) 17d 18d 19d

1.0

0.8

U/Ucentral

The DNS program was coded in C language. For parallel computing, the domain should be divided into suitable blocks which contain different segments of the total data and treated by separate processor. The necessary data exchanges among processors are realized by MPI functions through network. To evaluate the parallel speed-up and computational accuracy, one air jet flow was first simulated. The diameter of the central jet is 1 mm. The inlet velocity is 120 m/s surrounded with a low speed co-flow at 0.5 m/s. Additionally, the entire domain was initialized at room temperature with size of 8d  8d  20d. The meshes were assigned to be 160  160  390 (z  y  x) at three propagation directions. 48 to 64 CPUs were used to test the parallel efficiency. The machine used in the present study is an IBM blade cluster system equipped with Intel quad core 2.0 GHz CPUs and 10 GB/s Infiniband network and Linux operation system.

4×4×4

2×4×7

2.5

0.6

0.4

0.2

3.1. Parallel efficiency

3.2. Accuracy evaluation

0.0 0.0

0.8

1.6

2.4

3.2

4.0

4.8

L/R half Fig. 2. Radial profile of Favre-averaged jet velocity. Normalized using centerline velocity and half-width.

0.35

Exp. (Panchapakesan and Lumley) 17d 18d 19d

0.30

Uc

0.25

0.20

Ruu1/2/

The processors distribution in three directions has obvious influence on the parallel efficiency besides computational scale and CPU numbers. If processors are not topologically properly assigned, large proportion of running time could be wasted in communications between processors, of course depending on the configured network and processor, eventually undermining the parallel performance. In Fig. 1 speed-up generally shows an increasing tendency with the rising-up of CPU amount but grows at a slightly less than linear speed. In 60-CPU runs, when CPUs are oriented 2  3  10, the speed-up is only 0.61 greatly less than that of the 3  4  5 orientation with the same CPU numbers. It could be explained by the reason that more communications are required in this orientation, especially in x-dimension, so that more time is expensed to wait for sending and receiving data while computational capacity of local processors are not fully utilized. Based on the above test, we fixed the parallel processor arrangement to be 2  4  8 in the following applications.

0.15

0.10

The Favre-averaged jet velocity and Reynolds stress along the axis are first normalized and then distributed along the spanwise direction. Three sections are used to obtain the relevant statistics as shown in Figs. 2 and 3. By comparing the simulated results with experimental data by Panchapakesan and Lumley [39], the accuracy of the code was well demonstrated. A perfect self-similarity of the air jet flow is predicted after 16d downstream where turbulent flow is fully developed [40]. Reynolds stress, as a second-order variable, is correctly computed as well, further warranting the

0.05

0.00 0.0

0.8

1.6

2.4

3.2

4.0

L/R half Fig. 3. Radial profile of streamwise Reynolds stress Ruu ¼ U 02 : Normalized using DUc = (U1 + U2)/2 and half-width.

1386

Z. Wang et al. / Computers & Fluids 39 (2010) 1381–1389

v¼2

k

qc p

ðrZÞ2

ð43Þ

Low values of this variable represent good mixing. This variable can also indicate the vortex exertions on flame front. When flame element is compressed, mixture fraction gradient locally increases and the local scalar dissipation rate is intensified. 4.1. Flame structure overview

Fig. 4. DNS computational mesh.

applicability of this code. Therefore, this code can be further considered for reactive jet flow simulation including sophisticated chemical kinetics. 4. Methane/air diffusion flame A methane/air jet flame was accordingly simulated using more than 10 million meshes and the physical domain is shown in Fig. 4. The meshes in Y and Z directions are locally refined while maintaining uniformly in X direction. To get relatively finer meshes in central location, the following function [41] is adopted to stretch meshes in Y and Z direction:

yl ¼ 0:6

sin h½sðl  0:5Þ sin hð0:5sÞ

ð40Þ

Here the parameter s controls the stretching intensity; l refers to the original coordinate in coordinate system; and yl is the new coordinate in the same coordinate system. In the present work, we set s to be equal to 0.5, and then the refining effect has been well shown in Fig. 4. Mixture in the central jet consists of 20% methane and 80% nitrogen in volume at temperature about 400 K. A round nozzle is employed to import central main jet with 120 m/s along X-axis. The around hot air is introduced as 0.5 m/s low velocity co-flow. The air co-flow is at 1200 K. The inlet velocities are consistent with those used in above air-jet case. ‘‘Top-hat” smoothing strategy specifies the velocity profile by a hyperbolic tangent function:

U¼

r U1 þ U2 U1  U2 þ tan h 2h 2 2

Fig. 5 shows the 3D structure of flame front, defined by the isosurface of Z = 0.2; the contours of OH concentration are also presented at three cross-sections along streamwise. As vortex develop into more fine structure along streamwise, flame front is wrinkled more seriously. The flame structure is well resolved and the turbulence stretch effect on reaction zone is clearly observed on OH contours. Fig. 6a–c shows the contours of CH4, CO, OH concentrations at central section respectively and Fig. 6d shows the contours of heat release rate (HRR) at central section. It is evident that locations with higher OH concentrations generally experience relatively larger HRRs, which agree well that OH can be used as indicator for reaction zone [43,44]. Along with the flow development, reaction zone becomes thicker due to the intensified wrinkling effect of vortex. Quite a few small OH pools, found inside the flame enclosure, tend to be fast depleted, which stems from flame–flame interaction [45]. Thus, CO can distribute more homogeneously downstream because of the quenching of internal OH oxidizer layers. The subtle flame structure and scalar fields resolved here imply that the proposed DNS method and imposed boundary conditions can effectively guarantee the non-reflecting effect during the simulation and no computational contamination occurs. 4.2. Mixture fraction effect on HRR The conditional PDF between heat release rate (HRR) and mixture fraction is shown in Fig. 7. An ignition process is found from 0.033 ms to 0.130 ms due to species-diffusion development, and then the maximum of heat release rate shows a steady characteristic. The most reactive zone ZMR, which refers to the mixture fraction with relatively larger HRR, can be observed at all four moments, consistent with the results in literatures [11,13,46]. Viggiano and Magi [11] have found that ZMR corresponds to relatively leaner mixture at lower fuel temperature, which better explains the reason why the predicted ZMR in this DNS is located at around 0.2. Therefore, in the present study the flame front can be defined

ð41Þ

where U1, U2 are jet velocity and co-flow velocity respectively; r refers to the jet radius and h is the momentum thickness. Stanley et al. [42] have summarized the scale of r/h and it should be around 20 based on the current Reynolds number. Thus this value was used here. The mixture fraction is defined following as [25]:

Z¼

2ðY C  Y C;2 Þ=M C þ ðY H  Y H;2 Þ=2M H  ðY O  Y O;2 Þ=M O 2ðY C;1  Y C;2 Þ=M C þ ðY H;1  Y H;2 Þ=2M H  ðY O;1  Y O;2 Þ=M O ð42Þ

where subscript ‘‘1” and ‘‘2” denote fuel mixture in central jet and co-flow air respectively; subscript ‘‘C”, ‘‘H” and ‘‘O” refer to carbon, hydrogen and oxygen elements respectively. To quantify the species mixing, scalar dissipation rate is accordingly defined as [19]:

Fig. 5. 3D flame front profile and OH contours at three cross-sections.

1387

Z. Wang et al. / Computers & Fluids 39 (2010) 1381–1389

0.2 0.1 0

0.2

0.4

0.6

0.3 0.2 0.1 0

0.8

Mixture fraction

0.2

0.4

0.6

0.3 0.2 0.1 0

0.8

Mixture fraction

0.033ms

0.4

Heat Release Rate

0.3

0.4

Heat Release Rate

0.4

Heat Release Rate

Heat Release Rate

Fig. 6. Central cross-section snapshot of flame profile at 0.310 ms including CH4 volume fraction profile (a), CO volume fraction profile (b), OH volume fraction profile (c) and normalized heat release rates (d).

0.2

0.4

0.6

0.4 0.3 0.2 0.1 0

0.8

Mixture fraction

0.066ms

0.2

0.4

0.6

0.8

Mixture fraction

0.130ms

0.310ms

Fig. 7. Conditional PDF between heat release rate and mixture fraction at four different moments.

Linear Slope = 6.8E-6

0.1

5000

10000

15000

20000

0.3

0.2

Linear Slope=7.6E-6

0.1

0

5000 10000 15000 20000 25000 30000

Heat Release Rate

0.2

Heat Release Rate

Heat Release Rate

0.3

0

0.4

0.4

0.4

0.3 Linear Slope = 7.6E-6

0.2

0.1

0

5000

10000

15000

Scalar Dissipation Rate

Scalar Dissipation Rate

Scalar Dissipation Rate

(a)

(b)

(c)

Fig. 8. Conditional PDFs between heat release rate and scalar dissipation rate: (a) 0.066 ms, (b) 0.130 ms (c) 0.310 ms.

1388

Z. Wang et al. / Computers & Fluids 39 (2010) 1381–1389

Fig. 9. Visualized relationships between scalar dissipation rate and heat release rate: (a) 0.066 ms, (b) 0.130 ms (c) 0.310 ms.

OH concentrations (see Fig. 6c and d) which indicate faster oxidation rate of fuel. Based on the above analysis, it concludes that the effect of scalar dissipation rate on HHR actually varies during flow transition. The flame located in refined meshes is typical 0.3 mm thick, while the mesh sizes vary from 30 to 40 lm. Thus there are 7–10 meshes to resolute flame front, which tend to be a reasonable number compared with previous study [26]. Due to vortex induced stretch, some flame elements tend to be thickened as shown in Fig. 10 and the framed region in Fig. 9b is zoomed to show the grid resolution. It is evident that the spike-type profiles of minor radicals and reaction rate are well resolved.

5. Summary Fig. 10. Zoomed local region in Fig. 9b to show grid resolution.

as the iso-surface of Z = 0.2. The HRR fluctuation at each mixturefraction bin may result from the difference of local thermo-chemistry states or turbulence effect. 4.3. Scalar dissipation rate effect on HRR The conditional PDF between HRR and scalar dissipation rate DisR at flame front is shown in Fig. 8. The contours corresponding to these two variables are shown in Fig. 9 to offer a visualized demonstration. At 0.066 ms moment when species mixing is poor, compression effect of vortex on flame front can lead to a local higher scalar dissipation rate. Meanwhile, heat transfers are intensified and mixtures are preheated faster at those locations, leading to enhanced reactions and higher HRRs, which well represent the obviously positive correlation between HRR and scalar dissipation rate (see Fig. 8a). This mechanism is dominating at 0.066 ms, and gradually becomes weak as species mixing is further developed. From 0.130 ms to 0.310 ms, another branch with low scalar dissipation rate and high HRR is gradually formed (see Fig. 8b and c), because when species are well mixed, local thermo-chemical state become another critical factor to control chemistry, such as local temperature and detailed species components. Thus, it is reasonably observed that some sites with higher HRRs are located at low scalar dissipation rate – good mixing – region (see Fig. 9c) and higher

A DNS code used to investigate reactive multi-component jet flow was developed on the basis of parallel computation. The governing equations, non-dimensionalization, chemistry integration, numerical methods are presented in detail in this paper. A modified boundary condition, which can better quantify the characteristic waves through the boundary, was used to ensure nonreflecting effect. Afterwards, a pure air jet was simulated first to evaluate the accuracy and parallel efficiency of the code. Finally, a transient jet flow with combusting methane/air mixture was simulated using this DNS code and general features about this jet flame was analyzed. An ignition process along with the development of the jet flow is predicted and the most reactive zone is found. The transient characteristics of the effect of scalar dissipation rate on heat release rate during flow transition is initially identified and discussed.

Acknowledgement This study was implemented under the financial support from National Natural Science Foundation of China (50806066). We are so grateful for that. References [1] Launder BE, Spalding DB. Lectures in mathematical turbulence. London, England: Academic Press; 1972.

models

of

Z. Wang et al. / Computers & Fluids 39 (2010) 1381–1389 [2] Patankar SV, Spalding DB. A calculation procedure for heat, mass and momentum transfer in three-dimensional parabolic flows. Int J Heat Mass Transfer 1972;15:1787–806. [3] Launder BE, Reece GJ, Rodi W. Progress in the development of a Reynoldsstress turbulence closure. J Fluid Mech 1975;68:537–66. [4] Watanabe H, Otaka M. Numerical simulation of coal gasification in entrained flow coal gasifier. Fuel 2006;85:1935–43. [5] Stanmore BR, Visona SP. Prediction of NO emissions from a number of coalfired power station boilers. Fuel Process Technol 2000;64:25–46. [6] Rizk NK, Mongia HC. Three-dimensional analysis of gas turbine combustors. J Propul Power 1991;7:445–51. [7] Reitz R, Rutland CJ. Development and testing of diesel engine CFD models. Prog Energy Combust Sci 1995;21:173–96. [8] Wang P, Bai XS. Large eddy simulation of turbulent premixed flames using level-set G-equation. Proc Combust Inst 2005;30:583–91. [9] Lesieur M, Metais O, Compte P. Large-eddy simulations of turbulence. Cambridge, UK: Cambridge University Press; 2006. [10] Veynante D, Vervisch L. Turbulent combustion modeling. Prog Energy Combust Sci 2002;28:193–266. [11] Viggiana A, Magi V. A 2-D investigation of n-heptane autoignition by means of direct numerical simulation. Combust Flame 2004;137:432–43. [12] Doom J, Mahesh K. Direct numerical simulation of auto-ignition of a hydrogen vortex ring reacting with hot air. Combust Flame 2009;156:813–25. [13] Hibert R, Thévenin D. Autpignition of turbulent non-premixed flames investigated using direct numerical simulation. Combust Flame 2002;128:22–37. [14] Choi CW, Puri IK. Contribution of curvature to flame-stretch effects on premixed flames. Combust Flame 2001;126:1640–54. [15] Chakraborty N, Cant RS. Effects of strain rate and curvature on surface density function transport in turbulent premixed flames in the thin reaction zones regime. Phys Fluids 2005;17:065108 [1–15]. [16] Chen JH, Echekki T. Analysis of the contribution of curvature to premixed flame propagation. Combust Flame 1999;118:308–11. [17] Bédat B, Egolfopoulos FN, Poinsot T. Direct numerical simulation of heat release and NOx formation in turbulent nonpremixed flames. Combust Flame 1999;119:69–83. [18] Lignell DO, Chen JH, Smith PJ, Lu T, Law CK. The effect of flame structure on soot formation and transport in turbulent nonpremixed flames using direct numerical simulation. Combust Flame 2007;151:2–28. [19] Bharadwaj N, Safta C, Madnia CK. Flame-wall interaction for a non-premixed flame propelled by a vortex ring. Combust Theor Model 2007;11:1–19. [20] Poinsot TJ, Haworth DC, Bruneaux G. Direct simulation and modeling of flamewall interaction for premixed turbulent combustion. Combust Flame 1993;95:118–32. [21] Wang Z, Zhou J, Fan J, Cen K. Direct numerical simulation of ozone injection technology for NOx control in flue gas. Energy Fuels 2006;20:2432–8. [22] Poinsot TJ, Lele SK. Boundary conditions for direct simulations of compressible viscous flows. J Comput Phys 1992;101:104–29. [23] Popp P, Baum M. Analysis of wall heat fluxes, reaction mechanisms, and unburnt hydrocarbons during the head-on quenching of a laminar methane flame. Combust Flame 1997;108:327–48. [24] Kee RJ, Rupley FM, Miller JA, Coltrin ME, Grcar JF, Meeks E, et al. CHEMKIN collection, Release 3.5, Reaction Design, Inc., San Diego, CA 94551-0969; 1999.

1389

[25] James S, Jaberi FA. Large scale simulations of two-dimensional nonpremixed methane jet flames. Combust Flame 2000;123:465–87. [26] Jenkins KW, Klein M, Chakraborty N, Cant RS. Effects of strain rate and curvature on the propagation of a spherical flame kernel in the thin-reactionzones regime Combustion and Flame 2006;145:415–34. [27] Jenkins KW, Cant RS. Curvature effects on flame kernels in turbulent environment. Proc Combust Inst 2002;29:2023–9. [28] Thévenin D. Three-dimensional direct simulations and structure of expanding turbulent methane flames. Proc Combust Inst 2005;20:629–37. [29] Chakraborty N, Klein M, Cant RS. Stretch rate effects on displacement speed in turbulent premixed flame kernels in the thin reaction zones regime. Proc Combust Inst 2007;31:1385–92. [30] Klein M, Sadiki A, Janicka J. A digital filter based generation of inflow data for spatially developing direct numerical or large eddy simulations. J Comput Phys 2003;186:652–65. [31] Chang WC, Chen JY. Reduced mechanisms based on GRI_Mech 1.2. . [32] Carpenter MH, Kennedy CA. Fourth-order 2n-storage Runge–Kutta schemes. NASA Tech. Memorandum 1994:109112. [33] Kennedy CA, Carpenter MH. Several new numerical-methods for compressible shear-layer simulations. Appl Numer Math 1994;14:397–433. [34] Lele SK. Compact finite difference schemes with spectral-like resolution. J Comput Phys 1992;103:16–42. [35] Baum M, Poinsot T, Thévenin D. Accurate boundary conditions for multicomponent reactive flows. J Comput Phys 1993;116:247–61. [36] Lodato G, Domingo P, Vervisch L. Three-dimensional boundary conditions for direct and large-eddy simulation of compressible viscous flows. J Comput Phys 2008;227:5105–43. [37] Thompson KW. Time-dependent boundary conditions for hyperbolic systems. J Comput Phys 1987;68:1–24. [38] Yoo CS, Wang Y, Trouvé, IM HG. Characteristic boundary conditions for direct simulations of turbulent counter flow flames. Combust Theor Model 2005;9:617–46. [39] Panchapakesan NR, Lumley JL. Turbulence measurements in axisymmetric jets of air and helium. Part 1. Air jet. J Fluid Mech 1992;246:197–223. [40] Anders JW, Magi V, Abraham J. Large-eddy simulation in the near-field of a transient multi-component gas jet with density gradients. Comput Fluids 2007;36:1609–20. [41] Uzun A, Koutsavdis K, Blaisdell GA, Lyrintzis AS. Direct numerical simulation of 3-D jets using compact schemes and spatial filtering. AIAA paper 2001-2543; 2001. [42] Stanley SA, Sarkar S, Mellado JP. A Study of the flow-field evolution and mixing in a planar turbulent jet using direct numerical simulation. J Fluid Mech 2002;450:377–407. [43] Vervisch L, Poinsot T. Direct numerical simulation of non-premixed turbulent flames. Annu Rev Fluid Mech 1998;30:655–91. [44] Baum M, Poinsot TJ, Haworth DC, Darabiha N. Direct Numerical Simulation of H2/O2/N2 flames with complex chemistry in two-dimensional turbulent flows. J Fluid Mech 1994;281:1–32. [45] Hawkes ER, Chen JH. Direct numerical simulation of hydrogen-enriched lean premixed methane–air flames. Combust Flame 2004;138:242–58. [46] Sreedhara S, Lakshmisha KN. Direct numerical simulation of autoignition in a non-premixed, turbulent medium. Proc Combust Inst 2000;28:25–34.