open TELEMAC-MASCARET The mathematically superior suite of solvers

Hi all!

I am building a new model which will aim, among other, to simulate recreation and resting areas for fishes that are to be created on the bottom and analyze the induced head losses. These areas are typically composed of blocks of rather small dimensions compared with the mesh size (0,2 – 0,5 m).
I was then thinking of reproducing them by simulating their drag force, corrected by a factor = min(h,d)/h, d being the structure’s diameter, as suggested in “Hydrodynamics of free surface flows”.

I had a look at the dragfo subroutine, which presents a case of a drag force simulated on a quadrilateral. In my case I will need to model several places where drag force is to be modeled. Does someone have any advice on how to adapt the example given in order to add several areas (irregular polygons) or on how to set the drag force on given mesh points (such as for the marine turbines example)?

The factor min(h,d)/h can be set as a constant since water depth will be rather constant on the given areas. But I was also thinking of using it with the real value of h in order to model a pipe crossing the river on its bottom in order to avoid refining the mesh. Can this be programmed in dragfo? If yes, can someone give me some hints on how h should be set?

Thank you in advance for your help and advices!
PL

Hi PL,

In dragfo subroutine, as you have the statement "USE DECLARATIONS_TELEMAC2D", you can have access to all variables defined in declarations_telemac2d.f file, in particular the water depths H (at the new time step) and HN (at the old time step). Is it what you want with your final question?

If you want to define another polygon different from a quadrilater, you just have to change the number of vertices NSOM, the maximum size of the arrays XSOM and YSOM and the coordinates of the polygon XSOM(1), YSOM(1) and so on...
If you have more than one polygon, I think that you can do a loop on the number of polygons, change the coordinates of the polygon inside the loop and perhaps NSOM if you use different kinds of polygons.

Hope this helps,

Chi-Tuan

Hi Chi-Tuan,

I wrote a reply but it seems to have disappeared…
Anyway, I write it again.

Regarding the use of water depth, I assume this is HN that should be used. Should I write it in the form HN%R(I) if I want the depth to be taken into account at each mesh node?

Regarding the possibility to add other polygons, I still need to have a closer look at the subroutine, but I share here my firsts questions.
How are polygons declared? I don’t really understand the relation between SOM, NSOM, XSOM and YSOM.
Once a new polygon is declared, I think I can manage to add quite easily a loop within the IF (INPOLY…) THEN section. Just a question: what means the command INPOLY(…)? I guess it means that for mesh points X(I) and Y(I) located in a polygon defined by XSOM,YSOM,NSOM, then the following actions are done. Is that right?

Thanks in advance!
PL

Hi again!

I did some test and I think I managed to do something.
I used the dragfo validation case given with 6.2, and I have added a second polygon where to apply the drag force near the upstream end of the mesh (between x=-90 and -60). Polygons now have 5 points and have the same area.

I removed the porosity option, just for the test.

As you can see on the graph, I have 2 zones with similar head losses.

I join the modified subroutine in case you want to check if there are no programmning errors. I especially don't understand what is "SOM" which is used in the CALL OS command. What does "CALL OS" ?
Is "SOM" related in any way the the polygon arrays (since they also use SOM in their names)?

Thanks for your help!
PL

HI PL,

Just a "quick" answer to your requests. I will take more time sooner to complete.
You can use HN%R(I) if you want the water depth at the previous time step for node number I (HN is a structure, HN%R is a pointer on floats).

XSOM and YSOM are arrays containing abscissae and ordinates of the polygon you want to define and NSOM is the number of vertices of this polygon (FYI SOM is the abbreviation of "sommet" which means vertex in French, it can be helpful to know some French words to understand some implementations in TELEMAC ). So you can fill in XSOM and YSOM with a lot of values (e.g. 10) but if you declare NSOM = 4 and give this value to INPOLY, you will then only consider a quadrilateral.
If you read carefully the dragfo subroutine, you may see that SOM is equal to the inverse of the area (the drag force applies on an area, not only for one point).
INPOLY is a function defined in bief directory. OS (for Operation on Structures) is a subroutine defined in bief directory also.
You can read the description of these functions and subroutines in the heading.
FYI, as told in other posts, the description of every subroutine, function and module follow doxygen standards and by looking at docs.opentelemac.org/ website, you can have the automatically generated description of each of them + call graphs. It is very convenient when you do not know so much about TELEMAC implementation.
So I can say that you were right with INPOLY.

'CALL OS('X=CX ',X=FUDRAG,C=SOM)' means that after this operation, FUDRAG is replaced by the product of the constant SOM and FUDRAG calculated before the calling of OS subroutine (SOM=1/AIRE is a constant and FUDRAG a structure).

Hope this helps,

Chi-Tuan

Hi Chi-Tuan, and thank you for this "quick" answer

I will have a look at docs.opentelemac.org.
In the meanwhile, I am just wondering if the section for division by area should be changed when declaring several polygons... In my test (see txt attached in my previous post) I have only changed the treatment of FUDRAG and FVDRAG, assuming that the division by area will be applied to all polygons.
In other words, I don't really konw if SOM = 1/AIRE will be applied to all polygons or not. And what if polygons have different areas? Will AIRE be exactly the area of each polygon separately all the way through? Or should I define AIRE and SOM for each polygon?

Otherwise, I am French, that's why I was suspecting SOM, and XSOM,YSOM to be linked somehow :)

Thanks again for your help!
PL

Hello,

In your case, SOM must be calculated for every zone. Basically we need a stress which is a force divided by an area, in your case this must be done zone by zone, and you must be also careful with parallelism, this is the reason why there is a P_DSUM(AIRE), in case a zone is split into several sub-domains.

Regards,

Jean-Michel

Hello !

I understand the process but am not sure on how should it be translated into code...

So far I have done this, for 2 polygons:
...
!POLYGON 1
IF(INPOLY(X(I),Y(I),XSOM,YSOM,NSOM)) THEN
UNORM = SQRT(UN%R(I)**2+VN%R(I)**2)
FUDRAG%R(I) = - 0.5D0 * N * DIAM * CD * UNORM
FVDRAG%R(I) = - 0.5D0 * N * DIAM * CD * UNORM
AIRE = AIRE + T1%R(I)
ENDIF
!POLYGON2
IF(INPOLY(X(I),Y(I),XSON,YSON,NSON)) THEN
UNORM = SQRT(UN%R(I)**2+VN%R(I)**2)
FUDRAG%R(I) = - 0.5D0 * N * DIAM * CD * UNORM
FVDRAG%R(I) = - 0.5D0 * N * DIAM * CD * UNORM
AIRE = AIRE + T1%R(I)
ENDIF
...
and with:
...
SOM = 1.D0 / AIRE
...
CALL OS('X=CX ',X=FUDRAG,C=SOM)
CALL OS('X=CX ',X=FVDRAG,C=SOM)
...

In this way, are there 2 or only 1 value affected to AIRE (for XSON,YSON polygon, the last one processed)?
If only 1, can I solve this by declaring:
AIRE1, AIRE2, SOM1, SOM2... replacing them in respective parts, and then finishing by:

CALL OS('X=CX ',X=FUDRAG,C=SOM1)
CALL OS('X=CX ',X=FVDRAG,C=SOM1)
CALL OS('X=CX ',X=FUDRAG,C=SOM2)
CALL OS('X=CX ',X=FVDRAG,C=SOM2)

A last question: maybe T1 also needs to be adapted (I will try to look exactly what T1 is...).

Thanks for your help...
PL

Hello,

It is not exactly like that...

I think you have to compute AIRE for the first polygon in a separate loop:

!POLYGON 1
AIRE=0.D0
LOOP ON I:
IF(INPOLY(X(I),Y(I),XSOM,YSOM,NSOM)) THEN
AIRE = AIRE + T1%R(I)
ENDIF
IF(NCSIZE.GT.1) AIRE=P_DSUM(AIRE)

then do SOM=1.D0/AIRE, then

another loop on I:
!POLYGON 1
IF(INPOLY(X(I),Y(I),XSOM,YSOM,NSOM)) THEN
UNORM = SQRT(UN%R(I)**2+VN%R(I)**2)
FUDRAG%R(I) = (- 0.5D0 * N * DIAM * CD * UNORM)*SOM
FVDRAG%R(I) = (- 0.5D0 * N * DIAM * CD * UNORM)*SOM
ENDIF

and the same for polygon 2

and you do no CALL OS after.

Regards,

JMH

Hello,

Thanks a lot, this makes sense.
I will give it a try as soon as I can have some time for it...

I change topic now, but have you seen my question on tubes treatment in this post:
www.opentelemac.org/index.php/kunena/16-...tubes-data-file#6522

I am a bit surprised by the conditioning used in the validation case. But it is maybe because I am not used to formulate orifice laws with head loss coefficients. It would be great to have a feeback from you about it.

Thanks!
PL