what is efforts in code saturne

[César Vecchio]

Thank you, Mubashir, it works as you say. However, I'm having some issues:
The drag_lift.dat file is created and the first value is given, but then no other timestep is written despite the 'append' indication. The newer values are written in a 'fort.70' file but not in an appended way, the same line is updated with every new output. Anyway, this is a minor issue.
Also, and I think here comes the paralell thing, the first value is big (Cd=90.8) but then all the coefficients lower to 10e-300 and I don't see why. I changed the coeff formula in the following way because my chord is not 1 but 0.1524m:
 
XCOF(II) = 2.D0 * xfor(II) / (RO0(1) * UREF(1)**2 * 0.1524D0)
 
I also tried using ALMAX(1) with the same value but I get the very same results. Also, my compiler is gfortran and my CPU is a Sempron LE-1100 32bit and only one (modest) core :(. Should I remove the parrsm lines?
Thanks.

[Mubashir Ali]

hi Cesar,

If you r not doing parallel run then you should remove the following lines

if (irangp.ge.0) then
   call parrsm(ndim,xfor)
endif

In calculating Cd one divides drag by Dynamic pressure where Dynamic pressure is = 0.5*U*U*density  (U = free stream velocity)

while  UREF(IPHAS) is the characteristic flow velocity, used for the initialisation of the turbulence. Check if your UREF is same as U if not then replace it by free stream velocity U

one more thing you can check is dimensions of your grid. Code Saturne (CS) takes the grid im meters.   for example if you prepared the grid in mm it will be read as m in CS . i.e is to say if  chord length = 10 mm  then in Cs chrod length would be = 10 m

and lastly add these lines in your program to have appended file ..... at the very begining  as shown in my program above in this THREAD

Code: Select all

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))

[Yvan Fournier]

Hello,

You may remove the parrsm line when not running in parallel, but I do not recommend it: if you ever want to run the computation again in parallel (perhaps on another machine), you may forget that you removed it and have incorrect results. In any case, in serial mode, irangp is equal to -1, so parrsm is never called.

Which means there must be another cause.

To debug the issue, I would print (to the log file, using write(nfecra, *)) the first few xfor(ii) terms, to see how they change at different time steps. I would also print rho0(1) and uref(1) to make sure they were not overwritten.

Best regards,

  Yvan

[César Vecchio]

Ok, I think I got it. Apparently my gfortran is not very comfortable with some syntaxes, I had to change ACCESS='APPEND' to POSITION='APPEND', otherwise it would only write the first line and nothing else. Then, and this could be because I made my mesh with Gmsh and saved it directly as a MED file from there, the group name has to be enclosed by double quotation marks and then by the usual single quotation marks, for example in my case:
'Perfil' is NOT correct (the supersmall values I had before were actually zero, as they weren't reading any face) but '"Perfil"' works.
Also, and this is an issue with my compiler as well, the & symbol cannot be at the beggining of the line, otherwise I get a compile error, so it was matter only of moving those symbols to the end of the previous line. My gfortran is v4.4.3 and my Linux distro is Unity 32bit, in case you wonder. Anyway, small issues probabaly due to decaying support of F77 in favor of F90/95, just my guess.
What I'm getting now is that my values are like 50 to 70 times smaller than expected, for example my lift coefficient ( xcoff(2) ) should be something like 6.3e-1 but the usproj routine throws 9.07e-3. The "wrong" values are the efforts xfor; Rho and Uref values remain constant and equal to the freestream parameters all along the process, and my mesh and input is definetely in SI units (meters, mass Kg, etc.). The coefficient formula is being well computed, I checked it manually for some values, and I think it's all about some scaling factor, because the ratio lift/drag remains almost constant and with reasonable values, around 26.5 at 6º incidence.
If anyone wants, I can send the mesh, xml and usproj, it's a simple NACA 0012 airfoil and 1000 steady steps are enough for a decent convergence (one hour in my "monocore cluster" :D). Thank you all so far, and I guess I have no more excuses for not learning Fortran...

[Mubashir Ali]

send it to me at my email address  (Also send the listing file )
mubashir@iitk.ac.in
or
moobyknp@gmail.com
 
lets see whats the problem :D

[Alexandre Douce]

Hi César,
I am interested to, but I can not receive big files by email. Perhaps could you attach these files in a post or upload them in a file sharing web site like " sendspace ", " yousendit " or " mediafire " :)

[César Vecchio]

Ok, great, then there will be triple validation .

The mesh was created in Gmsh through an Octave script. As soon as I translate it I'll post it in these forums, it will be intended to mesh airfoils fully hexahedral. It's saved as a MED file so you can open it in Salome if you wish. The folders were created with:

cs create -s n0012 -c turb6,032

I hope Mubashir you're reading this because I forgot to add that in the mail I sent you

The compressed bz2 file attached includes the necessary files to run the simulation as intended, as well as the full listing file and the drag_lift.dat for the first 500 steps. You will notice the physical properties have somewhat odd pressure and density values, but it's simply because I used the data from a nitrogen cryogenic tunnel test and wanted to stick to the published data. You may find the full NASA report here: http://hdl.handle.net/2060/19880009181 (I don't attach it because it's 21MB). My files and mesh correspond to alpha=6.032deg from the first set of tables of the report, which roughly correspond to the Re=3e+6 published in Abbot's book (both experimental values differ a bit). Using only pressure information, my computed values are Cl=0.61 and Cd=0.017. The mesh has an y+=1 so using a k-e won't work well, I myself used SST (the hydraulic diameter value in the xml is ficticious and is only to initialize turbulence with a reasonable value). The groups defined for the mesh are the following:

• "Perfil" = airfoil surface (both upper and lower)
  "Entrada" = inlet
  "Salida" = outlet
  "Simetria" = lateral walls, for symmetry condition
  "Paredes" = top and bottom walls

In case you need any other information, don't hesitate in asking it. Thank you very much again.

EDIT
IMPORTANT: the parrsm lines are commented in my usproj file because I don't need them, so uncomment them if you will use more than one core.

[Mubashir Ali]

hi Cesar,
sorry i neva used GMSH :( .. but i tried to run your case but unfortunately med files r not supported in my CS installation.... idk y  but in the meanwhile hv u tried the suggestions i gave u ..... i am egarly w8ing  how different would be the reuslt :D

[César Vecchio]

Mubashir, your advices improved sensibly the results, specially for the extrados, so thank you for that :). I also sent you an I-DEAS file, perhaps it didn't arrive, don't worry.
However, usproj still refuses to compute correctly, the efforts summatory (xfor variable) is around 70 times less than expected in both components (Z component was already zero). Like I said before, if I compute the pressure efforts manually with a spreadsheet the results are good, which is what I've been doing so far, but accounting for viscous efforts should not modify those numbers greatly, and less reduce drag.
I tried making a basic debug. The adimensionalization by the dynamic pressure and the chord is well done, that's easy to check. Then I supposed that the routine was not cycling correctly around the airfoil contour, so I printed to the screen the number of elements, current face and dimension read, and everything seems to be in order, so I guess the routine does contour the airfoil correctly. The only thing that remains strange to me is then the RA(IFORBR + (IFAC-1)*NDIM + II-1) part. According to the formula, that should extract the effort components of each element we are cycling on. Could it be that it is extracting some other value? Or perhaps summing incorrectly because the projection on the normals is not well computed?
I'm just guessing, because looking at the pressures, velocities, streamlines, etc, everything suggests me the flow characteristhics have been well estimated, except for the efforts variable. In one of my previous posts on this thread I'm reporting an issue of Paraview displaying the efforts values, which perhaps is related.
For the record, I re-enabled the parrsm statements and it produces no effect, so it's working well.

[Yvan Fournier]

Hello,

I just checked the postprocessing of the eforts field, and the formula used there to obtain a value in Newtons is:

Code: Select all

surfbn = ra(isrfbn-1+ifac)
ra(iforbr+(ifac-1)*ndim)/surfbn
ra(iforbr+(ifac-1)*ndim + 1)/surfbn
ra(iforbr+(ifac-1)*ndim + 2)/surfbn

Confirming ra(iforbr+(ifac-1)*ndim) contains values already integrated over the face surface (which is consistent with what was written in the documentation and in previous posts here)

According to previous posts you postprocessed the efforts field, but the values seemed wrong. Did you use use filter -> extract blocks in ParaView to separate this part from the volume so as to make it usable.

Did you may also try this when moving usproj.f90 outside SRC, in case there is a bug there ? I noticed nothing wrong in the code, but sometimes debugging by bisection between what works and what doesn't may be the most efficient route (when all else fails).

Also, did you do any mesh sensitivity study, or simper, turbulence model sensitivity study for your calculation. If the mesh is too fine or too coarse and not well adapted to a turbulence model, viscous terms might change quite a bit.

Best regards,

  Yvan