what is efforts in code saturne

[Mubashir Ali]

hi,
 
I wish to know what is efforts when boundary post processing is switched on in code saturne.
 
does it give shear stresses or pressure or force ??
 
any info on efforts would be of use for me
 
thnx
mubashir

[Alexandre Douce]

Hi Mubashir,
In the efforts you have pressure forces, viscous and turbulent stresses. In the sources you can follow the computation of efforts with the keyword ineedf.

[Mubashir Ali]

hi thnx for the reply. Actually i want to calculate forces on the airfoil (b.c = wall)

If i can calculate force on the airfoil ..then i can also calculate drag co-efficient which is = force on airfoil /  dynamic pressure

dynamic pressure = 0.5*density*u*u  (where u = free stream velocity= 36 m/sec in my case)

UREF(1) = 36 m/sec

I have a code which calculate efforts on the wall..... i divide these efforts by dynamic pressure to obtain Cd. But i get some absurd values in oder of 10 to the  power of 5

for example for 100th iteration i get

Cd= 0.65936E+05

i think i am making some mistake.  code you please check my code and tell what am i doing wrong

Code: Select all

IDEBIA = IDBIA0
      IDEBRA = IDBRA0
C
C=======================================================================
C CALCULATE VALUES ON A BOUNDARY FACE:
C=======================================================================
C
C Open file
      IF(NTCABS.EQ.1) THEN
        OPEN(FILE='drag_lift.dat',UNIT=IMPUSR(1))
C Write headings to the dat file
C CX, CY, CZ are the x,y,z coefficients
        WRITE(IMPUSR(1),'(4(A,1X))')
     &   ' TIME STEP  ',
     &   '     CX     ',
     &   '     CY     ',
     &   '     CZ     '
C 
        CLOSE(UNIT=IMPUSR(1))
C
      ENDIF
C
C XFOR is forces in x,y,z direction
C XCOF is coefs...      
      DO II=1,NDIM
        XFOR(II) = 0.D0
        XCOF(II) = 0.D0
      ENDDO
C      
C open the file
      OPEN(FILE='drag_lift.dat',UNIT=IMPUSR(1),ACCESS='APPEND')
C      
C load up the boundary face
      CALL GETFBR('wall,NLELT,LSTELT)
C     ===========
C loop over the boundary faces we called up      
      DO ILELT = 1, NLELT
C        
C tmp variable IFAC is face we are on
        IFAC  = LSTELT(ILELT)
C
C get (pressure + stress) forces in x,y,z
        DO II=1,NDIM
          XFOR(II)   = XFOR(II) + RA(IFORBR + (IFAC-1)*NDIM + II-1)
        ENDDO
C
C out the face loop
      ENDDO
C
C if in parallel do a gloabal sum
        IF(IRANGP.GE.0) THEN
          CALL PARRSM(NDIM,XFOR)
        ENDIF
C
C calculate x,y,zz coefficents
        DO II=1,NDIM
          XCOF(II) = 2.D0*XFOR(II)/RO0(1)/UREF(1)/UREF(1)
        ENDDO
C print our coefficients for current timestep
        WRITE(IMPUSR(1),'(I12,1X,3(E12.5,1X))')
     &   NTCABS,XCOF(1),XCOF(2),XCOF(3)
C
C close
        CLOSE(UNIT=IMPUSR(1))
C 
      RETURN
C
      END
c@z

[James McNaughton]

I'm not sure but will

Code: Select all

CALL GETFBR('wall,NLELT,LSTELT)

work correctly? Surely you need to end the quotes?

Also that code is for a unit length aerofoil is this true of your mesh?

James

[Mubashir Ali]

i have the right call getfbr('wal', NLELT, LSTELT)
without closing the quotes i get compiling error .... so i had corrected in USPROJ.F   but here i missed it somehow
never the less for checking purpose i replaced wall with outlet and the pressure i was getting was correct but for wall its so absurd.  i am attaching a pressure coutour symetry plane then you would get idea for what kind of body i am trying to calculate Cd.

[César Vecchio]

Hi, I'm a bit confused on the efforts variable definition. Are they forces or stresses? I mean, their units would be Newtons or Pascals? I ask because I'm trying to validate Saturne for airfoil aerodynamics and the visualizations are quite according to what I expect, but when I check the numbers there's definitely something wrong (from my part surely).
I do the postprocess in Paraview and then export the data over the airfoil contour ignoring Z components (2D case) previously using extract surface, clip and slice filters. Then in OOCalc I compute Cp from the pressure field and the freestream dynamic pressure and the results are as expected, even Cl. So far I'm doing for a symmetric airfoil at no incidence, I haven't tried more because I don't know how to compute drag and Cl from the efforts:
If I simply sum up all the efforts components individually (that's in X and Y alone) and then divide it by the product of dynamic pressure and the chord, I end up with a quite high Cl (0.045 vs. 7e-5 using Cp) and Cd is way too big (Cd=111 :( ). So then I tried weighting each value over the cell length projected on X, and there I got more reasonable results. Is this last procedure the correct one? Should I weight according to the cell length without projecting (and if were a 3D case over the cell area)?
Thank you, and my apologies if the question seems silly.

[César Vecchio]

Hi again. I tried this time with alpha=5º; the pressure distributions and Cl agree very well with experimental data, but when I try to compute Cl and Cd from the efforts I cannot get any logical value no matter how I do it.
I also tried using usproj.f90 but despite it compiles well, I don't get any output anywhere and I did define the output file in the GUI. I'm attaching both a spreadsheet with the output and my usproj file. I know it's some error from my part, but I just cannot find where. Thank you very much.

[David Monfort]

Hi all,

First, sorry for this late answer...

Efforts in Code_Saturne are actually forces (so, unit is Newton). This is what we compute because this is the value that can be summed upon the different faces. For example, the pressure efforts are computed as follow in the code (assuming we have an homogeneous Neumann boundary condition for the pressure variable).

Code: Select all

! Loop on faces
do ifac = 1, nfac
  ! Loop on dimensions
  do idim = 1, 3
    xfor(idim) = xfor(idim) + pressure(ifabor(ifac))*surfbo(idim,ifac)
  enddo
enddo

In the previous code example, pressure is the pressure in the cells, ifabor the boundary-face-to-cell connectivity array and surfbo the boundary surface array. If you want the Cp coefficient, you have to calculate the global surface area, divide the efforts by its value and projet the resulting vector on the freestream velocity.

Also, please note that what we call total pressure in Code_Saturne is the calculated pressure plus the hydrostatic pressure (usually 10e5 bar)...

Hope this helps... if not let me know!

David

[César Vecchio]

David, please don't apologize, I noticed you answered on Sunday so I appreciate it greatly . Also, thank you very much for your hints, they were indeed useful and I was able not only to get reasonable values but also close to experimental.

However, I think there's then some kind of bug in the way the "efforts" variable is written for Ensight values or the way Paraview inteprets them. Visit 2.1 cannot read at all the efforts if I import the Ensight file; it shows the name among the posibilities to plot, but if I try to make a graph (any kind) then it says Efforts doesn't exist and the same happens with other variables (valocity and pressure work ok). So I use Paraview 3.8.0, but it shows the efforts with correct direction but wrong sense and worse magnitude. I calculated with a Calc spreadsheet the efforts using the procedure you show in the Fortran code (Total Pressure = Pressure + Reference pressure). My case is an airfoil in a nitrogen cryogenic windtunnel, so my reference pressure is simply the freestream static pressure, no hydrostatic considerations given its negligible value. In this way, my results are as expected, and the maximum efforts magnitudes per cell are in the order of 300N, but Paraview shows values of almost 9000N and pointing outwards the surface (negative total pressure, physical impossibility) and that leads to the oversized drag and lift values.

Now, regarding the code you showed me, you said those are the pressure efforts only. As there's viscosity some of its effects are translated into the pressure, but then I suppose skin friction is not accounted, is it? Do you compute it from the dU/dy gradient and the user defined viscosity?

And finally, I cannot get my usproj routine to export the data anywhere, even after defining the output file in the GUI. How should I use the WRITE command? (I never programmed in Fortran, so I simply copy-paste codelines from here and there with the hope they'll work). I attached my usproj routine in my previous post, if you have time to give it a look.

Thank you for your time and patience

[Mubashir Ali]

USE    this condition  IF(NTCABS.EQ.1) THEN instead of if (ineedf.eq.1) then

and also remove  this checking condition from your program

iutile = 0

if(iutile.eq.0) return

Check the temp_saturne directory .... there you will find your drag_lift.dat file

write(nfecra,'(A, F8.3)')'usproj: Drag Co-efficient =',XCOF(1)

Copy paste the above line JUST before the WRITE command ..... this will print Drag-coeffiecint in the listing file.... it is useful to debug the program

last thing .... are you doing serial run or parallel run ....???  but your program it seems you are doing parallel run.

I hope these suggestions will help you.