Page 1 of 2
calculating Drag over an airfoil
Posted: Wed Sep 22, 2010 3:44 pm
by 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
Re: calculating Drag over an airfoil
Posted: Wed Sep 22, 2010 10:58 pm
by 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
Re: calculating Drag over an airfoil
Posted: Thu Sep 23, 2010 2:24 pm
by 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
Re: calculating Drag over an airfoil
Posted: Thu Sep 23, 2010 3:39 pm
by 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
Re: calculating Drag over an airfoil
Posted: Fri Sep 24, 2010 12:54 pm
by 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
Re: calculating Drag over an airfoil
Posted: Fri Sep 24, 2010 3:11 pm
by 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
Re: calculating Drag over an airfoil
Posted: Sun Sep 26, 2010 9:24 am
by 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
Re: calculating Drag over an airfoil
Posted: Sun Sep 26, 2010 8:02 pm
by 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
Re: calculating Drag over an airfoil
Posted: Mon Sep 27, 2010 1:07 pm
by 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 ??
Re: calculating Drag over an airfoil
Posted: Mon Sep 27, 2010 1:59 pm
by 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