open TELEMAC-MASCARET The mathematically superior suite of solvers

Hello,

I would like to ask some questions about subroutine VC13PP. I need to use some values of the gradient alone in my own fortran code and what I do is using script(velocity gradient along x):

CALL VECTOR(UG1,'=','GRADF X',IELM3,
& 1.D0,U,SVIDE,SVIDE,SVIDE,SVIDE,SVIDE,
& MESH3,MSK,MASKEL)

I think in telemac3d linear prism structure it will call vc13pp.f, but I don't understand how this subroutine works:

1. What do these coefficients mean:
XS24 = XMUL/24.Do
XS144 = XMUL/144.D0

2. When calculate derivative along X for instance, why use real coordinates of Y but the distance of Z?

! REAL COORDINATES OF THE POINTS OF THE ELEMENT (ORIGIN IN 1)
!
! Y2 = Y(I2) - Y(I1)
! Y3 = Y(I3) - Y(I1)
Y2 = Y(IELEM,2)
Y3 = Y(IELEM,3)
!
Z2 = Z(I2) - Z(I1)
Z3 = Z(I3) - Z(I1)
Z4 = Z(I4) - Z(I1)
Z5 = Z(I5) - Z(I1)
Z6 = Z(I6) - Z(I1)

3. Can I have more details about the equations which are used in calculations of W1,W2,...,W6? Basically, I want to know how base functions are used here.

Thank you very much!
Best Regards,
Yue

Hello,

XMUL is a multiplication factor for the user, so no problem with this. The important thing to understand is that VECTOR computes the variational formulation of the gradient, i.e. the integral over the domain of the gradient multiplied by the test function. In VC13PP (so for prisms in 3D) this is done at prism level (and it will be assembled later). The X and Y here are relative to point 1, so that X2, X3, Y2, Y3 only are not 0. This is why you have a difference of treatment with Z. Note that the formulas with coefficients 1/24 or 1/144 have been computed exactly and written in Fortran by software like Maple. They are done in the reference prisms, with the requested multiplication by the Jacobian to account for the geometrical transformation, but the gradient introduces a division by the Jacobian and it simplifies.

After assembling all the element contributions on nodes, you do not have yet the gradient, you must then divide by the integral of the test function to get the true value of gradient at nodes. This integral is the "volume" around points, which is stored in arrays called VOLUN, VOLU (old and new time). In fact we do not divide by the volumes but multiply by the inverse of volumes, which are stored in UNSV3D. You can see in subroutine soukep.f at the beginning how we compute a number of gradients for the k-epsilon model, in a way that works also in parallel.

With best regards,

Jean-Michel Hervouet

Hi Jean-Michel.

Thanks for your quick reply. Because I am looking at the subroutine smago.f, does that mean actually,in the code,

2 * SIJ * SIJ = 2* (DU/DX *meshsize)**2 + 2* (DV/DY*meshsize)**2 +
(DU/DX*meshsize + DV/DX*meshsize)**2

I mean if:
g1 = gradF| x-direction, the surface is normal to x (let's call it S1)
g2 = gradF| y-direction, the surface is normal to y (let's call it S2)
g3 = gradF| z-direction, the surface is normal to z (let's call it S3)

Then the actual gradients are:

G1 = g1/S1
G2 = g2/S2
G3 = g3/S3

and when you compute

g1**2 + g2**2 + g3**2

it is different from

(G1**2 + G2**2 + G3**2) * meshsize **2

Is this an issue?

Best Regards
Yue

Hello,

Yes, you may be right, there could be an issue here. The 3D subroutine is copied on the 2D. In 2D it is correct as the extra area that we get from the variationnal formulation stands for the square of the mesh size. In 3D we get an extra volume... we have to think twice and review this to make sure, I'll do this as soon as possible.

With best regards,

Jean-Michel Hervouet

Jean-Michel

Thank you very much.

I am thinking if we could do like this:

2 * SIJ * SIJ = 2* (DU/DX * UNSV3D)**2 + 2* (DV/DY*UNSV3D)**2 +
(DU/DX*UNSV3D + DV/DX*UNSV3D)**2

Then

Nusmag = CS2 * (1/UNSV3D)**2 (2 * SIJ * SIJ)**(1/2)

Cheers
Yue

Hello,

I'll post a corrected version today. My idea was that compared to what is currently done, we need to multiply by unsv3d**(1/3). Dimensionnally we need to have M**2/S and Dij is 1/S in reality, but M**3/S as given by vector in 3D. unsv3d**(1/3) is 1/M (M: metre S: second).

With best regards,

Jean-Michel Hervouet

Hello,

Here are tentative corrections of smago.f and smago3d.f.

When there are large aspect ratios (horizontal mesh size much larger than vertical mesh size) this cannot be very good and in this case I would advise to take Smagorinski only on the horizontal (then it is smago.f), keeping a mixing-length model on the vertical. Anyway the coefficient CS can be tuned.

If we agree on this implementation I'll take it for the enxt release.

With best regards,

Jean-Michel Hervouet

Hi Jean-Michel,

Thanks a lot.

I will do some tests by this.

Best Regards
Yue

Hello Jean-Michel,

I did some test cases with corrected subroutines, modelling flow around a cylinder. It couldn't give better agreements comparing with experiment data. Then I change the code to calculate Sij by pure velocity gradient and time meshsize outside of square root. The results are better. But I am still struggling about the velocity profile close to the bottom. Does telemac use log law or any wall functions at layer above the bottom?

Best Regards
Yue

Hello,

Yes a log law is used near the bottom if you use a Nikuradse friction formula, or in smooth turbulence regime for the friction. You can look at subroutine tfond.f in library telemac3d.

With best regards,

Jean-Michel Hervouet