Calculation of a local Nusselt number

This function is called at the end of each time step, and has a very general purpose (i.e. anything that does not have another dedicated user subroutine)

cs_f_user_extra_operations subroutine

Local variables

! Local variables

integer iscal, ivar, iortho, iprev
integer inc, iccocg, ilelt, irangv
integer ifac, iel, npoint, iel1, irang1, ii, iun, impout, nlelt, neltg

double precision diipbx, diipby, diipbz
double precision xtbulk, xubulk, xyz(3), xtb, xub, tfac, lambda, xab

character*19 namfil

integer, allocatable, dimension(:) :: lstelt
double precision, dimension(:), pointer :: coefap, coefbp, cofafp
double precision, allocatable, dimension(:,:) :: grad
double precision, allocatable, dimension(:) :: treco,treglo,treloc
double precision, allocatable, dimension(:) :: xnusselt
double precision, allocatable, dimension(:) :: xabs, xabsg
double precision, dimension(:,:), pointer :: vel
double precision, dimension(:), pointer :: cvar, viscl

double precision :: height, prandtl, qwall
parameter(height = 1.0d0)
parameter(prandtl = 1.0d0)
parameter(qwall = 1.0d0)

Calculation of the Nusselt number

if (ntcabs.eq.ntmabs) then

  allocate(lstelt(max(ncel,nfac,nfabor)))
  call field_get_val_v(ivarfl(iu), vel)

  iscal = iscalt         ! temperature scalar number
  ivar =  isca(iscal)    ! temperature variable number
  call field_get_val_s(iviscl, viscl)
  call field_get_coefa_s(ivarfl(ivar), coefap)
  call field_get_coefb_s(ivarfl(ivar), coefbp)

  call field_get_val_s(ivarfl(ivar), cvar)

Compute value reconstructed at I' for boundary faces

  allocate(treco(nfabor))

  iortho = 0

General cas (for non-orthogonal meshes)

  if (iortho.eq.0) then

Allocate a work array for the gradient calculation

    allocate(grad(3,ncelet))

Compute gradient

    iprev = 0
    inc = 1
    iccocg = 1

    call field_gradient_scalar(ivarfl(ivar), iprev, 0, inc, iccocg, grad)

Compute reconstructed value in boundary cells

    do ifac = 1, nfabor
      iel = ifabor(ifac)
      diipbx = diipb(1,ifac)
      diipby = diipb(2,ifac)
      diipbz = diipb(3,ifac)
      treco(ifac) =   coefap(ifac) + coefbp(ifac)*(cvar(iel)       &
           + diipbx*grad(1,iel)  &
           + diipby*grad(2,iel)  &
           + diipbz*grad(3,iel))
    enddo

Free memory

    deallocate(grad)

Case of orthogonal meshes

  else

Compute reconstructed value

Note

Here, we assign the non-reconstructed value.

    do ifac = 1, nfabor
      iel = ifabor(ifac)
      treco(ifac) = coefap(ifac) + coefbp(ifac)*cvar(iel)
    enddo

  endif

  impout = impusr(1)
  if (irangp.le.0) then
    open(file="Nusselt.dat",unit=impout)
  endif

  call getfbr('normal[0,-1,0,0.1] and y < 0.01',nlelt,lstelt)

  neltg = nlelt
  if (irangp.ge.0) then
    call parcpt(neltg)
  endif

  allocate(xabs(nlelt))
  allocate(treloc(nlelt))
  allocate(xnusselt(neltg))
  if (irangp.ge.0) then
    allocate(xabsg(neltg))
    allocate(treglo(neltg))
  endif


  do ilelt = 1, nlelt
    ifac = lstelt(ilelt)
    xabs(ilelt) = cdgfbo(1,ifac)
    treloc(ilelt) = treco(ifac)
  enddo

  if (irangp.ge.0) then
    call paragv(nlelt, neltg, xabs, xabsg)
    call paragv(nlelt, neltg, treloc, treglo)
  endif

  do ilelt = 1,neltg

Calculation of the bulk temperature

    if (irangp.ge.0) then
      xab = xabsg(ilelt)
    else
      xab = xabs(ilelt)
    endif

    xtbulk = 0.0d0
    xubulk = 0.0d0

    npoint = 200
    iel1   = -999
    do ii = 1, npoint
      xyz(1) = xab
      xyz(2) = float(ii-1)/float(npoint-1)
      xyz(3) = 0.d0
      call findpt(ncelet, ncel, xyzcen, xyz(1), xyz(2), xyz(3), iel, irangv)
      !==========
      if ((iel.ne.iel1).or.(irangv.ne.irang1)) then
        iel1   = iel
        irang1 = irangv
        if (irangp.eq.irangv) then
          xtb = volume(iel)*cvar(iel)*vel(1,iel)
          xub = volume(iel)*vel(1,iel)
          xtbulk = xtbulk + xtb
          xubulk = xubulk + xub
          lambda = cp0 * viscl(iel) / prandtl
        endif
      endif
    enddo

    if (irangp.ge.0) then
      iun = 1
      call parall_bcast_r(irangv, lambda)
      call parsom(xtbulk)
      call parsom(xubulk)
    endif

    xtbulk = xtbulk/xubulk
    if (irangp.ge.0) then
      tfac = treglo(ilelt)
    else
      tfac = treloc(ilelt)
    endif

    xnusselt(ilelt) = qwall * 2.d0 * height / lambda / (tfac - xtbulk)

  enddo

  if (irangp.eq.-1) then
    do ii = 1, neltg
      write(impout,'(2E17.9)') xabs(ii)*10, xnusselt(ii)/(0.023d0*30000.d0**(0.8d0)*0.71d0**(0.4d0))
    enddo
  endif

  if (irangp.eq.0) then
    call sortc2(xabsg, xnusselt, neltg)
    do ii = 1, neltg
      write(impout,'(2E17.9)') xabsg(ii)*10, xnusselt(ii)/(0.023d0*30000.d0**(0.8d0)*0.71d0**(0.4d0))
    enddo
  endif

  close(impout)

  deallocate(treco)
  deallocate(xabs)
  deallocate(xnusselt)
  deallocate(lstelt)
  deallocate(treloc)
  if (irangp.ge.0) then
    deallocate(xabsg)
    deallocate(treglo)
  endif
endif

return
end subroutine cs_f_user_extra_operations

sortc2 subroutine

Local variables

implicit none
integer n
double precision tab1(n), tab2(n)

integer ns, ii, jj, kk
double precision tabmin, t2

Body

ns = 1
do ii = 1, n-1
  tabmin = 1.d20
  kk = 0
  do jj = ns,n
    if (tabmin.gt.tab1(jj)) then
      tabmin = tab1(jj)
      kk = jj
    endif
  enddo
  t2 = tab2(kk)
  tab1(kk) = tab1(ns)
  tab2(kk) = tab2(ns)
  tab2(ns) = t2
  tab1(ns) = tabmin
  ns = ns + 1
enddo

return
end subroutine sortc2