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

A generalized thermal/dynamic wall function: Part 1

Register Blogs Community New Posts Updated Threads Search

Rate this Entry

A generalized thermal/dynamic wall function: Part 1

Posted October 14, 2016 at 13:27 by sbaffini
Updated December 21, 2016 at 10:07 by sbaffini

In a previous post i wrote about an extension of the Reichardt law of the wall to pressure gradient effects. That was derived by assuming the Reichardt profile as a solution for the case without pressure gradient and using integration by parts. In particular, given the the Reichardt function (k is the Von Karman constant):

f^+\left(y^+\right) = \frac{1}{k}\log\left(1+ky^+\right) +A\left(1-e^{-\frac{y^+}{B}}-\frac{y^+}{B}e^{-\frac{y^+}{C}}\right)

with:
A=\frac{1}{k}\log\left(\frac{D}{k}\right)

B=11

C=\frac{2}{\frac{1}{B}+\frac{Bk}{A}}

D=11.25

the pressure gradient sensitized formula for the velocity is given by:

u^+\left(y^+\right)=f^+\left(y^+\right)\left(1+F_u^+y^+\right)-F_u^+ \int{f^+dy^+}

where F_u^+=\left(\mu dp/dx\right)/\left(\rho^2 u_\tau^3\right) is the nondimensional pressure gradient.

However, while useful for its purpose, this formulation turns out to have at least two drawbacks, both of which depending, somehow, on the fact that Reichardt didn't actually solved the underlying problem exactly. The first one is that the base profile (i.e., even without the pressure gradient effect) does not even reproduce the correct profile derivatives at the wall. Hence, as a main consequence, the formula might not be suitable as toy model for certain theoretical purposes. The second drawback, which is also another consequence of the same fact, is that it can be used as thermal wall function only for Pr/Pr_t \le 1, as higher values will tend to amplify the part of the profile which is uncorrect.

To overcome such issues, i follow here the same steps as before but with a different formulation for \mu_t/\mu, using one which can be integrated exactly. We start from the temperature equation, which is more general:

\frac{d}{dy}\left[C_p\left(\frac{\mu}{Pr}+\frac{\mu_t}{Pr_t}\right)\frac{dT}{dy}\right]=F_T

where C_p (constant pressure specific heat), \mu (dynamic viscosity), Pr (Prandtl number), Pr_t (turbulent Prandtl number) and F_T (defined in the following) are all constant. The fact that this equation correctly represents a large set of velocity/temperature viscous boundary conditions should be self-evident. For example, F_T=\mu_t=0 gives the classical linear relationship for the temperature (which becomes the one for the velocity when C_p=Pr=1). Also, C_p=Pr=Pr_t=1 with F_T=0 leads to the classical wall-function formulation (for certain forms of \mu_t/\mu). According to the specific equation form, F_T can have different meanings. For example, in the case of the velocity, it can include the wall-parallel pressure gradient (note that y is assumed to be the wall-normal coordinate, x the wall parallel one), the buoyancy term, etc. In general, it can include any of the terms that have been neglected in the LHS for the given equation. The only assumption here is that they are assumed constant along y (as will be shown in the following, it turns out that this can be actually relaxed).

The first step to proceed is integrating the equation above once. This produces (using \mu_t\left(y=0\right)=0):

C_p\left(\frac{\mu}{Pr}+\frac{\mu_t}{Pr_t}\right)\frac{dT}{dy}=F_T y +q_w

where q_w=\left(C_p \mu/Pr\right)\left(dT/dy\right)_w (the w underscore meaning that it is evaluated at the wall, i.e., for y=0). Multiplying by Pr/q_w and rearranging, we obtain:

\left[1+\left(\frac{Pr}{Pr_t}\right)\left(\frac{\mu_t}{\mu}\right)\right]\frac{dT^+}{dy^+}=Pr\left(1+F_T^+y^+\right)

where:

u_\tau = \sqrt{\frac{\mu}{\rho}\left(\frac{du}{dy}\right)_w}

y^+=\frac{\rho y u_\tau}{\mu}

F_T^+ = \frac{\mu F_T}{\rho^2 u_\tau^2 C_p T_\tau}

T_\tau = \frac{q_w}{\rho C_p u_\tau}

T^+ = \frac{T-T_w}{T_\tau}

The same equation, interpreted for the velocity, would read:

\left[1+\left(\frac{\mu_t}{\mu}\right)\right]\frac{du^+}{dy^+}=1+F_u^+y^+

where u^+=u/u_\tau. This confirms that the temperature solution can also be used for the velocity by a proper interpretation of the terms. Going back to the temperature equation, we then end up with the following equation for the nondimensional temperature derivative:

\frac{dT^+}{dy^+}=\frac{Pr\left(1+F_T^+y^+\right)}{\left[1+\left(\frac{Pr}{Pr_t}\right)\left(\frac{\mu_t}{\mu}\right)\right]}

Then, by a further integration:

T^+\left(y^+\right)=Pr \left\{\int_0^{y+}{\frac{1}{\left[1+\left(\frac{Pr}{Pr_t}\right)\left(\frac{\mu_t}{\mu}\right)\right]}dz^+}+\int_0^{y+}{\frac{F_T^+z^+}{\left[1+\left(\frac{Pr}{Pr_t}\right)\left(\frac{\mu_t}{\mu}\right)\right]}dz^+}\right\}

Note that i temporarily suspended the hypothesis that F_T^+ is constant. Now, let us also assume that we know the first integral, that is:

\int_0^{y+}{\frac{1}{\left[1+\left(\frac{Pr}{Pr_t}\right)\left(\frac{\mu_t}{\mu}\right)\right]}dz^+}=f^+\left(y^+\right)

is a known function. Then, the solution becomes:

T^+\left(y^+\right)=Pr \left\{f^+\left(y^+\right)+\int_0^{y+}{F_T^+z^+\frac{df^+}{dz^+}dz^+}\right\}

And after integration by parts:

T^+\left(y^+\right)=Pr f^+\left(y^+\right)\left[1+F_T^+\left(y^+\right)y^+\right]-Pr\int_0^{y+}{f^+\left(F_T^++z^+\frac{dF_T^+}{dz^+}\right)dz^+}

For the sake of conciseness, let us also assume:

\int_0^{y+}{f^+\left(F_T^++z^+\frac{dF_T^+}{dz^+}\right)dz^+}=g^+\left(y^+\right)

so that the full solution is:

T^+\left(y^+\right)=Pr f^+\left(y^+\right)\left[1+F_T^+\left(y^+\right)y^+\right]-Pr g^+\left(y^+\right)

While we have just moved the original problem to that of solving the two integrals for f^+ and g^+, those will actually turn out to be easier to solve. Note also that the parenthesis in g^+ is a polynomial of the same order of F_T^+ (assuming that it is approximated with a certain degree polynomial). Particularly relevant is the case with F_T^+ constant (as is typically assumed). However, note that ,for the f^+ that we are going to derive below, a linear approximation is also integrable, as well as any other polynomial of higher degree. Obviously, things rapidly become cumbersome, with the degree of the polynomial, so for the present purposes it is sufficient to go back to the case F_T^+ constant. In this case g^+ is simply the integral of f^+ times F_T^+. Note that Popovac and Hanjalic have shown (for the velocity) that this is actually relevant for several practical cases. Nonetheless, is worth remembering that this is not a costriction of the present approach (probably it wasn't neither when f^+ was based on the Reichardt profile).

So, at this point we are left with determining f^+. To this end, we need a model for \mu_t. Note that this is not different from the approach used by Van Driest to determine his velocity profile. It's just that the mixing length hypothesis he used led to a non integrable function. The same also happened to Reichardt, who then started modificating his function until something integrable came out.

To overcome this limit, here we abandon those formulations and use instead the one from Musker:

\frac{\mu_t}{\mu}=\frac{\left(ky^+\right)^3}{\left(ky^+\right)^2+\left(ka\right)^3-\left(ka\right)^2}

with a being, in its original formulation, a constant (see below). It is worth mentioning that this functional form respects both the ky^+ behavior for y^+ -> \infty and the \left(y^+\right)^3 behavior for y^+ -> 0.

Note also that the Spalart-Allmaras model has a simlar behavior near the wall (see Allmaras et al., ICCFD 2012):

\frac{\mu_t}{\mu}=\frac{\left(ky^+\right)^4}{\left(ky^+\right)^3+c_{v1}^3}

which then fails to recover the correct \left(y^+\right)^3 behavior for y^+ -> 0 (However, Spalart and Allmaras explicitly mention in their original 1994 work that thay didn't notice any meaningful difference deriving from using this specific form instead of a correct one, say, the mixing length). This just to say that the Musker form is not only meaningful, but also representative of the near wall behavior of the Spalart-Allmaras model.

To solve for f^+ with the Musker turbulent viscosity, we first consider the case Pr/Pr_t=1, as this already has a known solution, provided in its most general form by Monkewitz et al. (POF 20, 105102, 2008):

u^+=\frac{1}{k}\ln\left(\frac{y^++a}{a}\right)+
\frac{\alpha}{a+4\alpha}\left\{\left(a-4\alpha\right)\ln\left[\frac{a\left[\left(y^+-\alpha\right)^2+\beta^2\right]}{2\alpha\left(y^++a\right)^2}\right]+\gamma\left[\arctan\left(\frac{y^+-\alpha}{\beta}\right)+\arctan\left(\frac{\alpha}{\beta}\right)\right]\right\}

where:

2\alpha=a-\frac{1}{k}

\beta^2=2a\alpha-\alpha^2

\gamma = \frac{2\alpha\left(5a-4\alpha\right)}{\beta}

The constant a is found by calibration with the turbulent velocity profile. Monkewitz et al. use a=10.306. In contrast, i calibrated the constant using the mixing length model with a constant A^+=19, resulting in a=11.489.

Now, the whole point of the procedure is: how to extend such formulation to the case Pr/Pr_t \ne 1? See part 2
Views 2348 Comments 0 Edit Tags Email Blog Entry
« Prev     Main     Next »
Total Comments 0

Comments

 

All times are GMT -4. The time now is 16:40.