dT in wall function for momentum when iwalltemp/=2

[tomgrylls]

The wall function for momentum uses the temperature difference between the wall and the adjacent cell, dT. However, under neutral conditions or under non-neutral conditions when iwalltemp/=2 (zero flux or constant flux BCs used) this term is still calculated even though the wall temperature is not being modelled. This means that the wall function will have a stability term based off a temperature difference that is not being modelled.

This will occur:

1. in neutral conditions if there is a difference between the arbitrarily defined air temperature (typically thl=288) and the values assigned in Tfacinit.f90. This can be avoided by ensuring that these are equal and therefore dT=0 by definition (currently done in preprocessing but maybe this should be automated as it should not need to be considered when modelling neutral conditions). Perhaps have a switch for if (ltempeq) where Tfacinit.f90 is read.
2. When iwalltemp=1 the air temperature is non-uniform and dT will be non-zero. Potential solution would be to define a switch within each case in wfuno.f90 where dT=0 if iwalltemp=1. Or may be easier to just edit function unom and we just need Ribl0=0 under these conditions. @ivosuter - thoughts?

The consequence of this bug is that under these specific conditions we can expect changes to the applied wall shear stress due this stability term.

[ivosuter]

1. If you have a switch 'lneutral', it's easy to set the wall temperature equal to the fluid temperature
2. there is no real solution: maybe set the Ri or directly Cd according to the stability, i.e. Cd=0.001 for a T-flux towards the wall, Cd=0.01for T-flux from the wall; or you do iwallmom=1 and also fix the momentum flux?

Ivo

On 22 Jul 2019 19:38, tomgrylls [email protected] wrote:

The wall function for momentum uses the temperature difference between the wall and the adjacent cell, dT. However, under neutral conditions or under non-neutral conditions when iwalltemp/=2 (zero flux or constant flux BCs used) this term is still calculated even though the wall temperature is not being modelled. This means that the wall function will have a stability term based off a temperature difference that is not being modelled.

This will occur:

1. in neutral conditions if there is a difference between the arbitrarily defined air temperature (typically thl=288) and the values assigned in Tfacinit.f90. This can be avoided by ensuring that these are equal and therefore dT=0 by definition (currently done in preprocessing but maybe this should be automated as it should not need to be considered when modelling neutral conditions). Perhaps have a switch for if (ltempeq) where Tfacinit.f90 is read.
2. When iwalltemp=1 the air temperature is non-uniform and dT will be non-zero. Potential solution would be to define a switch within each case in wfuno.f90 where dT=0 if iwalltemp=1. Or may be easier to just edit function unom and we just need Riblo0=0 under these conditions. @ivosuterhttps://github.com/ivosuter - thoughts?

The consequence of this bug is that under these specific conditions we can expect changes to the applied wall shear stress due this stability term.

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHubhttps://github.com/uDALES/u-dales/issues/23?email_source=notifications&email_token=AMQO4V67O2QOFLIMEKT247LQAXWAXA5CNFSM4IF23KC2YY3PNVWWK3TUL52HS4DFUVEXG43VMWVGG33NNVSW45C7NFSM4HAWGUUA, or mute the threadhttps://github.com/notifications/unsubscribe-auth/AMQO4V3BLOMMGN7IWSN7HFLQAXWAXANCNFSM4IF23KCQ.

[tomgrylls]

@ivosuter are those proposed values of Cd estimates or is there some reference for setting those?

I don't think that we want to set iwallmom==1. The base/ easiest solution would be to neglect the thermal effects on the wall function for momentum completely such that the momentum flux is defined as if it was neutral - having a look in wfuno.f90 would Ribl0 = 0 suffice to do this? Or do you think this is too great a simplification? In these kind of idealised simulations it is likely that you have some facets with zero flux and others with a uniform flux.

Or I am trying to figure out if this wall function could be adapted seeing as we actually already "know" the heat flux instead of the temperature difference. In a way it should be easier with u*\theta* already known? E.g. if a flux Richardsons number was used. Looking at @ivosuter's thesis and Uno (1995) I think we would need some way to estimate the bulk Richardson's number from the heat flux and fluid velocity at the adjacent cell. This may be possible using the implicit relationship between the Obukhov length and the bulk Richardson's number but is not trivial.

[ivosuter]

If you want it to be 'neutral' skip the whole wfuno routine and set Cd directly based on wall roughness length? From subroutine unom it seems that in this case Cd = Ctm = fkar2/(logdz2), where logdz2 = log(z/z_0)^2 and fkar2 = \kappa^2

I don't think I would try to find the Ri based on the heat flux and then use that Ri to calculate the momentum flux

[tomgrylls]

I see Figure 3.3 now in your thesis and it is tempting to assign a general value for Cd for unstable and stable conditions. However it would mean that Cd is not a function of the roughness length. At least by assuming it is neutral, Cd is a function of z0 (albeit not the right magnitude due to neglecting the thermal effects).

In terms of adding this into the code, I can either:

1. add an additional switch in modibm so that everywhere that iwallmom==2, if (iwalltemp==1) then we apply the simplified neutral wall function instead of calling wf_uno. But would I have to set-up all of the different cases as is already done in wf_uno.f90?
2. I can have a switch within the momentum flux subroutine in wf_uno that just sets dT=0 if (iwalltemp==1). For me the latter would be easier to implement - but is there any reason I should not just do it like this?

[ivosuter]

see below

[ivosuter]

hmm. It seems it's still quite a bit of coding. So, I would recommend to make a wfmneutral: copy wf_uno.f90, rename delete cases 'x2' and function 'unoh' & 'unom' delete all lines in cases 'x1': logdzh = LOG(delta/z0h) logzh = LOG(z0/z0h) dT = ((Tcell(i, j, k) + Tcell(i, j, k - 1)) - (Twall + Twall))0.5 Ribl0 = gravdeltadT2/((Twall + Twall)*utangInt) !Eq. 6, guess initial Ri change 'ctm = ...' to: ctm = fkar2/(logdz**2) remove unused in/output done! call the same way as wfuno: else if (iwallmom == 3) then do n = 1, nxwall k = block(ixwall(n), 8) !west side call wfmneutral(ih, jh, kh, vp, wp, momfluxb, v0, w0, facz0(k), ixwall(n), 1, 11) k = block(ixwall(n), 9) !east side call wfmneutral(ih, jh, kh, vp, wp, momfluxb, v0, w0, facz0(k), ixwall(n), 1, 21) end do end if

[tomgrylls]

Okay thanks. Agreed that makes sense. I will create a new branch, implement those changes and then create a pull request back into the master that you can review.

We may also want some check at startup e.g. if ((ltempeq=.false.) .or. (iwalltemp==1)) then iwallmom=3.

[tomgrylls]

Still have unom called in cases x1 on lines starting dummy = ?

CORRECTION: not for all lines but for some. I assume on these lines I can replace unom(logdz, logdzh, logzh, sqdz, utangInt, dT, Ribl0, fkar2) with fkar2/(logdz**2)?

[ivosuter]

yes usually its 'ctm=unom(...)', followed by 'dummy=u^2ctm'; sometimes its directly 'dummy=u^2unom()'. it's pretty much the same thing.

wherever there is a unom() you can just replace it with 'fkar2/(logdz**2)'

"We may also want some check at startup e.g. if ((ltempeq=.false.) .or. (iwalltemp==1)) then iwallmom=3." I agree

[bss116]

We believe this is essentially solved and just need to confirm this. Moving to 0.2.0 milestone.

[dmey]

@bss116 do you know what's happening with this one?

[bss116]

I think this is solved, wfmneutral.F90 got implemented with #24 and the switch in modstartup.F90 exists. @tomgrylls can you confirm please?

[dmey]

@bss116 @samoliverowens this seems to be the last bug to check --would you be able to create an example/test case to test for this?

[samoliverowens]

@bss116 @samoliverowens this seems to be the last bug to check --would you be able to create an example/test case to test for this?

I think this is tested in 002, because wfmneutral will be used for neutral simulations and simulations not using the wall functions for temperature:

https://github.com/uDALES/u-dales/blob/9c5e8e308ae811cabc15e1f90d1a28cccb0f2970/src/modstartup.f90#L566-L570

[dmey]

Do we need further checks for this or can we close it?

[ivosuter]

I said I would have a last look at it, but that we believe it's done. So gonna check wfmneutral now

[ivosuter]

Hmm it appears to me there are indeed some mistakes in both, wfm_neutral and wf_uno in the edge cases. On the staggered grid the tangential velocities on a surface are at the vertices of the gridcell/surface. The edge cases (the edge of the building) therefore only experieces half a gridcell of solid -> we simply assume the flux is also halfed. Can someone double check that the indeces make sense and were not forgotten: (i'll post a list in a few minutes)

[ivosuter]

ui-1 ui ui+1 ui+2 o-----o-----o-----o |...................................| |...................................| |...................................| ui-1 and ui+2 only experience half a cell

[ivosuter]

wfm_neutral: do we need w bottom cases at all? since blocks are never floating the bottom is always the floor? If we do, shouldn't it be (on lines 142,251,352 and 461) k = block(n, 5)

On line 518/519 a 'i = block(n, 2) + 1' is missing

wfm_uno: same story should it be 'k = block(n, 5)' on lines 357,497,626,763

On line 835/836 a 'i = block(n, 2) + 1' is missing

[dmey]

@samoliverowens should we ask @tomgrylls to look into this? @ivosuter made some comments some time ago about this issue but I am not the right person to look into this.

[samoliverowens]

@samoliverowens should we ask @tomgrylls to look into this? @ivosuter made some comments some time ago about this issue but I am not the right person to look into this.

I'll have a look too, but yes @tomgrylls will be better with this.

[tomgrylls]

wfm_neutral: do we need w bottom cases at all? since blocks are never floating the bottom is always the floor? If we do, shouldn't it be (on lines 142,251,352 and 461) k = block(n, 5)

Agreed I think the best option is to edit these to k = block(n, 5). This loop will be used in situations where blocks start from non-zero values. We do not have them floating but in complex morphologies where blocks are positioned above other blocks this will occur.

On line 518/519 a 'i = block(n, 2) + 1' is missing

Good spot - yes agreed we should add this line.

wfm_uno: same story should it be 'k = block(n, 5)' on lines 357,497,626,763

On line 835/836 a 'i = block(n, 2) + 1' is missing

Both of these are equivalent to the points above. And therefore again I am in agreement of changing them. As shown from the discussion above wfmneutral is directly based upon wf_uno and therefore it makes sense that these are consistent across both.

I will create a pull request in the next couple of days and set @ivosuter as reviewer in order to have a final check on the above.