Large-eddy simulation with parabolized stability equations for turbulent transition using OpenFOAM

Computers and Fluids 189 (2019) 108–117

Contents lists available at ScienceDirect

Computers and Fluids journal homepage: www.elsevier.com/locate/compfluid

Benchmark solutions

Large-eddy simulation with parabolized stability equations for turbulent transition using OpenFOAM Minwoo Kim a, Jiseop Lim a, Seungtae Kim a, Solkeun Jee a,∗, Jaeyoung Park b, Donghun Park b a b

School of Mechanical Engineering, Gwangju Institute of Science and Technology, 123 Cheomdan-gwagi-ro, Buk-gu, Gwangju 61005, South Korea Department of Aerospace Engineering, Pusan National University, 2 Busandaehak-ro 63beon-gil, Geumjeoung-gu, Busan 46241, South Korea

a r t i c l e

i n f o

Article history: Received 25 November 2018 Revised 27 February 2019 Accepted 17 April 2019 Available online 17 May 2019 Keywords: Turbulent transition Large-eddy simulation Parabolized stability equations Openfoam

a b s t r a c t Laminar-to-turbulent transition is simulated with large-eddy simulation (LES) coupled with parabolized stability equations (PSE). A canonical transition on a flat plate is computed. The PSE-coupled approach accurately reproduces the subharmonic resonance and subsequent turbulent boundary layer using the solver OpenFOAM. Various LES inlet locations are simulated to reduce the computational domain. Nonlinear PSE provides appropriate instabilities at various inlet locations. The computational cost and domain are reduced by 30% with the LES inlet close to the final stage of the nonlinear transition region, compared to the baseline, upstream LES inlet. Current LES results are compared with relevant data in the literature.

1. Introduction Laminar-to-turbulent transition is a critical phenomenon in boundary layer flows. The ability to predict the boundary layer transition is significant in estimating the drag force exerted on a moving object. The prediction capability also allows us to determine the state of the boundary layer in a region where the flow may separate due to pressure gradient. While the extension of the laminar region can be beneficial to drag reduction for a streamlined body, turbulent boundary layer is sometimes necessary to avoid the laminar separation. For high-speed vehicles, the boundary layer transition is associated with heat transfer from the heated gas to the vehicle surface, and as such the transition location is one of key parameters to determine thermal protection systems. Likewise, it is significant to be able to predict the transition in designing a moving object. Several theoretical approaches have been established for the understanding of boundary layer transition. Such theoretical approaches are based on equations of disturbances which are obtained from the Navier–Stokes equations subtracted by the undisturbed (mean) equations. Linear stability theory (LST) [1,2] provides how disturbances evolve in the laminar region, particularly in the region where disturbances are linearly amplified. LST can indicate dominant instabilities which can trigger nonlinear interactions downstream once the magnitude of the instabilities is

∗

Corresponding author. E-mail address: [email protected] (S. Jee).

https://doi.org/10.1016/j.compfluid.2019.04.010 0045-7930/© 2019 Elsevier Ltd. All rights reserved.

© 2019 Elsevier Ltd. All rights reserved.

not negligible. Albeit LST is powerful and widely used for the stability analysis, it is limited to the linear-growth region. In addition, LST is based on the assumption of quasi-parallel flow, limiting the stability analysis. Parabolized stability equations (PSE) were developed with taking into account of non-parallel flow effects [3]. Disturbance growth in the boundary layer is better captured in PSE, rather than LST, because of the growth of the boundary layer itself [4]. PSE is also able to capture nonlinear interactions between instabilities, which allows the stability analysis in the nonlinear region [4–7]. A downside of PSE is that the computational cost grows exponentially as the number of significant instabilities increases, which inhibits to apply PSE into the subsequent turbulent region. Instead of disturbance equations in LST and PSE, the Navier– Stokes equations can be numerically solved to simulate transitional boundary layer. Direct-numerical simulation (DNS) can be used to resolve small-magnitude instabilities in the pre-transitional region and their interactions (linear and nonlinear) with high accuracy. Joslin et al. [5] conducted both DNS and nonlinear PSE, validating both the approaches for the canonical transition on the flat plate of Kachanov and Levchenko [8]. DNS was also carried for the transition on the flat plate in Fasel et al. [9]. Both the DNS studies [5,9] were focused on the transition region, not covering the downstream turbulent region. As the computational power becomes sufficient, a complete transition process has been simulated with the DNS approach for the canonical case [10,11]. A well-known limitation of DNS is its immense computational cost which increases further for complex flow. Computational cost in eddy-resolving simulation depends on the Reynolds number Re. While the cost of DNS is proportional

M. Kim, J. Lim and S. Kim et al. / Computers and Fluids 189 (2019) 108–117

to Re2.64 for turbulent boundary layer, it is reduced to Re1.85 with large-eddy simulation (LES) [12]. Taking the advantage of the reduced cost in the turbulent region, the LES approach has been investigated for a complete transitional process [11,13–15]. Early LES studies [13,14] provided instabilities with blowing and suction on a wall. Although this wall treatment is able to reproduce dominant instabilities in the experiment of the reference [8], it is not straightforward to apply their 5th-order polynomial function of the wall-normal velocity to other transition processes. This is probably because instabilities are generated via their own receptivity process from the wall blowing/suction. Recently, LES coupled with stability theory has been explored, introducing precise instabilities into the computational domain [11,15]. Jee et al. [15] validated an LES framework coupled with LST against the experiment [8] and the DNS study [10]. Lozano-Durán et al. [11] investigated the LES framework coupled with PSE for the same transition case of [8]. It is noted that both wall-resolved (WRLES) and wall-modeled LES (WMLES) approaches were studied in the reference [11]. While WRLES provides DNS-like prediction capabilities such as the overshoot of the skin friction near the end of the transition and the transitional length, WMLES was not able to reproduce such features [11]. It can be conjectured that WRLES is suitable to capture details in the late transition stage where the log layer, which is desirable for WMLES, is not yet fully developed in the flow simulation. One major goal of the current work is to investigate the LES framework coupled with stability theory for a complete transition process using OpenFOAM. It would be useful for the community to be able to use an open-source code for computations of boundary layer transition with DNS-like fidelity. The previous LES study [15] tested the LST-coupled LES in a proprietary flow solver which is not accessible in the public domain. Since PSE can handle nonparallelism of the mean flow and nonlinear interactions between instabilities, the PSE-coupled LES is pursued here. Following two previous LES studies [11,15], WRLES is used in the current study. The second goal is to reduce the computational domain by locating the LES inlet further downstream. Nonlinear PSE is useful here because it provides appropriate instabilities at various inlet locations for LES. Moving the LES inlet close to the end of the transition could significantly reduce the size of the computational domain, bringing down the computational cost. This approach is similar to Lozano-Durán et al. [11] where two LES inlet locations are used. WRLES is also used in the current investigation because the reference [11] reported that WRLES outperforms WMLES in capturing the transition and the subsequent turbulent boundary layer at relatively low Reynolds numbers. The paper is organized as follows. Numerical methods in the current study are presented in Section 2 including PSE and LES approaches, computational domains, and numerical schemes. Current results are discussed in Section 3, which is further divided in two subsections corresponding to the two main objects. 2. Numerical methods A framework of LES combined with PSE is investigated in this study. Main computations are conducted with the filtered incompressible Navier–Stokes equations, and the computational domain sufficiently extends to turbulent flow. PSE is used to provide instabilities at a given location in laminar boundary layer, which is the inlet for the LES domain. The following decompositions in Eqs. (1) and (2) to a total variable φˇ are used in this paper.

φˇ =  + φ

(1)

where  is the mean (undisturbed) and φ is the fluctuation, i.e., instability in this paper. The variable φˇ can be decomposed to the

109

spatially filtered φ¯ and the filtered-out, residual part φ  for LES:

φˇ = φ¯ + φ 

(2)

It should be noted that φˇ = φ¯ in LES of laminar flow, yielding no distinction between DNS and LES in the laminar region. The residual part φ  is a turbulence fluctuation, so φ  = 0 in laminar flow. Note that the instability in the laminar boundary layer is not turbulence. In this study, allvariables  are scaled with the boundary layer thickness δ˜0 at R = Rex = U∞ x/ν = 400 and the freestream velocity U˜∞ , unless notified otherwise. Dimensional variables are denoted by the tilde such as φ˜ . 2.1. Parabolized stability equations The parabolized stability equations (PSE) is an efficient method in describing linear and nonlinear evolution of disturbances [3]. PSE takes into account of the flow non-parallelism and thus capture the effect of mean flow growth on disturbances and vice versa. With the capability of incorporating nonlinear interactions between instabilities, PSE has been used for the transition region in both incompressible [5–7] and compressible [4] boundary layers. It is well recognized that PSE is an excellent method for a weakly nonlinear region where, for subharmonic resonance studied here, amplified subharmonics are still small in the amplitude so that the back influence of subharmonics on the fundamental wave is negligible [6,7,16]. In this study, PSE computations are used up to the early nonlinear stage, but not all the way through the late nonlinear region and final breakdown. LES is carried out for the downstream nonlinear and the subsequent turbulent regions. A formal approach of PSE is to decompose the disturbance φ into a fast-varying oscillatory-wave part and a slowly-varying shape function, using a Fourier expansion, as written in Eq. (3).

φ (x, y, z, t ) =

Nm 

Nn 

x φˆ (m,n) (x, y ) exp i α(m,n) (ξ )dξ

m=−Nm n=−Nn

+ inβ z − imωt



x0

(3)

where the streamwise, wall-normal, and spanwise directions are denoted by x, y, and z, respectively. The wave part is assumed periodic in the spanwise direction and time, and is governed by the complex wave number α (m,n) where the subscript m and n indicates the temporal and spanwise wave modes, respectively. The shape function φˆ (m,n ) is a complex function of both x and y. The spanwise wave number β and the angular frequency ω are real for spatially growing instabilities studied here. A set of partial differential equations for the shape functions with unknown variable α (m,n) is obtained, using the disturbance expression (Eq. (3)) in the disturbance governing equation which is derived from the Navier–Stokes equations for given fundamental frequency and spanwise wave number. These equations, then, are parabolized by neglecting the second derivatives (∂ 2 /∂ x2 ) of the shape functions based on the assumption that the change of shape function along the streamwise direction is of order of 1/R and the terms of O(1/R) are negligible. The PSE code for a generalized coordinate system of Park and Park [7,17] is used for the present study with sufficiently small Mach number Mach = 0.03. The PSE code can handle both linear and nonlinear PSE analyses. In this study, nonlinear PSE is carried out to examine nonlinear evolution of disturbances. The initial condition of the nonlinear PSE (i.e. for each of fundamental and subharmonic mode) is obtained from linear PSE analysis which starts from a further upstream location. The solution of linear PSE is chosen here to take into account the flow non-parallelism even in the initial condition and to minimize possible transient phenomena during downstream marching. It should be noted

110

M. Kim, J. Lim and S. Kim et al. / Computers and Fluids 189 (2019) 108–117

Fig. 1. Schematic diagram of the current computational domain.

that the solution from either LST or linear PSE can be used as the initial condition, since they are found to be almost identical at the initial location R = 400 considered in the present study. Additional nonlinear PSE with LST results as the initial condition was also conducted, and it was confirmed that no noticeable differences to the current nonlinear PSE with linear PSE as the initial condition. This implies that the flow non-parallelism is not significant enough to affect the stability result near R = 400. Nonlinear PSE is conducted in the domain of 400 ≤ R ≤ 700. The number of the temporal and spanwise modes are set to Nm = 6 and Nn = 3, respectively, keeping total 28 modes including the mean distortion, the (0,0) mode. The 4th-order central scheme is used for the wall-normal derivatives, and the 2nd-order backward scheme is employed to integrate the equations in the streamwise direction. Uniform grids with 107 points are used for the streamwise direction. At least 80 points are placed in the boundary layer with the total wall-normal points of 220 in the wall-normal domain size of 40δ 0 . Further details of disturbance equations and numerical approaches for PSE are well documented in the reference [7,17]. 2.2. Large-eddy simulation The filtered incompressible Navier–Stokes Eqs. (4) and (5) are numerically solved in the present LES computations using the code OpenFOAM v.4.1.

∂ u¯ i =0 ∂ xi

(4)

∂τirj 1 ∂ p¯ ∂ u¯ j ∂ 2 u¯ j ∂ + (u¯ i u¯ j ) = ν − − ∂t ∂ xi ∂ xi ∂ xi ∂ xi ρ ∂ x j

(5)

Variables with the overbar are filtered (resolved) quantities as described in Eq. (2). It is required to model the residual stress tensor τirj in the momentum Eq. (5). A common approach is to use the assumption that the stress tensor is proportional to the resolved strain rate S¯i j .

τirj = −2νt S¯i j S¯i j =

1 2



(6)

∂ u¯ i ∂ u¯ j + ∂ x j ∂ xi

 (7)

Then, the turbulence viscosity ν t needs to be modeled. The wall-adapting local eddy-viscosity (WALE) model [18] is selected in this study. The WALE model is given by Eq. (8).

νt = (Cw )2

(Sidj Sidj )3/2 (S¯i j S¯i j )5/2 + (Sidj Sidj )5/4

(8)

where Sidj is the traceless, symmetric tensor of the square of the velocity gradient tensor.

Sidj

1 = 2



∂ u¯ i ∂ u¯ k ∂ u¯ j ∂ u¯ k + ∂ xk ∂ x j ∂ xk ∂ xi



−

1 ∂ u¯ k ∂ u¯ k δi j 3 ∂ xk ∂ xk

(9)

The model coefficient is Cw = 0.325, following the reference [19]. Our study indicated that the turbulent transition is not sensitive to the exact value of the model coefficient at least in the range 0.325 ≤ Cw ≤ 0.5. The idea of the implicit filtering is adopted here with = ( x y z )1/3 . The WALE model is designed for the appropriate behavior of the turbulent viscosity near the wall, ν t ∼ O(y3 ). The WALE model provides almost zero turbulence viscosity for two-dimensional shear flow such as laminar boundary layer. As such, WALE would be a proper sub-grid-scale model for the current study. 2.3. Computational domain The current computational domain consists of two parts, as shown in Fig. 1. Laminar simulation is conducted on a twodimensional slice of the whole domain consisted of the upstream, the main LES and the outlet LES zones. Then, LES is carried out only in the three-dimensional LES zones. At the LES inlet, instabilities in the boundary layer are introduced from PSE computations. The main LES zone extends to R = 10 0 0, which is sufficiently large enough for fully turbulent boundary layer after the transition. The height of the domain is Ly = 230δ0 . The width of the domain is Lz = λz where λz is the wavelength of the spanwise disturbance. Three inserts in Fig. 1 indicate the three velocity components on the inlet plane, which are obtained from the superposition of disturbances into the laminar solution. The outlet LES zone with coarse grids is chosen here to minimize any unfavorable effect from the outlet boundary. Based on additional computations without the outlet zone, the computational results reported in this study are not visually sensitive to the presence of the outlet zone. For the validation of the PSE-coupled LES approach against the experiment [8], the main LES region starts at R = 400, close to the location of the vibrating ribbon in the experiment. Following the 2 × 106 = 124 test condition, the fundamental wave of F = 2π f˜ν˜ /U˜∞ and a pair of subharmonic wave are introduced at the LES inlet of R = 400. The current PSE computations provide the mode shape of instability modes, some of which are used for the LES inlet. For the inlet at R = 400, two modes, fundamental (2,0) and subharmonic (1,1) are used. The amplitude (rms) of the 2D Tollmien-Schlichting (TS) and the subharmonic oblique modes at R = 400 are determined from the reference [15]. The wave number of the 2D TS ˜ x = 0.677, and the spanwise wave numwave is α(2,0 ) = 2π δ˜0 /λ

M. Kim, J. Lim and S. Kim et al. / Computers and Fluids 189 (2019) 108–117 Table 1 Amplitude of instability modes at current LES inlets. √ R = Rex 400 500 641 A(0,0) A(2,0) A(1,1) A(4,0) A(3,1)

– 3.3 × 10−3 7.1 × 10−5 – –

– 5.2 × 10−3 8.5 × 10−5 – –

4.4 × 10−4 9.9 × 10−3 5.6 × 10−3 4.1 × 10−4 4.2 × 10−4

Table 2 Previous and current numerical studies of the subharmonic resonance of the experiment [8]. Case

Method

Spatial accuracy

Range of R

Fasel [9] Joslin et al. [5] Sayadi et al. [10] Lozano-Durán et al. [11] Huai et al. [13] Sayadi & Moin [14] Lozano-Durán et al. [11] Jee et al. [15]

DNS DNS DNS DNS WRLES WRLES WRLES WRLES

4th-O in x & y, spectral in z 4th-O in x & y, spectral in z 4th-O 2nd-O 4th-O 6th-O 2nd-O 2nd-O

Current

WRLES

2nd-O

Joslin et al. [5] Esfahanian et al. [6] Park & Park [7] Current

PSE PSE PSE PSE

2nd-O 2nd-O 2nd-O 2nd-O

378–928 426–700 316–1030 424–870 360–800 316–1030 632–870 40 0–10 0 0 40 0–10 0 0 50 0–10 0 0 641–10 0 0 426–670 400–720 393 - 680 400 - 700

in in in in

x, x, x, x,

4th-O 4th-O 4th-O 4th-O

in in in in

y y y y

˜ z = 0.6506. For the ber of the subharmonic mode is β = 2π δ˜0 /λ angular frequency of the TS wave, ω = 2π f˜δ˜0 /U˜∞ = 0.2477 is used. The influence of the inlet location on the current computations is also investigated in this study. Two additional downstream locations are tested: R = 500 and R = 641. In all locations, the fundamental and subharmonic modes are given by the nonlinear PSE. Additional modes of A > 10−4 are also introduced at the most downstream location because their magnitude is not negligible. The amplitude of instability modes at current LES inlets is listed in Table 1. The subharmonic resonance of Kachanov and Levchenko [8] has been reproduced in numerical studies in the past decades [3,5–7,9–11,13–15]. To model the vibrating ribbon in the test, a blowing-suction method on the wall has been explored [9,10,13,14]. The wall blowing and suction requires a fifth-order polynomial with fine-tuned coefficients for the polynomial function, which presumably hinders the general application of the blowing and suction method. A direct introduce of disturbances in the eddyresolving framework (LES or DNS) has been proposed based on the boundary layer stability analysis [5,11,15]. Imposing stabilityanalysis-based disturbances is relatively straightforward even if the disturbance location is changed, where the blowing and suction method needs to determine the empirical polynomial function. Selected previous numerical studies are listed in Table 2 with the turbulence treatment, the spatial accuracy and the streamwise domain size in each case. 2.4. Numerical schemes and grids LES Eqs. (4) and (5) are numerically solved using OpenFOAM v.4.1. The second-order Crank–Nicolson is used for the temporal discretization, and the second-order central difference scheme is used for the spatial discretization. The pressure-implicit with splitting of operators (PISO) algorithm [20] is used for the pressure solver. Unsteady boundary conditions are required for the LES inlet due to instabilities which are periodic in time. Mode shapes obtained from PSE at a selected location are represented with a 7th-order

111

Fourier series. Then, the inlet fluctuation is calculated at each inlet grid point at each time step. Finally, the time-dependent inlet condition is generated by superimposing the instability component to the laminar solution. All boundary conditions are listed in Table 3. The size of the current time step is determined, based on the period of the fundamental instability T(2,0) . Three values of the time step are tested: t = T(2,0 ) /512, T(2,0) /256, and T(2,0) /128. The smallest time step size is t + = 0.085 in the wall unit. Our study indicates that the t = T(2,0 ) /256 is small enough, providing no significant difference between the case of t = T(2,0 ) /512 and

t = T(2,0) /256. A typical CFL number for t = T(2,0) /256 is 2.5 on the fine grid. The mode shape for two disturbance modes at R = 400 is shown in Fig. 2. Both the amplitude and the phase profiles are plotted for the 2D TS and 3D subharmonic oblique waves. The maximum amplitude of the component u of the 2D TS wave appears near y/δ ∼ 0.26, whereas the maximum of the component v is around 0.7δ with about one third amplitude of the u amplitude. The phase difference between the u and the v modes is about π /2 for y/δ ࣠ 0.7. Near y/δ ∼ 0.7, the phase of the u disturbance is shifted by π due to the sign change. The subharmonic oblique mode includes three components u, v and w, as shown in Fig. 2(c) and (d). Similar to the fundamental mode, the streamwise u component is dominant for the subharmonic. The maximum amplitude of u and the second dominant component w are located around y/δ ࣃ 0.25. The wall-normal component v is the smallest one in the oblique mode. The phase difference between u and w is about π in the boundary layer, indicating that these two components are almost out of phase. The phase of the u component of the 2D TS wave is relatively constant within the boundary layer except near the wall, compared to the oblique wave. Grids used in the current LES are listed in Table 4. The wavelength of the disturbance is used to determine the spatial resolution. The streamwise wavelength λx of the TS wave helps to guide the streamwise grid size x , and the spanwise wavelength λz of the subharmonic oblique mode for the spanwise grid size

z . Choi and Moin [12] suggested the following requirement for WRLES: x+  100, z+  30, and the number of wall-normal grid points in a turbulent boundary layer Ny,TBL ࣡ 80. The current fine grid fully satisfies this resolution requirement. The medium grid marginally satisfies the requirement with both z and Ny,TBL near the criteria. Note that the first wall normal grid point is located in the same position in both the fine and medium grids. The first wall normal resolution y+ ∼ 1 would be required 1,max for LES to capture the transition and the subsequent turbulent boundary layer. It is expected that the fine grid allows us to conduct WRLES with sufficient grid points in all three directions. Computations on the fine grid for the longest domain R = 400 − 1000 are conducted using 512 cores of the Intel Xeon X5570 processor in Korea Institute of Science and Technology Information. One flow through time requires about 85 wall clock hours. Initial 1.3 flow through times are required to allow numerical artifacts from the initial laminar condition to leave completely the computational domain. Flow statistics is accumulated over four periods of the TS wave 4T(2,0) , which is about 0.15 flow through time. Statistics is not sensitive between the two time windows 4T(2,0) and 8T(2,0) . 3. Results and discussion 3.1. LES of subharmonic resonance One of the main goals of the current study is to propose the PSE-coupled LES method in the OpenFOAM framework. The canonical subharmonic resonance of Kachanov and Levchenko [8] is

112

M. Kim, J. Lim and S. Kim et al. / Computers and Fluids 189 (2019) 108–117 Table 3 Boundary conditions of current LES computations. The Einstein summation convention is used for the velocity ui . The direction n indicate the boundary-normal direction.

Velocity Pressure

Inlet

Bottom

Sides

Top

Outlet

Dirichlet ∇p = 0

no-slip wall ∂ p/∂ n = 0

periodic periodic

∂ ui /∂ n = 0

convective outlet p=0

p=0

Fig. 2. Amplitude and phase profiles of the 2D fundamental and the subharmonic oblique waves at R = 400.

Table 4 Two grids used in the current LES. Grid level

Fine

Medium

x

y1

z

x+

y+1,max

z+

λx /64 δ 0 /80 λz /64

λx /32 δ 0 /80 λz /32

14 1.3 15 80 160 3088 2976 112 256 64

28 1.3 30 42 82 1544 1488 56 128 32

Ny,LBL Ny,TBL Nx Nx in Nx in Ny Nz

at R = 400 at R = 10 0 0 main LES zone outlet LES zone

simulated in the current study. Two instabilities, a plane TollmienSchlichting wave and a pair of subharmonic oblique waves, are obtained in the current PSE computation. These two modes are superposed into the laminar solution at the LES inlet. Current wall-resolved LES computations are conducted with the WALE model for the sub-grid-scale stress tensor. Instantaneous vortical structures in the overall computational domain are shown in Fig. 3. The current domain is duplicated in

the span just for the visualization. The Q criteria, which is based on the second invariant of the velocity gradient tensor [21], is used for the visualization. Two dimensional TS waves are dominant in the upstream near the inlet because (1) the amplitude of the 2D TS is about 50 times as large as the subharmonic oblique wave, and (2) the oblique mode requires a certain distance to grow exponentially. It is observed that the nominal, two  dimensional shape in the span begins to undulate around R = Rex = 630. Then, -shaped structures appear between 670 ࣠ R ࣠ 720. Staggered  structures are associated to the subharmonic breakdown, which is often called either H-type [22] or N-type transition [23] in the literature. The staggered pattern is also observed in DNS studies [10,11]. Thus, current computations can resolve the secondary instability associated to subharmonic growth. Near R = 720,  structures break down to small eddies, and the boundary layer completely transits to turbulent flow in the downstream. Hairpin vortices are captured well in the turbulent boundary layer (see Fig. 3c where large Q values are used for the visualization of small-scale eddies). Laminar-to-turbulent transition can be clearly depicted in the skin friction on the wall, as shown in Fig. 4. Current LES on the fine grid is compared to DNS data [10,11] along with the Blasius andthe turbulent skin friction. For the Blasius skin friction, 2 C f = 0.664 Rex . The “exact” formulation C f = 0.455/ ln (0.06Rex ) from the asymptotic analysis is used for the turbulent skin

M. Kim, J. Lim and S. Kim et al. / Computers and Fluids 189 (2019) 108–117

113

Fig. 3. Instantaneous vortical structures visualized with the iso-volume of Q-criterion colored by normalized streamwise velocity: Q > 0 for (a) and (b); Q > 350 for (c). Fine grid.

Fig. 4. Skin friction in the current LES and two previous DNS studies [10,11].

friction (see Ch. 6.6 of White [24]). Current LES computations accurately reproduce the deviation from the laminar solution and also the overshoot of the skin friction in the late stage of the transition. Evolution of instabilities in the boundary layer is shown in Fig. 5. The amplitude of four selected modes are calculated from the current LES and PSE computations. The fundamental mode

Fig. 5. Amplitude of four selected instabilities in the current LES on the fine grid and in the current PSE.

is indicated by the subscript (2,0) and the subharmonic by (1,1). Two higher harmonics (3,1) and (4,0) are also plotted. The current PSE was already validated with the experiment in the reference [7]. Current LES results are matched well with the PSE analysis. Two plane modes F(2,0) and F(4,0) grow exponentially in the range of 450 ≤ R ≤ 600. Note that the linear growth of F(2,0) and F(4,0) in the logarithmic plot is associated to the exponential growth in R. The amplification of two oblique modes F(1,1) and F(3,1) is double-exponential (exponent-in-exponent). The amplitude of the subharmonic mode surpasses the amplitude of the fundamental mode around R = 650 where the spanwise undulation in vortical structures is clear in Fig. 3. The significant influence of amplified oblique modes on the plane modes occurs after R = 650. Staggered- structures are observed around R = 680 (see Fig. 3) where all instabilities grow exponentially in Fig. 5. The amplification of two given instabilities is compared to the experiment [8] and two DNS studies [10,11], as shown in Fig. 6. Current LES on the fine grid provides accurate growth of these two modes observed in the experiment. Note that instabilities are generated with a wall blowing and suction around R = 416 in the DNS study of Sayadi et al. [10], so the initial amplitude near R = 400 is different to the current and the experiment. It

Fig. 6. Amplitude of the fundamental and the subharmonic modes in the current LES, the previous experiment [8] and the two previous DNS studies [10,11].

114

M. Kim, J. Lim and S. Kim et al. / Computers and Fluids 189 (2019) 108–117

Fig. 7. Mean velocity U + at R = 837.

seems that the amplitude of the subharmonic mode yields insignificant impact on the subharmonic growth, based on the DNS of Lozano-Durán et al. [11] where the subharmonic amplitude is one order of the magnitude lower compared to the other data sets. The amplification of the experiment [8] and the DNS [10] is reported up to R ࣃ 650. The subharmonic amplification is similar to the experiment [8] and the DNS [10], which is consistent with the skin friction deviation near R = 700 that is more similar to the DNS [10] in Fig. 4, compared to the DNS [11]. Different methods to introduce perturbations in the flow, such as the vibrating ribbon [8], wall blowing and suction [10] and stability-analysis-based inlet condition ([11] and here), can lead to slightly different evolution of the flow given the sensitivity of transition to small changes. After the boundary layer transits to turbulent, the mean velocity profile U + (y+ ) is analyzed (see Fig. 7). The current LES accurately predicts the viscous sublayer and the log-law layer. The agreement to the DNS solution is excellent. It should be noted that the von Kármán constant κ is slightly less than the usual value κ ࣃ 0.38 for high-Re wall-bounded flow [25]. The value of κ is reported around 0.36 for Reτ ࣠ 10 0 0 [26]. Statistics of turbulent fluctuations from the current LES agree well with the previous DNS [10] as shown in Fig. 8. The velocity fluctuations are well captured, including the urms peak and the overall profiles of all the three components. The shear component of the Reynolds stress tensor u v also agrees well with the DNS data. Based on the mean velocity in Fig. 7 and the turbulent quantities in Fig. 8, current wall-resolved LES provide DNS-like fidelity in the turbulent boundary layer. The current LES is not significantly sensitive to the temporal and spatial resolution tested here (see Fig. 9). The current time

Fig. 9. Skin friction with varying temporal and spatial resolution in the current LES.

step is determined using the period of the fundamental mode T(2,0) , following the approach used in the previous LES study [15]. Two fine-level time steps, t = TT S /512 and TTS /256, provide almost the identical skin friction as shown in Fig. 9. It should be noted that a large time step such as t = TT S /128 causes the solver unstable, has not been pursued further. The most of the current LES results are reported in this paper with t = TT S /256 on the fine grid unless notified otherwise. The transition predicted on the medium grid is very similar to the case of the fine grid - the deviation from the laminar skin friction is almost identical. Due to the marginal resolution of the medium grid for wall-resolved LES (see discussion in Section 2.4), the turbulent skin friction from the medium grid is slightly lower than the fine grid. The lower skin friction is associated to under-resolved fluctuations in the turbulent boundary layer. Albeit the slight underestimation on the turbulent skin friction, the medium grid provides acceptable prediction in the turbulent transition. Coarser grids are also tested in this study (results are not shown); they are not able to reproduce either the transition location nor the turbulent skin friction. Similar behavior with coarse resolution is also observed in Bodart and Larsson [27]. The current medium grid would be close to the minimum grid resolution required for the current laminar-to-turbulent transition using the second-order code OpenFOAM. The turbulence viscosity is essentially zero in the laminar region on both fine and medium grids (see Fig. 10), providing an adequate environment for small-amplitude instabilities to grow without artificial viscosity. The medium grid also captures well the staggered  structures in the region near R = 700. The  structures on the medium grid seem to break down to small eddies

Fig. 8. Turbulent quantities at R = 837. Lines are from current LES on the fine grid, and symbols are DNS data [10].

M. Kim, J. Lim and S. Kim et al. / Computers and Fluids 189 (2019) 108–117

115

Fig. 12. Skin friction from current LES starting from three different inlet locations on the medium grid.

Fig. 10. Instantaneous Q-criterion (a and b) and the normalized eddy viscosity ν t /ν (c and d) on the plane at y/δ0 = 0.4 with the fine (a and c) and medium (b and d) grids.

earlier, compared to the fine grid, which is consistent to the slightly early transition on the medium grid shown in Fig. 9. 3.2. Variations in LES inlet locations The second main goal of the current study is to reduce the computational domain by locating the LES inlet in the downstream closer to the end stage of the transition. Two additional inlet locations are tested: R = 500 and R = 641. Instabilities forced at each LES inlet are listed Table 1. The component u of the 2D TS and the oblique modes is shown in Fig. 11 for the three locations of the current LES. In Fig. 11, each profile is normalized by the maximum of the inlet instability at R = 400. Note that the mode shape of the 2D TS wave remains almost the same in the range of 400 ≤ R ≤ 641 (see Fig. 11a). The subharmonic mode grows rapidly from R = 500 to R = 641 so that the subharmonic amplitude is comparable to the fundamental amplitude at the downstream location R = 641. LES computations with the variation in the inlet location are conducted on the medium grid (see Fig. 12). The current LES is not sensitive to the inlet location. All the three profiles of the skin friction deviate from the laminar skin friction around R = 690. The

transition length is also very similar to each other, regardless the inlet location. The amplitude growth of four instability modes is shown in Fig. 13. The two LES cases of downstream inlets provide very similar growth of the instabilities to the baseline case of the upstream inlet.  The streamwise length in the smallest domain R = Rex = 641 − 10 0 0 is reduced by 30%, compared to the baseline domain R = 400 − 1000. This domain reduction results in the same 30% reduction in the computational cost due to the uniform grid in x in this study. This reduction corresponds to about 40% reduction in the range of R. It was considered that locating the LES inlet close to the highly nonlinear region R ࣡ 650 is challenging in the DNS and LES study of Lozano-Durán et al. [11]. The current LES approach with inlet disturbances determined from nonlinear PSE is able to move the inlet location close to R ∼ 650 without altering the transition in the boundary layer. The WALE model is practically dormant until enough number of instabilities grow (see Fig. 14). The eddy viscosity is sufficiently small νt /ν < 10−2 before the nonlinear interaction is severe. In the advent of the  structures, the eddy viscosity is νt /ν ∼ O (10−2 ). The value of the eddy viscosity quickly increases as the  structures break down to smaller structures. It is expected that the sub-grid-scale model activates sizable eddy viscosity when the LES grid is not sufficiently small to resolve all the spectrum of the turbulence including the smallest scale (i.e., the Kolmogorov scale). The current computational study provides quality data which can

Fig. 11. Amplitude profiles at three LES inlet locations.

116

M. Kim, J. Lim and S. Kim et al. / Computers and Fluids 189 (2019) 108–117

Fig. 13. Amplitude of four selected instabilities in the current LES on the medium grid.

nonlinear PSE. Current LES computations with three inlet locations provided very similar results. The most downstream location reduced the domain by 30% from the baseline case, and, as a result, the computational cost was also reduced by 30%.

Acknowledgments

Fig. 14. Instantaneous, normalized eddy viscosity ν t /ν on the plane at y/δ0 = 0.4 with inlets located at R = 400 (a), R = 500 (b), and R = 641 (c).

be used for the assessment of the sub-grid-scale model in the transitional boundary layer flow.

This work is carried out as part of Space Core Technology Development Program. This work was supported by the National Research Foundation of Korea (NRF) grant funded by the Ministry of Science and ICT (MSIT) of Korea government (NRF-2017M1A3A3A02016810 and NRF-2017R1C1B5018300). This work is also supported by the National Institute of Supercomputing and Network/Korea Institute of Science and Technology Information with supercomputing resources including technical support (Project No. KSC-2018-C2-0031, KSC-2018-CRE-0113, KSC-2017-C3-0048, and KSC-2018-T1-0 0 05).

4. Conclusions The PSE-coupled LES approach was investigated for laminar-toturbulent transition using the open-source solver OpenFOAM. The canonical, subharmonic resonance on a flat plate was simulated. PSE provided instabilities required at the LES inlets. The nonlinear PSE also provided the evolution of instabilities of interest, validating the LES computation in the pre-transitional boundary layer. The current LES was compared well with the experiment and DNS data, including the skin friction, amplitude growth of instabilities, and turbulence statistics. The implementation of the PSE-coupled LES in OpenFOAM was successfully validated and can be used in other types of laminar-to-turbulent transition. The inlet for current LES was re-located in order to reduce the computational domain. The most challenging case was to locate the inlet significantly close to the final stage of the transition. Instabilities at various inlet locations were obtained from the

Appendix A. LES inlet condition based on PSE in OpenFOAM The LES inlet condition based on PSE analysis is implemented in OpenFOAM-4.1. Two files, U and ISstream.C, are modified. At the LES inlet, the current code requires (1) laminar solution and (2) instabilities. The laminar solution is provided from an independent computation. Instabilities are provided from a separate PSE simulation. The time-varying profiles of instability modes are superposed to the laminar solution at the inlet. The modified file U along with exemplary files are accessible in the web theory.gist.ac.kr. Check the news item 1, PSE-based LES Inlet, under the Community tab. It is required to use the passcode of seven letters LES+PSE (all capital) to unzip the archived file PSE-BASED_LES_INLET.ZIP. Exemplary files for the current LES on the medium grid starting at R = 400 are also included in the archived file.

M. Kim, J. Lim and S. Kim et al. / Computers and Fluids 189 (2019) 108–117

References [1] Mack LM. Boundary-layer linear stability theory. Special course on stability and transition of laminar flow. Advisory Group for Aerospace Research and Development (AGARD); 1984. [2] Reed HL, Saric WS, Arnal D. Linear stability theory applied to boundary layers. Annu Rev Fluid Mech 1996;28(1):389–428. doi:10.1146/annurev.fl.28.010196. 002133. [3] Herbert T. Parabolized stability equations. Annu Rev Fluid Mech 1997;29(1):245–83. doi:10.1146/annurev.fluid.29.1.245. [4] Chang C-L, Malik MR, Erlebacher G, Hussaini MY. Linear and nonlinear PSE for compressible boundary layers. Tech. Rep.. NASA Langley Research Center; 1993. [5] Joslin RD, Streett CL, Chang C-L. Spatial direct numerical simulation of boundary-layer transition mechanisms: validation of PSE theory. Theor Comput Fluid Dyn 1993;4(6):271–88. doi:10.10 07/BF0 0418777. [6] Esfahanian V, Hejranfar K, Sabetghadam F. Linear and nonlinear PSE for stability analysis of the Blasius boundary layer using compact scheme. J Fluids Eng 2001;123(3):545–50. [7] Park D, Park SO. Linear and non-linear stability analysis of incompressible boundary layer over a two-dimensional hump. Comput Fluids 2013;73:80–96. doi:10.1016/j.compfluid.2012.12.007. http://www.sciencedirect. com/science/article/pii/S0 0457930120 04586 [8] Kachanov YS, Levchenko VY. The resonant interaction of disturbances at laminar-turbulent transition in a boundary layer. J Fluid Mech 1984;138:209– 47. doi:10.1017/S0 0221120840 0 010 0. [9] Fasel H, Rist U, Konzelmann U. Numerical investigation of the three-dimensional development in boundary-layer transition. AIAA J 1990;28(1):29–37. [10] Sayadi T, Hamman CW, Moin P. Direct numerical simulation of complete Htype and K-type transitions with implications for the dynamics of turbulent boundary layers. J Fluid Mech 2013;724:480–509. doi:10.1017/jfm.2013.142. [11] Lozano-Durán A, Hack MJP, Moin P. Modeling boundary-layer transition in direct and large-eddy simulations using parabolized stability equations. Phys Rev Fluids 2018;3:023901. doi:10.1103/PhysRevFluids.3.023901. [12] Choi H, Moin P. Grid-point requirements for large eddy simulation: Chapman’s estimates revisited. Phys Fluids 2012;24(1):011702. doi:10.1063/1.3676783. [13] Huai X, Joslin R, Piomelli U. Large-eddy simulation of transition to turbulence in boundary layers. Theor Comput Fluid Dyn 1997;9(2):149–63. doi:10.1007/ s0 01620 050 037.

117

[14] Sayadi T, Moin P. Large eddy simulation of controlled transition to turbulence. Phys Fluids 2012;24(11):114103. doi:10.1063/1.4767537. [15] Jee S, Joo J, Lin R-S. Toward cost-effective boundary layer transition computations with large-eddy simulation. J Fluids Eng 2018;140(11):111201–111201–12. doi:10.1115/1.4039865. [16] Borodulin VI, Kachanov YS, Koptsev DB. Experimental study of resonant interactions of instability waves in a self-similar boundary layer with an adverse pressure gradient: i. Tuned resonances. J Turbul 2002;3:N62. doi:10. 1088/1468-5248/3/1/062. [17] Park D, Park SO. Study of effect of a smooth hump on hypersonic boundary layer instability. Theor Comput Fluid Dyn 2016;30(6):543–63. doi:10.10 07/s0 0162-016- 0396- 7. https://link.springer.com/article/10.1007% 2Fs00162-016- 0396- 7 [18] Nicoud F, Ducros F. Subgrid-scale stress modelling based on the square of the velocity gradient tensor. Flow Turbul Combust 1999;62(3):183–200. doi:10. 1023/A:10 099954260 01. [19] Qin S, Koochesfahani M, Jaberi F. Large eddy simulations of unsteady flows over a stationary airfoil. Comput Fluids 2018;161:155–70. doi:10.1016/ j.compfluid.2017.11.014. [20] Issa RI. Solution of the implicitly discretised fluid flow equations by operatorsplitting. J Comput Phys 1986;62(1):40–65. doi:10.1016/0021-9991(86) 90099-9. [21] Hunt JCR, Wray AA, Moin P. Eddies, streams, and convergence zones in turbulent flows. In: Summer program of center for turbulent research at Stanford University. Stanford, CA; 1988. p. 193–208. [22] Herbert T. Secondary instability of boundary layers. Annu Rev Fluid Mech 1988;20(1):487–526. doi:10.1146/annurev.fl.20.010188.002415. [23] Kachanov YS. Physical mechanisms of laminar-boundary-layer transition. Annu Rev Fluid Mech 1994;26(1):411–82. doi:10.1146/annurev.fl.26.010194.002211. [24] White FM. Viscous fluid flow. 3rd ed. McGraw-Hill; 2006. [25] Lee M, Moser RD. Direct numerical simulation of turbulent channel flow up to Reτ ≈ 5200. J Fluid Mech 2015;774:395–415. doi:10.1017/jfm.2015.268. [26] Mizuno Y, Jiménez J. Mean velocity and length-scales in the overlap region of wall-bounded turbulent flows. Phys Fluids 2011;23(8):085112. doi:10.1063/ 1.3626406. [27] Bodart J, Larsson J. Sensor-based computation of transitional flows using wall-modeled large eddy simulation. In: Center for turbulence research - Annual research briefs. Stanford University; 2012. p. 229–40.