RFLU_ModViscousFlux.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: Collection of HLLC flux routines.
!
! Description: None.
!
! Notes: None. 
!
! ******************************************************************************
!
! $Id: RFLU_ModViscousFlux.F90,v 1.10 2008/12/06 08:44:24 mtcampbe Exp $
!
! Copyright: (c) 2005 by the University of Illinois
!
! ******************************************************************************

MODULE rflu_modviscousflux

  USE moddatatypes
  USE modglobal, ONLY: t_global
  USE modparameters
  USE moderror
  USE moddatastruct, ONLY: t_region
  USE modbndpatch, ONLY: t_patch
  USE modgrid, ONLY: t_grid

  IMPLICIT NONE

  PRIVATE
  PUBLIC :: rflu_enforceheatflux, & 
            rflu_viscousfluxes, &  
            rflu_viscousfluxespatches
  
! ******************************************************************************
! Declarations and definitions
! ******************************************************************************  
      
  CHARACTER(CHRLEN) :: & 
    RCSIdentString = '$RCSfile: RFLU_ModViscousFlux.F90,v $ $Revision: 1.10 $' 
              
! *****************************************************************************
! Routines
! *****************************************************************************

  CONTAINS





! ******************************************************************************
!
! Purpose: Enforce heat flux.
!
! Description: None.
!
! Input: 
!  pRegion      Pointer to region
!  tv           Transport variables
!  tvIndxCond   Index to conductity entry in tv
!
! Output: None.
!
! Notes: None.
!
! ******************************************************************************

SUBROUTINE rflu_enforceheatflux(pRegion,tv,tvIndxCond)
   
  IMPLICIT NONE

! ******************************************************************************
! Declarations
! ******************************************************************************

! ==============================================================================  
! Arguments 
! ==============================================================================  

  INTEGER, INTENT(IN) :: tvindxcond
  REAL(RFREAL), DIMENSION(:,:) :: tv
  TYPE(t_region), POINTER :: pregion

! ==============================================================================  
! Locals
! ==============================================================================  

  INTEGER :: c1,distrib,ifg,ifgbeg,ifgend,ifl,ipatch
  REAL(RFREAL) :: cond,dtdn,dtdx,dtdy,dtdz,nx,ny,nz
  REAL(RFREAL) :: tvf(tvindxcond)
  TYPE(t_global), POINTER :: global
  TYPE(t_grid), POINTER :: pgrid
  TYPE(t_patch), POINTER :: ppatch  

! ******************************************************************************
! Start
! ******************************************************************************

  global => pregion%global

  CALL registerfunction(global,'RFLU_EnforceHeatFlux',&
  'RFLU_ModViscousFlux.F90')

! ******************************************************************************
! Set pointers and variables
! ******************************************************************************

  pgrid => pregion%grid

! ******************************************************************************
! Loop over patches
! ******************************************************************************

  DO ipatch = 1,pgrid%nPatches
    ppatch => pregion%patches(ipatch)

    distrib = ppatch%mixt%distrib

! ==============================================================================   
!   Select boundary type
! ==============================================================================   

    SELECT CASE ( ppatch%bcType )    

! ------------------------------------------------------------------------------
!     No-slip wall with imposed heat flux
! ------------------------------------------------------------------------------

      CASE ( bc_noslipwall_hflux )

! ----- Loop over faces and compute gradients ----------------------------------

        DO ifl = 1,ppatch%nBFaces
          c1 = ppatch%bf2c(ifl)

! ------- Get face geometry 
    
          nx = ppatch%fn(xcoord,ifl)
          ny = ppatch%fn(ycoord,ifl)
          nz = ppatch%fn(zcoord,ifl)       

! ------- Get face state 
   
          cond = tv(tvindxcond,c1)

! TEMPORARY
!          CALL RFLU_InterpCells2FacePatch(pRegion,pPatch,ifl, &
!               pRegion%mixt%tv(tvIndxCond:tvIndxCond,:), &
!               tvf(tvIndxCond:tvIndxCond))
! 
!          cond = tvf(tvIndxCond)
! END TEMPORARY
    
! ------- Compute and set gradient

          dtdn = -ppatch%mixt%vals(bcdat_noslip_q,distrib*ifl)/cond

          dtdx = dtdn*nx
          dtdy = dtdn*ny
          dtdz = dtdn*nz                  
    
! ------- Set gradient     

          ppatch%mixt%gradFace(xcoord,grbf_mixt_temp,ifl) = dtdx
          ppatch%mixt%gradFace(ycoord,grbf_mixt_temp,ifl) = dtdy
          ppatch%mixt%gradFace(zcoord,grbf_mixt_temp,ifl) = dtdz
        
        END DO ! ifl
        
! ------------------------------------------------------------------------------
!     Any other patches
! ------------------------------------------------------------------------------

      CASE default
        
    END SELECT ! pPatch%bcType
  END DO ! iPatches

! ******************************************************************************
! End
! ******************************************************************************

  CALL deregisterfunction(global)

END SUBROUTINE rflu_enforceheatflux







! ******************************************************************************
!
! Purpose: Compute viscous fluxes for actual faces.
!
! Description: None.
!
! Input: 
!  pRegion      Pointer to region
!  tv           Transport variables
!  tvIndxVisc   Index to viscosity entry in tv
!  tvIndxCond   Index to conductity entry in tv
!
! Output: None.
!
! Notes: None.
!
! ******************************************************************************

SUBROUTINE rflu_viscousfluxes(pRegion,tv,tvIndxVisc,tvIndxCond)
   
  USE rflu_modinterpolation, ONLY: rflu_interpcells2face 
   
  IMPLICIT NONE

! ******************************************************************************
! Declarations
! ******************************************************************************

! ==============================================================================  
! Arguments 
! ==============================================================================  

  INTEGER, INTENT(IN) :: tvindxcond,tvindxvisc
  REAL(RFREAL), DIMENSION(:,:) :: tv
  TYPE(t_region), POINTER :: pregion

! ==============================================================================  
! Locals
! ==============================================================================  

  INTEGER :: c1,c2,ifg
  REAL(RFREAL), PARAMETER :: two_thirds = 2.0_rfreal/3.0_rfreal
  REAL(RFREAL) :: beta,cond,divterm,dtdx,dtdy,dtdz,dudx,dudy,dudz,dvdx,dvdy, &
                  dvdz,dwdx,dwdy,dwdz,nm,nx,ny,nz,s11,s12,s13,s21,s22,s23, &
                  s31,s32,s33,u,v,visc,w
  REAL(RFREAL) :: fd(4)
  REAL(RFREAL) :: cvf(cv_mixt_xvel:cv_mixt_zvel),tvf(tvindxvisc:tvindxcond)
  REAL(RFREAL), DIMENSION(:,:), POINTER :: pcv,pdiss
  TYPE(t_global), POINTER :: global
  TYPE(t_grid), POINTER :: pgrid

! ******************************************************************************
! Start
! ******************************************************************************

  global => pregion%global

  CALL registerfunction(global,'RFLU_ViscousFluxes',&
  'RFLU_ModViscousFlux.F90')

! ******************************************************************************
! Set variables and pointers
! ******************************************************************************

  pgrid => pregion%grid
  pcv   => pregion%mixt%cv
  pdiss => pregion%mixt%diss
    
  beta = pregion%mixtInput%betrk(pregion%irkStep)  

! ******************************************************************************
! Check state of conserved state variable vector (defensive coding)
! ******************************************************************************

  IF ( pregion%mixt%cvState /= cv_mixt_state_duvwt ) THEN 
    CALL errorstop(global,err_cv_state_invalid,__line__)
  END IF ! pRegion%mixt%cvState
  
! ******************************************************************************  
! Loop over faces and compute viscous fluxes
! ******************************************************************************

  DO ifg = 1,pgrid%nFaces
    c1 = pgrid%f2c(1,ifg)
    c2 = pgrid%f2c(2,ifg)

! ==============================================================================   
!   Get face geometry
! ==============================================================================    
    
    nx = pgrid%fn(xcoord,ifg)
    ny = pgrid%fn(ycoord,ifg)
    nz = pgrid%fn(zcoord,ifg) 
    nm = pgrid%fn(xyzmag,ifg)        

! ==============================================================================   
!   Get face state
! ==============================================================================

    u = 0.5_rfreal*(pcv(cv_mixt_xvel,c1) + pcv(cv_mixt_xvel,c2))
    v = 0.5_rfreal*(pcv(cv_mixt_yvel,c1) + pcv(cv_mixt_yvel,c2))
    w = 0.5_rfreal*(pcv(cv_mixt_zvel,c1) + pcv(cv_mixt_zvel,c2))

    visc = 0.5_rfreal*(tv(tvindxvisc,c1) + tv(tvindxvisc,c2))    
    cond = 0.5_rfreal*(tv(tvindxcond,c1) + tv(tvindxcond,c2))

! TEMPORARY
!    CALL RFLU_InterpCells2Face(pRegion,ifg, &
!                               pRegion%mixt%cv(CV_MIXT_XVEL:CV_MIXT_ZVEL,:), &
!                               cvf(CV_MIXT_XVEL:CV_MIXT_ZVEL))
! 
!    CALL RFLU_InterpCells2Face(pRegion,ifg, &
!                               pRegion%mixt%tv(tvIndxVisc:tvIndxCond,:), &
!                               tvf(tvIndxVisc:tvIndxCond))
! 
!    u = cvf(CV_MIXT_XVEL)
!    v = cvf(CV_MIXT_YVEL)
!    w = cvf(CV_MIXT_ZVEL)
!    
!    visc = tvf(tvIndxVisc)
!    cond = tvf(tvIndxCond)
! END TEMPORARY

! ==============================================================================   
!   Get gradients
! ==============================================================================

    dudx = pregion%mixt%gradFace(xcoord,grf_mixt_xvel,ifg)
    dudy = pregion%mixt%gradFace(ycoord,grf_mixt_xvel,ifg)
    dudz = pregion%mixt%gradFace(zcoord,grf_mixt_xvel,ifg)
    
    dvdx = pregion%mixt%gradFace(xcoord,grf_mixt_yvel,ifg)
    dvdy = pregion%mixt%gradFace(ycoord,grf_mixt_yvel,ifg)
    dvdz = pregion%mixt%gradFace(zcoord,grf_mixt_yvel,ifg)

    dwdx = pregion%mixt%gradFace(xcoord,grf_mixt_zvel,ifg)
    dwdy = pregion%mixt%gradFace(ycoord,grf_mixt_zvel,ifg)
    dwdz = pregion%mixt%gradFace(zcoord,grf_mixt_zvel,ifg)

    dtdx = pregion%mixt%gradFace(xcoord,grf_mixt_temp,ifg)
    dtdy = pregion%mixt%gradFace(ycoord,grf_mixt_temp,ifg)
    dtdz = pregion%mixt%gradFace(zcoord,grf_mixt_temp,ifg)
  
! ==============================================================================   
!   Compute fluxes
! ==============================================================================
  
    divterm = two_thirds*(dudx + dvdy + dwdz)
  
    s11 = 2.0_rfreal*dudx - divterm
    s12 = dudy + dvdx
    s13 = dudz + dwdx 
  
    s21 = s12
    s22 = 2.0_rfreal*dvdy - divterm
    s23 = dvdz + dwdy  
  
    s31 = s13
    s32 = s23
    s33 = 2.0_rfreal*dwdz - divterm
  
    fd(1) = visc*(s11*nx + s12*ny + s13*nz)*nm
    fd(2) = visc*(s21*nx + s22*ny + s23*nz)*nm
    fd(3) = visc*(s31*nx + s32*ny + s33*nz)*nm
    
    fd(4) = u*fd(1) + v*fd(2) + w*fd(3) + cond*(dtdx*nx + dtdy*ny + dtdz*nz)*nm   

! ==============================================================================   
!   Accumulate into residual     
! ==============================================================================

    pdiss(cv_mixt_xmom,c1) = pdiss(cv_mixt_xmom,c1) + beta*fd(1)
    pdiss(cv_mixt_ymom,c1) = pdiss(cv_mixt_ymom,c1) + beta*fd(2)
    pdiss(cv_mixt_zmom,c1) = pdiss(cv_mixt_zmom,c1) + beta*fd(3)
    pdiss(cv_mixt_ener,c1) = pdiss(cv_mixt_ener,c1) + beta*fd(4)

    pdiss(cv_mixt_xmom,c2) = pdiss(cv_mixt_xmom,c2) - beta*fd(1)
    pdiss(cv_mixt_ymom,c2) = pdiss(cv_mixt_ymom,c2) - beta*fd(2)
    pdiss(cv_mixt_zmom,c2) = pdiss(cv_mixt_zmom,c2) - beta*fd(3)
    pdiss(cv_mixt_ener,c2) = pdiss(cv_mixt_ener,c2) - beta*fd(4)
  END DO ! ifg

! ******************************************************************************
! End
! ******************************************************************************

  CALL deregisterfunction(global)

END SUBROUTINE rflu_viscousfluxes









! ******************************************************************************
!
! Purpose: Compute viscous fluxes for actual faces on boundary patches.
!
! Description: None.
!
! Input: 
!  pRegion      Pointer to region
!  tv           Transport variables
!  tvIndxVisc   Index to viscosity entry in tv
!  tvIndxCond   Index to conductity entry in tv
!
! Output: None.
!
! Notes: None.
!
! ******************************************************************************

SUBROUTINE rflu_viscousfluxespatches(pRegion,tv,tvIndxVisc,tvIndxCond)
    
  USE rflu_modinterpolation, ONLY: rflu_interpcells2facepatch  
   
  IMPLICIT NONE

! ******************************************************************************
! Declarations
! ******************************************************************************

! ==============================================================================  
! Arguments 
! ==============================================================================  

  INTEGER, INTENT(IN) :: tvindxcond,tvindxvisc
  REAL(RFREAL), DIMENSION(:,:) :: tv
  TYPE(t_region), POINTER :: pregion

! ==============================================================================  
! Locals
! ==============================================================================  

  INTEGER :: c1,ifg,ifgbeg,ifgend,ifl,ipatch
  REAL(RFREAL), PARAMETER :: two_thirds = 2.0_rfreal/3.0_rfreal
  REAL(RFREAL) :: beta,cond,divterm,dtdx,dtdy,dtdz,dudx,dudy,dudz,dvdx,dvdy, &
                  dvdz,dwdx,dwdy,dwdz,icfref,ichref,nm,nx,ny,nz,rref,s11, &
                  s12,s13,s21,s22,s23,s31,s32,s33,visc,vref
  REAL(RFREAL) :: fd(4)
  REAL(RFREAL) :: tvf(tvindxvisc:tvindxcond)
  REAL(RFREAL), DIMENSION(:,:), POINTER :: pcv,pdiss
  TYPE(t_global), POINTER :: global
  TYPE(t_grid), POINTER :: pgrid
  TYPE(t_patch), POINTER :: ppatch

! ******************************************************************************
! Start
! ******************************************************************************

  global => pregion%global

  CALL registerfunction(global,'RFLU_ViscousFluxesPatches',&
  'RFLU_ModViscousFlux.F90')

! ******************************************************************************
! Set variables and pointers
! ******************************************************************************

  pgrid => pregion%grid
  pcv   => pregion%mixt%cv
  pdiss => pregion%mixt%diss
  
  beta = pregion%mixtInput%betrk(pregion%irkStep)  

  rref = global%refDensity
  vref = global%refVelocity
  
  icfref = 2.0_rfreal/(rref*vref*vref)
  ichref = 2.0_rfreal/(rref*vref*vref*vref)

! ******************************************************************************
! Check state of conserved state variable vector (defensive coding)
! ******************************************************************************

  IF ( pregion%mixt%cvState /= cv_mixt_state_duvwt ) THEN 
    CALL errorstop(global,err_cv_state_invalid,__line__)
  END IF ! pRegion%mixt%cvState
  
! ******************************************************************************  
! Loop over patches 
! ******************************************************************************

  DO ipatch = 1,pgrid%nPatches
    ppatch => pregion%patches(ipatch)

! ==============================================================================   
!   Select boundary type
! ==============================================================================   

    SELECT CASE ( ppatch%bcType )    

! ------------------------------------------------------------------------------
!     No-slip wall
! ------------------------------------------------------------------------------

      CASE ( bc_noslipwall_hflux,bc_noslipwall_temp )

! ----- Loop over faces and compute gradients ----------------------------------

        DO ifl = 1,ppatch%nBFaces
          c1 = ppatch%bf2c(ifl)

! ------- Get face geometry 
    
          nx = ppatch%fn(xcoord,ifl)
          ny = ppatch%fn(ycoord,ifl)
          nz = ppatch%fn(zcoord,ifl) 
          nm = ppatch%fn(xyzmag,ifl)        

! ------- Get face state 

          visc = tv(tvindxvisc,c1)   
          cond = tv(tvindxcond,c1)

! TEMPORARY
!          CALL RFLU_InterpCells2FacePatch(pRegion,pPatch,ifl, &
!               pRegion%mixt%tv(tvIndxVisc:tvIndxCond,:), &
!               tvf(tvIndxVisc:tvIndxCond))
! 
!          visc = tvf(tvIndxVisc)   
!          cond = tvf(tvIndxCond)
! END TEMPORARY

! ------- Get gradients 

          dudx = ppatch%mixt%gradFace(xcoord,grbf_mixt_xvel,ifl)
          dudy = ppatch%mixt%gradFace(ycoord,grbf_mixt_xvel,ifl)
          dudz = ppatch%mixt%gradFace(zcoord,grbf_mixt_xvel,ifl)

          dvdx = ppatch%mixt%gradFace(xcoord,grbf_mixt_yvel,ifl)
          dvdy = ppatch%mixt%gradFace(ycoord,grbf_mixt_yvel,ifl)
          dvdz = ppatch%mixt%gradFace(zcoord,grbf_mixt_yvel,ifl)

          dwdx = ppatch%mixt%gradFace(xcoord,grbf_mixt_zvel,ifl)
          dwdy = ppatch%mixt%gradFace(ycoord,grbf_mixt_zvel,ifl)
          dwdz = ppatch%mixt%gradFace(zcoord,grbf_mixt_zvel,ifl)

          dtdx = ppatch%mixt%gradFace(xcoord,grbf_mixt_temp,ifl)
          dtdy = ppatch%mixt%gradFace(ycoord,grbf_mixt_temp,ifl)
          dtdz = ppatch%mixt%gradFace(zcoord,grbf_mixt_temp,ifl)
   
! ------- Compute fluxes 
  
          divterm = two_thirds*(dudx + dvdy + dwdz)

          s11 = 2.0_rfreal*dudx - divterm
          s12 = dudy + dvdx
          s13 = dudz + dwdx 

          s21 = s12
          s22 = 2.0_rfreal*dvdy - divterm
          s23 = dvdz + dwdy  

          s31 = s13
          s32 = s23
          s33 = 2.0_rfreal*dwdz - divterm

          fd(1) = visc*(s11 *nx + s12 *ny + s13 *nz)
          fd(2) = visc*(s21 *nx + s22 *ny + s23 *nz)
          fd(3) = visc*(s31 *nx + s32 *ny + s33 *nz)
          fd(4) = cond*(dtdx*nx + dtdy*ny + dtdz*nz)
    
! ------- Set friction and heat-transfer coefficients

          ppatch%cf(xcoord,ifl) = -icfref*fd(1)
          ppatch%cf(ycoord,ifl) = -icfref*fd(2)
          ppatch%cf(zcoord,ifl) = -icfref*fd(3)       
          
          ppatch%ch(ifl) = ichref*fd(4)               
    
! ------- Accumulate into residual     

          pdiss(cv_mixt_xmom,c1) = pdiss(cv_mixt_xmom,c1) + beta*fd(1)*nm
          pdiss(cv_mixt_ymom,c1) = pdiss(cv_mixt_ymom,c1) + beta*fd(2)*nm
          pdiss(cv_mixt_zmom,c1) = pdiss(cv_mixt_zmom,c1) + beta*fd(3)*nm
          pdiss(cv_mixt_ener,c1) = pdiss(cv_mixt_ener,c1) + beta*fd(4)*nm      
        
        END DO ! ifl
        
! ------------------------------------------------------------------------------
!     Any other patches
! ------------------------------------------------------------------------------

      CASE default
        
    END SELECT ! pPatch%bcType
  END DO ! iPatches

! ******************************************************************************
! End
! ******************************************************************************

  CALL deregisterfunction(global)

END SUBROUTINE rflu_viscousfluxespatches








  
! ******************************************************************************
! End
! ******************************************************************************
  
END MODULE rflu_modviscousflux


! ******************************************************************************
!
! RCS Revision history:
!
! $Log: RFLU_ModViscousFlux.F90,v $
! Revision 1.10  2008/12/06 08:44:24  mtcampbe
! Updated license.
!
! Revision 1.9  2008/11/19 22:17:35  mtcampbe
! Added Illinois Open Source License/Copyright
!
! Revision 1.8  2007/08/07 15:18:30  rfiedler
! Bug fix: Wrong loop limit in RFLU_EnforceHeatFlux.
!
! Revision 1.7  2006/08/19 15:39:19  mparmar
! Renamed bGradFace, removed bf2bg, used GRBF_ for boundary grad arrays
!
! Revision 1.6  2006/04/07 15:19:21  haselbac
! Removed tabs
!
! Revision 1.5  2005/10/16 18:03:37  haselbac
! Bug fix: Missing definition of pGrid
!
! Revision 1.4  2005/10/16 17:16:42  haselbac
! Bug fix: distrib set in wrong place
!
! Revision 1.3  2005/10/14 14:08:41  haselbac
! Added RFLU_EnforceHeatFlux - temporary until have proper constr reconstr
!
! Revision 1.2  2005/10/05 14:14:04  haselbac
! No longer distinguish between isothermal and adiabatic walls
!
! Revision 1.1  2005/05/16 20:36:29  haselbac
! Initial revision
!
! ******************************************************************************