c
TU Delft, The Netherlands, 2006
APPLICATION OF LOCAL DEFECT CORRECTION IN A
3D TURBULENT CHANNEL FLOW
J. de Hoogh∗, J.G.M. Kuerten†
∗Technische Universiteit Eindhoven,Faculty W
P.O box 513, 5600 MB Eindhoven, The Netherlands e-mail:j.d.hoogh@tue.nl
†Technische Universiteit Eindhoven,Faculty W
P.O box 513, 5600 MB Eindhoven, The Netherlands e-mail:j.g.m.kuerten@tue.nl
Key words: Turbulence, Local Defect Correction, Passive Scalar
Abstract. To investigate the behavior of the concentration of a passive scalar in a tur-bulent flow with realistic high Schmidt numbers, one needs a very fine grid to capture all length-scales that are present in the concentration. An efficient method to increase the spatial resolution without a drastic raise of the computational costs is to apply a Local De-fect Correction method. To compute the velocity field of a turbulent channel flow, a Direct Numerical Simulation is used that consists of a pseudo-spectral method in two periodic directions and a Chebyshev expansion in the wall-normal direction. The concentration equation for the passive scalar is solved using a finite volume method. The main concern with an LDC method in a hyperbolic time-dependent problem is the interpolation of the concentration when the refinement area is moved. The interpolated solution has to be smooth and continuous to ensure numerical stability. Using the mean averaged radius and the total surface area of the concentration, the behaviour of several numerical methods can be investigated.
1 Introduction
fluid, the governing equation for particle concentration is a convection-diffusion equation. The convection velocity equals the local fluid velocity and the diffusion coefficient depends on the properties of particles and fluid molecules. In many applications involving small particles in water or air, the diffusivity is orders of magnitudes smaller than the kinematic fluid viscosity. This implies that structures in the particle concentration are much smaller than the structures in the fluid velocity and a very fine grid is needed to accurately capture all scales. If the particles are distributed over the total computational domain this goes beyond presently available computer capacities. However, if the particles are mostly located in a confined region, a very fine local grid could be superposed on a global coarse grid to capture all relevant structures. This situation resembles that of a stationary fine-scaled flame in a coarser stationary fluid flow, which has been studied by Anthonissen 1. His work was based on previous work by Ferket 2 and uses local defect
correction (LDC). In LDC first a global coarse grid solution is calculated, which is used as a boundary condition for a local fine grid solution in the region where the small structures are located. The difference between the fine and coarse grid solutions on the fine grid region is a measure for the error on the coarse grid. Subsequently the coarse grid solution is improved by a defect term in the governing equation based on this error. The present problem differs from the original setting in which LDC was applied, since the present problem is governed by a time-dependent hyperbolic equation and is three dimensional. In a time-dependent problem the local fine grid region should move with the fluid flow. Minero 3 already applied LDC to a two-dimensional convection-diffusion equation. The
main problem in extension to three spatial dimensions is that the number of grid points needed to maintain stability is too large in three dimensions. In this paper it will be investigated how the LDC method can be applied to a three-dimensional convection-diffusion equation, describing the concentration of a passive particle in turbulent channel flow.
In section 2 the governing equations for fluid and particles and their discretization will be presented, as well as several numerical schemes that are used to test the LDC scheme. In section 3 the implemented numerical methods are visually and quantitatively compared. Section 4 will describe the basic principles of the LDC method and section 5 shows some results that are computed with LDC.
2 Governing equations and numerical method
In the DNS the three-dimensional Navier-Stokes equation for incompressible flow is solved in the vorticity formulation. The equation reads:
∂u
∂t + ω × u + 1
ρ∇p = ν∆u (1)
Here, P denotes the total pressure, P = p+21ρu2, u the fluid velocity, ω the fluid vorticity,
In the DNS the computational domain consists of a channel of a finite length equal to 4πH and finite width equal to 2πH, where H is half the height of the channel. In both directions periodic boundary conditions are applied. This makes a spectral method with a Fourier-Galerkin approach in the two periodic directions the natural choice. In the wall-normal direction a Chebyshev-collocation method is applied.
For integration in time a second-order accurate method is chosen, which is a combi-nation of the explicit Adams-Bashforth method for the nonlinear terms and the implicit Crank-Nicolson method for the linear terms. In the first step the nonlinear terms are calculated pseudo-spectrally by fast Fourier transform, where the 3/2-rule is applied to prevent aliasing errors. In the second step the pressure is calculated and in the last step the viscous terms are treated implicitly. Since the Fourier modes are completely decou-pled, only a one-dimensional system of equations needs to be solved for each Fourier mode. The walls of the channel act as no-slip walls. In order to satisfy the continuity equation the velocity field is made divergence free within machine accuracy following the approach proposed by Kleiser and Schumann 4 applied to the collocation approximation 5. The
mean axial pressure gradient is chosen in such a way that the volume flow, and hence the Reynolds number based on the bulk velocity, remains constant. The simulations are started from an arbitrary initial solution. After a large number of time steps a state of statistically stationary turbulent flow is reached. All results shown have been calculated at a Reynolds number based on the bulk velocity, Re = UbH/ν = 2800 where Ub is the
bulk velocity. This corresponds to a Reynolds number based on the friction velocity uτ
of 180,
uτ = s
νdux
dz |z=H (2)
where x, y, and z, denote the stream wise, span wise and wall normal coordinate. At this Reynolds number results can be compared with DNS results by Moser et al.6. As
an example the fluctuations of the three velocity components are compared in figure 1. It can be seen that the agreement is good. The same conclusion can be drawn for other quantities, such as the terms in the balance equation for the turbulent kinetic energy. 2.1 Concentration equation for a passive scalar
The DNS code is extended with a convection-diffusion equation for the concentration of a passive scalar. The governing equation reads,
∂c
∂t + u · ∇c = D·c (3)
where D is the diffusivity. The equations of fluid and particle concentration are non-dimensionalized with half the channel height and the friction velocity. When these scaling quantities are applied to Eq. (3), this results in,
∂c
∂t + u · ∇c = 1
0
50
100
150
0
0.5
1
1.5
2
2.5
3
y
+rms(u)
wall−normal
spanwise
streamwise
Figure 1: Velocity fluctuations as functions of the wall-normal coordinate wall units. Lines: present results, symbols: Moser et al6
With Re the Reynolds number and Sc the Schmidt number that is defined as the ratio of the kinematic viscosity ν and the diffusivity D. The numerical method that is used for the fluid velocity, is not well suited for the concentration. The velocity has many points located near the wall, but less point in the center of the channel. Due to the initial condition, the concentration will only be unequal to zero in the center of the channel and therefore has no need for many points in the wall region. Secondly, a spectral method is not suited for grid refinement techniques, because the solution inside the refinement area is not periodic.
For the computation of the particle concentration a finite volume is implemented on a Cartesian mesh. The convection-diffusion equation is integrated over small control volumes that are formed by connecting the midpoints of the grid cells closest to a grid point. Since the Cartesian mesh for the concentration is different from the grid of the fluid velocity, the latter has to be converted towards the Cartesian mesh on every time step. This is done using a fourth-order Hermite-Lagrange interpolation method 7.
To allow the use of the finite volume method, Eq. (4) is rewritten, using Gauss’s theorem, into ZZZ V ∂c ∂tdV + ZZ A(u · ncdA) = 1 ReSc ZZ A ∂c ∂ndA, (5)
unit-normal on the surface. In discrete form at time level n Eq. (5) is written as, (cn j,k,l− cn−1j,k,l)Vj,k,l ∆t = 6 X i=1 ( 1 ReScQd n j,k,l− Qc n j,k,l) (6) with Qdn
j,k,l and Qcnj,k,l the fluxes for the convection and diffusion terms resp. with the
indices j, k and l denoting the stream wise, wall-normal and span wise point indices of the control volumes. The total flux that leaves or enters a control volume is computed through summation over i of the partial fluxes through each surface of the control volume. On a Cartesian grid, the surface area is, in each direction, independent of the position and can be taken as a constant.
2.2 Selection of the numerical method
The flux terms of the concentration have to be evaluated at the surface of each control volume and not at the center where the concentration is known. Several numerical schemes are available to obtain the concentration at the surface. To investigate which numerical method gives the best results, a comparison between different schemes is made.
To keep the comparison between all available numerical schemes manageable, three different methods are implemented: central differencing, Roe’s first order upwind method and a 3rd order accurate upwind method. Central differencing is the most obvious choice
because it is second order accurate and fast. However, due to the absence of numerical diffusion, oscillations might occur. These oscillations are usually due to the large grid size and can be prevented using an upwind method. Note that adding numerical diffusion and increasing the Schmidt number at the same time is not very meaningful. Increasing the Schmidt number will cause smaller scaled structures to appear and the added numerical diffusion will smoothen the solution, removing some of these details.
The choice for these three specific methods is based on the properties of LDC and the added numerical diffusion. An expensive, high order accurate numerical scheme is not needed, because LDC improves the accuracy and stability of the solution through the applied defect correction.
2.2.1 Diffusive flux
In the remaining part of the paper only one direction of the flow is treated and the indices k and l are dropped to improve readability. The diffusive flux is independent of the implemented method, because it only requires the gradient on the surface of a control volume and becomes,
Qdj = Aj(
cj+1− 2cj + cj−1
∆x ), (7)
2.2.2 Convective flux
The convective part of the flux is more difficult to compute and requires a different approach. For this purpose, following a method found in 8, Qc
j is first split up into,
Qcj = Aj(Q1|j+1 2 − Q1|j− 1 2 + Q2|j− 1 2 − Q2|j+ 1 2), (8)
where the index Q1|j+12 denotes the flux Q1 evaluated at the surface between point j and
j + 1. For central differencing, the terms from Eq. (8) become, Q1|j+12 =
1
2(uj+1cj+1+ ujcj) (9)
Q2|j+1
2 = 0.
To add numerical diffusion to the numerical scheme, Roe’s first order upwind method can be used. The coefficients are almost identical to the central difference method, with the exception that a 2nd order derivative is added, that depends on the flow direction. The flux terms for Roe’s first order upwind method become,
Q1|j+1 2 = 1 2(uj+1cj+1+ ujcj) (10) Q2|j+12 = 1 4{|uj+1+ uj|(cj − cj+1)}.
To check the total amount of numerical diffusion that is actually added, using Roe’s first order upwind method, a third order scheme is implemented as well,
Q1|j+12 =
1
12(−uj+2cj+2+ 7uj+1cj+1+ 7ujcj− uj−1cj−1) (11) Q2|j+12 =
1
12(−|uj+2|cj+2+ 3|uj+1|cj+1− 3|uj|cj+ |uj−1|cj−1).
The last method is more accurate then the two other schemes, but does require a larger stencil and is therefore more expensive to compute. On the boundaries of the computa-tional domain and the edge of the refinement area extra effort has to be taken to supply the boundary condition. A suitable approach for this is the usage of the central difference method.
Time integration is performed using a standard first order accurate Euler forward method.
3 Numerical experiments
quantitatively compared. The initial condition is given by a Gaussian blob of concentra-tion, c0j,k,l= 1 σx √ 2πe −(x−x0)2 σ2 x −(y−y0)2 σ2 y −(z−z0)2 σ2 z (12)
where σx, σy and σz are the standard deviations of the Gaussian distribution in
wall-normal, span wise and stream wise direction respectively. The parameters x0, y0 and z0
are the offset of the blob and are chosen in such a way that the concentration at the boundaries of the computational domain is approximately equal to zero. The values of the standard deviations are set to σx = σy = 5 · 10−2 and σz = 5 · 10−1. The choice for a
larger standard deviation in stream wise direction is unavoidable because of the amount of points that would otherwise be needed in that direction. To investigate the influence of the Schmidt number and the numerical methods used to compute the convective flux, the initial shape is not important. The grid size of the computational domain is set to ∆x = 641, ∆y = π
128 and ∆z = π
256. Note that a semi grid-refinement technique is applied
to obtain the small grid size in stream wise direction. The concentration is only computed over a length of 2π and the grid is shifted as the concentration moves along with the fluid velocity.
All experiments are done on the same system and using the same initial setup for the concentration and the velocity. The Reynolds number is kept constant to remove its influence on the development of the scalar. In figure 2 and 3 the three different numerical methods are visually compared using a contour plot of concentration in the XY-plane through the center of the channel, using a Schmidt number of Sc = 1.0 and Sc = 0.1. To give a more quantitative analysis the total surface area of the concentration is computed using a method that is developed by Geurts9. With this method it is possible to compute
the surface area of a specific value of the concentration. An increase in surface area tells something about the size of the structures that are present within the solution and about the amount of numerical diffusion. Figure 4 shows the development of the surface area as a function of time.
3.1 Conclusion on the numerical schemes
From figures 2,3 and 4, three conclusion about the implemented numerical methods can be drawn.
• The first order upwind method adds too much numerical diffusion and causes some of the important details to be removed. This method is not suitable for computation with Schmidt numbers that are larger than unity.
• The central differencing and the 3rd order accurate upwind method give almost
identical results, with only the small details in the 3rdorder upwind method removed
Y X −0.5 0 0.5 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 5.5 Y X −0.5 0 0.5 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 5.5 Y X −0.5 0 0.5 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 5.5
Figure 2: Contour plots of concentration through the center of the blob with Sc = 1.0 at t = 0.14H/uτ.
Left: central differencing; Middle: 3rd order upwind; Right: Roe’s first order upwind method.
Y X −0.5 0 0.5 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 5.5 Y X −0.5 0 0.5 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 5.5 Y X −0.5 0 0.5 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 5.5
Figure 3: Contour plots of concentration through the center of the blob with Sc = 0.1, at t = 0.14H/uτ.
Left: central differencing; Middle: 3rd order upwind; Right: Roe’s first order upwind method.
0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18 0.2 10 15 20 25 30 35 Time (H/u tau) Area
Time development of the surface area
Roe first order upwind, Sc=1.0 Central difference, Sc=1.0 3rd order upwind, Sc=1.0 Roe first order upwind, Sc=0.1 Central difference, Sc=0.1 3rd order upwind, Sc=0.1
Figure 4: The development of the surface area at two different Schmidt numbers, using all 3 available computational methods. The top lines represent the iso-surface at c = 1 · 10−2 and the bottom lines are
at c = 1 · 10−1.
numerical diffusion present in the upwind method has no influence on the solution. Based on these results, the central difference method gives the best representation of the physics of the problem and is therefore implemented into the LDC method. The principle behind LDC will be explained in the next section.
4 General principle of LDC
The principle of LDC 1,3 is based on the assumption that a solution computed on a
the fine grid.
1. Compute the coarse grid approximation using the implicit Crank-Nicolson scheme, Cn j − C n−1 j ∆T + Qn j 2Vj +Q n−1 j 2Vj = 0 (13) where Qn
j contains both the convective and the diffusive flux terms from Eq. (6).
2. Compute the fine grid approximation using Euler forward time integration. cn j − cn−1j ∆t + qjn−1 vj = 0 (14)
Because of the difference in grid size, the time step applied to the coarse grid is larger than the time step of the fine grid. To be able to compute the best possible solution that is needed in step 3 and the defect term from step 5, both the coarse and fine approximations need to be at the same position in time. On the fine grid multiple steps have to be taken to reach the next coarse time level. For each time step on the fine grid, the boundary conditions are interpolated in time using linear interpolation. Note that the time interpolation also needs to be applied to the fluid velocity to evaluate the flux terms on the corresponding time level.
3. Assemble the best possible coarse approximation for the concentration and the flux by combining the coarse and the fine grid approximations. The best possible con-centration can be obtained by taking,
Cbnj (
Cn
j Outside the refinement area
cn
k Inside the refinement area
(15) where k denotes those points where the fine and coarse grid coincide. The best flux approximation Qbn
j can be computed by inserting the new best solution for the
concentration into Eq. (7) and Eq. (8). Figure 5 shows a schematic overview of a refinement area with the corresponding coarse grid.
4. Determine the defect terms by taking the difference between the coarse approxima-tion and the best possible soluapproxima-tion for both the concentraapproxima-tion and the flux. This results in the defect terms to become
Figure 5: Schematic overview of a refinement area with the corresponding coarse grid. The coarse grid points are denoted with an (×) , the fine grid points using a (·) The fluxes through the surfaces of the control volumes are visualized using different sized arrows.
5. Correct the global approximation by applying the defect, as a source term, to Eq. (13), (Cn j − C n−1 j ) − Dcnj ∆T + Qn j − Dqjn 2Vj + Q n−1 j − Dqjn−1 2Vj = 0 (18)
and solve the updated equation to obtain the corrected coarse solution Ccjn.
6. Update the boundary conditions of the refinement area and recompute the local solution.
7. Iterate from point 3 until convergence has been reached and move towards the next time level when convergence is reached.
In a time-dependent problem, the total number of iteration steps that are needed to reach convergence is small. This is due to the fact that the solution at the previous time step is a good initial condition for the iteration process. Numerical experiments have shown that one or two LDC-iterations are sufficient to reach convergence. If many iterations would be needed, this would severely reduce the overall efficiency of LDC in a time dependent problem.
Fur-thermore, it keeps the difference in grid size between two subsequent levels of refinement small. The latter is important for a proper estimate of the updated boundary conditions and the interpolation of the concentration that is needed when the refinement area is moved. Figure 6 shows a schematic drawing for one LDC time step containing 3 levels of refinement, using one LDC iteration. After each correction the underlying fine grid has to be recomputed using the new boundary conditions from the corrected coarse con-centration. These last steps are straightforward and have been left out of figure 6 to maintain readability. A side note towards time integration, when using multiple levels
n+2/3 t n+1/3 t ttn+1 tn n+2/9 t n+1/9 t Top level (1) Intermediate level (2) (3) (4) (5) (8) (9) (10) (13) (14) (15) (6) (2)
Intermediate level correction (11) Top level Correction (17)
(12) (16)
Finest level
Figure 6: The time levels that are needed in a complete LDC cycle.
of refinement, is that only on the finest level an explicit time integration scheme can be applied. On every coarser grid, an implicit scheme is necessary to be able to apply the defect correction.
5 LDC in practice
In the previous section it was shown that the central difference method gives the best representation of the physics of the problem. Using LDC, the grid size on the finest level of refinement, in stream wise direction, can be much smaller (∆z = 1.25 ·10−3H) than the
grid size that is practically reachable using the standard approach (∆z = 1.25 · 10−2H).
The decrease in grid size has several advantages. The smaller structures, that appear when the Schmidt number is increased, can be captured in more detail. Secondly, due to the correction step applied to the coarser grids and the small grid size on the finest levels, LDC is very stable. A third advantage is that the method is very flexible. The total amount of refinement levels and the refinement factor can be chosen freely. The refinement factor is the ratio between the grid size on the coarse and fine grid. Despite the smaller grid size, the computational time that is needed to perform one LDC iteration is much less then the time that would be needed to compute a conventional method that uses the same grid size everywhere. The total amount of points in 3D that are needed to satisfy the requirements for a Total Variation Diminishing scheme for the central difference method would be larger then 1 · 1010. Using LDC, the actual amount of grid points on
decrease in the amount of points that need to be treated. However, the computation and the application of the defect term require additional computational effort, but this is small compared to the time needed to get the same accuracy using a conventional approach. A side note to this large reduction in grid points is that the refinement area needs to be small compared to the whole computational domain.
5.1 Numerical experiments
To check the application of LDC to a three-dimensional convection-diffusion equation, a sphere of concentration is released in the same turbulent flow. The initial concentration field is chosen to be spherical with a standard deviation of σ = 5 · 10−2 in each direction.
This initial condition results in much steeper gradients in stream wise direction that were impossible to compute using a regular method. The Schmidt number is set to unity and the time step of the coarsest grid is equal to ∆t = 8.5 · 10−5H/uτ. A total of four stacked
levels of refinement are used to reach a final grid size of approximately 2.5 · 10−3H in
wall-normal and span wise direction and 1.25 · 10−3H in stream wise direction. The time
step on each level decreases accordingly and matches the grid size. Figure 7 shows contour plots of concentration at the two topmost levels of refinement at six different moments in time. From figure 7 it can be concluded that it is possible to compute the development
3.6 3.8 4 4.2 4.4 −0.2 −0.1 0 0.1 0.2 (1) 4.2 4.4 4.6 4.8 5 −0.2 −0.1 0 0.1 0.2 (2) 4.8 5 5.2 5.4 5.6 −0.2 −0.1 0 0.1 0.2 (3) 5.2 5.4 5.6 5.8 6 −0.2 −0.1 0 0.1 0.2 (4) 5.8 6 6.2 6.4 6.6 −0.2 −0.1 0 0.1 0.2 (5) 6.2 6.4 6.6 6.8 7 7.2 −0.3 −0.2 −0.1 0 0.1 0.2 (6)
Figure 7: Contour plot of the two topmost levels of refinement area at six different moments in time. The time sequences goes from left to right and top to bottom, ranging from t = 0 to t = 0.3H/utau in
steps of ∆t = 0.06H/utau
further research. Future work will focus on increasing the overall efficiency of the LDC algorithm and increasing the Schmidt number even further.
6 Acknowledgment
This work was sponsored by the National Computing Facilities Foundation, NCF for the use of supercomputer facilities, with financial support from the Netherlands Orga-nization for Scientific Research, NWO. The author would like to thank the Netherlands Organization of Scientific Research (NWO) for funding this research, Computational Sci-ence grant 635.000.002.
REFERENCES
[1] Anthonissen, M.J.H., Local defect correction techniques: analysis and application to combustion, Ph.D. Thesis, Eindhoven University of Technology, 2001.
[2] Ferket, P.J.J., Solving boundary value problems on compsite grids with an application to combustion, PhD Thesis, Eindhoven University of Technology, 1996
[3] Minero, R., Anthonissen, M.J.H., Mattheij, R.M.M, Solveing parabolic problems us-ing local defect correction in combination with the finite volume method, Eindhoven University of Technology, 2005.
[4] Kleiser, L. and Schumann, U., ”Treatment of incompressibility and boundary con-ditions in 3-D numerical spectral simulations of plane channel flows,” in Proceeding of the Third GAMM-Conference on Numerical Methods in Fluid Mechanics, edited by E.H. Hirschel (Vieweg-Verlag, Braunschweig, 1980), pp. 165–173.
[5] Canuto, C., Hussaini, M.Y., Quarteroni, A. and Zang, T.A., Spectral methods in fluid dynamics (Springer-Verlag, Berlin,1988).
[6] Moser, R.D.,Kim, J., and Mansour, N.N., ”Direct numerical simulation of turbulent channel flow up to Reτ = 590,” Phys. Fluids 11, 943-945 (1999).
[7] Veenman, M.P.B., Statistical Analysis of turbulent pipe flow: a numerical approach, Ph.D Thesis, Eindhoven University of Technology, 2004.
[8] Kuerten J.G.M., Computational Fluid Dynamics III, Lecture notes, University of Twente, 1999.
