New, scaleselective discretization for LES (and RANS)
Hi,
A proposed, new approach for discretizing the convection terms for the NavierStokes equations is given in the paper below based on implementation of the new SSD scheme in the OpenFOAM environment together with the RK4 scheme and projection method: Vuorinen V., Schlatter P., Boersma B., Larmi M., and Fuchs L., A ScaleSelective, Low Dissipative Discretization Scheme for the NavierStokes Equation, to appear in Computers and Fluids. The paper shows an "easy trick" to reduce numerical diffusion due to flux limiting or upwinding. Normally, in unstructured codes, the diffusive error is O(dx) and it makes the overall accuracy also O(dx) for turbulent flows. The SSD scheme is able to maintain O(dx^2) as the overall accuracy and it makes the diffusive error O(dx^3). This is shown in the paper. In case you want more info or want to implement the scheme in e.g. PISO solver you can ask me here for more info. Best, Ville %%%%%%%%%%%%%%%%%%%%%% Brief instructions for implementation at toplevel code are provided below. We use it mostly to stabilize implicit LES computations but it is generally a good candidate for explicit LES simulations as it is very nondissipative. About implementing SSDscheme (Vuorinen et al., to appear in Comp.&Fluids, 2012) to OpenFOAM toplevel code 1) A Laplacian highpass filter can be defined using the function Up = (Delta/(2*pi))^2 laplacian(U), where cutoff wavenumber is kcut = 2pi/Delta. Setting Delta = Nyqvist wave gives kcut = pi/dx which is recommended for SSD. We implement the filter here in conservative/self adjoint form (see the publication). Create fields for a filter (e.g. in createFields.H) surfaceScalarField DxDyDz ( IOobject ( "DxDyDz", runTime.timeName(), mesh, IOobject::NO_READ, IOobject::NO_WRITE ), mesh, dimensionSet(0,2,0,0,0,0,0) ); DxDyDz = Foam::pow(mesh.surfaceInterpolation::deltaCoeffs() ,2.0); scalar pi = 3.141592653589793; DxDyDz = (filterScope*filterScope)*DxDyDz/(4.0*pi*pi); %%%%%%%%%%%%%%%%%%%%%%%%%% %additionally implement highpass filtered fields for %the quantities you are solving e.g. velocity and set the BCs %in /0/ directory to be e.g. 0 for the fluctuation (I used 0 typically) % For example for the momentum volVectorField rhoUprime ( IOobject ( "rhoUprime", runTime.timeName(), mesh, IOobject::READ_IF_PRESENT, IOobject::AUTO_WRITE ), rho*U ); rhoUprime.correctBoundaryConditions(); 2) Now things are settled up and the simulation begins normally. You have to compute the fluctuation before solving the NSequations. For example, what I do in an explicit density based RK4 code: rhoprime = fvc::laplacian(DxDyDz,rho); rhoprime.correctBoundaryConditions(); rhoUprime = fvc::laplacian(DxDyDz,rhoU); rhoUprime.correctBoundaryConditions(); rhoEprime = fvc::laplacian(DxDyDz,rhoE); rhoEprime.correctBoundaryConditions(); pprime = fvc::laplacian(DxDyDz,p); pprime.correctBoundaryConditions(); 3) The SSD scheme is designed and validated for FULLY explicit time integration schemes. Normally the convection term in e.g. momentum eq would be written as fvc::div(phiv, rhoU). The convection term is rewritten as fvc::div(phiv, (rhoU+rhoUprime))  fvc::div(phiv, (rhoUprime)) Note that if you discretize both terms with the same scheme you are not doing anything different from fvc::div(phiv, rhoU) i.e. fvc::div(phiv, rhoU) = fvc::div(phiv, (rhoU+rhoUprime))  fvc::div(phiv, (rhoUprime)) The trick is to choose "linear" for the first term (which is smooth because rhoU+rhoUprime is smooth; the fluctuation is taken away; i.e. we low pass filter by writing rhoU+rhoUprime) and "blended" for the latter term which adds upwinding to it. That's it! 4) FOR IMPLICIT formulation we tried it as well with good experiences. Now you have probably terms such as fvm::div(phiv, rhoU) . Do not do anything to them BUT discretize them with "linear". You have to add an EXPLICIT source term to the RHS of the equations which is fvc::div(phiv, rhoUprimeLinear)  fvc::div(phiv, rhoUprimeBlended). Here rhoUprimeLinear and rhoUprimeBlended are the same field computed in the same way as rhoUprime i.e. with the HPfilter. You just need to copy that to two fields to separate the terms. As you might have guessed the "Linear" is discretized with "linear" and "Blended" with "blended". Perhaps you need to change the signs + >  but that works, we tested it. What effectively happens is that you get a sink of momentum on the RHS of your eqs which is essentially of form ~Uf''''(x) i.e. contains fourth derivatives. You see this easily in the case of linear convection eq from Taylor series expansion and the facts do not change overhere. More info in the paper. 
Many thanks, I just downloaded the paper for a reading

Thanks, very interesting. Filtering and numerics have been seldom taken into account simultaneously in the past, at least to my knowledge. I recall a similar attempt in the paper
Bose, Moin & You, "Gridindependent largeeddy simulation using explicit filtering", Physics of Fluids, 2010 although their approach is somewhat simpler. They explicitly filter the velocity field (using a commutative filter) and then just take care that the resolved scales are not affected by the numerical error (simply checking the modified wavenumber diagram, for a given scheme). Is there any chance that your new solver / scheme would be included within a future OpenFOAM release? 
SSD scheme and RK4
Hi,
Just a link to the final version of the article: http://dx.doi.org/10.1016/j.compfluid.2012.09.022 Thanks for the hint on the previous paper as well. The scheme is so easy to implement in the toplevel code that I mostly see it as two steps: 0) add the lines shown above which create a surfaceScalarField that is used in the hpfilter to take into account grid spacing variation 1) add a .H file where you evaluate the hpfiltered fields 2) add a source term to the governing equations. That's why I'd see it as a 5 minute modification of any existing solver to implement. From my behalf the RK4solver could well become a part of the OF distribution (extended) in the future. Would you have some instructions of the proper actions to take? Best, Ville 
Dear Ville,
thanks for your reply, I'll try to follow your instructions and modify an existing OF solver. Concerning the OFext distribution, unfortunately I have no idea about the procedure to follow, but I'm pretty sure that other members could help you, especially in the OpenFOAM subsection of the forum. Regards, Francesco 
I just had a look to the paper, what I see is that the definition of the fluctuation u' is actually the recovering only of the resolved frequency smoothed by a tophat like filter. I mean, in no way you can define a fluctuation velocity possessing wavenumbers outside the Nyquist frequency. You just use a deconvolutionlike approach to define the subgrid scale resolution.
It is correct? 
Hi,
I'm not doing deconvolution as we are limited to the grid scales i.e. k<=kNyq. So the argument is that we are not resolving well the waves at the grid scale ~2*dx. Thereby, it makes sense to split the resolved scales into smooth and fluctuating parts using the filter. The upwinding acts a dissipation model for the resolved part of the fluctuation (i.e. models the unresolved effects). Regarding what is really physical is a harder question. Here, the stabilization was carried out so that the diffusive error becomes of "higher order" i.e. acts as a high order filter which works technically and stabilizes enough. As the paper showed, the SSD scheme is quite nondissipative and it is even less dissipative than UW3. Any further comments? It is interesting to here the opinions here! Ville 
Quote:
Maybe I have not understood ... it seems to me that you are using the classical definition of a differential filter (tophat), that is you integrate the Taylor expansion around a volume that if (and only if) is centered around x leads to u_bar(x)=u(x) + (h^2/24) Lap(u(x)) +.... = (I + A^2 + ... )u(x) Then, you truncate and defines u'(x) =  (h^2/24) Lap(u(x)) This is a well known regularization procedure (Leray) that is a deconvolution procedure since you are formally inverting the truncated tophat operator (I + A^2 + ... )^1 = (IA^2 + ....). This is if I have well understood.... I have many references on such procedure 
Hi,
Seems like we are not talking about the same thing although I agree with you that a Laplacianterm is indeed one of the first terms in the expansion when inverting a filter.. I'd say that using a Laplacian filter in CFD and OpenFOAM is quite a standard way of separating scales. I prefer the formalism where the Laplacian filter is defined using (1/kcut)^2 prefactor because it relates the filter with the wave length very clearly. See the cited reference in Phys.Fluids by Stoltz et al. for a similar approach. One could use some other prefactor(cutoff) equally well and even argue a better choice for sure. Here it does not really matter too much because the truncation error is eventually controllable with the blending factor. So as long as one separates the scales in a meaningful way one can control the dissipation with alpha. Are you familiar with so called implicit LES? I have mostly in mind ILES when writing the paper. NSeqs are of form NS(u,p) = 0 and the normal LES filtered eqs of for NS(uf,pf) = SGSterms, where uf and pf represent the resolved numerical solution. In ILES the SGSterms are modeled using a truncation error of the convection term. Here I showed how the diffusive error can be made Odiff(dx) > Odiff(dx^3) using this numerical trick. Deconvolution approaches belong to their own category. I don't really prefer them because one tries to invert a filter using the van Cittert procedure and "gains" SGS terms represented on the grid scale. So what are the SGS terms represented on the grid scale is another question. Anyways, even in deconvolution methods damping terms are added and I think there is not very much difference with the SSD scheme. V 
Ok, I think that the issues could be seen in a unified way...
For example, see Eq.s (6) to (9) in http://onlinelibrary.wiley.com/doi/1...O;29/abstract as well as Eq.s (6) to (9) in http://onlinelibrary.wiley.com/doi/1...d.278/abstract I think that is possible to give a common framework Regards Filippo 
Thanks for the links, I'll have a look at them on monday as I don't have access to
them right now. Best, Ville 
All times are GMT 4. The time now is 03:12. 