Heat flux calculation

[Axel Cablé]

Hi,
I would like to compute the heat fluxes on some surfaces in my calculation domain. Which routine should I use to do that?
Thanks
 
 

[Matthieu Guillaud]

Hello,
 
to compute heat fluxes you have to use the routine usproj.f90.
 

[Axel Cablé]

Thanks. I tried to specify the face I want to calculate the heat flux on by using the following lines in usproj.F:

         CALL GETFBR('X=0
     .and Y<1.1 and Y>0.0 and Z<2.19
     .and Z>1.81',NLELT,LSTELT)

However, I get the following error message in listing:

Error parsing expression:
X=0 and Y<1.1 and Y>0.0 and Z<2.19 and Z>1.81
 ^ Expected operator instead of operand.

Changing < to <= doesn't change anything.

What did I do wrong?

[Yvan Fournier]

Hello,
The X = syntax is not allowed (only inequalities are, as an equality for real numbers is not recommended, so 'x <0.0001 and x > -0.0001' for example should work).
In any case, I prefer the plane[a, b, c, d, epsilon] operator syntax (such that abs(ax+by+cz+d) < epsilon), to the inequality syntax, as the latter can lead the user to believe we have a true equation interpreter, which is not the case, and than be surprised when "x + y < 0" for example doesn't work).
Search for "selection" in the user manual for more examples and the list of all operators.
Best regards,
  Yvan

[Axel Cablé]

Hi,

Thank you for your answer Yvan.

I modified usproj in order to get the flux on one face (see usproj.f at the end of this message).

However, if I use k-omega SST or k-epsilon std with scalable wall function or 2 scale wall function, the calculation keeps running but does not go further than 2 iterations, without giving any error message though(the "derive" value of omega reaches 10^12 after the first iteration though).

Any idea where the problem could come from?

usproj.F:
 
c BILAN SUR LES PAROIS OPAQUES
c PAROI OPAQUE HAUT
      CALL GETFBR('X <=0.0000 and X >=0.00
     .and Y <= 2.500 and Y >= 2.250 and Z <= 4.00
       .and Z >= 0.00',NLELT,LSTELT)

      DO ILELT = 1, NLELT
C
        IFAC = LSTELT(ILELT)
C
C ---   Element de bord
C
        IEL    = IFABOR(IFAC)
C
C ---   Variables geometriques
C
        SURFBN = RA(ISRFBN-1+IFAC)
        DISTBR = RA(IDISTB-1+IFAC)
C
C ---   Variables physiques
C
        VISCT  = PROPCE(IEL,IPCVST)
        FLUMAB = PROPFB(IFAC,IFLMAB)
C
        IF(IPCCP.GT.0) THEN
          XCP = PROPCE(IEL,IPCCP)
        ELSE
          XCP    = CP0IPH
        ENDIF
C
        IF(IPCVSL.GT.0) THEN
          XVSL = PROPCE(IEL,IPCVSL)
        ELSE
          XVSL = VISLS0(ISCAL)
        ENDIF
C
C ---   Calcul de la contribution au flux sur la facette courante
C         (flux de diffusion et de convection, negatif si entrant)
C
          XFLUXF =          SURFBN * DT(IEL) * XCP *
     &     (XVSL+VISCT/SIGMAS(ISCAL))/DISTBR *
     &     (COEFA(IFAC,ICLVAR)+(COEFB(IFAC,ICLVAR)-1.D0)
     &                                               *RA(ITRECO+IFAC-1))
     &                    - FLUMAB * DT(IEL) * XCP *
     &     (COEFA(IFAC,ICLVAR)+ COEFB(IFAC,ICLVAR)
     &                                               *RA(ITRECO+IFAC-1))
C
          XOPAQUEHAUT = XOPAQUEHAUT + XFLUXF
C
      ENDDO
C
C     - Somme des grandeurs sur tous les processeurs (calculs paralleles)
C
      IF (IRANGP.GE.0) THEN
    CALL PARSOM (XOPAQUEHAUT)
      ENDIF

c Ecriture de la valeur renvoyee

         OPEN(unit = 11, file = 'Fluxopaquehaut.hst', status = 'unknown',
     &   position = 'append')       WRITE(11,5600) XOPAQUEHAUT
 5600 FORMAT(' ',E14.5)

[Yvan Fournier]

Hello,
I would recommend closing the file you use in usproj each time, or keeping it open (using a counter to check if you are running the first pass or not; search for "ipass" in user subroutine examples to see how this may be done).
Otherwise, the crash may be due to a diverging calculation. What is the turbulence model/wall model combination which dos work ? What are the values of y+ you obtain (you may visualize them if you activate post-processing of the boundary mesh) ?
Best regards,
  Yvan

[Axel Cablé]

Hello,
It seems the problem occured only when running the calculation on multiple processors. Anyway, I now used the usproj heat flux calculation routine by only printing the results in the listing file and it works on multiple processors.
I have a question on the results I obtain though. From what I can read in usproj.F, the results obtained are given in Joules. However, I have a surface of 1.8m² in my domain which emits 69.4W/m² (-69.4 set as scalar boundary condition), but the computed value is 8.84, so it's probably not the real heat flux value in W (or J). Any advice on how I could get the heat flux? (the turbulence model I use is k-omega SST with a two-scale wall law).
Regards,
Axel

[Axel Cablé]

Would anyone have an idea on the problem I'm facing?
Thanks