open TELEMAC-MASCARET The mathematically superior suite of solvers

Dear all,

I am simulating an open channel with Telemac3d. The dimensions of the channel are 60 m (depth) * 500 m (width) * 2km (length).
I am using k-epsilon on both horizontal and vertical directions.
My mesh is refined up to 1m in the middle of the channel (see the figure mesh.jpg). This refinement has been determined so that the velocity converges in the wake of a turbine (the turbine implementation is presented in Thiebot et al., 2016, Oc Eng). There are 30 horizontal planes.
The case is enclosed to the message (the turbine has been removed).

My problem is that the values of k near the bed (at the 2nd plane for instance) are strongly affected by the mesh size. The zone where the cell size is 1 m acts as a sink term on k (see the figure k.jpg). The mesh has a much more limited effect on the velocity.

Do you have any idea to explain the sensitivity of k (and epsilon) to the size of the mesh ? Is there a minimal cell size (or a maximal growth rate to use when building the mesh) ?

Any help would be appreciated.

Thank you

Best regards

Jerome

Dear all,
I made additional tests and it seems that the problem comes from the velocities rather than from the turbulence model.
Simulations with PSI convection scheme on velocities give much better results (the results are not mesh dependent anymore).
Best regards.
Jérôme

Hello,

Good to hear that, what scheme were you using previously for advection? Have you tried the characteristics on the velocity?
Also you should be careful that the k-epsilon turbulence model contains a modelling of the near-wall turbulence, which assumes that your first plane is located around yu_*/nu = 30 with y the distance to the bed, nu the molecular viscosity and u_* the friction velocity (can be calculated based on the log profile, the bulk velocity and the roughness of the bed). You should not place the first plane much closer or further to the bed than this.

I hope this helps,
Best regards

Agnès

Hello Agnès,

Thank you for your reply.

I have tried both PSI and characteristics on the velocities.
When the turbulence model is disconnected from the momentum equation (VISCI = 0 in subroutine viscke), the PSI gives much better results than the characteristics as the results are not mesh dependant (see the enclosed figure (up) charact ; down) PSI)). However, when I reconnect the k-epsilon model to the momentum equation, the results are stil mesh dependent. I foresee that something is wrong with my implementation of the turbulence model...

As regards the vertical mesh, I use 30 planes. The first plane is thus located 2m above the seabed (the water depth is 60 m). Therefore yu_*/nu is much greater than 30.
However, unless I am wrong, when using a Strickler friction formulation, ustar (which is used as input of the wall function for k) is calculated from the mean velocity (and not from the velocity at the 2nd plane using the Reichard law). It is thus not necessary to use a fine mesh (yu_*/nu = 30) to have a proper estimate of ustar (and thus a proper estimate of k at the 1st plane). Could you confirm if it is right?

Thank you.

Best regards

Jérôme

Hello,

Are you using the same convection scheme for the velocity and the turbulent quantities?

I would say that even when using a Strickler friction model, your second plane should be at y+ around 30, calculating y+ with the u* given by the Strickler model (0.12(zeta/h)^1/6 if I remember well, check it though), but maybe I'm wrong? Also, you could try using a Nikuradse friction law instead. You have an empirical equivalence between the Strickler coefficient and the roughness zeta used in the Nikuradse model.

Best regards,
Agnès

Hello Agnès,

Thank you for your reply.

I use PSI for both the velocity and the turbulence.

I am using Nikuradse for regional modelling (with the roughness determined from the type of seadbed) and Strickler for numerical tests on tidal turbines.
I will check again the implementation of the wall function to see why the first plane should be located at y+ = 30 when using a Strickler formulation. Unless I am wrong, this condition is only required when a no-slip conditions is used at the bottom (as done in Thiebot et al., 2016, Oc Eng where y+ = 30).

I made additonnal tests to assess the sensitivity of the results to the mesh size. Those tests consist in simulating a flow in a simple canal (60 m x 500 m x 2000 m) with uniform mesh size (1m, 5m, 10m and 20m). The results indicate that the vertical profiles of k are mesh dependent (the differences reach up to 15% near the bottom). Although it is difficult to discriminate the different source of diffusion, I foresee that the differences are imputed to numerical diffusion on the turbulent quantities (when using the greatest cells).

Thank you for your help.

Best regards,

Jerome

Hi,

The PSI scheme is more diffusive than the characteristics, you could try using the latter for the advection of both the velocities and the turbulent quantities. Yes you're right the condition y+=30 is only required if you have a no-slip boundary condition, which seems to me the correct condition if you want to represent the turbulence effects due to the bottom...

Otherwise, you could try to tune the diffusion operator on k and epsilon:
ACCURACY FOR DIFFUSION OF K-EPSILON : 1e-8 or 1e-10
SOLVER FOR DIFFUSION OF K-EPSILON : 7
(7 corresponds to the GMRES solver).

Finally, in your cas file there is a very high value of prescribed flow rate, is this normal?

Best regards
Agnès

Hi Agnès,

Thank you for your reply.

My requirement is to find a model configuration suited to simulate both the turbine and the regional scales. I am using PSI because I didn't manage to run the regional model with the characteristics (the model rapidly crashes).
The requirement to simulate both scales is also the reason why I do not use no-slip condition (Strickler or Nikuradse formulation are more appropriate than no-slip conditions at the regional scale).

Thank you for your suggestion as regards the diffusion. I will try to tune the diffusion operator and see if it reduces the sensibility of the results to the mesh size.

The flow rate corresponds to a depth-averaged velocity of 2.3 m/s multiplied by the cross-sectionnal area (60 m * 500 m).

Thank you for your time.

With best regards,

Jérôme

Hello,

to support the suggestion from Agnès,
1.
indeed in some our cases in order to get steady state conditions the increasing of the accuracy -> ACCURACY FOR DIFFUSION OF K-EPSILON = 1.E-10 or 1.E-12 (default = 1.E-6) was the solution.
2.
Sometimes the default values for the limiters KMIN and EMIN created spurious oscillations at the inlet zone and there was no way to get rid of them. Reason is (at least my hypothesis), that with the Telemac default KMIN and EMIN values the calculated turb. viscosity is higher than the laminar viscosity which can cause some negative side effects. By changing the values to KMIN = 1.E-12 and EMIN = 1.E-14 the oscillations vanished (I have "stolen" these numbers from another software).
3.
Depending on the Telemac version we had also problems with the MOC scheme, in that the simulation crashed and we still don't know why. So we used for both the velocities and the k-e transport equations the N-SCHEME FOR TIDAL FLATS (14).
4.
By default Telemac-3D implicitly should apply a no-slip boundary condition when using the rough or smooth Nikuradse wall law or the Strickler law?! In the CFD community it is commonly accepted that the y+ value should be in the range of 30 to 300 in order to assure that the grid point is located in the logarithmic layer. This is conventional theory but maybe I have a misunderstanding of your case. In that case forget the last point.

Hope this helps a little bit.
Clemens

Hello Clemens,

Thank you for your feedback. Sorry for the delay in responding to you.

1. Following your suggestion (and the suggestion of Agnès), I made additional tests with different tuning of the diffusion. The results do not change significantly wich suggests that the mesh dependency can not be imputed to the numerical scheme for the diffusion.

2. Idem. The test with the new values of KMIN and EMIN do not significantly change the results.

3. I have to test this scheme on both the regional model and turbine model.
For information, I am using v7p1 and I do not have tidal flat in my regional model.

4. To my point of view, smooth regime and Nikuradse (rough regime) require a y+ value in the range 30 - 300 because ustar is calculated from the velocity of the 2nd plane (1st plane above the seabed) which should be located in the log layer. When using Chezy of Strickler, there is no requirement as regards the position of the 2nd plane because ustar is calculated from the depth-averaged velocity. Could you confirm it is right?

Thank you.

Best regards.

Jérôme