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, imrgra, 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