A dynamic closure modeling framework for large eddy simulation using approximate deconvolution: Burgers equation

Abstract We put forth a dynamic closure modeling framework for the large eddy simulations of the Burgers equation based upon the use of the approximate deconvolution (AD) procedure to compute the Smagorinsky constant self-adaptively from the resolved flow quantities. In our proposed framework, the test filtering process of the standard dynamic model is replaced by the AD procedure. The robustness of the model has been tested considering the Burgers equation in its conservative and skew-symmetric forms. Our numerical assessments for solving the single-mode sine wave and the decaying Burgers turbulence problems show that the present framework effectively damps grid-to-grid oscillations and yields an improved shock capturing property for central numerical schemes as underlying discretizations.


Introduction
Turbulent flows are encountered in a variety of engineering and geophysical systems involving a wide range of spatial and temporal scales. In a direct numerical simulation (DNS), the full spectra of

PUBLIC INTEREST STATEMENT
Multiple physical processes in the engineering and environmental sciences are tightly coupled with turbulent flows where the accuracy and generalizability of closure models play a key role. Large eddy simulation (LES) is a popular technique for simulating turbulent flows. In practice, direct numerical simulation (DNS), which resolves every scale of the solution, is prohibitively expensive for nearly all systems with complex geometry or flow configurations. To achieve the same order of resolved physical accuracy as DNS, however, LES requires to correctly treat the well-known closure problem: the effect of the small scales on the large ones needs to be modeled. Therefore, there has been a substantial effort in developing these LES closure models in the past few decades. The present study investigates the feasibility of a dynamic eddy viscosity model for solving the decaying Burgers turbulence problem with quadratic nonlinearity, which can be considered a simplified prototype integrable system for more complex flow dynamics.
turbulence should be resolved down to the Kolmogorov scale where the smallest feature of flow motion is captured. The resolution requirements of small scale turbulence, however, are computationally prohibitive to fully resolve for all associated scales. Large eddy simulation (LES) aims to reduce this computational complexity and has been proven to be a promising approach for calculations of complex turbulent flows (Boris, Grinstein, Oran, & Kolbe, 1992;Galperin & Orszag, 1993;Lesieur & Metais, 1996;Piomelli, 1999;Meneveau & Katz, 2000). LES allows for much coarser spatial meshes as it is designed to resolve the most energetic large scales of the turbulent motion while modeling small scales.
The LES equations are derived formally by applying a low-pass spatial filter to the governing equations of turbulent motion and can be considered a weak form or regularized form of the governing equations to remove the resolution requirement of small-scale turbulence. The nonlinear nature of the governing equations leads to the well-known closure problem in LES: the requirement of subgrid-scale (SGS) parameterization where interactions between the small scales and the large ones need to be modeled. In the past few decades, there has been a substantial effort on developing LES closure models using physical or mathematical arguments (Berselli, Iliescu, & Layton, 2006;Sagaut, 2006).
The eddy viscosity (EV) approach is the foundation of many functional turbulent closure models. This approach yields one of the most celebrated closure models -the Smagorinsky model (1963)which assumes that the eddy viscosity is computed from the resolved strain rate magnitude and a characteristic mixing length scale which is assumed to be proportional to the filter width via a modeling parameter usually called the Smagorinsky constant. Applications of the Smagorinsky model to various problems have revealed that the constant is not single-valued and varies depending on resolution and flow characteristics (Galperin & Orszag, 1993). A major advance took place with the development of a dynamic model proposed by Germano, Piomelli, Moin, Cabot (1991) and its modified form by Lilly (1992) in which the Smagorinsky constant was self-adaptively determined from the simulation using a low-pass spatial test filter. This has made the eddy viscosity LES closure models more widely applicable in many fields (e.g. see Piomelli (1999) and Meneveau and Katz (2000)). Layton (2016) recently analyzed the dissipative behavior of the Smagorinsky model and addressed the dynamic parameter selection as an important open problem, which is the topic of focus in the present study.
We emphasize that the iterative AD process can successfully deconvolve the subfilter-scale (SFS) interactions on a given numerical grid. Additional penalty terms (Stolz, Adams, & Kleiser, 2001) or relaxation filtering techniques (Mathew, Lechner, Foysi, Sesterhenn, & Friedrich, 2003;Mathew et al. 2006) are used to account for the effective SGS interactions (i.e. scales beyond the grid cut-off scale where they cannot be deconvolved from the filtered quantities on resolved flow and must be modeled). A family of selective filters has been studied to eliminate grid-to-grid oscillations (Berland, Lafon, Daude, Crouzet, Bogey, & Bailly, 2011;Bogey & Bailly, 2004;Bogey & Bailly, 2006a;Bogey & Bailly, 2006b). We also note that the LES closure term (i.e. a representation of total residual interactions) consists of the SFS and effective SGS interactions Gullbrand and Chow (2003), where the summation of these interactions is traditionally called the SGS stress term in LES literature. We highlight that the SFS and SGS stresses are different quantities and the AD procedure recovers an approximation of the truncated LES field, not of the DNS. Therefore, an additional penalty or stabilizer term should be used to model subgrid scales in structural modeling AD frameworks.
In this note, we put forth a dynamic approach for functional closure modeling of Burgers equation based upon using the AD process to estimate the eddy viscosity coefficient in a self-adaptive manner during simulations. The approach presented here is fundamentally different from the mixed methods where the AD model is coupled with the eddy viscosity models by adding a dynamically computed source term to increase the stability of the AD model (Gullbrand & Chow 2003;Habisreutinger, Bouffanais, Leriche, & Deville 2003), and it is also different from the dynamic mixed scale-similarity models Sarghini, Piomelli, & Balaras, 1999;Zang, Street, & Koseff, 1993). In this framework, the test filtering process of the standard dynamic model is replaced by the AD procedure, and it is shown that the order of accuracy increases by increasing the number of the Van Cittert iterations (i.e. damping grid-to-grid oscillations more effectively). It could be also interpreted that the deconvolved SFS interactions are utilized to compute total eddy viscosity interactions minimizing the error between SFS and SGS interactions. This procedure provides an additional stabilization and yields a good approximation when resolved SFS interactions are pronounced. Since the resolved SFS contributions diminish in pseudospectral simulations, our framework will be limited to finite domain discretization methodologies (i.e. finite differences). In our previous studies, the performance of the proposed methodology has been assessed for solving incompressible flow problems without any shocks and discontinuities (Maulik & San, 2017a, 2017b, 2017c. This paper focuses on the numerical assessments of dissipation mechanism for the stationary and moving shock problems. Indeed, our assessments for solving the Burgers equation numerically prove the concept of using the AD procedure in a dynamic eddy viscosity model as a proposed alternative to the classical test filtering approach. Although the Burgers turbulence problem is a limiting case, it has been used to perform assessments for turbulence closure models (Adams & Stolz, 2002;De Stefano & Vasilyev, 2002;Love, 1980), since it retains some important properties of the Navier-Stokes equations (Falkovich & Sreenivasan, 2006).

Modeling framework
In this section, we present the proposed dynamic model for the Burgers equation, a simple prototype for fluid flows having a similar quadratic nonlinearity. The viscous Burgers equation in the conservative form reads as: where u is the velocity, and is the kinematic viscosity. To clarify the presentation in our analysis, we rewrite the governing equation as where we note that its extension to the Navier-Stokes equations is straightforward. Before moving on to the derivation of the proposed framework, it would be appropriate to define the formal LES framework. In computational studies the governing equations are solved numerically on a grid. We consider LES on a grid of size h with the truncation of the DNS field u on this coarse grid denoted by ũ. For LES computations, the governing equations are filtered in space and solved numerically on a grid, hence the filtered equations of motion can be obtained by performing a low-pass filtering in the following form where ū denotes the filtered velocity on an LES grid. The effective low-pass filtered LES equation can be rewritten as where g is the SGS stress term to account for the effects of the truncated small scales on the resolved scales, and it can be formally written as The nonlinear term can be explicitly filtered in order to avoid aliasing (Carati, Winckelmans, & Jeanmart, 2001;Fauconnier, De Langhe, & Dick, 2009;Gullbrand & Chow, 2003;Orszag, 1971) where In this note, we follow Equation (4) since it only alters the original governing equation through the definition of a source term g. In the present notation, Equation (5) can be further decomposed into where the first term represents the "resolved" SFS interactions and the second term is for "effective" SGS term. The SFS term can be partially deconvolved using the approximate deconvolution (AD) procedure; however, the second term must be modeled. If ũ Q is assumed to be an approximately recovered value for the unfiltered grid variable ũ, the Van Cittert iterations can be used to estimate ũ Q from the filtered grid variable ū : where G is the low-pass spatial filtering operator (i.e. G * ũ =ū), and this approach yields increasingly accurate reconstructions with increasing Q, the order of the Van Cittert iteration (Dunca & Epshteyn, 2006). Here, is the over-relaxation parameter, and it can be shown that 0 < ≤ 2 when the transfer function of the spatial filter 0 ≤ T(k) ≤ 1 (i.e. see Section 2.4 for further details). The interested reader is referred to Biemond et al. (1990) for convergence characteristics of the Van Cittert iterations and AD procedure. We have set = 2 in our computations. In an AD framework, Equation (8) can be rewritten as where we denote ũ ≈ũ Q as Qth order approximation for the recovered grid value (i.e. partially deconvolved velocity field on an LES grid obtained using the iterative inverse filtering AD procedure). (4) x 2 +g,

Eddy viscosity modeling for the Burgers equation
Assuming the dissipation of kinetic energy at subgrid scales to be analogous to molecular diffusion, Smagorinsky's eddy viscosity model approximates the total SGS term as where e is the eddy viscosity which is given by in which C S is the Smagorinsky constant and ̄ is formally defined by the filter width and usually set to the representative LES mesh size (e.g. ̄= h). We note that ū denotes the filtered velocity field on the coarse numerical grid (i.e. LES numerical resolution) which is readily available from the LES computations.

Dynamic eddy viscosity modeling for the Burgers equation
Application of the Smagorinsky model to various flow problems has revealed that the C S constant is not single-valued Smagorinsky (1993) and needs for an ad-hoc specification depending on the flow characteristics. Following the classical dynamic model procedure (Germano, Piomelli, Moin, & Cabot, 1991;Lilly, 1992), the modeling constant (C S̄) 2 can be computed self-adaptively by defining a test filter. For the dynamic implementation, using Equations (4) and (5) 1 2 Factoring the model coefficient (i.e. a practical frozen assumption), Equation (17) can be recast in the following form where in which =̂∕̄ is the filter ratio (e.g. = 2 in the present study). The coefficient (C S̄) 2 in Equation (18) can be obtained by minimizing the average square error ⟨E 2 ⟩, where the error is given by In the present study, the averaging operator is expressed by for any quantity f. On differentiating the mean square error with respect to the model parameter , we get where we assume that the Smagorinsky constant can be factored out of the averaging operators (Lilly, 1992). Then, we can minimize the square of the error when where H and M are given by Equations (19) and (20), respectively. This completes the dynamic Smagorinsky (D-SMA) model given by Equations (4), (11) and (12).

Derivation of the proposed dynamic AD model for the Burgers equation
The combination of the structural and functional closure models has been successfully used in the LES literature (e.g. see mixed dynamic models Sarghini, Piomelli, & Balaras, 1999;Zang, Street, & Koseff, 1993). A coupled AD and eddy viscosity model has been also proposed by Habisreutinger et al. (2007) for large eddy simulation of incompressible flows. Here, we purpose a different approach to determine the modeling constant based upon using the AD procedure. In the present notation, the SGS field is modeled functionally using the eddy viscosity hypothesis where we use Equations (11) and (12). It can be also approximated structurally using the AD procedure where we take into account for only the SFS interactions (i.e. we ignore the second term in the righthand side of Equation (8)). What is new in the proposed model here is the way of computing (C S̄) 2 using the following approximation between Equations (24) and (25) where we use the SFS portion of Equation (8) in order to approximate the free modeling parameter of the eddy viscosity model dynamically. Consequently, Equation (26) can be recast in the following form where Similar to the procedure prescribed in Section 2.2 , the coefficient (C S̄) 2 in Equation (27)

Low-pass filtering operator
The definition of the spatial filter G completes the formulation of the proposed dynamic model. The filter shape and associated filtering parameters are free parameters in LES computations and various forms of the low-pass filters have been utilized in the LES literature (e.g. see Germano (1986), Jordan and Ragab (1996), Najjar and Tafti (1996), Vasilyev et al. (1998), Sagaut and Grohens (1999), Mullen and Fischer (1999), Pruett and Adams (2000), Najafi-Yazdi et al. (2015)). Here, we consider the following -parameter sixth-order compact Padé filter (Lele, 1992): where f j represents the filtered value of a discrete quantity f j . Here, the filtering coefficients are given by The free parameter, , which is in the range | | < 0.5, determines the filtering properties, with high values of yielding less dissipative results. It is easy to verify that the value of = 0.5 yields back the unfiltered quantity of interest. A modified wavenumber analysis yields a transfer function, T(k), which correlates the Fourier coefficients of the filtered variable to those of the unfiltered variable in the Fourier domain, as follows (San, 2016): where k is the wavenumber and h is the grid spacing. The transfer function of the Padé filter for different values of are plotted in Figure 1, showing that the Padé filter remains close to unity over a wide range of wavenumbers, resulting in minimal dissipation over resolved inertial scales. Also, another characteristic of the compact Padé filter is that the maximum attenuation happens for the largest wavenumber (i.e. T(k) goes zero when k becomes the grid cut-off wavenumber). As a special case in Equation (30), where the free parameter = 0, we obtain the following discrete filter: We also note that Equation (33) is also called sixth-order binomial smoothing filter and its transfer function can be written as Jahne (1997) Unless otherwise stated, we use Equation (33) in our numerical assessments.

Numerical schemes
In this study, we use the method of lines procedure by utilizing an explicit Runge-Kutta scheme for time integration. Specifically, a sixth-order compact central difference and a third-order optimal Runge-Kutta schemes are used for the spatial and temporal discretizations, respectively. In the compact central difference scheme, the first derivative of any function f on a collocated finite difference grid can be computed as follows (Lele, 1992): where h is the uniform grid spacing and superscript prime denotes the first derivative (i.e. f � = f x ). Here, Equation (35) provides a sixth-order tridiagonal scheme which can be solved efficiently using the Thomas algorithm (1992). Similarly, the compact sixth-order centered scheme for the second derivative is given by where superscript double prime denotes the second derivative (i.e. f �� = 2 f x 2 ).
(32) T(k) = a 0 + a 1 cos(hk) + a 2 cos(2hk) + a 3 cos (3hk) 1 + 2 cos(hk) , In the method of lines, a system of semi-discrete ordinary differential equations is obtained after the spatial discretization of the governing equations. For example, the discrete version of Equation (4) can be written as where £(ū j ) is the discrete operator for the nonlinear, linear, and SGS terms. We assume that the numerical approximation for time level n is known, and we seek the numerical approximation for time level n + 1, after the time step Δt. The optimal third-order accurate total variation diminishing Runge-Kutta scheme is then given as follows (Gottlieb & Shu, 1998):

Results
In this section, we present our results for two test cases: (i) the shock formation problem initiated by a single-mode sine wave, u(x, 0) = sin(x), and (ii) the decaying Burgers turbulence problem initiated by a specified energy spectrum considering the 64 sample initial conditions associated with randomly generated phases. Further details on the initial energy spectrum can be found in San (2016).
Both problems are considered in a domain of x ∈ [0, 2 ] with periodic boundary conditions using = 5 × 10 −4 . Time integration for our discrete system is carried out using a third-order accurate total variation diminishing Runge-Kutta scheme with Δt = 1 × 10 −5 , and the sixth-order central compact difference schemes are used for the spatial discretizations to ensure a minimal numerical error.  Figure  2(d-f), the proposed dynamic AD (D-AD) model yields substantially better results with the damping of grid-to-grid oscillations by increasing Q, the order of Van Cittert iteration, which can also be seen from Figure 3 demonstrating the total dissipation rate −dE∕dt where E = ∫ E(k)dk is the total kinetic energy. Figure 3 also demonstrates that the peak value of the dissipation rate agrees well with the DNS data.

Single-mode sine wave simulations
The nonlinear term given in Equation (1) can be also written in the skew-symmetric form, Ritchmyer and Morton (1967), Jameson (2008), leading to It has been well known that the discrete version of these formulations yields slightly different results.
To analyze the effect of the formulation, we perform the same analysis using the skew-symmetric form for the nonlinear term. Having parameters identical to those used in Figures 2 and 3 for the 1 2 conservative form, Figures 4 and 5 illustrate the results obtained by the skew-symmetric form. Further comparisons are shown in Figure 6 for energy spectra and in Figure 7 for velocity fields. Although similar conclusions can be made from these observations regarding the closure modeling performance of the proposed framework, we can easily see that the skew-symmetric formulation yields slightly better results due to its discrete conservation properties. A first conclusion from these  results suggests that the proposed D-AD model is more dissipative than the D-SMA model and reduces grid-to-grid oscillations even for considerably coarse grids. It is also evident on examinations of Figure 7 that the shock capturing property of the proposed D-AD model is superior to the D-SMA model.

Decaying Burgers turbulence simulations
Next, we present our results for the Burgers turbulence problem in a periodic domain of x ∈ [0, 2 ].
We solve 64 different cases generated from the following initial energy spectrum  where k 0 = 10. Figure 8 illustrates the solution fields for three of them obtained using DNS with a resolution of N = 32768. We note that each simulation is obtained using the same initial energy spectrum with randomly assigned phase distributions (i.e. in order to generate different initial conditions in physical space). Figure 9(a) illustrates the time evolution of the decaying Burgers turbulence based upon the ensemble average of 64 sample fields, obtained by DNS with a resolution of N = 32768. Figure 9(b) shows the energy spectra at t = 0.05 for DNS and UNS computations for different resolutions demonstrating the ideal k −2 scaling at inertial range. It is clear from figure that the resolution of N = 32768 yields DNS data resolving all associated shock scales at = 5 × 10 −4 . One must note that t = 0.05 corresponds approximately to the time for the maximum dissipation rate of the problem. The energy pile-up close to grid cut-off scale is also evident at lower resolutions.  Note: The conservative formulation is used for discretization of the nonlinear term. In Figure 10, we present our LES analysis using the coarsest resolution of N = 512 illustrating the time evolution of the total dissipation rate and energy spectra at time t = 0.05. Results obtained by DNS, UNS, and the dynamic Smagorinsky (D-SMA) model are also included for the purpose of comparison. It is shown that the proposed dynamic model yields a dissipation rate with considerably higher accuracy and the energy pile-up can be effectively eliminated by increasing Q. It is also evident on careful examination of Figure 10 that the overall accuracy of the proposed scheme is fairly good over the entire simulation time. However, the UNS and the dynamic Smagorinsky (D-SMA) models cannot capture the physics beyond t = 0.05 due to accumulation of grid-to-grid oscillations.
In a manner similar to the previous test case, increasing Q yields an improved approximation to the representation of the SFS interactions. This implies that more dissipative behavior can be obtained using a higher order Van Cittert iterative procedure, which would be crucial for solving LES equations on very coarse grid resolutions. We also found that results are negligibly improved by increasing the order of Van Cittert iterations beyond Q = 5.
To analyze the discrete formulation effect, we perform the same analysis using the skew-symmetric form for the nonlinear term. Utilizing parameters and initial conditions identical to those used for Figures 9-12 illustrate the results obtained by the skew-symmetric form. It is noted that the same resolving performance is obtained for the finest resolution, which can be considered as a verification of its DNS characteristics at N = 32768. As discussed earlier, we can easily see that the skew-symmetric formulation yields slightly better results due to its discrete conservation properties. In the Note: The skew-symmetric formulation is used for discretization of the nonlinear term. skew-symmetric formulation, all coarse grid simulations yield stable numerical solutions for all time evolutions for the decaying Burgers turbulence cases. However, it is clear that the D-AD model predicts the peak dissipation rate better.
It is also shown that the proposed model yields more dissipative results. It follows a similar trend that we observed in the case of the single-mode sine wave experiments in which we obtain more dissipative results for higher order Q. This can be explained in such a way that an increase in the order Q yields an increase of the estimated C S value. An examination of energy spectra plots given in Figure 12 reveals that the D-SMA yields a longer inertial range compared to the D-AD. This can be explained by the use of relatively dissipative behavior of the low-pass filter given by Equation (33). Therefore, a sensitivity analysis with respect to the free filtering parameter is also presented in Figure 13. As shown here less dissipative results can be obtained using a less dissipative filter (i.e. increasing ). Here, we note that larger inertial range can be obtained by increasing value in the D-AD model. It can also be seen that UNS data can be recovered when → 0.5. We emphasize that a-posteriori results of D-SMA and D-AD models depend on the free filtering parameters, underlying numerical scheme, convective term formulation, and computational resolution. Using the same filtering kernel and the same discrete method, our numerical assessments show that the D-SMA performs better for capturing the larger inertial range since the D-AD model acts more dissipative. On the other hand, the D-AD performs better to eliminate grid-to-grid oscillations and provide more stable results for coarser grids due to the its more dissipative behavior. Furthermore, the level of dissipation can be controlled by the free filtering parameter.

Conclusion
We propose a simplified dynamic framework for SGS modeling of LES of turbulent flows based upon the use of AD process to compute the Smagorinsky constant self-adaptively from the resolved flow quantities. Our proposed model replaces the test filtering procedure of the standard dynamic model with the AD procedure. Although the presentation of the model is given using the Burgers equation, its generalization to the Navier-Stokes equations is straightforward. Our numerical assessments for solving the single-mode sine wave and the Burgers turbulence problems show that the proposed approach could represent a viable tool for subgrid scale parametrization of turbulent flows by combining the AD process with an eddy viscosity model. Since we consider only the SFS portion of the residual interactions when computing the free modeling parameter, the AD procedure provides for a more dissipative eddy viscosity model than the classical test filtering approach and results in a more stable behavior reducing the grid-to-grid oscillations in coarse grid central scheme-based LES computations.