CFD Online Logo CFD Online URL
Home > Forums

Closing on wall functions - part 8: coupled/thin wall boundary conditions

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

Rate this Entry

Closing on wall functions - part 8: coupled/thin wall boundary conditions

Posted April 27, 2022 at 08:37 by sbaffini
Updated May 16, 2022 at 11:51 by sbaffini

It might happen that wall functions for temperature or scalars are needed when the assigned boundary condition is not simply the assigned flux or temperature/scalar, but rather a more complex one. One example is when the wall is a coupled one, with either a fluid or a solid on the other side, or it has some thickness, possibly with source terms in it and, say, radiation or convective boundary conditions assigned on the other side, or maybe some other combination of the above.

In all such cases, the formulas presented before are still valid but need a slight rearrangement in order to fit the new conditions. Nothing of this is really new, the concpet is as old as the book of Patankar (probably older) and this is just one of the latest additions. Also, major commercial CFD codes have been offering this for decades. Still, I am not aware of full formulas available in the more general case. Everything I write here for the temperature just straightforwardly applies to other scalars with similar equations. Of course, as per the original wall function ODE, the assumption is that of steady state. That is, the solved problem might or not be steady, but the boundary conditions are, in fact, derived by solving a steady state problem (this was, as a matter of fact, true also for the wall function ODE, even if the unsteady term was considered among the non equilibrium ones).

We start by noting that the wall heat flux formula provided here can actually be rewritten as follows (also taking into account that the original derivation used a wrong sign for ease of exposition):

q_w = -\frac{\left\{\frac{\left[T\left(y\right) - T_w\right]}{y} \left(\frac{\mu C_p}{Pr}\right) - y\sum_{i=0}^{N}\frac{F_T^i}{i+1}\left(\frac{y}{y_p}\right)^i\left(\frac{{s_T^i}^+}{{y^+}^{i+2}}\right) \right\}}{\left(\frac{{s_T^{-1}}^+}{y^+}\right)} = -\frac{\left(T_p - T_w\right)}{y_p} k_p C_{WF} + \dot{Q}_p

Where C_{WF}=1/\left(\frac{{s_T^{-1}}^+}{y^+}\right) and it is recognized that, being independent, at the first order, from the flow conditions and being an integral that grows with wall distance, the non equilibrium terms are, indeed, just an explicit source term \dot{Q}_p for the near wall cell. In practice, the source term is also assimilable to a non orthogonal correction thus, in the following, we will simply consider the point p, with temperature T_p and distance from the wall y_p to be along the normal to the wall. A similar reasoning can also be done for the viscous dissipation that, as presented here, is independent from the temperature distribution and can be absorbed by the same source term as well.

We want to extend the formula above to the case where neither T_w nor q_w are actually given. Also, we assume that T_w=T_n and q_w=q_n, that is, they are only known as the values at the extreme of an n layers thin wall. So we have T_n, T_{n-1}, ..., T_1 and T_0, and the same for q. They are are the temperatures and fluxes at the intefraces between the n layers of the thin wall. Each one of these n layers will have its own thickness \Delta x_i, thermal conductivity k_i and possibly a source term \dot{Q}_i. Finally, we want to consider two possible boundary conditions, either T_0 directly given or q_0 given as:

q_0 = h_{\infty}\left(T_{\infty}-T_0\right) + \epsilon_R \sigma \left(T_R^4-T_0^4\right) = h_{\infty}\left(T_{\infty}-T_0\right) + \widetilde{h}_R\left(T_R-T_0\right)

with \widetilde{h}_R = \epsilon_R \sigma \left(T_R+T_0\right)\left(T_R^2+T_0^2\right) being non-linearly dependent from T_0 and where we assumed that the q_i are positive if entering the domain.

In order to formalize a solution, we need to complement the fundamental conservation statement that holds for a layer in the thin wall q_i=q_{i-1}+\dot{Q}_i \Delta x_i with a similar relation for the temperature jump across the layer, T_i-T_{i-1}. Considering that, in our model, the heat conduction equation that holds in each layer has the form:

\frac{d^2T}{dx^2} + \frac{\dot{Q}}{k} = 0

solving it with proper boundary conditions one easily obtains that:

T_i-T_{i-1} = -\frac{\Delta x_i}{k_i}\left(q_{i-1}+\frac{\dot{Q}_i \Delta x_i}{2}\right)

The two jump relations for single layers above can then be used to obtain jump relations across the whole thin wall of n layers. For the flux it is just, again, a simple conservation statement:

q_n = q_0 + \sum_{i=1}^n \dot{Q}_i \Delta x_i

For the temperatures one obtains:

T_n - T_0 = \sum_{i=1}^n \left(T_i -T_{i-1}\right) = -q_0 \sum_{i=1}^{n}\frac{\Delta x_i}{k_i} - \sum_{i=1}^{n} \frac{\Delta x_i}{k_i} \left(\frac{\dot{Q}_i \Delta x_i}{2} + \sum_{j=1}^{i-1}\dot{Q}_j \Delta x_j\right)


h_p = \frac{k_p C_{WF}}{y_p}

R_T = \sum_{i=1}^{n}\frac{\Delta x_i}{k_i}

\alpha = \sum_{i=1}^{n} \frac{\Delta x_i}{k_i} \left(\frac{\dot{Q}_i \Delta x_i}{2} + \sum_{j=1}^{i-1}\dot{Q}_j \Delta x_j\right)

\beta = \sum_{i=1}^n \dot{Q}_i \Delta x_i

where the latter 3 can all be pre-computed and stored, then it is a matter of simple manipulation (yet a quite long one, which is omitted here), to show that our initial formula can now be generally expressed as follows:

q_w = \left(T_w - T_p\right)h_p + \dot{Q}_p = \frac{\left(T^*-T_p\right)+\left(R^*\beta-\alpha\right)}{R^*+\frac{1}{h_p}} + \dot{Q}_p

where T^* and R^* depend from the specific boundary condition in use. More specifically, for T_0 directly assigned, one has T^*=T_0 and R^*=R_T. For the general convection/radiation boundary condition (or just one of the two, by zeroing the coefficient of the other) one has:

R^* = R_T + \frac{1}{h_{\infty}+\widetilde{h}_R}

T^* = \frac{h_{\infty}T_{\infty}+\widetilde{h}_R T_R}{h_{\infty}+\widetilde{h}_R}

Finally, the coupled case is simply obtained by using \widetilde{h}_R = 0, h_{\infty} = h_{pc} and T_{\infty}=T_{pc} where the subscript pc refers to the quantities taken from the coupled side.

The corresponding value of T_w is instead given by:

T_w = \frac{T^*+R^*\left(\beta+h_pT_p\right)-\alpha}{1+h_pR^*}

The last step missing is the determination of \widetilde{h}_R, which depends from T_0. A relation for the latter can be obtained, but itself, of course, will depend from \widetilde{h}_R. This relation (whose derivation is again simple but long, so it is omitted here) can then be used in tandem with the one for \widetilde{h}_R iteratively:

\widetilde{h}_R = \epsilon_R \sigma \left(T_R+T_0\right)\left(T_R^2+T_0^2\right)

T_0 = \frac{\left(h_{\infty}T_{\infty}+\widetilde{h}_R T_R\right)\left(R_T+\frac{1}{h_p}\right)+T_p+\alpha+\frac{\beta}{h_p}}{\left(h_{\infty}+\widetilde{h}_R\right)\left(R_T+\frac{1}{h_p}\right)+1}

I have found that, starting from T_0=T_p no more than 10 iterations are necessary to converge on \widetilde{h}_R and T_0.

One thing which is worth highlighting here in the more general context of wall functions is that, as a matter of fact, non equilibrium/viscous dissipation terms (from both sides of the thin wall, if a coupled bc is used) do not directly enter the modifications presented above. That is, they still appear as in the original wall function formulation. This might very easily go unnoticed if a thermal wall function is just presented as a single relation for T^+ without distinguishing the roles of the terms.

Another thing worth noting is that, if one decided to solve the original ODE as it was (say, with a tridiagonal algorithm as in one of the scripts provided here), instead of directly integrating it as done here, it would have been impossible to let the exact form above emerge, leading to 3 major consequences: 1) inability to separate the non-equilibrium/viscous dissipation part from the rest (that is, in the ODE solution one only has the wall values, not their dependence), 2) either the need to solve the 1D problem also in the thin wall and/or coupled side or the need to iterate the solution on both sides exchanging wall values at each iteration and 3) which is a consequences of 1 and 2, the impossibility to set up the problem for an implicit implementation in the coupled case, which might have major consequences on convergence in certain cases.

Finally, I want to mention that all the above developments are also relevant for the Musker-modified Spalart-Allmaras model for which, given the proper conditions (equilibrium for the velocity, steady state in general), they represent a full analytical solution, which now is thus extended to the present more general boundary conditions.
Posted in Uncategorized
Views 716 Comments 0 Edit Tags Email Blog Entry
« Prev     Main     Next »
Total Comments 0



All times are GMT -4. The time now is 01:31.