TURB_floWlmUpdateLoglay.F90

Go to the documentation of this file.

! *********************************************************************
! * Rocstar Simulation Suite                                          *
! * Copyright@2015, Illinois Rocstar LLC. All rights reserved.        *
! *                                                                   *
! * Illinois Rocstar LLC                                              *
! * Champaign, IL                                                     *
! * www.illinoisrocstar.com                                           *
! * sales@illinoisrocstar.com                                         *
! *                                                                   *
! * License: See LICENSE file in top level of distribution package or *
! * http://opensource.org/licenses/NCSA                               *
! *********************************************************************
! *********************************************************************
! * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,   *
! * EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES   *
! * OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND          *
! * NONINFRINGEMENT.  IN NO EVENT SHALL THE CONTRIBUTORS OR           *
! * COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER       *
! * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE,   *
! * Arising FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE    *
! * USE OR OTHER DEALINGS WITH THE SOFTWARE.                          *
! *********************************************************************
!******************************************************************************
!
! Purpose: Model wall stresses based on logarithmic assumption at constant p.
!
! Description: Wall stresses are derived from modeling terms in BL equations.
!              Only viscous term is modeled based on log layer assumption at 
!              zero pressure gradient. BL convective term, pressure gradient
!              and time derivative term are not considered. Surface roughnes
!              is however taken into account. The method is inspired by paper
!              of Hoffman and Benocci, "Approximate wall BC for LES", 5th
!              Advances in Turbulence, Siena, Italy 1994. Surface roughness
!              were not considered in their model.
!
! Input: region  = data of current region.
!        patch   = current patch.
!
! Output: total and wall parallel components of wall stresses in body fitted 
!         coordinate.
!
! Notes: the pressure gradient and unsteady terms are modeled if BNDLAY 
!        (boundary layer) model is selected and added to the model viscous 
!        term computed here.
!
!******************************************************************************
!
! $Id: TURB_floWlmUpdateLoglay.F90,v 1.5 2008/12/06 08:44:44 mtcampbe Exp $
!
! Copyright: (c) 2001 by the University of Illinois
!
!******************************************************************************

SUBROUTINE turb_flowlmupdateloglay( region,patch )

  USE moddatatypes
  USE modbndpatch, ONLY   : t_patch
  USE moddatastruct, ONLY : t_region
  USE modglobal, ONLY     : t_global
  USE modinterfaces, ONLY : rflo_getpatchindices, rflo_getpatchdirection, &
                            rflo_getcelloffset
  USE moderror
  USE modparameters
  USE turb_modparameters
  IMPLICIT NONE

#include "Indexing.h"

! ... parameters
  TYPE(t_region) :: region
  TYPE(t_patch)  :: patch

! ... loop variables
  INTEGER :: i, j, k

! ... local variables
  CHARACTER(CHRLEN)       :: rcsidentstring
  TYPE(t_global), POINTER :: global

  INTEGER :: ilev, lbound, idir, jdir, kdir, iref, jref, kref, nref
  INTEGER :: ibeg, iend, jbeg, jend, kbeg, kend, icoff, ijcoff
  INTEGER :: n1, n2, ioff, ijbeg, ijend, errorflag
  INTEGER :: ijkc, ijkcip, ijkcim,ijkcjp, ijkcjm,ijkckp, ijkckm, ijkval

  REAL(RFREAL)              :: dens, wdist, rvonk, consb, ksize, abvel, utau
  REAL(RFREAL)              :: ypls, kspls, ksplsmax, func, dfunc, tauwall
  REAL(RFREAL)              :: minut, maxut, relaxcoef 
  REAL(RFREAL), POINTER     :: cv(:,:), dv(:,:), tv(:,:), vals(:,:)
  REAL(RFREAL), ALLOCATABLE :: mul(:)

!******************************************************************************

  rcsidentstring = '$RCSfile: TURB_floWlmUpdateLoglay.F90,v $ $Revision: 1.5 $'

  global => region%global
  CALL registerfunction( global,'TURB_FloWlmUpdateLoglay',&
  'TURB_floWlmUpdateLoglay.F90' )

! get dimensions and parameters -----------------------------------------------

  ilev   =  region%currLevel
  lbound =  patch%lbound
  nref   =  patch%valBola%switches(wlm_input_refpoint) - 1

  CALL rflo_getpatchindices( region,patch,ilev,ibeg,iend,jbeg,jend,kbeg,kend )
  CALL rflo_getpatchdirection( patch,idir,jdir,kdir )
  CALL rflo_getcelloffset( region,ilev,icoff,ijcoff )

  relaxcoef = 0.2_rfreal   ! coef. for utau iterat. in log-law (usr inp. later)

! get pointers and allocate temporary workspace -------------------------------

  cv    => region%levels(ilev)%mixt%cv
  dv    => region%levels(ilev)%mixt%dv
  tv    => region%levels(ilev)%mixt%tv
  vals  => patch%valBola%vals 

  n1    = abs(patch%l1end-patch%l1beg)
  n2    = abs(patch%l2end-patch%l2beg)
  ioff  = n1 + 1
  ijbeg = indij( 0, 0,ioff)
  ijend = indij(n1,n2,ioff)

  ALLOCATE( mul(ijbeg:ijend) )

! begin and end indices for i, j and k directions

! IF (lbound==1 .OR. lbound==2) THEN
!   indxb = jbeg
!   indxe = jend
!   jndxb = kbeg
!   jndxe = kend
! ELSEIF (lbound==3 .OR. lbound==4) THEN
!   indxb = kbeg
!   indxe = kend
!   jndxb = ibeg
!   jndxe = iend
! ELSE
!   indxb = ibeg
!   indxe = iend
!   jndxb = jbeg
!   jndxe = jend
! ENDIF

! Note that indx correspond to zeta-coordinate and jndx to xi.
! Get p gradient in xi and zeta direction and velocities depending on lbound.

  IF (lbound==1 .OR. lbound==2) THEN

! - j-dir == indx-dir == zeta-dir;  k-dir == jndx-dir == xi-dir
    DO k=kbeg,kend
      DO j=jbeg,jend
        DO i=ibeg,iend

          iref  = i+idir*nref
          ijkc  = indijk(iref   ,j      ,k      ,icoff,ijcoff)
          ijkcjp= indijk(iref   ,j+1    ,k      ,icoff,ijcoff)
          ijkcjm= indijk(iref   ,j-1    ,k      ,icoff,ijcoff)
          ijkckp= indijk(iref   ,j      ,k+1    ,icoff,ijcoff)
          ijkckm= indijk(iref   ,j      ,k-1    ,icoff,ijcoff)
          ijkval= indij(j-jbeg ,k-kbeg ,jend-jbeg+1)

          vals(ijkval,wlm_vals_dpdxi) = &
                          (dv(dv_mixt_pres,ijkckp)-dv(dv_mixt_pres,ijkckm))* &
                          0.5_rfreal/vals(ijkval,wlm_vals_dxi)
          vals(ijkval,wlm_vals_dpdzt) = & 
                          (dv(dv_mixt_pres,ijkcjp)-dv(dv_mixt_pres,ijkcjm))* &
                          0.5_rfreal/vals(ijkval,wlm_vals_dzt)
          vals(ijkval,wlm_vals_xiv) = &
                          dv(dv_mixt_uvel,ijkc)*vals(ijkval,wlm_vals_xix) + &
                          dv(dv_mixt_vvel,ijkc)*vals(ijkval,wlm_vals_xiy) + &
                          dv(dv_mixt_wvel,ijkc)*vals(ijkval,wlm_vals_xiz) 
          vals(ijkval,wlm_vals_etv) = &
                          dv(dv_mixt_uvel,ijkc)*vals(ijkval,wlm_vals_etx) + &
                          dv(dv_mixt_vvel,ijkc)*vals(ijkval,wlm_vals_ety) + &
                          dv(dv_mixt_wvel,ijkc)*vals(ijkval,wlm_vals_etz) 
          vals(ijkval,wlm_vals_ztv) = &
                          dv(dv_mixt_uvel,ijkc)*vals(ijkval,wlm_vals_ztx) + &
                          dv(dv_mixt_vvel,ijkc)*vals(ijkval,wlm_vals_zty) + &
                          dv(dv_mixt_wvel,ijkc)*vals(ijkval,wlm_vals_ztz) 
          vals(ijkval,wlm_vals_dens) = cv(cv_mixt_dens,ijkc)
          mul(ijkval)   = tv(tv_mixt_muel,ijkc)
        ENDDO
      ENDDO
    ENDDO

  ELSEIF (lbound==3 .OR. lbound==4) THEN

! - k-dir == indx-dir == zeta-dir;  i-dir == jndx-dir == xi-dir
    DO k=kbeg,kend
      DO j=jbeg,jend
        DO i=ibeg,iend

          jref  = j+jdir*nref
          ijkc  = indijk(i      ,jref   ,k      ,icoff,ijcoff)
          ijkcip= indijk(i+1    ,jref   ,k      ,icoff,ijcoff)
          ijkcim= indijk(i-1    ,jref   ,k      ,icoff,ijcoff)
          ijkckp= indijk(i      ,jref   ,k+1    ,icoff,ijcoff)
          ijkckm= indijk(i      ,jref   ,k-1    ,icoff,ijcoff)
          ijkval= indij(k-kbeg ,i-ibeg ,kend-kbeg+1)

          vals(ijkval,wlm_vals_dpdxi) = &
                          (dv(dv_mixt_pres,ijkcip)-dv(dv_mixt_pres,ijkcim))* &
                          0.5_rfreal/vals(ijkval,wlm_vals_dxi)
          vals(ijkval,wlm_vals_dpdzt) = & 
                          (dv(dv_mixt_pres,ijkckp)-dv(dv_mixt_pres,ijkckm))* &
                          0.5_rfreal/vals(ijkval,wlm_vals_dzt)
          vals(ijkval,wlm_vals_xiv) = &
                          dv(dv_mixt_uvel,ijkc)*vals(ijkval,wlm_vals_xix) + &
                          dv(dv_mixt_vvel,ijkc)*vals(ijkval,wlm_vals_xiy) + &
                          dv(dv_mixt_wvel,ijkc)*vals(ijkval,wlm_vals_xiz) 
          vals(ijkval,wlm_vals_etv) = &
                          dv(dv_mixt_uvel,ijkc)*vals(ijkval,wlm_vals_etx) + &
                          dv(dv_mixt_vvel,ijkc)*vals(ijkval,wlm_vals_ety) + &
                          dv(dv_mixt_wvel,ijkc)*vals(ijkval,wlm_vals_etz) 
          vals(ijkval,wlm_vals_ztv) = &
                          dv(dv_mixt_uvel,ijkc)*vals(ijkval,wlm_vals_ztx) + &
                          dv(dv_mixt_vvel,ijkc)*vals(ijkval,wlm_vals_zty) + &
                          dv(dv_mixt_wvel,ijkc)*vals(ijkval,wlm_vals_ztz) 
          vals(ijkval,wlm_vals_dens) = cv(cv_mixt_dens,ijkc)
          mul(ijkval)   = tv(tv_mixt_muel,ijkc)
        ENDDO
      ENDDO
    ENDDO

  ELSE

! - i-dir == indx-dir == zeta-dir;  j-dir == jndx-dir == xi-dir
    DO k=kbeg,kend
      DO j=jbeg,jend
        DO i=ibeg,iend

          kref  = k+kdir*nref
          ijkc  = indijk(i      ,j      ,kref   ,icoff,ijcoff)
          ijkcip= indijk(i+1    ,j      ,kref   ,icoff,ijcoff)
          ijkcim= indijk(i-1    ,j      ,kref   ,icoff,ijcoff)
          ijkcjp= indijk(i      ,j+1    ,kref   ,icoff,ijcoff)
          ijkcjm= indijk(i      ,j-1    ,kref   ,icoff,ijcoff)
          ijkval= indij(i-ibeg ,j-jbeg ,iend-ibeg+1)

          vals(ijkval,wlm_vals_dpdxi) = &
                          (dv(dv_mixt_pres,ijkcjp)-dv(dv_mixt_pres,ijkcjm))* &
                          0.5_rfreal/vals(ijkval,wlm_vals_dxi)
          vals(ijkval,wlm_vals_dpdzt) = & 
                          (dv(dv_mixt_pres,ijkcip)-dv(dv_mixt_pres,ijkcim))* &
                          0.5_rfreal/vals(ijkval,wlm_vals_dzt)
          vals(ijkval,wlm_vals_xiv) = &
                          dv(dv_mixt_uvel,ijkc)*vals(ijkval,wlm_vals_xix) + &
                          dv(dv_mixt_vvel,ijkc)*vals(ijkval,wlm_vals_xiy) + &
                          dv(dv_mixt_wvel,ijkc)*vals(ijkval,wlm_vals_xiz) 
          vals(ijkval,wlm_vals_etv) = &
                          dv(dv_mixt_uvel,ijkc)*vals(ijkval,wlm_vals_etx) + &
                          dv(dv_mixt_vvel,ijkc)*vals(ijkval,wlm_vals_ety) + &
                          dv(dv_mixt_wvel,ijkc)*vals(ijkval,wlm_vals_etz) 
          vals(ijkval,wlm_vals_ztv) = &
                          dv(dv_mixt_uvel,ijkc)*vals(ijkval,wlm_vals_ztx) + &
                          dv(dv_mixt_vvel,ijkc)*vals(ijkval,wlm_vals_zty) + &
                          dv(dv_mixt_wvel,ijkc)*vals(ijkval,wlm_vals_ztz) 
          vals(ijkval,wlm_vals_dens) = cv(cv_mixt_dens,ijkc)
          mul(ijkval)   = tv(tv_mixt_muel,ijkc)
        ENDDO
      ENDDO
    ENDDO

  ENDIF

! compute tau-wall from total contributions

  rvonk    = 2.5_rfreal
  consb    = 8.5_rfreal
  ksplsmax = 70._rfreal
  vals(:,wlm_vals_tauux:wlm_vals_tauwz) = 0._rfreal

  DO ijkval=ijbeg,ijend

    wdist = vals(ijkval,wlm_vals_wdist)
    ksize = vals(ijkval,wlm_vals_rough)
    dens  = vals(ijkval,wlm_vals_dens)
    abvel = sqrt( vals(ijkval,wlm_vals_xiv)**2+vals(ijkval,wlm_vals_ztv)**2 )

    IF (ksize > real_small) THEN

! --- surface roughness may be large enough to be taken into account

      utau  = rvonk*log( wdist/ksize ) + consb
      utau  = abvel/utau   

! --- utau at LHS above is large for high ksize, but turns to negative when
!     ksize becomes even larger. Theoretically there is maximum value of utau, 
!     i.e. when ksPls reaches value 70. This contrain is applied here.

      maxut = ksplsmax*mul(ijkval)/(ksize*dens) 
      utau  = min( utau, maxut )          ! for large ks but utau still > 0.
      IF (utau < 0._rfreal) utau = maxut  ! for large ks that utau < 0.

      kspls = utau*ksize*dens/mul(ijkval)
      IF (kspls >= 5._rfreal) THEN
!wlmCheckprobe-----------------------------------------
!        yPls  = ksPls*wdist/ksize
!        write(*,*)1,ijkVal,ksPls,yPls,utau,utau/abVel
!------------------------------------------------------
        goto 30
      ELSE
        goto 20
      ENDIF
    ENDIF

20  CONTINUE

! - surface hidraulically smooth, find utau using Newton-Raphson iteration

!    minUt = 1.E-8_RFREAL*abVel   ! constrain for low utau (near separation)
!    utau  = MAX( minUt,vals(ijkVal,WLM_VALS_UTAU) ) ! initial estimate

!    func = 1._RFREAL
!    DO WHILE (ABS( func ) > STOP_VALUE )
!      yPls  = utau*wdist*dens/mul(ijkVal)
!      func  = rVonk*LOG( yPls ) + 5.5_RFREAL - abVel/utau
!      dfunc = (rVonk*yPls + abVel)/(utau*utau)
!      utau  = MAX( minUt , (utau  - relaxCoef*func/dfunc) )
!wlmCheckprobe---------------------
!      write(*,*) ijkVal,func,yPls
!----------------------------------
!    ENDDO 

! - or using rough estimate (replaced by Dean`s correlation later)
    
    utau = 0.05_rfreal*abvel
    ypls = utau*wdist*dens/mul(ijkval)
!wlmCheckprobe---------------------
!    write(*,*) ijkVal,utau,wdist,dens,mul(ijkVal),yPls
!----------------------------------

! - check roughness height in wall unit

    kspls = utau*ksize*dens/mul(ijkval)
!wlmCheckprobe-------------------------------------
!    write(*,*)2,ijkVal,ksPls,yPls,utau,utau/abVel
!--------------------------------------------------
30  CONTINUE

! - wall stress in xi and zeta direction and store utau for next stage

    tauwall = dens*utau*utau
    vals(ijkval,wlm_vals_tauuy) = vals(ijkval,wlm_vals_xiv)/abvel*tauwall
    vals(ijkval,wlm_vals_tauwy) = vals(ijkval,wlm_vals_ztv)/abvel*tauwall

    vals(ijkval,wlm_vals_tauvx) = vals(ijkval,wlm_vals_tauuy)     
    vals(ijkval,wlm_vals_tauvz) = vals(ijkval,wlm_vals_tauwy)      

    vals(ijkval,wlm_vals_utau)  = utau

! - store total wall stress in heat-flux array for heat transfer comp. later
    vals(ijkval,wlm_vals_hflux) = tauwall
!checkprobe--------------------------------------------------------------------
!    write(*,*) region%procId,patch%lbound,ijkVal,vals(ijkVal,WLM_VALS_TAUUY), &
!               vals(ijkVal,WLM_VALS_TAUWY),utau/abVel
!------------------------------------------------------------------------------   

  ENDDO       ! ijkVal

! deallocate temporary arrays

  DEALLOCATE( mul )

! finalize --------------------------------------------------------------------

  CALL deregisterfunction( global )

END SUBROUTINE turb_flowlmupdateloglay

!******************************************************************************
!
! RCS Revision history:
! 
! $Log: TURB_floWlmUpdateLoglay.F90,v $
! Revision 1.5  2008/12/06 08:44:44  mtcampbe
! Updated license.
!
! Revision 1.4  2008/11/19 22:17:56  mtcampbe
! Added Illinois Open Source License/Copyright
!
! Revision 1.3  2004/03/19 02:53:24  wasistho
! prepared for RFLU
!
! Revision 1.2  2004/03/12 02:55:36  wasistho
! changed rocturb routine names
!
! Revision 1.1  2004/03/08 23:35:46  wasistho
! changed turb nomenclature
!
! Revision 1.4  2004/02/14 03:43:30  wasistho
! added new WLM parameter: reference point
!
!
!******************************************************************************