GitLab

• The results presented in the journal paper for the Reynolds stress tensor of Moser et al.'s (1999) smooth-wall plane channel flow with Re_\tau = 395 (Fig. 4) could not be reproduced with or without CFD despite all our efforts. (We are more than happy to be proven incorrect in our attempts.)
  • In the release notes of v1606+ (link), a high level of resemblance to the journal paper results as well as the benchmark results were reported.​
    • These results were unfortunately revealed to be involving a case setup error: the integral-length scale input was incorrect in those tests; more specifically, they were almost zero.
    • I think the reason why we had obtained fairly good results in the release notes is due to the nCellPerEddy term (which doesn't exist in the original journal paper) - that dominated the integral-length scale computation for relatively low-valued integral-length scale input sets.
  • Similar observations had been made by various CFD developers/users as well. For example, Kopper (2017) has stated in Implementation and Verification of a Synthetic Eddy Method (SEM) in the Eagle3d Compressible Flow Solver (SEM) in the Eagle3d Compressible Flow Solver for failed attempts of a DFSEM implementation into the solver Eagle3D:

While his published results are impressive and OpenFOAM Ltd announced the implementation of this method in the current OpenFOAM-plus release (OpenCFD Ltd, 2016), Poletto’s results could not be reproduced using the currently available versions of OpenFOAM-plus and CodeSaturne.

• See the following OpenFOAM results obtained from a replication of the journal-paper test case using the DNS dataset provided by Moser et al. (1999) / /data/MKM/chan395 (i.e. verificationAndValidation/turbulentInflow/oneCellThickPlaneChannel) (and compare it with the original Fig. 4). The input wall-normal profile of integral-length scales were computed by using the x-correlations (correlations) as well as k^1.5/epsilon (balances) provided by the DNS dataset (both approaches resulted similar scale values). The input Reynolds-stress tensor was obtained from chan395.reystress. The Reynolds-stress tensor results were sampled on the patch of the computational domain. Sampling has been performed with different sampling methods available in OpenFOAM to avoid any potential pitfalls existing in one of the sampling methods. Also note that the patch and entire computational domain were preserved in a single processor to prevent any parallelisation effects.

• Reynolds stress tensor results sampled on the patch:

• Reynolds stress tensor results sampled on the first-inlet cell centres:

• Mean streamwise velocity component sampled on the patch:

• Using the complete channel case also could not help to reproduce the Reynolds stress tensor results of the journal paper (Eqs. 13-15). See the test case and results in !454 (merged).

• One of the potential culprits we have identified is the coefficient C1 in Eq. 11 being used in Eq. 10, which is at the heart of the DFSEM.

  • C1 in Eq. 11 has units of m^(1/2), yet it must be dimensionless for Eqs. 9-10 to be dimensionally consistent.
    • Checked also the thesis, whereat C1 expression is identical to that in the journal paper.
  • This sounds a minor issue (which should easily have been picked by the reviewers of the journal paper, I believe), but Eq. 11's Vo term affects the instantaneous velocity field as well as the Reynolds stress tensor - refer to #1744 and #1004.
  • Changing the sqrt to cbrt, like suggested by Michael, can alleviate velocity-field issues, but some local tests indicated that the change could cause scaling issues for the Reynolds stress tensor this time.
  • No explanation was given on how C1 was derived.

• Another culprit might be the C2 in Eq. 12 which should also be dimensionless. No detailed derivation or elaboration was given on C2 in the journal paper or in the thesis; hampering the reverse engineering.

• One problem is related to the average length scale expressed in Eq. 14. Assuming \Delta {x,y,z} are the local mesh size, this equation returns the local mesh size as the average length scale of the flow for any input length scales larger than the local mesh size. For example, L=0.1, \kappa \delta = 0.41, and \Delta {x,y,z} = 0.01 - this expression returns the length scale as \sigma=0.01 rather than the input L.

• Another problem we observed in the DFSEM is that the reproduction of input Reynolds stress tensors on a patch is severely affected by input choices for integral-length scales without any deterministic relation between the two.

• Also, we have observed that, due to the (theoretical) limiters being applied on input integral-length scales, the scales generated on a patch differ from the user input. For example, if the user measures autocorrelation function on the patch with Taylor's hypothesis and calculates scales from the autocorrelation, the (physical) input scales and calculated scales will likely be different. For this reason, this situation leaves us with a method which could not reproduce user inputs of the Reynolds stress tensor and integral length scales on a given computational domain patch (i.e. no CFD effects) and inside the domain.

• DFSEM method does not allow any input length scales larger than kappa*delta where kappa is the von Karman constant and delta is the characteristic length of the domain. The rationale of this is the premise that a length scale of an eddy cannot be larger than the domain length.

  • However, this premise is only valid for lateral directions in a wall-bounded flow and is not valid for freestream flows. In a wall-bounded flow, the streamwise-direction components can be larger than kappa*delta. This was even showed in one of the figures in the original work - the linearly-increasing line is computed from Moser et al.'s (1999) DNS computation:

• The method is reportedly divergence-free. However, to what extent the solenoidal field persists inside a CFD domain is unquantified. Also, it is not clear whether the pressure fluctuations were reduced due to the divergence-freeness of the method or the ad-hoc application of the mass flux correction.

• Various solution suggestions:

  • Discard DFSEM, or limit its usage for isotropic turbulence cases.
  • Replace DFSEM with the other SEM variants being developed by Uni. of Manchester, since to our observation, the active developments and usage on DFSEM in the originator institute has been ceased as well in favour of other SEM variants being developed thereat.