CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums

Closing on wall functions - part 5: testing scripts

Register Blogs Members List Search Today's Posts Mark Forums Read

Rate this Entry

Closing on wall functions - part 5: testing scripts

Posted April 23, 2022 at 21:26 by sbaffini
Updated June 3, 2022 at 08:02 by sbaffini

I provide here a set of MATLAB scripts to test all the claims made in the first 4 parts.

The first group of scripts is actually made of functions, that you are not supposed to directly call or modify:
  • muskersp.m: returns \left(\frac{{s_{U,T}^i}^+}{{y^+}^{i+2}}\right), \left(\frac{{p^i}^+}{{y^+}^{i+2}}\right) and {q^i}^+ as shown here. It only works for N up to 0 (constant non equilibrium terms) EDIT: There is an apparently innocuous mistake in the limiting behavior of s, as it should not use the factor prprt for the thermal case. It seems innocuous as it only affects the first order term, which probably is already negligible, and only for the temperature case (so in the scripts it is only used for plotting purposes, and probably never at those low y^+).
  • standardsp.m: returns the same quantities but for the standard wall function presented here. It works for arbitrary N but the scripts below only test it up to 0 (constant non equilibrium terms)
  • iteryv.m: returns the value of y_T^+ needed by standardsp.m for Pr>Pr_t
  • numericalsp.m: a generic wall function that takes as input a function handle for the turbulent viscosity ratio and provides the required integrals by numerical integration. It works for arbitrary N (but, again, it is only tested for N up to 0). It only works for turbulent viscosities going to 0 at least as fast as {y^+}^2.
  • standard.m: the actual standard wall function adapted to the present framework for comparisons. Only works for equilibrium cases (N=-1)
  • spiter.m: the function that performs the iterations on the wall functions to find \tau_w following what is presented here
  • sa1d.m: solves the original ODE of the problem (actually a slightly more general one) using a second order cell centered finite volume discretization, the Spalart-Allmaras turbulence model (actually a generalized version that also works with the Musker profile) and the tridiagonal solver

The second group of scripts is the one actually making the tests and the comparisons:
  • check.m: for given flow conditions, it tests that a) the analytical musker profile in muskersp is consistent with the one obtained trough numerical integration with numericalsp, b) that the numerical musker profile (and, by point a above, also the analytical one) is consistent with the Spalart-Allmaras direct solution of the ODE with a properly modified fv1 function c) that the numericalsp routine (already tested in a and b above) is consistent with the non modified Spalart-Allmaras solution of the ODE when the relative turbulent viscosity profile is used. Note that the velocity wall function can only match the Spalart-Allmaras solution in the equilibrium case (dpdx=0), but in the same case the temperature equation can have non equilibrium terms and the solutions will match.
  • wfiter.m: for a given wall function (user selectable) and flow conditions plots the iterations to find \tau_w as produced by itersp.m
  • wfcompare.m: a straightforward comparison of multiple wall functions for given flow conditions
  • wfmap.m: for the musker profile (but easily adaptable to other profiles) it maps the solution \theta_p of the iterative procedure to find \tau_w in a certain range of Re_p and {F_p}^0 parameters.
  • apr.m: a comparison of the variations of the constant a in the Musker profile and y_T^+ in the standard profile as a function of Pr/Pr_t giving a sensation of how the variation of the Musker constant with Pr/Pr_t is, indeed, a sort of Jayatilleke term.

Few comments on the use of the scripts above and some claims made in previous posts.

In the second post of this series I mentioned that, when using the Spalart-Allmaras model, one could avoid iterations to find \tau_w by using y_p^+ = \frac{\rho \widetilde{\nu}_p}{\mu \kappa}. This is confirmed in check.m where the Spalart-Allmaras boundary conditions for the solution of the 1D ODE are obtained from the relation above using as input the y_p^+ value obtained by the underlying wall function.

In the same post it is also mentioned that one could precompute the wall function for a given range of parameters and later interpolate from it. Such precomputation is what the script wfmap.m does, and you can also see from the commented parts how it could be adapted to directly use the 1D ODE solution instead of the musker profile actually used by the script.

Still in the same post, I mentioned how, for non equilibrium cases, the iterative computation of \tau_w might become cumbersome as multiple solution branches might exist for non favourable pressure gradients. This can be easily visualized by using wfiter.m.

It has been mentioned here that y_T^+ in the standard profile can be obtained by few iterations of the Halley's method. The specific implementation is in iteryv.m.

In the last post I have mentioned that the van Driest damped mixing length is, indeed, mostly equivalent to the Musker profile for the incompressible formulation used here. This is clearly shown in wfcompare.m. While I don't directly provide the 1D ODE solution with the mixing length (just its wall function form trough numericalsp), I do for two different turbulent viscosity formulations trough the Spalart-Allmaras model in check.m and show that they correspond to the analogous formulation used as wall function.

All the scripts are in a zipped folder for convenience (as a maximum of 5 files can be uploaded here) and have been tested with MATLAB R2022a.
Attached Files
File Type: zip wallfunction.zip (14.2 KB, 644 views)
Posted in Uncategorized
Views 784 Comments 0 Edit Tags Email Blog Entry
« Prev     Main     Next »
Total Comments 0

Comments

 

All times are GMT -4. The time now is 09:06.