Fix wall treatment for mixing length model in TELEMAC-2D
Clemens Dorfmann (Flow Engineering) has suggested a revision for the mixing length turbulence model in TELEMAC-2D. See some details of a few e-mails below.
1st e-mail: It had some minor mistakes in the wall treatment of the viscosity to be computed since it didn't account for liquid boundaries. Hence, it treated liquid boundary nodes as solid boundary nodes which could produce wrong viscosity values at the liquid boundary nodes.
I tested and implemented the appraoch by Jia et al., see the enclosed paper, especially Figure 1. It is a reasonable and elegant approach. For the wall distance I used the wall distance (WDIST) computed as in the Spalart-Almaras model. For the boundary nodes themselfes as wall distance I used the MeshDisbor variable * 0.33 as in other subroutines.
Furthermore, I distinguish if the wall is fully slip (wall roughness = 0) which corresponds to a symmetry boundary condition or if there is a wall roughness. With symmetry conditions we should no have damping at all. Maybe one should revise the other turbulence models also?
2nd e-mail:
especially damping and differentiation between solid and liquid boundary, which I considered as a bug before.
[...]
Enclosed you will find the modified mixlength subroutine and the paper by Jia et al. Furthermore, the memory allocation in "point_telemac2d.f" as well as the call to the WDIST subroutine for the mixing length model in "telemac2d_init.F" have to be changed.
I tested it on some simple test cases and it works well. Maybe my coding style in the subroutine would need some optimizations.
Even if in large scale models the corrections don't make a huge difference, still I believe that we should always try to get the boundary conditions from a theoretical point of view right.
3rd e-mail:
0.2113 is the approximate solution of the exact solution (-sqrt(3)+3)/6. The exact solution results from the assumption, that at a relative distance d/h the parabolic profile is equal to the depth averaged value:
vt = 1/6 x kappa x h* x U* = kappa x h x U* x d/h (1-d/h)
Final decision, as no opposition: exact solution for (-sqrt(3)+3)/6 rather than 0.2113 direct computed solution for 0.3245 (Cardan's formulae), double precision: 0.32444233594882926 (no changes in vegetation results file if values are truncated to 4 significant digits or with double precision accuracy.
Enclosed, Jia's paper and program to solve the 2nd value.