Variational Multiscale Closures for Finite Element Discretizations Using the Mori-Zwanzig Approach

TL;DR

This paper introduces a novel coarse-grained modeling approach combining Variational Multiscale decomposition with the Mori-Zwanzig formalism to improve finite element simulations of multiscale problems, adaptively determining memory effects.

Contribution

It develops a parameter-free, adaptive sub-scale model for multiscale finite element discretizations using Mori-Zwanzig formalism, applicable to turbulence and Burgers equation.

Findings

Model effectively captures unresolved scale effects.

Adaptive memory length improves simulation accuracy.

Applicable to turbulence and nonlinear problems.

Abstract

Simulation of multiscale problems remains a challenge due to the disparate range of spatial and temporal scales and the complex interaction between the resolved and unresolved scales. This work develops a coarse-grained modeling approach for the Continuous Galerkin discretizations by combining the Variational Multiscale decomposition and the Mori-Zwanzig (M-Z) formalism. An appeal of the M-Z formalism is that - akin to Greens functions for linear problems - the impact of unresolved dynamics on resolved scales can be formally represented as a convolution (or memory) integral in a non-linear setting. To ensure tractable and efficient models, Markovian closures are developed for the M-Z memory integral. The resulting sub-scale model has some similarities to adjoint stabilization and orthogonal subscale models. The model is made parameter free by adaptively determining the memory length during the simulation. To illustrate the generalizablity of this model, it is employed in coarse-grained simulations for the one-dimensional Burgers equation and in incompressible turbulence problems.

Tables

Table 1. Table 1 : Simulation parameters for DNS and LES of the Burgers Equation for an initial sine profile.

Case	Domain Size $L$	Degrees of Freedom $N$	Grid Size $dx$	Time Step $dt$	Viscosity $\nu$	Memory Length $\tau$
DNS (Spectral)	$2\pi$	4096 modes	$7.67\times {10}^{-4}$	$3.83\times {10}^{-4}$	${10}^{-4}$	-
Dynamic $\tau$	$2\pi$	32 elements	$1.96\times {10}^{-1}$	$1.96\times {10}^{-2}$	${10}^{-4}$	Dynamic
FM $\tau$=0.01	$2\pi$	32 elements	$1.96\times {10}^{-1}$	$1.96\times {10}^{-2}$	${10}^{-4}$	0.01
FM $\tau$=0.11	$2\pi$	32 elements	$1.96\times {10}^{-1}$	$1.96\times {10}^{-2}$	${10}^{-4}$	0.11
FM $\tau$=0.23	$2\pi$	32 elements	$1.96\times {10}^{-1}$	$1.96\times {10}^{-2}$	${10}^{-4}$	0.23
OSS [15]	$2\pi$	32 elements	$1.96\times {10}^{-1}$	$1.96\times {10}^{-2}$	${10}^{-4}$	-
No-Model	$2\pi$	32 elements	$1.96\times {10}^{-1}$	$1.96\times {10}^{-2}$	${10}^{-4}$	0

Table 2. Table 2 : Simulation parameters for DNS and LES of the Burgers Equation for the Burger turbulence problem.

Case	$L$	$N$	$dx$	$dt$	$\nu$	$U^{\ast }$	$K_c$	FM $\tau$’s
Case A (DNS)	$2\pi$	1024 modes	$3.06\times {10}^{-3}$	$2.33\times {10}^{-4}$	${10}^{-2}$	1	8	-
Case A (LES)	$2\pi$	32 elements	$1.96\times {10}^{-1}$	$8.5\times {10}^{-3}$	${10}^{-2}$	1	8	0.01, 0.1 and 0.4
Case B (DNS)	$2\pi$	4096 modes	$7.66\times {10}^{-4}$	$3.41\times {10}^{-6}$	$5\times {10}^{-4}$	10	32	-
Case B (LES)	$2\pi$	64 elements	$9.81\times {10}^{-2}$	$4.67\times {10}^{-4}$	$5\times {10}^{-4}$	10	32	0.0001, 0.001 and 0.01

Table 3. Table 3 : Simulation parameters for DNS and LES of the Taylor Green Vortex problem.

Case	$DOFs$	$dx$	$dt$	$\nu$	$U_0$	$L$	FM $\tau$’s
DNS-400	${64}^3$	$9.81\times {10}^{-2}$	$2\times {10}^{-2}$	$2.5\times {10}^{-3}$	1	1	-
LES-FM-400	${32}^3$	$1.96\times {10}^{-1}$	$1.96\times {10}^{-2}$	$2.5\times {10}^{-3}$	1	1	0.01 and 0.002
LES-DY-400	${32}^3$	$1.96\times {10}^{-1}$	$1.96\times {10}^{-2}$	$2.5\times {10}^{-3}$	1	1	Dynamic
DNS-800	${128}^3$	$4.90\times {10}^{-2}$	$2\times {10}^{-2}$	$1.25\times {10}^{-3}$	1	1	-
LES-FM-800	${32}^3,{48}^3$	$1.96\times {10}^{-1}$, $1.31\times {10}^{-1}$	$1.96\times {10}^{-2}$	$1.25\times {10}^{-3}$	1	1	0.01 and 0.002
LES-DY-800	${32}^3,{48}^3$	$1.96\times {10}^{-1}$, $1.31\times {10}^{-1}$	$1.96\times {10}^{-2}$	$1.25\times {10}^{-3}$	1	1	Dynamic
DNS-1600	${256}^3$	$2.45\times {10}^{-2}$	$5\times {10}^{-3}$	$6.25\times {10}^{-4}$	1	1	-
LES-FM-1600	${32}^3,{64}^3$	$1.96\times {10}^{-1}$, $9.81\times {10}^{-2}$	$1.96\times {10}^{-2}$, ${10}^{-2}$	$6.25\times {10}^{-4}$	1	1	0.01 and 0.002
LES-DY-1600	${32}^3,{64}^3$	$1.96\times {10}^{-1}$, $9.81\times {10}^{-2}$	${10}^{-2}$	$6.25\times {10}^{-4}$	1	1	Dynamic

Equations

$$\frac{dx}{dt}=A_{11}x+A_{12}y$$

$$\frac{dy}{dt}=A_{21}x+A_{22}y,$$

$$\frac{dx}{dt}=A_{11}x+F(x).$$

$$\frac{dx}{dt}=A_{11}x+A_{12}{\int }_0^t{e}^{A_{22}s}{A}_{21}x(t-s)ds+A_{12}{e}^{A_{22}t}y(0).$$

$$\frac{d\varphi }{dt}=R(\varphi ),$$

$$\frac{\partial }{\partial t}u({\varphi }_0,t)=\mathcal{L}u({\varphi }_0,t),$$

$$\mathcal{L}=\sum _{k=1}^N{R}_k({\varphi }_0)\frac{\partial }{\partial {\varphi }_{0k}},$$

$$u({\varphi }_0,t)=g(\varphi ({\varphi }_0,t)).$$

$$u({\varphi }_0,t)=e^{t\mathcal{L}}g(\varphi ({\varphi }_0,0))=g(\varphi ({\varphi }_0,t)),$$

$$\frac{\partial }{\partial t}{e}^{t\mathcal{L}}{\varphi }_{0j}=\mathcal{L}{e}^{t\mathcal{L}}{\varphi }_{0j}=e^{t\mathcal{L}}\mathcal{L}{\varphi }_{0j},$$

$$\frac{\partial }{\partial t}{e}^{t\mathcal{L}}{\varphi }_{0j}=e^{t\mathcal{L}}\mathcal{P}\mathcal{L}{\varphi }_{0j}+e^{t\mathcal{L}}\mathcal{Q}\mathcal{L}{\varphi }_{0j},$$

$$e^{t\mathcal{L}}=e^{tQ\mathcal{L}}+{\int }_0^t{e}^{(t-s)\mathcal{L}}\mathcal{P}\mathcal{L}{e}^{s\mathcal{Q}\mathcal{L}}ds,$$

$$\frac{\partial }{\partial t}{e}^{t\mathcal{L}}{\varphi }_{0j}=e^{t\mathcal{L}}\mathcal{P}\mathcal{L}{\varphi }_{0j}+e^{t\mathcal{Q}\mathcal{L}}\mathcal{Q}\mathcal{L}{\varphi }_{0j}+{\int }_0^t{e}^{(t-s)\mathcal{L}}\mathcal{P}\mathcal{L}{e}^{s\mathcal{Q}\mathcal{L}}\mathcal{Q}\mathcal{L}{\varphi }_{0j}ds.$$

$$\frac{\partial }{\partial t}{\mathcal{F}}_j({\varphi }_0,t)=\mathcal{Q}\mathcal{L}{\mathcal{F}}_j({\varphi }_0,t).$$

$$\frac{\partial }{\partial t}{e}^{t\mathcal{L}}{\hat{\varphi }}_{0j}=e^{t\mathcal{L}}\mathcal{P}\mathcal{L}{\hat{\varphi }}_{0j}+{\int }_0^t{e}^{(t-s)\mathcal{L}}\mathcal{P}\mathcal{L}{e}^{s\mathcal{Q}\mathcal{L}}\mathcal{Q}\mathcal{L}{\hat{\varphi }}_{0j}ds.$$

$${\int }_0^t{e}^{(t-s)\mathcal{L}}\mathcal{P}\mathcal{L}{e}^{s\mathcal{Q}\mathcal{L}}\mathcal{Q}\mathcal{L}{\hat{\varphi }}_{0j}ds\approx \tau e^{t\mathcal{L}}\mathcal{P}\mathcal{L}\mathcal{Q}\mathcal{L}{\hat{\varphi }}_{0j}.$$

$$\frac{\partial u}{\partial t}+\mathcal{R}(u)-f=0,$$

$$(\frac{\partial u}{\partial t},w)+(\mathcal{R}(u),w)=(f,w),$$

$$\mathcal{V}=\tilde{\mathcal{V}}\oplus {\mathcal{V}}^{\prime },$$

$$u=\tilde{u}+u^{\prime },$$

$$w=\tilde{w}+w^{\prime },$$

$$(\frac{\partial (\tilde{u}+u^{\prime })}{\partial t},\tilde{w}+w^{\prime })+(\mathcal{R}(\tilde{u}+u^{\prime }),\tilde{w}+w^{\prime })=(f,\tilde{w}+w^{\prime }).$$

$$(\frac{\partial (\tilde{u}+u^{\prime })}{\partial t},\tilde{w})+(\mathcal{R}(\tilde{u}+u^{\prime }),\tilde{w})=(f,\tilde{w}),$$

$$(\frac{\partial (\tilde{u}+u^{\prime })}{\partial t},w^{\prime })+(\mathcal{R}(\tilde{u}+u^{\prime }),w^{\prime })=(f,w^{\prime }).$$

$$(\frac{\partial \tilde{u}}{\partial t},\tilde{w})+(\mathcal{R}(\tilde{u}),\tilde{w})-(f,\tilde{w})=-(\frac{\partial u^{\prime }}{\partial t},\tilde{w})-(\mathcal{R}(\tilde{u}+u^{\prime })-\mathcal{R}(\tilde{u}),\tilde{w}).$$

$$u^{\prime }(y)=-{\int }_{\Omega }{g}^{\prime }(x,y)(\mathcal{R}(\tilde{u})-f)(x)d{\Omega }_x.$$

$$u^{\prime }=-\tau (\mathcal{R}(\tilde{u})-f).$$

$$\frac{\partial u}{\partial t}+\mathcal{R}(u)-f=0,$$

$$\bigg{(}\frac{\partial u}{\partial t},w\bigg{)}_{\Omega}+(R(u),w)_{\Omega}+(b(u),w)_{\Gamma}=(f,w)_{\Omega}\quad\forall{w}\in{\mathcal{V}},$$

$${\Omega }^{\prime }=\bigcup _{i=1}^{n_{el}}{\Omega }^e,$$

Full text

††journal: Computer Methods in Applied Mechanics and Engineering

1 Introduction

Numerical simulation of multi-scale phenomena requires the development of coarse-grained models, which attempts to resolve a sub-set of the scales while providing a model for unresolved features. As an example, a popular approach employed in the simulation of turbulent flows is large-eddy simulation [1, 2, 3, 4, 5, 6], which filters the flow field to resolve the largest energy containing scales and providing a model for the scales smaller than the filter length referred to as the sub-grid scales (SGS). The success of these sub-scale models largely depends on the validity of assumptions at simulated flow conditions. For example, the Smagorinsky model [1], which is one of the most commonly employed SGS models, is based on the assumption that modeled rate for turbulence kinetic energy transfer from large to small scale balances dissipation [7]. This assumption is clearly not valid for all turbulent flows. An alternate approach to SGS modeling without employing phenomenological assumption, which we pursue, is to derive sub-grid models directly from the structure of the PDE and the numerical discretization.

In addition to physics-based sub-grid scale models [1, 3, 2, 8] presented above, models based on the variational multiscale method have proven to be quite successful for both linear [9, 10, 11, 12, 13] and non-linear problems [14, 15, 16, 17, 18, 19]. The variational multiscale method is based on a similar idea of decomposing the flow field into resolved (coarse-scale) and un-resolved (fine-scale) variables. The fine-scales are then approximated using simple algebraic operators acting on the residual of the coarse scales which give rise to additional stabilization terms in the standard Galerkin procedure [10, 11, 12]. A link between the stabilization terms and implicit sub-grid models [9] has been established. Different stabilization techniques such as the Galerkin least square (GLS) [10], the streamwise upwind Petrov Galerkin (SUPG) [11, 13], the adjoint-stabilized methods [12, 20, 15], and the orthogonal sub-scale stabilization [15] can be derived based on the type of algebraic model for the fine-scales. In the early days of the development of the VMS techniques, these methods were formulated for linear problems, and their application to the non-linear multiscale problems such as the Navier-Stokes equations though successful, depended on constructs such as transformations to linear problems such as the Oseen equations at every non-linear iteration. Recently, however, several attempts have been extended to develop non-linear VMS closure models by Codina et al. [14, 15], Bazilevs et al. [16] and many others [17, 18, 19]. These methods utilize a stabilization parameter $\tau$ which is an approximation to the inverse of the differential operator of the governing equation. This model parameter $\tau$ is typically defined in terms of a local length-scale, elemental Reynolds and Courant numbers [20] or derived from the Fourier analysis of the fine scale equation [15]. Further, it is chosen to aid optimal convergence and stability of the method.

In this work, we aim to develop a general coarse-graining approach in the context of the continuous Galerkin method, that is: (i) built using a non-linear model reduction strategy akin to Greens functions for linear problems; (ii) capable of generating a fine-scale description directly from the structure of the PDE and the underlying numerics; and (iii) model parameters are adaptive to the resolution, and are dynamically determined.

The VMS decomposition of a PDE leads to a set of coupled equations which govern the coarse-scales (resolved) and the fine-scales (un-resolved) respectively. However, the fine-scale closure problem still persists. In our approach, the dependence of the fine-scale variables on the coarse-scale variable is removed by using the optimal prediction framework developed by Chorin [21, 22]. This framework, originally developed in the context of non-equilibrium statistical mechanics, enables the higher dimensional non-linear Markovian dynamical system to be written into an exactly equivalent lower dimensional non-Markovian dynamical system [23, 24, 25]. The advantage is that the evolution of any observable in time can be represented solely in terms of the resolved variables. The cost of evaluating the resulting closure term, however, is enormous. The possible simplifications will be discussed later in the paper. Similar ideas have been put forward by Stinis [26], Parish and Duraisamy [23, 24] in context of spectral methods and discontinuous Galerkin (DG) [27] methods. The MZ-VMS philosophy presented herein closely follows the approach presented in [27, 28]. However, the derivation of the closure terms in a continuous Galerkin setting requires specific approximations. As will be discussed in this article, the contributing term to the final closure model in CG is different from that in DG. The main contribution of this work is to extend these dynamic closure models to the continuous Galerkin (CG) method previously not explored by Parish and Duraisamy [27, 28] and testing them on canonical turbulent flow problems.

The outline of the paper is as follows: We introduce the M-Z formalism in Section 2 and the VMS method in Section 3. In Section 4, we develop the VMS-MZ method in the context of a continuous Galerkin (CG) discretization. In Section 4, we derive a dynamic model for the estimation of the memory length of the convolution integral to provide a parameter free closure to the model. In the final part, we discuss results for canonical turbulence cases in Section 5. Finally, we conclude our work in Section 6.

2 The Mori-Zwanzig (M-Z) formalism

In this section, we introduce the general principles of the M-Z formalism [21]. The concept of M-Z was first introduced in the context of statistical mechanics [29, 30] but was later extended by Chorin [21] to more generalized systems. To demonstrate the basic idea of M-Z, we introduce it in a simple linear dynamical system with two degrees of freedom. Following this, we present the generalization of this concept to non-linear systems via the Generalized Langevin Equations (GLEs).

2.1 Linear Dynamical System - An Example

Consider a dynamical system with two degrees of freedom given by

[TABLE]

[TABLE]

where $x\in\mathbb{R}$ and $y\in\mathbb{R}$ are the state space, and time $t\in(0,T]$ with initial conditions: $x(0)$ and $y(0)$ provided. Our aim is to write an exact evolution equation for just one variables, say $x$, i.e.

[TABLE]

Using the appropriate integration factor and integrating Equation (2) we have the following equation in $x$:

[TABLE]

Equation (4) has three terms: (i.) the first term represents the Markovian term containing the resolved variable; (ii.) the second term is the memory integral; and (iii.) the third term represents the dependence on the initial condition of $y$ on $x$. Although Equation (4) represents the evolution of $x$ without any dependence on the second variable $y$, the flow of $x$ at any point of time not only depends on the current values of $x$ but also on its history weighted by an exponential factor. The two variable Markovian system is now converted to a one variable non-Markovian system without loss of accuracy. For coarse-grained model development for $x$, Equation (4) requires the inclusion of closure for the memory integral term.

2.2 The Generalized Langevin Equation

Although the discussion in the previous section was limited to a linear system, the M-Z formalism can be generalized to non-linear problems as well. To this end, consider the spatial discretization of a space-time problem leading to the following set of $N$ coupled ODEs in time:

[TABLE]

where $\phi=\{\hat{\phi},\tilde{\phi}\}$, $\hat{\phi}\in\mathbb{R}^{M}$ and $\tilde{\phi}\in\mathbb{R}^{N-M}$ are the modes we want to resolve and model respectively. The choice of spatial discretization can be of non-tailored basis such as spectral methods [23, 24], continuous, and the discontinuous Finite Element (FE) basis functions [27] or tailored basis obtained from purely data driven techniques such as the proper orthogonal decomposition (POD) [31]. By assuming the initial condition of the problem to be $\phi_{0}$, we aim to solve for $\hat{\phi}$ without solving for the un-resolved modes $\tilde{\phi}$ to reduce the computation and cost. However, unlike the linear problem discussed in the previous section, non-linearity restricts us from using the integration factor approach. The Mori-Zwanzig approach [21, 22, 32] allows us to cast the above non-linear problem (Equation (5)) in the form of a linear PDE in the space of initial condition variables $\phi_{0}$ and time $t$ as follows

[TABLE]

where the Liouville operator $\mathcal{L}$ corresponding to Equation (5) is defined as

[TABLE]

and initial conditions $u(\phi_{0},0)=g(\phi(\phi_{0},0))$, where $g:\mathbb{R}^{N}\rightarrow\mathbb{R}$ is a scalar valued observable. Equation (6) can be shown to have the following solution [21, 22]:

[TABLE]

Next, the semi-group notation is introduced:

[TABLE]

where $e^{t\mathcal{L}}$ is called the Koopman operator, which is an infinite dimensional linear operator which when applied to an observable $g(\phi_{0},0)$ evolves it in time $t$. As a special case, the observables are chosen to be the same as the initial states $g(\phi_{0})=\phi_{0j}$. As a consequence, Equation (6) along with Equation (9) results in the following

[TABLE]

where the last equality is a result of commutative property [22] between $\mathcal{L}$ and $e^{t\mathcal{L}}$. We decompose the right hand side of Equation (10) into spaces of resolved initial conditions and un-resolved initial conditions as follows:

[TABLE]

where $\mathcal{P}:L^{2}\rightarrow\hat{L}^{2}$ is the projection operator, where the spaces formed by all initial conditions and resolved initial conditions are denoted by $L^{2}$ and $\hat{L}^{2}$ respectively and $\mathcal{Q}=I-\mathcal{P}$. Different forms of projectors $\mathcal{P}$ can be used [23, 24, 25]. In the present work, we use a truncation projector [22, 27], i.e. the application of the projector $\mathcal{P}$ to the function $f(\hat{\phi_{0}},\tilde{\phi_{0}})$ results in the truncation of the unresolved initial conditions $f(\hat{\phi_{0}},0)$. The projector $\mathcal{P}$ acts on the space formed by initial conditions and is different from the $L^{2}$-projectors commonly used to project onto finite dimensional spaces. By applying Duhamel’s formula [22],

[TABLE]

in Equation (11), which is equivalent to the integration factor approach for linear systems, we obtain the generalized Langevin equation (GLE) [21, 22, 32] also known as the Mori-Zwanzig identity:

[TABLE]

An important observation is that Equation (13) has a similar structure to Equation (4). The first term is the Markovian term, the second term is the noise due to uncertainty in the initial condition and the last term is called the memory integral. The noise term given by $\mathcal{F}_{j}(\phi_{0},t)=e^{t\mathcal{QL}}\mathcal{Q}\mathcal{L}\phi_{0j}$ is precisely the solution to the orthogonal dynamics [21, 33] equation given by

[TABLE]

It can also be shown [22] that $\mathcal{F}_{j}(\phi_{0},t)$ lies in the null space of the projector $\mathcal{P}$ i.e. $\mathcal{PF}_{j}(\phi_{0},t)=0$. As a result, application of $\mathcal{P}$ on Equation (13) results in the following simplification

[TABLE]

Equation (15) is exact and governs the evolution of the resolved modes $\hat{\phi}$ without any dependence on the unresolved modes. However, it does not lead to reduction in the overall computational cost as it requires the solution of the orthogonal dynamics equation (Equation (14)) which is a high-dimensional PDE and solving it is intractable. However, this marks the starting point for deriving coarse-grained models based on different approximations [22, 32, 34, 35, 24, 33] to the memory term. In this paper, we use the fixed memory type model which assumes that the memory integral is correlated to its integrand at $s=0$ and has finite support in time i.e.

[TABLE]

where $\tau$ is called the memory length. Examples on application of this formalism for non-linear model reduction of toy problems can be found in [21].

3 The Variational Multiscale Method

We now present a brief overview of the variational multiscale (VMS) method, which was originally formalized by Hughes et al. [9] . Consider the following PDE on an open and bounded domain $\Omega\subset\mathbb{R}^{d}$, where $d\geq 1$ is the dimension of the problem, with a smooth boundary $\Gamma=\partial\Omega$:

[TABLE]

where the operator $\mathcal{R}:\mathbb{R}^{d}\rightarrow\mathbb{R}^{d}$ can be both linear or non-linear, the function $f:\Omega\rightarrow\mathbb{R}$ and the time varying from $t\in(0,T]$. Let, $\mathcal{V}\equiv\mathcal{H}^{1}(\Omega)$ denote the Sobolev space containing square integral functions with square integral derivatives. We define the variational problem as follows:

[TABLE]

find $u\in\mathcal{V}$ for all $w\in\mathcal{V}$, where $(\cdot,\cdot)$ denotes the $L_{2}$ inner product. The solution and weighting space are decomposed as follows:

[TABLE]

where $\oplus$ represents a direct sum of $\tilde{\mathcal{V}}$ and $\mathcal{V}^{\prime}$. From the perspective of a numerical method, $\tilde{\mathcal{V}}$ is the resolved finite dimensional space and $\mathcal{V}^{\prime}$ represents the space of functions which is not resolved. This leads to a decomposition for $u$ and $w$:

[TABLE]

[TABLE]

where $\tilde{u},\tilde{w}\in\tilde{\mathcal{V}}$ and $u^{\prime},w^{\prime}\in\mathcal{V}^{\prime}$. By substituting Equations (20) and (21) into the variational problem given by equation (18) the following is obtained,

[TABLE]

Due to the linear independency of $w^{\prime}$ and $\tilde{w}$, equation (22) is separated into the coarse scale and fine-scale equations, respectively:

[TABLE]

[TABLE]

When the coarse-scale Equation (23) is rearranged, the following is obtained,

[TABLE]

The LHS of Equation (25) contains terms present in the standard Galerkin procedure. However, it is also depends on the solution to the fine-scale equation which can be considered as the error in the coarse scale approximation. The goal of VMS sub-grid modelling is to approximate the fine-scale solution using Equation (24) and substitute it in Equation (23). Different closures can be obtained for the fine-scales depending on the type of approximation, especially when the problem is non-linear [14, 15, 16, 17, 18, 19]. However, for linear problems Hughes et al. [9] demonstrated that the fine scale solution $u^{\prime}$ is related to the Green’s function $g^{\prime}(x,y)$ of the adjoint operator $\mathcal{R}^{*}$ and the coarse scale residual as follows:

[TABLE]

The simplest approximation to which is given by

[TABLE]

Although, Equation (27) and Equation (26) have been derived for linear problems, the idea that the coarse-scale residual can be linked to the fine-scale solution, remains the basis for developing non-linear VMS models [14, 15, 16, 17, 18, 19] as well.

4 The CG-MZ-VMS Framework

In this section, we combine methods from sections 2 and 3 to derive a coarse grained model. Figure 1 shows two kinds of errors in coarse grain modeling: (i) projection; and (ii) model error. A perfect model will give us the exact projection (green line), based on some optimality condition, of the full order solution on the subspace we are approximating our solution in. For a given set of FE basis functions, the projection error can never be reduced i.e. the high-dimensional full order solution cannot be represented using a small number of FE basis functions. In the present work, we seek to develop a model to accurately predict a low-dimensional projected solution rather than the high-dimensional solution itself. We begin with the governing equation in the domain $\Omega\subset\mathbb{R}^{d}$ with the boundary $\Gamma=\partial\Omega$, where $d\geq 1$ is the dimension of the problem as follows,

[TABLE]

where $u=g$ at the boundary $\Gamma$ and time $t\in(0,T]$. The weak form of the above PDE, obtained after integration by parts, can be written as follows,

[TABLE]

where $u\in\mathcal{V}$. The Sobolev space of functions $\mathcal{V}\equiv\mathcal{H}^{1}(\Omega)$ and first derivatives are square integrable. The functional space $\mathcal{V}$ is an infinite dimensional and must be approximated by a finite dimensional approximation $\tilde{\mathcal{V}}$. We consider the tessellation of $\Omega$ into non-overlapping finite elements. The domain and boundary of an element marked by $\Omega_{e}$ and $\Gamma_{e}$ respectively. Also consider the following notations:

[TABLE]

[TABLE]

where $\Omega^{\prime}$, $\Gamma^{\prime}$ denote the interior and boundaries of all the elements respectively. Let $\tilde{\mathcal{V}}\subset C^{0}\cap\mathcal{H}^{1}(\Omega)$ denote our finite dimensional FE approximation space containing basis functions having $C^{0}$ continuity everywhere including element boundaries. Approximating $w$ by $\tilde{w}$ and $u$ by $\tilde{u}$ in Equation (29), leads to the standard Galerkin procedure given by

[TABLE]

where $\tilde{u}\in\tilde{\mathcal{V}}$. The above method, although directly applicable to diffusion dominated problems, encounters stability issues when applied to convection dominated problems. The VMS procedure provides a solution to this problem by elegantly accounting for the sub-grid scale effects. By splitting the space of the solution $u=\tilde{u}+u^{\prime}$ and the weighting function $w=\tilde{w}+w^{\prime}$, and substituting it into Equation (32), we obtain the following integral equations for the coarse and fine scales, respectively,

[TABLE]

[TABLE]

where $u^{\prime}$ and $w^{\prime}$ lie in a space orthogonal to $\tilde{u}$ and $\tilde{w}$ i.e. ${u^{\prime}},{w^{\prime}}\in{\mathcal{V^{\prime}}}$ and $f\in\tilde{\mathcal{V}}$. This idea of decomposing the full space into orthogonal spaces has previously been pursued, for instance the Orthogonal Sub-Scale (OSS) method by Codina [15] and Parish and Duraisamy [28, 24]. By substituting the resolved and un-resolved variables in terms of their modal coefficients and their basis function as follows,

[TABLE]

[TABLE]

we get the following ODE systems for modal coefficients of the coarse and fine scales,

[TABLE]

[TABLE]

where the mass matrices for resolved scales and un-resolved orthogonal scales can be written as

[TABLE]

[TABLE]

The RHS of equation By utilizing the Mori-Zwanzig procedure to integrate out variables in Equation (38) from (37), we get the following system:

[TABLE]

where the additional term to the RHS is due to the memory effects. Different ways to model the memory term have been explored in the literature [22, 32, 34, 35, 24, 33]. In the present formulation, we will use the fixed memory model $\int_{0}^{t}K(\mathbf{\tilde{a}}(t-s),s)ds\approx\tau K(\mathbf{\tilde{a}}(t),0)$ which results in the following simplification:

[TABLE]

where $\tau$ is the memory length. The memory kernel at $s=0$ is given by

[TABLE]

First, we apply $\mathcal{L}$ on $\mathbf{\mathbf{\tilde{a}}_{0}}$, resulting in the RHS of Equation (37) given by

[TABLE]

Second, we apply the projection $\mathcal{Q=I-P}$ to Equation (44) which results in the following expression,

[TABLE]

Third, we apply the Liouville operator $\mathcal{L}$ to obtain $e^{\mathcal{L}t}\mathcal{LQL}\mathbf{\mathbf{\tilde{a}}_{0}}$. The effect of application of the Liouville operator $\mathcal{L}$ to any scalar function results in the Frechet derivative evaluated in the direction of the RHS i.e. $\mathcal{L}\mathbf{a}_{0}$. For example, for a scalar function g we have the following:

[TABLE]

Recognising that $u=\mathbf{w}^{T}{\mathbf{a}}$ and applying chain rules we get

[TABLE]

Finally $e^{\mathcal{L}t}\mathcal{LQL}\mathbf{\mathbf{\tilde{a}}_{0}}$ is obtained by linearising w.r.t to $u$ and evaluating the result model at RHS of Equation (37) and (38) as follows:

[TABLE]

Finally, we apply the projector $\mathcal{P}$ which removes the dependence on un-resolved variables $\mathbf{a}^{\prime}$ and results in

[TABLE]

Equation (49) can be compactly written as,

[TABLE]

where $\Pi^{\prime}$ is the orthogonal projector onto the space of the the fine scales i.e,

[TABLE]

An important aspect of the sub-scales $u^{\prime}$ is that it is dynamic [15, 14] in nature. When these sub-scales are approximated with the inverse of the spatial operator, the resulting VMS models are non-Markovian [15, 14] and the sub-scales $u^{\prime}$ have to be tracked in time. The M-Z formalism on the other hand, precisely integrates out the time dependency of these sub-scales, making the final formulation Markovian on the resolved variables only. To make this formulation computationally tractable, we assume that the memory integral is correlated to its integrand at $s=0$ and has finite support in time. This is the main reason an approximation to the inverse of the differential operator - as is popularly used to derive VMS models - is not required here.

The final form of closure given by Equation (50) is valid for both CG and DG methods. For a simpler derivation in case of smooth orthogonal basis, readers are encouraged to read Appendix A of [27]. As shown by Parish and Duraisamy [27, 28], the main contributing term in Equation (50) to the final closure model in DG is Term 4. However, in this formulation, we assume that the fine-scales vanish at the element boundaries analogous to the concept of bubble functions [36, 37, 38, 39]. This approximation has also been used in the OSS model [15, 14]. By using this approximation, Term 2 and Term 4 in Equation (50) are neglected, resulting in the following equation:

[TABLE]

The scale separation by projection of the residual on the fine scale can be computed as follows,

[TABLE]

where $\tilde{\Pi}$ is again the $L_{2}$ projector on the finite dimensional space spanned by $\mathbf{\tilde{w}}$. This concludes the derivation of CG-MZ-VMS framework for the fixed memory type model.

5 Dynamic Memory Estimation

While the constant memory length model provides a closure to the memory term in the M-Z expression, the parameter $\tau$ should adapt to the evolving resolution and not necessarily remain constant. Another approach is to allow the parameter $\tau$ to dynamically vary in time to attempt to represent the variations of the effects of the fine-scale quantities on the coarse scales. To facilitate the stabilization of our method with fewer parameters and account for the temporal variations of the memory length, we seek a dynamic memory length model using the variational counterpart of the Germano’s identity [4, 40, 41, 42]. A similar dynamic procedure has been previously used by Oberai et al. [41] and Akkerman et al. [42] to estimate model coefficients. We begin by applying a zero-variance phase space projector with a fully resolved initial condition with the large-scale equation (Eqn (33)) to obtain an exact solution to the closure problem as following:

[TABLE]

By assuming that the memory term has a finite support we obtain

[TABLE]

Similarly, for a separate coarser mesh with weighting function $\mathbf{\hat{w}}\in\hat{\mathcal{V}}$, where $\hat{\cdot}$ signifies a coarser mesh than $\tilde{\cdot}$, the memory terms can be written as

[TABLE]

We choose $\mathbf{\tilde{w}}$ such that it spans the weighting function on the coarser mesh $\mathbf{\hat{w}}$ i.e. $\hat{\mathcal{V}}\subset\tilde{\mathcal{V}}$, which results in the following equation:

[TABLE]

where $\mathbf{G}$ is a matrix which transforms $\mathbf{\tilde{w}}$ to $\mathbf{\hat{w}}$ given by

[TABLE]

If the finer grid is obtained by element wise refinement of the coarse grid, the fine grid $\mathbf{\tilde{w}}$ basis functions span all the weighting functions on the coarser mesh $\mathbf{\hat{w}}$. By subtracting Equations (56) and (57) we obtain

[TABLE]

To obtain $\hat{u}$ project the $\tilde{u}$

[TABLE]

where $\hat{\Pi}:\tilde{L}_{2}\rightarrow\hat{L}_{2}$ is the $L_{2}$ projector on the coarse grid. Figure 2 shows the projection of a fully-resolved simulation onto $\tilde{u}$ and $\hat{u}$. This is similar to test filtering in the dynamic Smagorinsky model [4] employed in LES. Here, we assume a scaling law similar to the one proposed by Parish and Duraisamy [24] relating the memory lengths $\tau$ at two different levels of coarsening as following

[TABLE]

where $\Delta_{1}$ and $\Delta_{1}$ denote the element sizes at the fine and coarse mesh. An important observation is Equation (50) cannot be satisfied for all $\mathbf{\hat{w}}$ with a single value of $\tau$, but is true only in the average sense. To satisfy this condition, three different possibilities are considered here:

Dynamic-$\tau$-AVG: Scale the modes with their respective modal values

[TABLE]

which gives the following final form for the dynamic memory length,

[TABLE] 2. 2.

Dynamic-$\tau$-LS: Solve the overdetermined system based on some optimality condition,

[TABLE]

where L and R are given by

[TABLE]

[TABLE]

The above system can be solved using the least-squares approach which is commonly used with the DSM [4] LES model, resulting in the following expression:

[TABLE] 3. 3.

Dynamic-$\tau$-$l_{2}$: Approximate $\tau$ based on the following equation:

[TABLE]

where $||\cdot||$ denotes any kind of norm. In the present work, we have used $l_{2}$ or Euclidean norm for all our calculations. By using this averaging procedure, we obtain a value of $\tau$ that is (i) always positive; and (ii) free from division errors.

Although the steps involved in derivation of the above formulation closely follow that of the DSM [4, 5], our approach is valid for general PDEs in that the functional form of the model is not chosen based on the underlying physical phenomena. Unlike other traditional LES SGS models such as DSM [4, 6, 5], WALE [3], VREMEN [2], and Sigma [8] which are derived exclusively for the Navier-Stokes equation or other scalar transport equations, our model is not equation specific.

6 One-dimensional viscous Burgers equation

As a first step towards deriving coarse-grained models for the Navier-Stokes equation, we apply our framework to a 1-D non-linear PDE exhibiting multiscale features. To this end, let $\mathcal{V}\equiv\mathcal{H}^{1}(\Omega)$ denote the Sobolev space where our solution $u$ and weighting functions $w$ exist. The viscous Burgers equation in the domain $\Omega\subset\mathbb{R}$ is given by the following equation:

[TABLE]

with periodic boundary conditions and the time varying from $t\in(0,T]$. The weak form of Equation (69) translates into a problem of finding ${u}\in{\mathcal{V}}$ such that

[TABLE]

Using integration by parts we obtain [43],

[TABLE]

where, $(a,b)_{\Gamma^{\prime}}=\sum_{k}\int_{\Gamma^{\prime}_{k}}ab\,d\Gamma^{\prime}$ and $J(u)=\nu n_{1}.\nabla u_{1}+\nu n_{2}.\nabla u_{2}$ (where subscripts 1 and 2 denote adjacent elements sharing a boundary). By utilizing the present coarse graining procedure to Equation 70 we get

[TABLE]

where the memory term $\mathbf{\tilde{M}}K(\mathbf{\tilde{a}}(t),0)$ is given by,

[TABLE]

and $R^{\prime}$ denotes the linearization of the non-linear operator about $\tilde{u}$. Using integration by parts and neglecting the sub-scale contributions at the elemental boundaries we have

[TABLE]

where $R^{*}$ is the adjoint of the linearized operator $R^{\prime}$. The integrand is computed as follows:

[TABLE]

The resulting closure is very similar to the adjoint stabilization method [12, 9, 20] except $R(\tilde{u}(x))-\tilde{\Pi}(R(\tilde{u}(x)))$ is present instead of $R(\tilde{u}(x))$. Here, $R^{*}$ denotes the adjoint of the linearized operator $R^{\prime}$ and not the operator $R$ itself. The adjoint operator $R^{*}$ is given by

[TABLE]

Substitution of Equation (74) into Equation (72) results in the following problem for the coarse scales ${\tilde{u}}\in{\tilde{\mathcal{V}}}$:

[TABLE]

Similarly the following coarse grained model can be derived for the linear advection-diffusion equation:

[TABLE]

Equations (77) and (78) are first discretized in time using the $\theta$ family of methods [44]. Equation (77) is then linearized using the standard Picard algorithm.

6.1 Steepening of sine wave

To benchmark our coarse-grained model, the solution to the viscous Burgers equation is computed at $T=3.0$ for an initial sine profile [17, 28] on a periodic domain of length $2\pi$. To discretize in time, we use the $\theta$ family of methods [44] with $\theta=0.5$ (Crank-Nicolson). The simulation parameters for DNS and coarse grained simulations are summarized in Table 1. When sufficient resolution is available (i.e. DNS limit), the viscosity is responsible for dissipating energy at the shock. However, when the resolution is insufficient, sub-grid models are responsible for dissipating the energy. To ensure that the viscous dissipation due to large scales is negligible, viscosity $\nu$ has been set to a small value of $10^{-4}$. As a consequence, the primary contribution to the total dissipation comes from the sub-grid model.

It can be observed in Figure 3, the coarse-grained model solution approaches the projected DNS. Figures 4  and  5  show the time evolution of resolved KE and its rate of dissipation. These results indicate that the performance of all the models are comparable except the case when no sub-grid model was used or a fixed $\tau=0.01$ was used. A comparison between the solutions on the space-time diagram obtained using our dynamic-$\tau$ model, OSS, no-model and projected DNS has been presented in Figure 7. Among our Finite Memory (FM) models, the case with $\tau=0.01$ performs the worst, as can be seen in Figures 3, 4 and 5. The solution at $T=3$ improves when $\tau$ is increased to $0.11$ and becomes worse when further increased to $\tau=0.23$, which suggests the existence of an optimum $\tau$ value. The uncertainty in choosing the value of $\tau$ close to its optimum value can be reduced with a dynamic model. To this end, methods described in Section 5 are used to compute $\tau$ dynamically in Figure 6. Results indicate that the Dynamic-$\tau$-AVG, the Dynamic-$\tau$-LS and the Dynamic-$\tau$-$l_{2}$ models predict a similar magnitude of $\tau$ for the period of time considered. However, the Dynamic-$\tau$-$l_{2}$ model, which ensures positivity of the $\tau$, was found to be most stable and was used for all the following calculations in the paper. Although it is possible to use Method 1 and Method 2 by clipping $\tau$ above zero, they were found to be unstable for the TGV problem which will be discussed later in Section 7.

At $T=3$, the Dynamic-$\tau$-$l_{2}$ predicts $\tau\approx 0.05$, which supports our argument that an optimum $\tau$ exists in the range of $0.01$ and $0.23$. The $t$-model which assumes $\tau=t$ predicts a $\tau$ which does not perform well in this case and becomes unstable. However, as noted by Stinis [45], the t-model needs to be re-normalized with a coefficient for the correct prediction of the memory length i.e. $\tau=C_{N}t$. When renormalization is used, $\tau=0.014t$ is the correct representative of the memory length with $C_{N}=0.014$ as shown in Figure 6.

6.2 Burgers turbulence

To further assess the performance of our coarse-grained model for turbulence, we use it to study the Burgers turbulence problem [23, 25]. The solution to the Burgers equations exhibits some similarities to realistic turbulence, both having an inertial and a dissipation range [17]. However, the solution to Burgers turbulence is non-chaotic, unlike physically realistic turbulence obtained though the Navier-Stokes equations. To obtain the initial flow field satisfying a given energy spectrum, the following initial condition has been used

[TABLE]

where, the phase $\beta$ is randomly set from $[-\pi,\pi]$ and the energy spectra $E(k)$ is set to $5^{-5/3}$ for $k$ = 1 to 5 and $k^{-5/3}$ thereafter. Two different test cases are considered here: (i.) a high viscosity case A with $U^{*}=1$, $K_{c}=8$ and $\nu=0.01$, and (ii.) a low viscosity case B with $U^{*}=10$, $K_{c}=32$ and $\nu=0.0005$. Simulation parameters are summarized in Table 2. The two cases are considered to demonstrate the effect of the sub-grid model on moderately and highly under-resolved simulations respectively.

In the first case, when no sub-grid model is employed, the time variation of resolved kinetic energy is close to the DNS solution. However, for the low viscosity case, a sub-grid model becomes necessary. For comparison, DNS using the Fourier-Galerkin method is performed using 1024 and 4096 modes for case A and case B, respectively. The de-aliasing of the non-linear terms for the Fourier-Galerkin method is conducted by zero-padding (3/2-rule). The LES is conducted using the present coarse grained model with just 32 and 64 linear elements, respectively. For each case, results from the FM model and dynamic-$\tau$ are compared to results obtained without using a sub-grid model, the OSS model and the DNS. For the FM model, different values of $\tau$ are considered in the range where our simulations are stable. Figure 8  and  9  show the time evolution of resolved KE and its rate of dissipation for all these cases. Both these figures indicate that both the dynamic-$\tau$ model and OSS model accurately predict the time-evolution of the resolved kinetic energy in comparison DNS. Figure 10 shows the variation of $\tau$ obtained from our dynamic model which for both case A and B, predict a large variation of $\tau$ in time, suggesting the importance of the adaptive selection of $\tau$. Figure 11 also shows the energy spectra at the final time. It can be observed for Case A, that all the models perform similarly and the resolved modes are able to capture most of energy. For Case B, the present dynamic-$\tau$ model performs better than other models especially at lower wavenumbers (large scales).

7 Application to the Navier-Stokes Equations

In this section, the coarse grained model is extended to the incompressible Navier-Stokes equations. Let $\mathcal{V}_{d}\equiv(\mathcal{H}^{1}(\Omega))^{d}$ and $\mathcal{K}\equiv{L}^{2}(\Omega)$ denote the Sobolev and Lebesque spaces where our solution and weighting functions exist and $d\geq 2$ in general. The weak form of the Navier-Stokes equations consists of finding $\mathbf{u}:[0,T]\rightarrow\mathcal{V}_{d}$, $p:[0,T]\rightarrow\mathcal{K}$ such that

[TABLE]

[TABLE]

for all $[\mathbf{w},k]\in\mathcal{V}_{d}\times\mathcal{K}$.

We start by assessing how our coarse-grained model stabilizes the forced viscous Burgers equation in higher dimensions. The weak form for the Burgers equations in higher dimensions can be written as

[TABLE]

By applying integration by parts to the viscous term we get

[TABLE]

where $\mathbf{J}(\mathbf{u})=\nu\mathbf{n}_{1}\cdot\nabla\mathbf{u}_{1}+\nu\mathbf{n}_{2}\cdot\nabla\mathbf{u}_{2}$ is the diffusive flux from adjacent elements sharing a boundary. Equations (82) and (83) can be equivalently written as,

[TABLE]

By decomposing the spaces into $\mathcal{V}_{d}=\mathcal{\tilde{V}}_{d}\oplus\mathcal{V}_{d}^{\prime}$ and $\mathcal{K}=\mathcal{\tilde{K}}\oplus\mathcal{K}^{\prime}$, and applying our finite memory based framework leads to the following formulation for the coarse scales $\tilde{\mathbf{u}}\in\tilde{\mathcal{V}_{d}}$:

[TABLE]

where $\tilde{\mathcal{V}_{d}}$ is our FE approximation space, and $\mathbf{R^{\prime}}$ and $\mathbf{b^{\prime}}$ represent the linearizations of $\mathbf{R}$ and $\mathbf{b}$ with respect to $\mathbf{\tilde{u}}$. The fine-scale variable $\mathbf{q}\in{\mathcal{V}_{d}}^{\prime}$ involving projection of the residuals on the fine-scales is given by

[TABLE]

In the above equation, we assume that the fine-scales vanish at elemental boundaries [36, 37, 38, 39], a thus neglect the second term. The quantity $\mathbf{q}$ is approximated as follows

[TABLE]

By further simplifying the memory term we obtain the following:

[TABLE]

Where the second term is simplified using the Green’s identity and calculated as follows

[TABLE]

From an implementation perspective, all the above terms are computed using numerical integration at the quadrature points. This results in the following problem for the coarse scales $\tilde{\mathbf{u}}\in\tilde{\mathcal{V}_{d}}$:

[TABLE]

[TABLE]

The role of pressure herein is to impose the divergence free condition on the velocity field. In this formulation, only the velocity sub-scales have been accounted for, and the pressure terms arising from standard Galerkin procedure are retained and treated like a forcing function. This leads to additional stabilization terms to the standard Galerkin procedure given by,

[TABLE]

[TABLE]

Although closure terms were obtained for the momentum equations in Equation (92), the effect of the velocity sub-scales on the continuity equation should also be accounted for. Hence, an approximate form of the velocity sub-scales is required. To this end, consider Equation (85) in a re-arranged form:

[TABLE]

where the operators $\mathbf{R}$ and $\mathbf{J}$ are non-linear and linear respectively and $\mathbf{R^{\prime}}$ is the linearization of $\mathbf{R}$ about $\mathbf{\tilde{u}}$. For small sub-scale $\mathbf{u^{\prime}}$ [16] approximation, we have,

[TABLE]

Consequently, we can express velocity sub-scales approximately as

[TABLE]

A similar form of sub-scales was also obtained by Wang et al. [17] by writing an asymptotic series in terms of residual [16]. Finally, the effect of the sub-scales on the continuity equation is taken into consideration as follows:

[TABLE]

[TABLE]

By applying integration by parts and using the fact that sub-scales vanish at the elemental boundaries, we have the following formulation for the continuity equation:

[TABLE]

The next step is to discretize the above equation in time using the $\theta$ family of methods. This resulting variational problem at each time step is to find $\mathbf{\tilde{u}}^{n+\theta}\in\tilde{\mathcal{V}}_{d}$ and $\tilde{p}^{n+\theta}\in\tilde{\mathcal{K}}$ such that

[TABLE]

[TABLE]

[TABLE]

One way to linearize the above set of non-linear equations is by using Picard iteration based technique given by

[TABLE]

where $i+1$ and $i$ denote the present and previous iteration respectively. It can be noted that ${\mathbf{q}^{a}}$ and ${\mathbf{q}^{n+\theta,i}}$ are defined differently. This has been done so that $({\mathbf{\tilde{u}}^{n+\theta,i}}\cdot\nabla{\mathbf{\tilde{u}}^{n+\theta,i+1}},{\mathbf{\tilde{w}}})_{\Omega^{\prime}}$ and $\tau({\mathbf{q}^{n+\theta,i}}\cdot\nabla{\mathbf{\tilde{u}}^{n+\theta,i+1}},\mathbf{\tilde{w}})_{\Omega^{\prime}}$ can be merged together. This is possible because ${{\mathbf{q}}^{n+\theta,i}}$ is calculated from previous iteration variables. This is similar to Codina’s procedure [15] of adding sub-scales to the convective velocity, and a direct consequence of retaining the non-linearity in the VMS formulation [14]. Defining ${\mathbf{q}^{a}}$ in this manner allows for an implicit calculation of the memory terms which is similar to the stabilization term in [15, 11, 10, 13] as follows:

[TABLE]

[TABLE]

Equations (102), (103), (104) and (105) are iterated until convergence of the relative norm of the solution vector between two consecutive iterations is achieved.

7.1 Homogeneous Isotropic Turbulence (HIT)

In this section, we present results for decaying homogeneous isotropic turbulence (HIT) and compare it to DNS. We choose the OSS model as a reference, as it has been shown to be a good VMS closure for turbulence [46]. The HIT problem has been extensively studied in literature both numerically [47, 48, 49] and experimentally [50]. This problem is well defined in a 3-D periodic box and the initialization of the initial velocity field for DNS is done using Rogallo’s procedure [51] which assumes the following energy spectrum at initial time:

[TABLE]

where $k_{p}$ is the wavenumber at which the energy spectra peaks and A is defined as $\int_{0}^{\infty}{k^{\sigma}exp(-\sigma k^{2}/2)}dk$. The velocity in spectral space is given by

[TABLE]

where $k$ denotes the magnitude of the wavenumber vector and $\alpha$ and $\beta$ are defined as follows:

[TABLE]

where $\phi$, $\theta_{1}$, $\theta_{2}$ are uniformly distributed random numbers from 0 to $2\pi$. In all the simulations, $k_{p}=3$, $q^{2}=3$ and $\sigma=4$. Although, the initial velocity field satisfies the divergence free condition, it does not represent a physical homogeneous isotropic turbulent velocity field. To achieve this state, the field is allowed to decay to a lower $Re_{\lambda}$ where the field will resemble a more realistic velocity field due to redistribution of energy [47]. Three different initial $Re_{\lambda}$ have been considered here: $Re_{\lambda}\approx$ 65, 75 and 164 where $Re_{\lambda}$ is defined as the Reynolds number based on the Taylor microscale $\lambda$ as follows:

[TABLE]

where $u^{\prime}$ is the velocity fluctuation/root mean square (rms) of the velocity field defined as $\sqrt{2k/3}$. The initial conditions for all these cases are generated from DNS simulations by starting at a higher $Re_{\lambda}$ and allowing it to reach our target $Re_{\lambda}$ of 65, 75 and 164 respectively. The kinematic viscosity for the three different cases are set to $\nu=$ 0.001, 0.0005 and 0.0001 respectively. The LES simulations utilize $64^{3}$ linear elements for all the three $Re_{\lambda}$ cases, the results of which are presented in Figures 12,13 and 14 respectively. This allows us to study the effects of increasing $Re_{\lambda}$ by retaining the same resolution.

At a relatively lower Reynolds number of $Re_{\lambda}\approx$ 65, all the models perform fairly well in predicting the time history of the resolved kinetic energy except the fixed memory models where an arbitrary choice of $\tau$ is used. As can be seen from the kinetic energy decay plots in Figures 12(b), 13(b) and 14(b), for all the three different Reynolds numbers, the OSS model is not stable and can be seen to oscillate initially. This indicates that for the initial time period, the OSS model incorrectly forces the turbulence. We compare the energy spectra at two different times of T=2.0 and T=4.0 in Figure 13(c) and 13(d) respectively and observe that all models perform very well at the lower wavenumber modes. However, the OSS is clearly more dissipative at higher wave-numbers where it predicts a lower energy content compared to the present model and DNS. At the higher $Re_{\lambda}$ = 75 case, we find that both the Dynamic-$\tau$ and the OSS model predict reasonably the evolution of kinetic energy and rate of kinetic energy decay. If we compare the energy spectra at T=2.0 in Figure 13(c), all the models predict a higher energy content across different wave-numbers compared to the DNS solution with the dynamic-$\tau$ again performing better at higher wave-numbers again. However, at T=4.0, when it decays to a lower $Re_{\lambda}$, the performance of all the models improve, as can be seen in Figure 13(d). At the highest $Re_{\lambda}$ case, $Re_{\lambda}$ = 174, we find that the performance of all the model becomes worse compared to the DNS results. The energy spectra for this case in presented is Figure 14(c) where only the lower wavenumber modes are resolved accurately in comparison to DNS. One possible reason for the deterioration of performance at high $Re_{\lambda}$ cases is that the current VMS models are efficient in modeling the cross-stress terms and not the Reynolds stress terms [17] which dominate at higher $Re_{\lambda}$ values.

The time variation of the predicted dynamic memory length for all the three cases is shown in Figure 15. From the plots, it can be observed that the memory length increases almost linearly with time similar to the viscous Burgers equation. This is consistent with Stinis [45] of re-normalizing the t-model for stability and accuracy. Also, the predicted value for $\tau$ by our dynamic model is higher in comparison to the randomly chosen $\tau$ values for our fixed memory model. The differences in the temporal evolution between the dynamically-selected $\tau$ and the imposed $\tau$ is a further indicate that the dynamic model is necessary for calculating the memory length.

7.2 Taylor Green Vortex (TGV)

The next step in understanding the applicability of the proposed method is to employ the model on a turbulent flow that undergoes complex dynamics such the Taylor-Green vortex. This problem involves transition to turbulence-like flow, as well as decay. Models such as the Smagorinsky which have been derived based on assumptions of homogeneity, isotropy and balance between sub-grid production and dissipation [7] might not optimally preform in such flows where there is non-homogenity and transition to turbulence. Similar to HIT, this problem is well-defined on a 3-D periodic box with smooth initial conditions which are given as follows:

[TABLE]

where $u,v,w$ denote the velocity in $x,y$ and $z$ directions respectively and $x,y,z$ $\in[-\pi L,\pi L]$. To study this problem, three different Reynolds numbers are considered: $Re={U_{o}L\over\nu}=$ 400, 800 and 1600. The values for $L$ and $U_{o}$ are unity and the Re is changed solely by varying the kinematic viscosity $\nu$. The initial conditions for the velocity field are kept the same for all the cases.

The profiles for the resolved kinetic energy for $Re=$ 400, 800 and 1600 are shown in Figures 16(a), 17(a) and 18(a) respectively. The evolution of the kinetic energy indicates that both the OSS and Dynamic-$\tau$ perform well with only $32^{3}$ degrees of freedom. Similar trends are also observed for the rate of KE energy decay for $Re=$ 400, 800 and 1600 in Figures 16(b),17(b) and 18(b) respectively, where the fixed $\tau$ models fail to accurately predict the correct results in comparison to DNS. When $48^{3}$ and $64^{3}$ degrees of freedoms are used for $Re=$ 800 and 1600 respectively, there is an overall improvement in the results for the fixed memory model. The present dynamic model and the OSS model perform well at finer resolutions.

Figures 16(c), 17(c) and 18(c), and Figures 16(d), 17(d) and 18(d) show the energy spectra of the resolved velocity fields at two time instants $T=5.0$ and $T=10.0$ respectively. At $T=5.0$, all the models are in agreement with DNS at the low wavenumber modes even with just $32^{3}$ degrees of freedom. However, the constant $\tau$ models produces a build-up of energy which grows with $Re$ at high-wavenumber modes. This suggests that either an incorrect value of $\tau$ is used for the FM models or the assumption of constant memory length throughout the simulation is not very accurate. As a result, model with constant memory length is not capable of producing enough dissipation and the energy increases at the high wavenumber modes. At a later time $T=10.0$, a similar trend is also observed with the constant $\tau$ model where energy increases at high-wave numbers. On the other hand, the dynamic $\tau$-model and OSS do not result in energy increase at high wavenumbers. Although our dynamic model and OSS model perform closely for the $32^{3}$ cases, at higher resolutions OSS is clearly more dissipative wherein lower energy is present at high wavenumber especially for the high Reynolds number case. In spite of the OSS model and the dynamic-$\tau$ model performing closely, the stabilization parameter in OSS and the memory length in dynamic-$\tau$ is computed differently.

8 Conclusion

The Variational Multiscale method and the Mori-Zwanzig formalism are combined within the Continuous Galerkin method to develop coarse grained models for multiscale PDEs. This approach utilizes the Variational Multiscale method to separate scales with the capability of the Mori-Zwanzig formalism to represent the impact of unresolved dynamics on the resolved dynamics. This approach - in a similar spirit to existing non-linear VMS models [14, 15, 16, 17, 18, 19] - is developed to provide sub-grid scale models without phenomenological assumptions. This procedure is generalizable and can potentially be applied to arbitrarily complex non-linear multiscale PDEs. In context of turbulent flows, this approach provides a general framework for large-eddy simulation that eschews assumptions such as those based on energy balance between scale. Sub-grid scale models are developed for the Burgers equation and the Navier-Stokes using the proposed approach. The sub-grid scale models include a parameter called the memory length, $\tau$, which represents the the time correlation of unresolved dynamics, and controls the stabilization. We impose different memory lengths $\tau$ and observe that there is an optimum memory length $\tau$ which provides results comparable to the full order solution [28]. To avoid the imposition of an adhoc memory length, and recognizing that the model should adapt to the instantaneous level of resolution, we derived a dynamic-$\tau$ model and found that it can accurately predict the temporal evolution of $\tau$. The predicted value of $\tau$ was observed to be generally linear in time, as conceptualized by the renormalized t-model [45]. In general, for the range of problems that were investigated, the proposed technique performs favorably in comparison to existing counterparts.

This work was focused on fixed memory type models leading to Markovian closures, however, other approaches leading to non-Markovian type closures can be implemented [28]. Alternatively, different approximations to the orthogonal dynamics [33, 25] can be used to construct models. The single memory length model was - in part - successful because of the nature of the problems investigated herein. Future extensions to highly anisotropic and inhomogeneous problems will require the development of local definitions for the memory length.

Acknowledgement

This research was funded by the AFOSR under the project LES Modeling of Non-local effects using Statistical Coarse-graining, grant number FA9550-16-1-0309. The authors also thank Dr. Eric Parish for providing DNS results for comparison for the HIT and TGV cases, and Dr. Daniel Foti for his valuable suggestions.