We present enhanced physics-based finite element schemes for two families of turbulence models, the NS-πœ” models and the Stolz-Adams approximate deconvolution models. These schemes are delicate extensions of a method created for the Navier-Stokes equations in Rebholz (2007), that achieve high physical fidelity by admitting balances of both energy and helicity that match the true physics. The schemes' development requires carefully chosen discrete curl, discrete Laplacian, and discrete filtering operators, in order to permit the necessary differential operator commutations.

1. Introduction

Energy and helicity play fundamental physical roles in the development of turbulent flow. They are conserved in inviscid flow, are delicately balanced when viscous and body forces are present [1, 2], are believed critical in flow structure development [3], and are cascaded jointly through the inertial range until they are removed by viscous forces [4]. Furthermore, their dissipation is closely tied to turbulence decay, and helicity has the topological interpretation that it is nonzero if and only if the velocity field is not reflectionally symmetric, that is, turbulent [3]. Therefore, the accurate treatment of these quantities is paramount for a computation to accurately predict fluid motion.

Unfortunately, often one or both of these quantities is not accurately treated in most models (in their continuous form), and thus cannot be correctly handled in their discretizations. Only a very few models correctly balance both energy and helicity, the most prominent being the Navier-Stokes equations (NSE, the true physical model), NS-𝛼 and related models [5, 6], the Stolz-Adams family of approximate deconvolution models (ADM) [7–11], and the recently introduced NS-πœ” family of regularization models [12, 13]. Although these models yield, in a sense, the correct balances of energy and helicity, standard discretizations for these models do not preserve this important feature. For example, the usual Crank-Nicolson method for any of these schemes nonphysically creates helicity via the nonlinearity, thus altering the balance and removing physical meaning from solutions. Correctly accounting for both of these fundamental quantities is difficult; however it is very desirable since incorrect treatment of these quantities can lead to their nonphysical accumulations or degradations, and thus computed solutions that lack physical relevance.

A solution to this problem was first introduced in 2004 by Liu and Wang in [14], where they proposed and tested the EHPS, a finite difference scheme for axisymmetric flow (derived from the vorticity-stream formulation of the NSE) that more accurately treats both energy and helicity than standard schemes, and found excellent results. The idea was later developed for the full 3d NSE in [15], where the main ideas were to use the rotational form with a projected vorticity in a finite element scheme. Since direct computation of the NSE is limited to lower Reynolds numbers, this scheme was extended for use with NS-𝛼 type models [16], where again excellent numerical results were reported. Herein we show how, through carefully chosen definitions of the discrete filter, discrete Laplacian and discrete curl operators, the enhanced-physics based scheme of [15], can be extended for use on the NS-πœ” and ADM families of models. The NS-πœ” family of models, in addition to having a sound mathematical foundation [13], provides several advantages over the other classes of energy and helicity preserving models, including a more physical treatment of mass conservation and allowing for more efficient computations [12]. Therefore the development of an enhanced-physics based scheme for this attractive family of models is warranted. For completeness we also develop an enhanced-physics scheme for the ADM, so that each class of models that accurately balances energy and helicity has such a scheme developed for it.

2. The NS-πœ” Family of Models

Given a domain Ξ© and endtime 𝑇, the NS-πœ” model is given by

π‘£π‘‘ξ€·βˆ’π‘£Γ—βˆ‡Γ—π‘£ξ€Έ]+βˆ‡π‘žβˆ’πœˆΞ”π‘£=𝑓inΩ×(0,𝑇,(2.1)βˆ‡β‹…π‘£=0in[]Ω×0,𝑇,(2.2)𝑣(0)=𝑣0inΞ©,(2.3) where 𝑣 denotes velocity, 𝑣0 the initial velocity, π‘ž pressure, 𝑓 a forcing, 𝜈 the kinematic viscosity, and overbar the Helmholtz filter with filtering radius 𝛼: ξ€·πœ™βˆΆ=πΉπœ™βˆΆ=βˆ’π›Ό2ξ€ΈΞ”+πΌβˆ’1πœ™.(2.4) Extensions of this model to include van Cittert approximate deconvolution, to achieve higher-order consistency error to the NSE, are discussed in [12]. By defining the deconvolution operator π·π‘βˆΆ=𝑁𝑛=0(πΌβˆ’πΉ)𝑛(2.5) and replacing the filtered term 𝑣 in (2.1) by 𝐷𝑁𝑣, where 𝑁 is a positive integer typically chose 1≀𝑁≀7, we have the NS-πœ” deconvolution models: π‘£π‘‘ξ€·βˆ’π‘£Γ—βˆ‡Γ—π·π‘π‘£ξ€Έ],+βˆ‡π‘žβˆ’πœˆΞ”π‘£=𝑓inΩ×(0,π‘‡βˆ‡β‹…π‘£=0in[],Ω×0,𝑇𝑣(0)=𝑣0inΞ©.(2.6) Since the deconvolution operator acts as an approximate inverse to the filter [17], its use creates a more formally accurate class of models. Note also that the 𝑁=0 model is the NS-πœ” model (2.1)-(2.2) since 𝐷0=𝐼.

NS-πœ” must be coupled to initial and boundary conditions in order to be a predictor of turbulent flow. As our interest is the accurate balance of energy and helicity, we consider for simplicity a domain Ξ© to be a 3d rectangular box with periodic boundary conditions. These boundary conditions are not likely to be physically realized, however their use here isolates the fundamental numerical difficulties from the additional complications that arise from more physically relevant boundary conditions. The work herein can be extended to most other boundary conditions, but would need to be done on an individual basis since boundary condition themselves can (and usually do) contribute to the energy and helicity balances.

Recall the balances of energy and helicity for NS-πœ” from [13]. Denote the 𝐿2 norm and inner product by β€–β‹…β€– and (β‹…,β‹…), respectively.

Lemma 2.1. Given initial condition 𝑣0∈𝐿2(Ξ©), endtime 𝑇, and forcing π‘“βˆˆπΏ2(0,𝑇;𝐿2(Ξ©)), the NS-πœ” model admits the following exact balances for energy and helicity 1Energy∢2‖𝑣(𝑇)β€–2ξ€œ+πœˆπ‘‡01β€–βˆ‡Γ—π‘£(𝑑)‖𝑑𝑑=2‖‖𝑣0β€–β€–2+ξ€œπ‘‡0ξ€·(𝑓(𝑑),𝑣(𝑑))𝑑𝑑,(2.7)HelicityβˆΆπ‘£(𝑇),βˆ‡Γ—ξ€Έπ‘£(𝑇)+𝛼2ξ€·βˆ‡Γ—π‘£(𝑇),(βˆ‡Γ—)2ξ€Έξ€œπ‘£(𝑇)+2πœˆπ‘‡0ξ€·ξ€·βˆ‡Γ—π‘£(𝑑),(βˆ‡Γ—)2𝑣+ξ€·(𝑑)(βˆ‡Γ—)2𝑣(𝑑),(βˆ‡Γ—)3𝑣=ξ€·(𝑑)𝑑𝑑𝑣0,βˆ‡Γ—π‘£0ξ€Έ+𝛼2ξ€·βˆ‡Γ—π‘£0,(βˆ‡Γ—)2𝑣0ξ€Έξ€œ+2𝑇0𝑓(𝑑),βˆ‡Γ—ξ€Έπ‘£(𝑑)𝑑𝑑.(2.8)

Remark 2.2. While the energy balance of NS-πœ” is the same as for the NSE, the helicity balance is the same as for the spatially filtered NSE. It is proven in [13] that up to a filtering radius dependent length scale, energy and helicity are cascaded jointly through the inertial range, just as in true fluid flow [4].

Remark 2.3. Balances for NS-πœ” with deconvolution can be easily derived by multiplying NS-πœ”-deconvolution by 𝑣 and βˆ‡Γ—π·π‘π‘£, for energy and helicity respectively, then integrating by parts. The energy balance will be identical to the NS-πœ” case, while an analogous result for the helicity balance will be obtained where the 𝑣 terms in (2.8) are replaced with 𝐷𝑁1/2𝑣.

3. An Enhanced-Physics Based Finite Scheme for NS-πœ”

The main obstacles in extending the scheme of [15] to NS-πœ” is the need for the discrete filter and curl operators to commute, and that a form of the vorticity is discretely-divergence free. These difficulties were not present in the extension to NS-𝛼 in [16], and hence this development is fundamentally more difficult. The extension is also quite delicate, but was achieved through the use of a new discrete implementation of the Helmholtz filter that is defined in terms of a discrete curl, which is discretely-divergence free via a suitably chosen projection in its definition. We present the scheme and show its energy and helicity balances after some necessary notation.

Denote the zero mean periodic subspaces of 𝐻1(Ξ©) and 𝐿2(Ξ©), respectively, by 𝐻1#(Ξ©) and 𝐿2#(Ξ©), let (π‘‹β„Ž,π‘„β„Ž)βŠ‚(𝐻1#,𝐿2#) be inf-sup stable finite element spaces (e.g., Taylor-Hood elements [18–20]) on a regular, conforming mesh of Ξ©. Let π‘‰β„Ž denote the discretely-divergence free subspace of π‘‹β„Ž: π‘‰β„Ž=ξ€½π‘£β„Žβˆˆπ‘‹β„Ž,ξ€·βˆ‡β‹…π‘£β„Ž,π‘žβ„Žξ€Έ=0βˆ€π‘žβ„Žβˆˆπ‘‹β„Žξ€Ύ.(3.1) We now define the discrete curl operator, which will be enforced to be discretely-divergence free by its construction.

Definition 3.1. Define the discrete curl operator to be the 𝐿2 projection of the curl operation into π‘‰β„Ž: Given πœ™βˆˆπ»1(Ξ©), curlβ„Žcurlβ„Žπœ™ is the unique element in π‘‰β„Ž satisfying ξ€·curlβ„Žπœ™,π‘£β„Žξ€Έ=ξ€·βˆ‡Γ—πœ™,π‘£β„Žξ€Έβˆ€π‘£β„Žβˆˆπ‘‰β„Ž.(3.2)
Note it is obvious that the linear operator curlβ„Ž is well defined, and also that it is self-adjoint since the usual curl operator is in this setting.
Since the Laplacian operator applied to a divergence free function πœ’ can be written as
Ξ”πœ’=βˆ’βˆ‡Γ—(βˆ‡Γ—πœ’)+βˆ‡(βˆ‡β‹…πœ’)=βˆ’βˆ‡Γ—(βˆ‡Γ—πœ’),(3.3) it is natural to define the discrete Laplacian in terms of the discrete curl operator, as follows.

Definition 3.2. Define the discrete Laplacian Ξ”β„ŽβˆΆπ»1(Ξ©)β†’π‘‰β„Ž by Ξ”β„ŽβˆΆ=βˆ’curlβ„Žcurlβ„Ž.(3.4)
We define the discrete filter πΉβ„Ž in terms of the discrete Laplacian:
πΉβ„Žπœ™βˆΆ=πœ™β„Žξ€·βˆΆ=βˆ’π›Ό2Ξ”β„Žξ€Έ+πΌβˆ’1πœ™.(3.5) Note that, by construction, the filter commutes with the discretely-divergence free projected curl operator and thus also the discrete Laplacian; recall these were the obstacles needed overcome for the development of an enhanced-physics based scheme for NS-πœ”, which we present now.

Algorithm 1 (an enhanced-physics based scheme for NS-πœ”). Given an endtime 𝑇, timestep Δ𝑑>0, filtering radius 𝛼>0, kinematic viscosity 𝜈>0, initial condition 𝑣0β„Žβˆˆπ‘‰β„Ž, forcing π‘“βˆˆπΏ2(0,𝑇;𝐿2(Ξ©)), set π‘€βˆΆ=𝑇/Δ𝑑 and for 𝑛=1,2,…,𝑀 find (π‘£π‘›β„Ž,π‘π‘›β„Ž)∈(π‘‹β„Ž,π‘„β„Ž) satisfying forall(πœ’β„Ž,π‘žβ„Ž)∈(π‘‹β„Ž,π‘„β„Ž), 1ξ€·π‘£Ξ”π‘‘β„Žπ‘›+1βˆ’π‘£π‘›β„Ž,πœ’β„Žξ€Έβˆ’ξ‚΅π‘£β„Žπ‘›+1/2Γ—curlβ„Žπ‘£β„Žβ„Žπ‘›+1/2,πœ’β„Žξ‚Άβˆ’ξ€·π‘β„Žπ‘›+1/2,βˆ‡β‹…πœ’β„Žξ€Έξ€·Ξ”βˆ’πœˆβ„Žπ‘£β„Žπ‘›+1/2,πœ’β„Žξ€Έ=𝑓𝑑𝑛+1/2ξ€Έ,πœ’β„Žξ€Έ,ξ€·(3.6)βˆ‡β‹…π‘£β„Žπ‘›+1,π‘žβ„Žξ€Έ=0.(3.7)

Remark 3.3. That solutions exist for Algorithm 1 follows analogous to the case for NS-𝛼 given in [16]. It is straight forward, following similar to [16, 19] to prove optimal convergence to a smooth, weak NSE solution for Taylor-Hood elements.

We now prove the balances of energy and helicity by the scheme. Note the energy balance is, in a sense, a stability estimate.

Theorem 3.4. The scheme of Algorithm 1 for NS-πœ” admits the following exact energy and helicity balances 1Energy∢2β€–β€–π‘£π‘€β„Žβ€–β€–2+πœˆΞ”π‘‘π‘€βˆ’1𝑛=0β€–β€–curlβ„Žπ‘£β„Žπ‘›+1/2β€–β€–2=12‖‖𝑣0β„Žβ€–β€–2+Ξ”π‘‘π‘€βˆ’1𝑛=0𝑓𝑑𝑛+1/2ξ€Έ,π‘£β„Žπ‘›+1/2ξ€Έ,ξ‚΅(3.8)HelicityβˆΆπ‘£π‘€β„Žβ„Ž,curlβ„Žπ‘£π‘€β„Žβ„Žξ‚Ά+𝛼2ξ‚΅curlβ„Žπ‘£π‘€β„Žβ„Ž,curl2β„Žπ‘£π‘€β„Žβ„Žξ‚Ά+2πœˆΞ”π‘‘π‘€βˆ’1𝑛=0βŽ›βŽœβŽœβŽœβŽξ‚΅curlβ„Žπ‘£β„Žβ„Žπ‘›+1/2,curl2β„Žπ‘£β„Žβ„Žπ‘›+1/2ξ‚Ά+𝛼2βŽ›βŽœβŽœβŽœβŽcurl2β„Žπ‘£π‘›+12β„Žβ„Ž,curl3β„Žπ‘£β„Žβ„Žπ‘›+1/2⎞⎟⎟⎟⎠⎞⎟⎟⎟⎠=𝑣0β„Žβ„Ž,curlβ„Žπ‘£0β„Žβ„Žξ‚Ά+𝛼2ξ‚΅curlβ„Žπ‘£0β„Žβ„Ž,curl2β„Žπ‘£0β„Žβ„Žξ‚Ά+2Ξ”π‘‘π‘€βˆ’1𝑛=0𝑓𝑑𝑛+1/2ξ€Έ,curlβ„Žπ‘£β„Žβ„Žπ‘›+1/2ξ‚Ά.(3.9)

Remark 3.5. The energy and helicity balances of the theorem are discrete analogs to the continuous NS-πœ” model's balances from Lemma 2.1, and thus are also discrete analogs of the energy balance of the NSE and helicity balance of the spatially filtered NSE respectively.

Proof. The energy balance follows from choosing πœ’β„Ž=π‘£β„Žπ‘›+1/2 in (3.6), π‘žβ„Ž=π‘β„Žπ‘›+1/2 in (3.7), adding the equations then summing over timesteps.
For the helicity balance, begin by choosing πœ’β„Ž=curlβ„Žπ‘£β„Žβ„Žπ‘›+1/2. This will immediately vanish the nonlinear term since the cross of two vectors is perpendicular to each of them, and will eliminate the pressure term since the discrete curl is discretely-divergence free. This gives us the equation
1ξ‚΅π‘£Ξ”π‘‘β„Žπ‘›+1βˆ’π‘£π‘›β„Ž,curlβ„Žπ‘£β„Žβ„Žπ‘›+1/2ξ‚Άξ‚΅Ξ”βˆ’πœˆβ„Žπ‘£β„Žπ‘›+1/2,curlβ„Žπ‘£β„Žβ„Žπ‘›+1/2ξ‚Ά=𝑓𝑑𝑛+1/2ξ€Έ,curlβ„Žπ‘£β„Žβ„Žπ‘›+1/2ξ‚Ά.(3.10) Since the discrete curl and discrete filter commute and are self-adjoint, the first term can be expanded and reduced via 1ξ‚΅π‘£Ξ”π‘‘β„Žπ‘›+1βˆ’π‘£π‘›β„Ž,curlβ„Žπ‘£β„Žβ„Žπ‘›+1/2ξ‚Ά=12Ξ”π‘‘ξ‚΅ξ‚΅π‘£β„Žβ„Žπ‘›+1,curlβ„Žπ‘£β„Žπ‘›+1ξ‚Άβˆ’ξ‚€π‘£π‘›β„Žβ„Ž,curlβ„Žπ‘£β„Žπ‘›+1+ξ‚΅π‘£β„Žβ„Žπ‘›+1,curlβ„Žπ‘£π‘›β„Žξ‚Άβˆ’ξ‚€π‘£π‘›β„Žβ„Ž,curlβ„Žπ‘£π‘›β„Žξ‚ξ‚Ά=12Ξ”π‘‘ξ‚΅ξ‚΅π‘£β„Žβ„Žπ‘›+1,curlβ„Žπ‘£β„Žπ‘›+1ξ‚Άβˆ’ξ‚€π‘£π‘›β„Žβ„Ž,curlβ„Žπ‘£π‘›β„Žξ‚ξ‚Ά=12Ξ”π‘‘ξ‚΅ξ‚΅π‘£β„Žβ„Žπ‘›+1,curlβ„Žπ‘£β„Žβ„Žπ‘›+1ξ‚Ά+𝛼2ξ‚΅curlβ„Žπ‘£β„Žβ„Žπ‘›+1,curl2β„Žπ‘£β„Žβ„Žπ‘›+1ξ‚Άβˆ’ξ‚€π‘£π‘›β„Žβ„Ž,curlβ„Žπ‘£π‘›β„Žβ„Žξ‚βˆ’π›Ό2ξ‚€curlβ„Žπ‘£π‘›β„Žβ„Ž,curl2β„Žπ‘£π‘›β„Žβ„Ž.(3.11) Similarly for the viscous term, we have that ξ‚΅Ξ”βˆ’πœˆβ„Žπ‘£β„Žπ‘›+1/2,curlβ„Žπ‘£β„Žβ„Žπ‘›+1/2ξ‚Άξ‚΅=𝜈curl2β„Žπ‘£β„Žπ‘›+1/2,curlβ„Žπ‘£β„Žβ„Žπ‘›+1/2ξ‚Άξ‚΅=𝜈curl2β„Žπ‘£β„Žβ„Žπ‘›+1/2,curlβ„Žπ‘£β„Žβ„Žπ‘›+1/2ξ‚Άβˆ’π›Ό2πœˆξ‚΅curl2β„ŽΞ”β„Žπ‘£β„Žβ„Žπ‘›+1/2,curlβ„Žπ‘£β„Žβ„Žπ‘›+1/2ξ‚Άξ‚΅=𝜈curl2β„Žπ‘£β„Žβ„Žπ‘›+1/2,curlβ„Žπ‘£β„Žβ„Žπ‘›+1/2ξ‚Ά+𝛼2πœˆξ‚΅curl3β„Žπ‘£β„Žβ„Žπ‘›+1/2,curl2β„Žπ‘£β„Žβ„Žπ‘›+1/2ξ‚Ά.(3.12) Summing over timesteps complete the proof.

Remark 3.6. From the proof, the need for the intricate design of the enhanced-physics scheme becomes more clear. Without the use of the discrete curl in the scheme, a test function could not have been chosen to vanish the nonlinearity in the helicity balance, which in turn would cause the nonphysical creation and dissipation of helicity instead of cascading it from the large to small scales. The derivation of differences in helicities at successive timesteps in (3.11) was made possible because the discrete curl and discrete filter commute. Commutation of these operators with the discrete Laplacian was necessary in the analysis of the viscous term.

3.1. Extension to the NS-πœ” Family of Deconvolution Models

The development of the discrete curl, Laplacian and filter operators allows for an easy extension to the deconvolution family of models (2.4)–(2.6) provided the natural definition for the discrete deconvolution operator is employed: π·β„Žπ‘βˆΆ=𝑁𝑛=0ξ€·πΌβˆ’πΉβ„Žξ€Έπ‘›.(3.13) Thus the discrete deconvolution can be seen to be positive and a polynomial in the discrete filter, and as such that will be self-adjoint and commute with the discrete curl, discrete Laplacian, and discrete filter. This leads us to the following enhanced-physics based algorithm for NS-πœ” family of deconvolution models.

Algorithm 2 (an enhanced-physics based scheme for NS-πœ”-deconvolution). Given an endtime 𝑇, timestep Δ𝑑>0, filtering radius 𝛼>0, deconvolution order 𝑁β‰₯0, kinematic viscosity 𝜈>0, initial condition 𝑣0β„Žβˆˆπ‘‰β„Ž, forcing π‘“βˆˆπΏ2(0,𝑇;𝐿2(Ξ©)), set π‘€βˆΆ=𝑇/Δ𝑑 and for 𝑛=1,2,…,𝑀 find (π‘£π‘›β„Ž,π‘π‘›β„Ž)∈(π‘‹β„Ž,π‘„β„Ž) satisfying forall(πœ’β„Ž,π‘žβ„Ž)∈(π‘‹β„Ž,π‘„β„Ž), 1ξ€·π‘£Ξ”π‘‘β„Žπ‘›+1βˆ’π‘£π‘›β„Ž,πœ’β„Žξ€Έβˆ’ξ‚΅π‘£β„Žπ‘›+1/2Γ—curlβ„Žπ·β„Žπ‘π‘£β„Žβ„Žπ‘›+1/2,πœ’β„Žξ‚Άβˆ’ξ€·π‘β„Žπ‘›+1/2,βˆ‡β‹…πœ’β„Žξ€Έξ€·Ξ”βˆ’πœˆβ„Žπ‘£β„Žπ‘›+1/2,πœ’β„Žξ€Έ=𝑓𝑑𝑛+1/2ξ€Έ,πœ’β„Žξ€Έ,ξ€·(3.14)βˆ‡β‹…π‘£β„Žπ‘›+1,π‘žβ„Žξ€Έ=0.(3.15)

Remark 3.7. As in the case of no deconvolution, proofs of solution existence and convergence follow similarly to work in [16, 21].

We now provide the energy and helicity balances for this algorithm.

Theorem 3.8. The scheme of Algorithm 2 for NS-πœ”-deconvolution admits the following exact energy and helicity balances 1Energy∢2β€–β€–π‘£π‘€β„Žβ€–β€–2+πœˆΞ”π‘‘π‘€βˆ’1𝑛=0β€–β€–curlβ„Žπ‘£β„Žπ‘›+1/2β€–β€–2=12‖‖𝑣0β„Žβ€–β€–2+Ξ”π‘‘π‘€βˆ’1𝑛=0𝑓𝑑𝑛+1/2ξ€Έ,π‘£β„Žπ‘›+1/2ξ€Έ.(3.16)𝐷HelicityβˆΆβ„Žπ‘ξ€Έ1/2π‘£π‘€β„Žβ„Ž,curlβ„Žξ€·π·β„Žπ‘ξ€Έ1/2π‘£π‘€β„Žβ„Žξ‚Ά+𝛼2ξ‚΅curlβ„Žξ€·π·β„Žπ‘ξ€Έ1/2π‘£π‘€β„Žβ„Ž,curl2β„Žξ€·π·β„Žπ‘ξ€Έ1/2π‘£π‘€β„Žβ„Žξ‚Ά+2πœˆΞ”π‘‘π‘€βˆ’1𝑛=0ξ‚΅ξ‚΅curlβ„Žξ€·π·β„Žπ‘ξ€Έ1/2π‘£β„Žβ„Žπ‘›+1/2,curl2β„Žξ€·π·β„Žπ‘ξ€Έ1/2π‘£β„Žβ„Žπ‘›+1/2ξ‚Άξ‚Ά+2πœˆπ›Ό2Ξ”π‘‘π‘€βˆ’1𝑛=0ξ‚΅ξ‚΅curl2β„Žξ€·π·β„Žπ‘ξ€Έ1/2π‘£β„Žβ„Žπ‘›+1/2,curl3β„Žξ€·π·β„Žπ‘ξ€Έ1/2π‘£β„Žβ„Žπ‘›+1/2=ξ‚΅ξ€·π·ξ‚Άξ‚Άβ„Žπ‘ξ€Έ1/2𝑣0β„Žβ„Ž,curlβ„Žξ€·π·β„Žπ‘ξ€Έ1/2𝑣0β„Žβ„Žξ‚Ά+𝛼2ξ‚΅curlβ„Žξ€·π·β„Žπ‘ξ€Έ1/2𝑣0β„Žβ„Ž,curl2β„Žξ€·π·β„Žπ‘ξ€Έ1/2𝑣0β„Žβ„Žξ‚Ά+2Ξ”π‘‘π‘€βˆ’1𝑛=0𝑓𝑑𝑛+1/2ξ€Έ,curlβ„Žξ€·π·β„Žπ‘ξ€Έ1/2π‘£β„Žβ„Žπ‘›+1/2ξ‚Ά.(3.17)

Proof. The energy balance follows exactly as in the case without deconvolution. The helicity balance follows similar to the case of NS-πœ”; choose πœ’β„Ž=π‘π‘’π‘Ÿπ‘™β„Žπ·β„Žπ‘ to vanish the nonlinearity and pressure, then use commutativity of the differential operators and that the deconvolution operator is positive and self-adjoint to obtain the result.

4. Extension to the Stolz-Adams ADM

The Stolz-Adams ADM family take a similar mathematical form to the NS-πœ” models, and the tools used to derive an enhanced-physics based scheme are sufficient to devise such a scheme for the ADM. For simplicity, we consider the ADM secondary regularization, often called time relaxation [8, 22]. This term does not correctly balance energy or helicity, however it is negligable away from the smallest resolved scales [22], and so it has little total effect on the overall balances.The ADM can be written in rotational form as

π‘£π‘‘βˆ’π·π‘π‘£β‹…βˆ‡π·π‘π‘£+βˆ‡π‘žβˆ’πœˆΞ”π‘£=𝑓,(4.1)βˆ‡β‹…π·π‘π‘£=0,(4.2)𝑣(0)=𝑣0,(4.3) where filtering and deconvolution are defined exactly as for NS-πœ”. Its balances for energy and helicity are derived in [10], and follow from multiplying (4.1) by (βˆ’π›Ό2Ξ”+𝐼)𝐷𝑁𝑣 and (βˆ’π›Ό2Ξ”+𝐼)(βˆ‡Γ—π·π‘π‘£), respectively, then integrating over the domain, integrating by parts, and reducing. These balances are as follows


Remark 4.1. Technical arguments in [11] discuss how these balances are equivalent, in a sense, to the energy and helicity balances of the spatially filtered NSE. This is not a coincidence, since the ADM can be derived as an approximation of the spatially filtered NSE [17].

Though perhaps not obvious, the discrete operators which permitted the development of an enhanced-physics scheme for the NS-πœ” family above are sufficient to devise such a scheme for the ADM. We now present this numerical scheme, which provides discrete analogs for the balances (4.4) and (4.5).

Algorithm 3 (an enhanced-physics based scheme for the ADM). Given an endtime 𝑇, timestep Δ𝑑>0, filtering radius 𝛼>0, deconvolution order 𝑁β‰₯0, kinematic viscosity 𝜈>0, initial condition 𝑣0β„Žβˆˆπ‘‰β„Ž, forcing π‘“βˆˆπΏ2(0,𝑇;𝐿2(Ξ©)), set π‘€βˆΆ=𝑇/Δ𝑑 and for 𝑛=1,2,…,𝑀 find (π‘£π‘›β„Ž,π‘žπ‘›β„Ž)∈(π‘‹β„Ž,π‘„β„Ž) satisfying forall(πœ’β„Ž,π‘Ÿβ„Ž)∈(π‘‹β„Ž,π‘„β„Ž), 1ξ€·π‘£Ξ”π‘‘β„Žπ‘›+1βˆ’π‘£π‘›β„Ž,πœ’β„Žξ€Έβˆ’ξ‚΅π·β„Žπ‘π‘£β„Žπ‘›+1/2Γ—curlβ„Žπ·β„Žπ‘π‘£β„Žβ„Žπ‘›+1/2,πœ’β„Žξ‚Άβˆ’ξ‚΅π‘žβ„Žβ„Žπ‘›+1/2,βˆ‡β‹…πœ’β„Žξ‚Άξ€·Ξ”βˆ’πœˆβ„Žπ‘£β„Žπ‘›+1/2,πœ’β„Žξ€Έ=𝑓𝑑𝑛+1/2ξ€Έ,πœ’β„Žξ€Έ,ξ€·(4.6)βˆ‡β‹…π‘£β„Žπ‘›+1,π‘Ÿβ„Žξ€Έ=0.(4.7)

Remark 4.2. Since by definition the discrete curl is discretely-divergence free, the discrete Laplacian must be also. Hence filtered elements in π‘‰β„Ž remain in π‘‰β„Ž, as do deconvolved ones. Thus, (4.7) is sufficient as a discrete enforcement of (4.2).

Remark 4.3. Similar to the NS-πœ” algorithms, it is straight forward and follows similar to work in [16, 23] that this algorithm admits stable solutions that convergence to weak NSE solutions or spatially filtered NSE solutions.

Theorem 4.4. The scheme of Algorithm 3 for the ADM admits the following exact energy and helicity balances, which are discrete analogs the energy and helicity balances of the continuous ADM 1Energy∢2ξ‚΅β€–β€–ξ€·π·β„Žπ‘ξ€Έ1/2π‘£π‘€β„Žβ€–β€–2+𝛼2β€–β€–curlβ„Žξ€·π·β„Žπ‘ξ€Έ1/2π‘£π‘€β„Žβ€–β€–2ξ‚Ά+Ξ”π‘‘πœˆπ‘€βˆ’1𝑛=0ξ‚΅β€–β€–curlβ„Žξ€·π·β„Žπ‘ξ€Έ1/2π‘£β„Žπ‘›+1/2β€–β€–2+𝛼2β€–β€–Ξ”β„Žξ€·π·β„Žπ‘ξ€Έ1/2π‘£β„Žπ‘›+1/2β€–β€–2ξ‚Ά=12ξ‚΅β€–β€–ξ€·π·β„Žπ‘ξ€Έ1/2𝑣0β„Žβ€–β€–2+𝛼2β€–β€–curlβ„Žξ€·π·β„Žπ‘ξ€Έ1/2𝑣0β„Žβ€–β€–2ξ‚Ά+Ξ”π‘‘π‘€βˆ’1𝑛=0𝑓𝑑𝑛+1/2ξ€Έ,π·β„Žπ‘π‘£β„Žπ‘›+1/2ξ€Έ+𝛼2ξ€·curlβ„Žπ‘“ξ€·π‘‘π‘›+1/2ξ€Έ,curlβ„Žπ·β„Žπ‘π‘£β„Žπ‘›+1/2,𝐷(4.8)HelicityβˆΆβ„Žπ‘ξ€Έ1/2π‘£π‘€β„Ž,curlβ„Žξ€·π·β„Žπ‘ξ€Έ1/2π‘£π‘€β„Žξ‚+𝛼2ξ‚€curlβ„Žξ€·π·β„Žπ‘ξ€Έ1/2π‘£π‘€β„Ž,curl2β„Žξ€·π·β„Žπ‘ξ€Έ1/2π‘£π‘€β„Žξ‚+2πœˆΞ”π‘‘π‘€βˆ’1𝑛=0ξ‚€ξ‚€curlβ„Žξ€·π·β„Žπ‘ξ€Έ1/2π‘£β„Žπ‘›+1/2,curl2β„Žξ€·π·β„Žπ‘ξ€Έ1/2π‘£β„Žπ‘›+1/2+2πœˆπ›Ό2π‘€βˆ’1𝑛=0ξ‚€ξ‚€curl2β„Žξ€·π·β„Žπ‘ξ€Έ1/2π‘£β„Žπ‘›+1/2,curl3β„Žξ€·π·β„Žπ‘ξ€Έ1/2π‘£β„Žπ‘›+1/2=ξ‚€ξ€·π·ξ‚ξ‚β„Žπ‘ξ€Έ1/2𝑣0β„Ž,curlβ„Žξ€·π·β„Žπ‘ξ€Έ1/2𝑣0β„Žξ‚+𝛼2ξ‚€curlβ„Žξ€·π·β„Žπ‘ξ€Έ1/2𝑣0β„Ž,curl2β„Žξ€·π·β„Žπ‘ξ€Έ1/2𝑣0β„Žξ‚+2Ξ”π‘‘π‘€βˆ’1𝑛=0𝑓𝑑𝑛+1/2ξ€Έ,curlβ„Žπ·β„Žπ‘π‘£β„Žπ‘›+1/2ξ€Έ+𝛼2ξ€·curlβ„Žπ‘“ξ€·π‘‘π‘›+1/2ξ€Έ,curl2β„Žπ·β„Žπ‘π‘£β„Žπ‘›+1/2.ξ€Έξ€Έ(4.9)

Proof. For the energy balance, choose πœ’β„Ž=πΉβ„Žβˆ’1π·β„Žπ‘π‘£β„Žπ‘›+1/2. This vanishes the nonlinear term, since after the πΉβ„Ž and πΉβ„Žβˆ’1 cancellation of, we have that the cross of two vectors is perpendicular to each of them. The pressure term is also zero, since π·π‘π‘£β„Žπ‘›+1/2 is discretely-divergence free. This leaves 1ξ€·π‘£Ξ”π‘‘β„Žπ‘›+1βˆ’π‘£π‘›β„Ž,πΉβ„Žβˆ’1π·β„Žπ‘π‘£β„Žπ‘›+1/2ξ€Έξ€·Ξ”βˆ’πœˆβ„Žπ‘£β„Žπ‘›+1/2,πΉβ„Žβˆ’1π·β„Žπ‘π‘£β„Žπ‘›+1/2ξ€Έ=𝑓𝑑𝑛+1/2ξ€Έ,πΉβ„Žβˆ’1π·β„Žπ‘π‘£β„Žπ‘›+1/2ξ€Έ.(4.10) Now using the definition of πΉβ„Žβˆ’1 we get 1ξ€·π‘£Ξ”π‘‘β„Žπ‘›+1βˆ’π‘£π‘›β„Ž,π·β„Žπ‘π‘£β„Žπ‘›+1/2ξ€Έβˆ’π›Ό2ξ€·π‘£Ξ”π‘‘β„Žπ‘›+1βˆ’π‘£π‘›β„Ž,Ξ”β„Žπ·β„Žπ‘π‘£β„Žπ‘›+1/2ξ€Έξ€·Ξ”βˆ’πœˆβ„Žπ‘£β„Žπ‘›+1/2,π·β„Žπ‘π‘£β„Žπ‘›+1/2ξ€Έ+𝛼2πœˆξ€·Ξ”β„Žπ‘£β„Žπ‘›+1/2,Ξ”β„Žπ·β„Žπ‘π‘£β„Žπ‘›+1/2ξ€Έ=𝑓𝑑𝑛+1/2ξ€Έ,π·β„Žπ‘π‘£β„Žπ‘›+1/2ξ€Έβˆ’π›Ό2𝑓𝑑𝑛+1/2ξ€Έ,Ξ”β„Žπ·β„Žπ‘π‘£β„Žπ‘›+1/2ξ€Έ,(4.11) and with the definition of Ξ”β„Ž coupled with commutativity of operators, 1‖‖𝐷2Ξ”π‘‘β„Žπ‘ξ€Έ1/2π‘£β„Žπ‘›+1β€–β€–2+𝛼2β€–β€–curlβ„Žξ€·π·β„Žπ‘ξ€Έ1/2π‘£β„Žπ‘›+1β€–β€–2βˆ’β€–β€–ξ€·π·β„Žπ‘ξ€Έ1/2π‘£π‘›β„Žβ€–β€–2βˆ’π›Ό2β€–β€–curlβ„Žξ€·π·β„Žπ‘ξ€Έ1/2π‘£π‘›β„Žβ€–β€–2ξ‚Άβ€–β€–+𝜈curlβ„Žξ€·π·β„Žπ‘ξ€Έ1/2π‘£β„Žπ‘›+1/2β€–β€–2+𝛼2πœˆβ€–β€–Ξ”β„Žξ€·π·β„Žπ‘ξ€Έ1/2π‘£β„Žπ‘›+1/2β€–β€–2=𝑓𝑑𝑛+1/2ξ€Έ,π·β„Žπ‘π‘£β„Žπ‘›+1/2ξ€Έ+𝛼2ξ€·curlβ„Žπ‘“ξ€·π‘‘π‘›+1/2ξ€Έ,curlβ„Žπ·β„Žπ‘π‘£β„Žπ‘›+1/2ξ€Έ.(4.12) Multiplying through by Δ𝑑 and summing over timesteps now gives the stated energy balance.
For the helicity balance, choose πœ’β„Ž=πΉβ„Žβˆ’1curlβ„Žπ·β„Žπ‘π‘£β„Žπ‘›+1/2. This vanishes the nonlinear term since after the filter and its inverse cancellation of, the cross of two vectors is perpendicular to each of them. The pressure term also vanishes, since the curl is discretely-divergence free. This leaves
1ξ€·π‘£Ξ”π‘‘β„Žπ‘›+1βˆ’π‘£π‘›β„Ž,πΉβ„Žβˆ’1curlβ„Žπ·β„Žπ‘π‘£β„Žπ‘›+1/2ξ€Έξ€·Ξ”βˆ’πœˆβ„Žπ‘£β„Žπ‘›+1/2,πΉβ„Žβˆ’1curlβ„Žπ·β„Žπ‘π‘£β„Žπ‘›+1/2ξ€Έ=𝑓𝑑𝑛+1/2ξ€Έ,πΉβ„Žβˆ’1curlβ„Žπ·β„Žπ‘π‘£β„Žπ‘›+1/2ξ€Έ.(4.13) Next we decompose πΉβ„Žβˆ’1, similar to the case of the energy balance. After simplifying, we have 1𝐷2Ξ”π‘‘ξ‚€ξ‚€β„Žπ‘ξ€Έ1/2π‘£β„Žπ‘›+1,curlβ„Žξ€·π·β„Žπ‘ξ€Έ1/2π‘£β„Žπ‘›+1+𝛼2ξ‚€curlβ„Žξ€·π·β„Žπ‘ξ€Έ1/2π‘£β„Žπ‘›+1,curl2β„Žξ€·π·β„Žπ‘ξ€Έ1/2π‘£β„Žπ‘›+1βˆ’1𝐷2Ξ”π‘‘ξ‚€ξ‚€β„Žπ‘ξ€Έ1/2π‘£π‘›β„Ž,curlβ„Žξ€·π·β„Žπ‘ξ€Έ1/2π‘£π‘›β„Žξ‚+𝛼2ξ‚€curlβ„Žξ€·π·β„Žπ‘ξ€Έ1/2π‘£π‘›β„Ž,curl2β„Žξ€·π·β„Žπ‘ξ€Έ1/2π‘£π‘›β„Žξ‚€ξ‚ξ‚+𝜈curlβ„Žξ€·π·β„Žπ‘ξ€Έ1/2π‘£β„Žπ‘›+1/2,curl2β„Žξ€·π·β„Žπ‘ξ€Έ1/2π‘£β„Žπ‘›+1/2+𝛼2πœˆξ‚€curl2β„Žξ€·π·β„Žπ‘ξ€Έ1/2π‘£β„Žπ‘›+1/2,curl3β„Žξ€·π·β„Žπ‘ξ€Έ1/2π‘£β„Žπ‘›+1/2=𝑓𝑑𝑛+1/2ξ€Έ,curlβ„Žπ·β„Žπ‘π‘£β„Žπ‘›+1/2ξ€Έ+𝛼2ξ€·curlβ„Žπ‘“ξ€·π‘‘π‘›+1/2ξ€Έ,curl2β„Žπ·β„Žπ‘π‘£β„Žπ‘›+1/2ξ€Έ.(4.14) We complete the proof by multiplying through by 2Δ𝑑 and summing over timesteps.

5. Conclusions

We have developed enhanced-physics based finite element schemes for the NS-πœ” and Stolz-Adams ADM families of turbulence models that provides physically accurate treatment of discrete energy and helicity. Development of a discretely-divergence-free discrete curl operator, a discrete Laplacian defined in terms of the discrete curl, a discrete filter defined in terms of the discrete Laplacian, and a discrete deconvolution operator defined in terms of the discrete filter, all of which commute, are the key technical tools which allow the discrete schemes to accurately balance two fundamental physical quantities. Future directions of this work will be extensions to commonly used boundary conditions, and efficient implementations of these schemes.