open TELEMAC-MASCARET The mathematically superior suite of solvers

Dear TELEMACers,

I'm unable to run a TELEMAC-3D simulation with sediment transport enabled when Nikuradse law is selected to define bed roughness. See the attached files for an example. I gets lots of lines like this one:
GRACJG (BIEF) : EXCEEDING MAXIMUM ITERATIONS: 200 ABSOLUTE PRECISION: NaN
but it also happens that my computer crashes with a blue screen or simply powers off. The model's behaviour is the same on other computers. I don't have these problems when specifying Strickler law.

I varied turbulence model, advection, boundary conditions, sidewall friction and many other parameters, but the behaviour remains the same. The problem comes from the combination of TELEMAC-3D, SISYPHE and Nikuradse.

The behaviour is the same on v6p2 and v7p0. I'm using Telemac on Windows 7 64-bit and it was compiled with Intel v11. I modified the distribution of horizontal layers' thickness to: 0.5, 10, 30, 50, 70, 90, 100% of water depth (see condim.f). Using equal cell height vertically makes no difference and the simulation still crashes with only two horizontal layers. If you wish to run this simulation without recompiling your TELEMAC system, you will need to remove the keyword MESH TRANSFORMATION. If you want to change the number of cells, you will need to set COMPUTATION CONTINUED to NO.

Any help would be greatly appreciated.

Yannick

Hello Yannick,

We'll try to look at this and reproduce this 'interesting' behaviour, but "we" is not a lot in August here, so please allow some time.

With best regards,

Jean-Michel Hervouet

Hello,

I am not a specialist of sediment transport (so sediment specialists are welcome to correct me), but I had a quick look at your case. It runs about 40 iterations before crashing, so it is rather a problem of numerical stability. I looked at the parameters, with a time step of 0.05 s amazingly small. So you probably are dealing with a physical model. You should however be able to run it with a higher time step.

Those who simulate physical models are generally plagued by a keyword :

THRESHOLD FOR VISCOSITY CORRECTION ON TIDAL FLATS

which is set to 0.2 by default and lets you little chance in your case. You can set it to 0.

I am also surprised by the very small height of the first layer. 0.005 is smaller than the bed roughness 0.0282. As we suppose that the first layer is in the logarithmic zone (this is used with the Nikuradse law) there seems to be a contradiction. Anyway using the Nikuradse law is always more unstable than other laws, because we have to find U* assuming the logarithmic law, so it is more sensitive to mesh issues on the vertical, whereas Strickler is a depth-averaged approach which is much less sensitive.

Such small layers are also a problem with aspect ratio, if you think of the horizontal size of elements. I am always worried that sedimentologists want to have such small layers on the bottom. Think also that very small layers may trigger tidal flats algorithms to work, as a matter of fact we have difficulties to design algorithms able to work both in the Pacific with depths of 5000 m (where a threshold of 0.2 is nothing) and in physical models with depths of 4 cm.

I hope this is helpful, with best regards,

Jean-Michel Hervouet

Hello Jean-Michel,

Thank you for your quick answer, especially for answering on a Friday afternoon during summer.

If I understand well, the coefficient that goes with Nikuradse law is the roughness height (and not Strickler or Manning coefficient). In that case, I was clearly using a coefficient that is too large (N=0.0262), considering the shallow water depth (4.14 cm). After adding the keyword THRESHOLD FOR VISCOSITY CORRECTION ON TIDAL FLATS = 0., setting bed roughness to 3*d50=0.00129 (using (26N)^6=9.99cm does not make much sense here), the sediment simulation works just fine in v7p0. In v6p2, the calculation does not even start and I obtain the following error message: LIT: REAR ERROR, ONE INTENDED TO READ A RECORD OF 1 VALUES OF TYPE I ON LOGICAL UNIT 1. In v6p2, I was able to run another simulation in a similar channel but I had to reduce the number of horizontal layers from 7 to 3 and use equal cell size vertically. But in v7p0, everything looks fine, which means that my problem is solved.

I’m conscious that using a near-bed horizontal layer (0.5% of water depth = 0.02cm) that is thinner than roughness elements (0.129cm) can cause numerical instability. However, v7p0 seems to be able to deal well with this situation (which was not the case with v6p2).

Regarding the time step, I'm usually selecting the larger value that the simulation can handle without crashing, as long as the Courant number is below unity. In my simulations (I'm studying a few cases like this one), the sediment simulation will not complete normally if the time step is larger than 0.1 (sometimes as low as 0.025) with the meshes that I’m using, depending on the exact parameter set.

Thank you again.

Yannick

Hello Yannick,

OK, I understand that the misunderstanding was the roughness coefficient, yes for Nikuradse law it is the grain size, the unit changes with the law.

I'll discuss with specialists (hopefully at least one...) tomorrow what we can do when the near-bed horizontal layer is thinner than the grain size, there is certainly some limitation here, but perhaps it is already done. Anyway this always happen on tidal flats.

I'm glad version 7.0 works better than 6.2, though your read error in 6.2 is on the geometry file and does not seem very normal.

For the rest, I agree with "selecting the larger value that the simulation can handle without crashing, as long as the Courant number is below unity".

With best regards,

Jean-Michel

Hello Jean-Michel,

I also encountered another issue in TELEMAC-3D v7p0 with the same case: I'm unable to select the k-epsilon turbulence model. This is why Smagorinski and mixed length models were selected, respectively for the horizontal and vertical turbulence. This happens even on a fixed bed (i.e. no sediment transport) and with equal horizontal layer thickness (7 horizontal layers) (see the attached file). The outcome is the same with reduced time step. It is possible that the problem originates from an incongruity in my case files. However, the same simulation with k-epsilon model works fine in v6p2 (if using Strickler law for friction).

Yannick