how to calculate vorticity in CS

[Mubashir Ali]

hi,
 
I want to calculate vorticity for my simulation. I know to calculate vorticity i need to write code in USPROJ.F  or USVPST.F  isnt it ??
But I dont know how to find curl of velocity vector  which is actually vorticity.
 
your help will be of great use for me.
 
thanking in advance  :d
 
Mubashir
 
 

[César Vecchio]

In Paraview when you use streamlines there's an offshoot value called "velocity rotation" or something similar. Perhaps you can check there first and avoid writing code.

[Alexandre Douce]

Code_Saturne does not propose the vorticity by default. If the solution proposed by César works it is the simplest way to have what you want.

If you are interested by computing this vector, you have to write the curl operator, with a combination of gradients. Perhaps you can work directly in usvpst.f90 as you said, with several calls of the GRDCEL function. Another possibility is to work in usproj.f90, storing the computations in the array RTUSER, and then still post-process them in usvpst.f90.

[César Vecchio]

Are you trying to find circulation? Maybe this other can be an option too (I never tried it but I know something similar is done in tunnel tests, there's even a special probe):
I will asume you're working with a (pseudo)2D case, let's say an airfoil (I suppose you can try it too for 3D). In Paraview use the clip tool and extract the flow around your airfoil with a sphere shape:
CelldataToPoindata
Clip filter-->sphere--> center on midchord and radius=chord (can be smaller or bigger, as long as you enclose the body).
Now apply the filter "extract surfaces". Then make a slice paralell to your airfoil plane and at 50% of your mesh depth.
Well, now you sould see only two planar lines, one being the airfoil contour and the other a circle. Get rid of the airfoil contour by applying another clip filter (this time you can use sphere or box) and apply it so that only the circle remains.
So, that should be it, you have a closed path to integrate speed along it. Don't do it around the airfoil contour as speed is null on the surface, that's why I suggest getting rid of it. If you ever do this procedure, tell us how it went.

[Mubashir Ali]

@Cesar nd alexandre  for your info
 
@Cesar.... i could compute vorticity by  gradient of unstructured data .. thr is an option of computing vorticity.... I also tried your last suggest of calculating circulation but i was lost after having only the CIRCLE at the end.... and yes with streamlines also there is option of calculating vorticity.   But i still want to write code for vorticity in CS as i prefer to check vorticity at every time step without clicking
too much in PARAVIEW ...
 
@Alexandre ..... in the user document i could find only one line for GRDCEL i.e. calculation of the gradient of A (use of grdcel)   
Could you please tell me where i can find more INFO about GRDCEL function.
 
Thanx again to all who replied :D

[David Monfort]

Hi Mubashir,
You have an example of gradient comptation in the usproj.f90 subroutine (temperature gradient IIRC). If it isn't enough, let me know so that I can provide you with more details.
Also, you can use some Fortran 90 features like dynamic memory allocation for temporary (local) arrays. This is not used in the 2.0 yet, but nothing prevents you from using this feature ;)
David

[Mubashir Ali]

i wish to know if this approach is write or wrong ??/
curl of x-component = (delW/del Y) - (delV/delZ)
if i write the above equation in finite difference i get
(delW/del Y) = (Wi+1 - Wi) / (Yi+1-Yi)
(delV/del Z) = (Vi+1 - Vi) / (Zi+1-Zi)
 the above two eq. can be written in CS as follows
(delW/del Y)= [RTP(iel, IW(1)) - RTP(iel-1, IW(1)] / [xyzcen(2,iel) - xyzcen(2,iel-1)]
(delV/del Z) = [RTP(iel, IV(1)) - RTP(iel-1, IV(1)] / [xyzcen(3,iel) - xyzcen(3,iel-1)]
 
similar  equations can be written for Y and Z direction. when are 3 direction are combined i can get curl of Velocity vertor.
 
@David :  It  would be great if you can provide more info on use gradcel ... thnx in advace :D

[Mubashir Ali]

this the final code i wrote to calculate vorticity .... idk how right values i am getting yet

Code: Select all

  iphas = 1
  if (ipart .eq. -1) then                  ! Initialize variable name 
    do ii = 1, 32
      namevr(ii:ii) = ' '
    enddo

    ! Variable name
    namevr = 'Vorticity'
    ! Variable dimension (1: scalar, 3: vector, 6/9: symm/non-symm tensor)
    idimt = 3
    ! Values are interlaced.
    ientla = 1
    iphas = 1   ! We only consider phase 1 in this exampl

    do iloc = 1, ncelps
      iel = lstcel(iloc)
      tracel(1+(iloc-1)*idimt)=((RTP(iel,IW(1))-RTP(iel-1,IW(1)))      &
                    / (xyzcen(2,iel)-xyzcen(2,iel-1)))      &
                    -((RTP(iel,IV(1))-RTP(iel-1,IV(1)))      &
                     /(xyzcen(3,iel)- xyzcen(3,iel-1)))
      tracel(2 +(iloc-1)*idimt)=((RTP(iel,IU(1))-RTP(iel-1,IU(1)))      &
                    / (xyzcen(3,iel)-xyzcen(3,iel-1)))      &
                    -((RTP(iel,IW(1))-RTP(iel-1,IW(1)))      &
                     /(xyzcen(1,iel)- xyzcen(1,iel-1)))
      tracel(3 +(iloc-1)*idimt)=((RTP(iel,IV(1))-RTP(iel-1,IV(1)))      &
                    /(xyzcen(1,iel)-xyzcen(1,iel-1)))      &
                    -((RTP(iel,IU(1))-RTP(iel-1,IU(1)))      &
                     /(xyzcen(2,iel)- xyzcen(2,iel-1)))
    enddo       ! Values are defined on the work array, not on the parent.

    ivarpr = 0    ! Output values
    call psteva(ipart, namevr, idimt, ientla, ivarpr,          &
          ntcabs, ttcabs, tracel, trafac, trafbr)
  ENDIF

[Yvan Fournier]

Hello,

For an example of usage of grdcel (so as to use a better gradient approximation, especially in case of a non-orthogonal mesh), you may look in users/base/usproj.f90 (temperature gradient in that example).

Best regards,

  Yvan