calculating Drag over an airfoil

[Mubashir Ali]

I am trying to calculate drag across an airfoil using the concept of momentum difference at inlet and exit
For example: inlet velocity is u1 and outlet velocity is u2
momentum at inlet (minlet)= density*u1*u1 integrated over the inlet boundary face momentum at outlet(moutlet) = density*u2*u2 integrated over the outlet boundary face
Drag = minlet - moutlet The above logic i am applying to code saturne, which is as follows using GETFBR to evaluate values at outlet boundaries

Code: Select all

         CALL GETFBR('outflow_rear',NLELT,LSTELT)
         DO ILELT = 1, NLELT
           IFAC = LSTELT(ILELT)
           IEL = IFABOR(IFAC)
           ux = ux+propfb(ifac,RTP(IEL,IU(IPHAS)))
           moutlet= den*ux*ux
         ENDDO

         ! Evaluating momentum at inlet boundary
         CALL GETFBR('vinlet',NLELT,LSTELT)
         DO ILELT = 1, NLELT
           IFAC = LSTELT(ILELT)
           IEL = IFABOR(IFAC)
           uxinlet= uxinlet+propfb(ifac,RTP(IEL,IU(IPHAS)))
           minlet= den*uxinlet*uxinlet
         ENDDO

         Drag = minlet- moutlet

Is this procedure right???? I am not sure if propfb(ifac,RTP(IEL,IU(IPHAS))) will provide me velocity at boundary face or not as velocity in code saturne is calculated at cell center And also i wish to know if there is possibility to apply free stream turbulence model to inlet boundary in code sature. If yes then how ? waiting for ur replies Regards Mubashir

[Jean Deschodt]

Hi,

RTP(IEL,IU(IPHAS)) is the value of the fisrt component of the velocity and not an index, so the expression propfb(ifac,RTP(IEL,IU(IPHAS))) is wrong, because an array waits for indexes.

Please look at this thread for additional information:  https://code-saturne.info/products/code-saturne/forums/general-usage/778113955

Otherwise I am not sure that the concept of momentum difference if the right way to do a drag evaluation with an incompressible algorithm...
 
Regards

[Mubashir Ali]

Previously Jean Deschodt wrote:

Hi,

RTP(IEL,IU(IPHAS)) is the value of the fisrt component of the velocity and not an index, so the expression propfb(ifac,RTP(IEL,IU(IPHAS))) is wrong, because an array waits for indexes.

Please look at this thread for additional information:  https://code-saturne.info/products/code-saturne/forums/general-usage/778113955

Otherwise I am not sure that the concept of momentum difference if the right way to do a drag evaluation with an incompressible algorithm...

 
Regards

Hi ,

yes you are absolutely right. I have finally got how to evaluate velocity in x-direction on inlet boundary face. And i am getting absolutely right answer. I write the code below

Code: Select all

 IPHAS= 1 
 CALL GETFBR('vinlet',NLELT,LSTELT)
 ICLVAR  = ICLRTP(IU(IPHAS),ICOEF)
 DO ILELT = 1, NLELT
   IFAC = LSTELT(ILELT)
   IEL = IFABOR(IFAC)
   uinlet=COEFA(IFAC,ICLVAR)+COEFB(IFAC,ICLVAR)*RTP(IEL,IU(IPHAS))
   minlet= den*uinlet**2
ENDDO 

 
where uinlet is the calculated inlet velocity and minlet is the momentum.
 
but now i want to calculate velocity on a plane parallel to my inlet boundary face...... this plane can be after a gap of few cells from inlet bounday face .

could you please help me with this
 
thnx in advance regards Mubashir
 

[Yvan Fournier]

Hello,

To select a plane "inside" the mesh, use "getfac" (instead of getfbr) with the plane[] operator (see user documentation for Syntax).

Caution: if you select interior cells, you will need to use ifacel instead of ifabor (with one cell on each side, so you may want to compute averages), and you do not have coefa or coefb arrays for interior faces.

A simpler way is probably directly to select cells using "getcel" with with plane[] operator (with a tolerance) or with the box[] operator. You can the use "rtp(iel(iu(iphas))" directly for the selected cells.

In any case, cell or face selection is based on the center of gravity, and Code_Saturne does not compute a "cut" the way a post-processing tool would. If your cells are aligned with the plane in the area, you want to select, things should b fine, but if this is not the case, using a post-processing tool such as EnSight, ParaView, or VisIt is probably a better idea (Code_Saturne's selection was designed to be used a a pre-selection to limit postprocessing output volume when using the usdpst and usvpst user subroutines).

Regarding the inlet turbulence, several models are available for the turbulence variable inlet conditions, so you may want to check the documentation regarding this (one is hydraulic diameter based, the other is turbulent intensity based if I remember correctly).
 
Regards,

  Yvan

[Mubashir Ali]

hi,

Below is the code i am using to calculate velocity in x direction at a distance of 400 mm from the origin parallel to yz plane.

Code: Select all

        call getcel('plane[400,0,0,0,inside]', nlelt, lstelt) 
   
         do ilelt = 1, nlelt 
         iel = lstelt(ilelt) 
   
         ux = ux+ rtp(iel,IU(iphas)) 
         mrear= den*ux**2   
         enddo 

but here i have some doubts .... firstly, in plane command  what does  "inside" means.  Secondly, how would i know where is the origin .. is it the same origin which i took while making the geometry and the grid or code Saturne defines its own origin. I used gambit for making the grid.
 
Thnx in advance
 
Regards
 
Mubashir

[Yvan Fournier]

Hello,

You probably do not want to use "inside" for the lan operator, as it includes all the half-space behind the plane (with the normal facing "outside"). You may use "epsilon=" to add a thickness to the plane, or use the "box" operator instead.

The origin is not changed by Code_Saturne, so it should be the one used by Gambit, unless Gambit defines a specific coordinate system which we might ignore. In any case, you should use a visualization tool to check how your mesh is imported under Code_Saturne, and the coordinates as they appear in EnSight or ParaView are the ones used by Code_Saturne.

Regards,

Yvan

[Mubashir Ali]

hi,
 
i tired out your suggestion and used the following code to calculate velocity at a plane parallel to yz plane and of thickness  0.01 mm. Below is the code:
 

Code: Select all

 call getcel('plane[4000,0,0,0,.01]', nlelt, lstelt) 
         do ilelt = 1, nlelt 
         iel = lstelt(ilelt) 
         uxplane = uxplane+ rtp(iel,IU(iphas)) 
      ENDDO 

but the value on velocity on this plane which i call uxplane is coming out to be much much larger but i expect it to be around 30-50 m/sec. Please check my code if i am doing something wrong.
 
And also i checked in paraview ...... the origin is taken where i had defined in gambit. So, this query is solved
 
Regards
 
Mubashir
 

[Yvan Fournier]

Hello,

I presume you are summing a velocity value over cells to compute an average, but then are you also summing a denominator ?

The following code would seem better:

Code: Select all

uxplane = 0
denom = 0

call getcel('plane[4000,0,0,0,.01]', nlelt, lstelt)

do ilelt = 1, nlelt
    iel = lstelt(ilelt)
    uxplane = uxplane + rtp(iel,IU(iphas))*volume(iel)
    denom = denom + volume(iel)
enddo

uxplane = uxplane / denom

You may also want to define an extra post-processing mesh in usdpst using the same selection criteria and load it under ParaView to make sure your selection works as expected.

Regards,

Yvan

[Mubashir Ali]

hi Mr Yvan,

i tried the above code for calculating velocity at a place defined using plane command.  but i am not getting correct values of velocity which i was expecting. Only thing i changed was i am having a place at 160 mm from origin not 4000 and also epsilon = 0.001
 
i think if i can get velocity at face of the plane it would be much better for my calculation. So, i have written the code below to calculate velocity at face of a plane at 160 mm from origin and epsilon = 0.001

Code: Select all

uxface= 0 
surface=0 
call getfac('plane[160,0,0,0,.0001]', nlelt, lstelt)
do ifac = 1, nfac
  iel = ifacel(2,ifac)
  uxface = uxface+ rtp(iel,IU(iphas))*ra(isrfan + ifac-1)
  surface = surface + ra(isrfan + ifac-1)
enddo
uxface= uxface/surface 

but i run code saturne i get an error which is pasted below

usp:velo at plane=   0.000
SIGSEGV signal (forbidden memory area access) intercepted! Call stack:
   1: 0x2ae8f9293a94 <fvm_selector_get_list+0x174>    (libfvm.so.0)
   2: 0x425d3c     <csgfac_+0xad>                   (cs13.exe)
   3: 0x411186     <getfac_+0x1e>                   (cs13.exe)
   4: 0x410b4c     <usproj_+0x294>                  (cs13.exe)
   5: 0x454205     <caltri_+0x5f85>                 (cs13.exe)
   6: 0x41172a     <main+0x54a>                     (cs13.exe)
   7: 0x3db7c1d994 <__libc_start_main+0xf4>         (libc.so.6)
   8: 0x410809     ?                                (?) End of stack

 
Could you please tell what am i doing wrong ??

[Yvan Fournier]

Hello,

You should have:

do iel = 1, nlelt
ifac = lstelt(iel)

And not:
do ifac = 1, nfac

If you want to use your selection.

Still, it seems it is the selection operation itself which is crashing, and this should not happen (looks like some list is not allocated correctly).

Also, I just noticed you are using millimeters instead of meters. I would recommend switching your mesh (and thus selection criteria) to meters. Otherwise, you'll need to adjust all the physical constants in the code to account for the fact that you are not using the standard units, which is quite risky.
 
If you still have the crash, you may want to post your case here (if it is both very small and non-confidential), or send to to Code_Saturne support via a large file send service (or put it on an ftp site) so that we may check why the crash occurs, but I would recommend adjustign your units first.

Regards,

 Yvan