TURB_fluLesBMij.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: Compute the mij stress appearing in the dynamic LES models.
!
! Description: We determine mij as appearing in the dynamic formulation:
!              M_{ij}=-test[bar(rho)](kappa \Delta)^2 |S(v)|S_{ij(v)
!                     +test[bar(rho)\Delta^{2}|S(u)|S_{ij}(u)]
!              in which v=test[bar(rho ui)]/test[bar(rho)] where test denotes
!              the test-filter with twice the filter-width. M_ij is computed
!              on the cell faces. S_ij(u) and Sij(test[u]) are already 
!              available, respectively, in mSij for mixture field and fSij for 
!              filtered field at cell faces. Moreover, we use the fact that we 
!              know 1/test[bar(rho)] (rhoFTb) and 1/[bar(rho)] (rhoF) at faces.
!
! Input: region = data of current region
!        ijk    = ijk-face is being treated
!
! Output: mij = resulting stress tensor in dynamic LES models
!
! Notes: This routine relevant only for physical boundaries of FLU regions.
!
!******************************************************************************
!
! $Id: TURB_fluLesBMij.F90,v 1.3 2008/12/06 08:44:44 mtcampbe Exp $
!
! Copyright: (c) 2001 by the University of Illinois
!
!******************************************************************************

SUBROUTINE turb_flulesbmij( region,ijk )

  USE moddatatypes
  USE moddatastruct, ONLY : t_region
  USE modglobal, ONLY     : t_global
  USE modbndpatch, ONLY   : t_patch
!  USE RFLU_ModInterpolation, ONLY : RFLU_InterpFaces2FacesPatches
  USE moderror
  USE turb_modparameters
  IMPLICIT NONE

! ... parameters
  TYPE(t_region), POINTER :: region
  INTEGER :: ijk

! ... loop variables
  INTEGER :: i, j, k, ijkn, ijkc0, ijkc1, ipatch

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

  INTEGER      :: ibn,ien,idbeg,idend,npatches,nbfaces,tndel(diri:dirk)
  REAL(RFREAL) :: kappa2,rkappa2,two3rd
  REAL(RFREAL) :: delfac2,delfkap2,forterm
  REAL(RFREAL), POINTER :: bfvol(:),bfsij(:,:),bfvar(:,:),bffvar(:,:)
  REAL(RFREAL), POINTER :: bsij(:,:),bmij(:,:)
  REAL(RFREAL), POINTER :: field(:,:),tfield(:,:)
  REAL(RFREAL), ALLOCATABLE :: modstrain(:)

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

  rcsidentstring = '$RCSfile: TURB_fluLesBMij.F90,v $'

  global => region%global
  CALL registerfunction( global,'TURB_FluLesBMij',&
  'TURB_fluLesBMij.F90' )

! get indices and pointers ----------------------------------------------------

  npatches = region%grid%nPatches
  nbfaces  = 0

  DO ipatch = 1,npatches
    patch   => region%patches(ipatch)
    nbfaces =  nbfaces + patch%nBTris + patch%nBQuads
  END DO ! iPatch
  ibn = 1
  ien = nbfaces

  bfvar  => region%turb%bfVar
  bffvar => region%turb%bffVar
  bmij   => region%turb%bMij
  bfvol  => region%turb%bfVolI
  bsij   => region%turb%bmIsij
  bfsij  => region%turb%bfISij

! allocate work arrays used within this scope

  ALLOCATE( modstrain(ibn:ien) )
  ALLOCATE( field(e11:e33,ibn:ien),tfield(e11:e33,ibn:ien) )

! test filter width is twice bar filter width

  tndel(diri) = 2*region%turbInput%filterWidth(diri)

! filterwidth on the fg-level (test-bar) is sqrt(5) times larger than 
! on the f-level (bar), i.e. kappa=sqrt(5)

  kappa2 = 5._rfreal
  rkappa2= 1._rfreal/kappa2
  two3rd = 2._rfreal/3._rfreal

! determine the bar filter width delta which appears in the dynamic model; 
! this factor in square is defined as (kappa*delfac)^2 and stored in delFKap2 

  delfac2  = region%turbInput%delFac2
  delfkap2 = delfac2*kappa2

! first calculate test[bar(rho)](\kappa \Delta)^2|S(v)|Sij(v) at cell faces; 
! we obtain Sij(v) from array fSij; 

  DO ijkn=ibn,ien

! - generate the 6-th component of fSij using the tracelessness of S_{ij} 

    bfsij(e33,ijkn)=-(bfsij(e11,ijkn)+bfsij(e22,ijkn))

! - get modulus of Sij(v): |S(v)|
    modstrain(ijkn)=sqrt(bfsij(e11,ijkn)*bfsij(e11,ijkn) &
                        +bfsij(e22,ijkn)*bfsij(e22,ijkn) &
                        +bfsij(e11,ijkn)*bfsij(e22,ijkn) &
                        +bfsij(e12,ijkn)*bfsij(e12,ijkn) &
                        +bfsij(e13,ijkn)*bfsij(e13,ijkn) &
                        +bfsij(e23,ijkn)*bfsij(e23,ijkn))

! - in turn, (kappa*delta)^2 in fg-level is stored in tField

    tfield(e11,ijkn) = delfkap2*bfvol(ijkn)**two3rd

! - fg-level contribution to M_{ij} can now be computed since
!   we already have test(bar[rho]) available at cell faces,
!   in fact we have 1/test(bar[rho]), so we devide instead of multiply;
!   install the first term of mij; the space of fSij in used in the process

    forterm       = -tfield(e11,ijkn)/bffvar(cv_turb_dens,ijkn)* &
                     modstrain(ijkn)

    bmij(e11,ijkn) = forterm*bfsij(e11,ijkn)
    bmij(e12,ijkn) = forterm*bfsij(e12,ijkn)
    bmij(e13,ijkn) = forterm*bfsij(e13,ijkn)
    bmij(e22,ijkn) = forterm*bfsij(e22,ijkn)
    bmij(e23,ijkn) = forterm*bfsij(e23,ijkn)
    bmij(e33,ijkn) = forterm*bfsij(e33,ijkn)

! - as final step we determine contribution of test-filtered f-level term;
!   note we have bar[rho] available at cell faces bRho stored in fVar; 
!   in fact we have 1/bRho in fVar(CV_TURB_DENS,:), so we devide by bRho
!   instead of multiply and compute modulus of Sij: |S(u)|

    modstrain(ijkn) = sqrt(bsij(e11,ijkn)*bsij(e11,ijkn) &
                          +bsij(e22,ijkn)*bsij(e22,ijkn) &
                          +bsij(e11,ijkn)*bsij(e22,ijkn) &
                          +bsij(e12,ijkn)*bsij(e12,ijkn) &
                          +bsij(e13,ijkn)*bsij(e13,ijkn) &
                          +bsij(e23,ijkn)*bsij(e23,ijkn))

! - finish computation of mij; filterwidth^2 at f-level is tField/kappa2;
!   build terms need to be filtered

    forterm = tfield(e11,ijkn)/bfvar(cv_turb_dens,ijkn)*rkappa2* &
              modstrain(ijkn)

    field(e11,ijkn) =  forterm*bsij(e11,ijkn)
    field(e12,ijkn) =  forterm*bsij(e12,ijkn)
    field(e13,ijkn) =  forterm*bsij(e13,ijkn)
    field(e22,ijkn) =  forterm*bsij(e22,ijkn)
    field(e23,ijkn) =  forterm*bsij(e23,ijkn)
    field(e33,ijkn) = -forterm*(bsij(e11,ijkn)+bsij(e22,ijkn))
  ENDDO

! perform test filtering:

  idbeg = e11
  idend = e33

!  CALL RFLU_InterpFaces2FacesPatches( region,field,tField )

! and complete components of mij:

  DO ijkn=ibn,ien
    bmij(e11,ijkn) = bmij(e11,ijkn) + tfield(e11,ijkn)
    bmij(e12,ijkn) = bmij(e12,ijkn) + tfield(e12,ijkn)
    bmij(e13,ijkn) = bmij(e13,ijkn) + tfield(e13,ijkn)
    bmij(e22,ijkn) = bmij(e22,ijkn) + tfield(e22,ijkn)
    bmij(e23,ijkn) = bmij(e23,ijkn) + tfield(e23,ijkn)
    bmij(e33,ijkn) = bmij(e33,ijkn) + tfield(e33,ijkn)
  ENDDO

! deallocate temporary arrays

  DEALLOCATE( modstrain,field,tfield )

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

  CALL deregisterfunction( global )

END SUBROUTINE turb_flulesbmij

!******************************************************************************
!
! RCS Revision history:
!
! $Log: TURB_fluLesBMij.F90,v $
! Revision 1.3  2008/12/06 08:44:44  mtcampbe
! Updated license.
!
! Revision 1.2  2008/11/19 22:17:56  mtcampbe
! Added Illinois Open Source License/Copyright
!
! Revision 1.1  2004/05/28 02:03:58  wasistho
! update unstructured grid LES
!
!
!
!******************************************************************************