open TELEMAC-MASCARET The mathematically superior suite of solvers

Hi all!

I am modelling a debris boom in a dam reservoir. The structure is modelled as a drag force applied on a raw of nodes (in horizontal plane, 2D) and on the 3 highest planes which are located between the water surface and 0,5 m below it (the boom is 0,5 m deep).

The drag force is implemented as an implicite source term in source.f, with following expression (with guidance from JMH from an earlier thread):

S1U%R(N)=0.5D0 * CD * VELOCITY * COEFF

With:
CD = 1,9 (drag coefficient for a long thin plate perpendicular to the current)
VELOCITY = reference current velocity read from a run without the boom (this is this velocity that should be taken into account for computing the drag force according to theory)

The real expression of the drag force per unit surface is

F = 0,5 * CD * RHO * V**2

If I understand right and when modelling this force as a source term (ie. S = 0,5 * CD * V), one can only define one of the two V terms, the other one being used is computed by T3D for the given node at the current time step. This means that even when reading V from a reference run, we can't really compute exactely V**2 but only VREF*V. And as V is lower than VREF, the drag force is thus underestimated.
To counteract this problem, at least partially and in a simple way, I have added in a second run a coefficient COEFF being equal to VREF / (V with drag force read from a first run).

My first question is: at which step is the source multiplied by (RHO * V) and is it then possible to change V into VREF for the relevant nodes?


My second question concerns the obtained velocity field. With the method described above, the boom induces a local head loss and the local velocity is lowered by approx. 50%. The vertical velocity around the boom does not show any plunging effect (I exepected the current to be diverted under the boom), instead I observe a positive vertical velocity at the water surface, corresponding to the decrease of velocity head and the corresponding increase of water level.

To make things even clearer, I have run another case with an artifically increased drag coefficient (* 10) to force the velocity to be approx 0 m/s at the boom. In that situation, I expected to observe a clear negative vertical velocity around the boom (plunging current). As you can see on the pictures below, the vertical velocity is at best nil or positive on the 3 highest planes, and slightly negative on planes 5 and 6 below water surface upstream of the boom. I am also surprised that the velocity vectors at the water surface are not parallel to it.

So my second and last question is: How is continuity achieved in such a case and why can't we obtain a negative vertical velocity around the region where the drag force is applied?

The 3D model is run in non-hydrostatic mode.

Picture 1: vertical section perpendicular to the boom which is located where the velocity is almost nil.



Picture 2: same section zommed in around the boom and with velocity vectors (no distortion effects there!).



Picture 3: same that no. 2 but with vertical velocity in colormap.



Any insight would be much appreciated!

Thank you in advance!

Best regards
PL

Hello,

All your statements before question 1 are correct. So our procedure assumes that the local velocity is not too different from the infinite upwind velocity, otherwise you must correct this with a coefficient VREF/V, with V taken at the previous time step, taking care of course of divisions by 0. Actually velocity V in the formula is implicit, so you cannot substitute VREF at some level in the algorithm. We always treat friction and head loss in that implicit way to avoid that a friction force inverts the velocity (for an explicit friction, there is always a value of the time step that will invert the velocity).

Regarding question 2, in the hydrostatic assumption, the vertical velocity is only a "cosmetic" value that is not used in the computation and which is calculated only with the continuity equation. It would be better to try the non hydrostatic option (in that case also however, the real velocity taken into account for advection is the W* of the transformed mesh, that ensures exactly the continuity, but there is a pressure step to compute the divergence free velocity, without putting all the stress on the vertical velocity). Be also aware that there may be a distorsion factor when you look at 2D vertical views, so that the velocities are not parallel to the free surface. For the velocity at the free surface, there is also a keyword : DYNAMIC BOUNDARY CONDITION : YES that enforces strongly the free surface dynamic boundary condition.

With best regards,

Jean-Michel Hervouet

Hello Jean-Michel!

Thank you for your comments on question 1.

Regarding question 2, I did run the model in non-hydrostatic mode precisely because vertical velocities are an important result. But even with that, I must say that the velocities look rather "cosmetic".
I don't have any distortion factor issues, there, I used 1.0 in postel.

What is done exactely during the pressure step, to obtain W from W*? Can it be that there are some methodological limitations when it comes to compute vertical velocities?

I will try with DYNAMIC BOUNDARY CONDITION = YES and see what's happening at the free surface. Hopefully I observe a plunging stream underneath the boom and not a "jumping" stream above it!


Best regards
PL

Hi again,

I have now the result from the simulation with DYNAMIC BOUNDARY CONDITION = YES, and it indeed looks much better at the free surface! We can also observe a slightly increased plunging current under the boom, not as much as expected though.

Can it be that there are some methodological limitations when it comes to compute vertical velocities during the pressure step (to obtain W from W*)?

I attach some pictures of the new results below.

Picture 1



Picture 2: vertical velocity




Best regards
PL

Hello,

What if you increase the drag coefficient ? I am not sure that only the classical drag coefficient alone can fully divert the flow. The methodological limitation here is that we give only a head loss. What is missing is some porosity coefficient stating that there is an obstacle. This is what we did in 2D for vertical structures : head loss and porosity, that locally change the velocity.

With best regards (from Western Cape in RSA),

Jean-Michel Hervouet

Hello,

Thank you for your comments.

The results I posted previously already correspond to a "fictive" case with CD * 10. This was done to make sure the velocity behind the boom could be lowered to approx. 0 m/s and to observe what will be the vertical velocity.

We also did a test by changing the turbulence model from Smagorinski + Mixing length to k-epsilon. The test was done with a real drag coefficient, ie. with CD = 1,9.

The images show that with k-epsilon the velocity around the boom is significantly more reduced than with Smagorinski + Mixing length. The velocity gradient in the vertical direction just below the boom is stronger (the area influenced by the boom being smaller), which seems to be a more realistic result (no calibration/measurement data is available to confirm this though).

The vertical velocity is still rather low around the boom, where only a slight plunging current being observed.

It seems that k-espilon is more suitable for modelling head losses. Is it a theoritically expected result?

Is it possible to define porosity in 3D?

Picture 1: Smagorinski + Mixing length



Picture 2: k-epsilon



Picture 3: k-epsilon (zoom with vectors)




Bets regards
PL

Hello,

No porosity in 3D yet...

The fact that the k-epsilon works better is quite normal, in the mixing length model there is an assumption that the velocity profile is logarithmic, which is not the case with the head loss, or with wind...

With best regards,

Jean-Michel Hervouet

Thank you for your reply!

Best regards
PL