BcondInflowPerf.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: set inflow boundary condition for one cell.
!
! Description: The subsonic boundary condition is based on the extrapolation of 
!   the Riemann invariant from the interior field (Holmes, D.G.: Inviscid 2D 
!   Solutions on Unstructured, Adaptive Grids. VKI Lecture Series 1989-06, 1989). 
!   The supersonic inflow boundary condition computes the conservative variables
!   from given velocity components, density and pressure.
!
! Input: bcOptType  = boundary treatment: subsonic, supersonic, or mixed
!        bcOptFixed = whether _computed_ inflow angle should be fixed or not
!        ptot       = given total pressure
!        ttot       = given total temperature
!        betah      = given inlet angle wrp. to y-axis
!        betav      = given inlet angle wrp. to z-axis
!        sx/y/zn    = components of normalized face vector (outward facing)
!        cpgas      = specific heat at constant pressure (boundary cell)
!        mm         = molecular mass at boundary cell
!        rl         = given density
!        ru/v/wl    = given velocity components
!
! Output: rr      = density at boundary
!         ru/v/wr = density * velocity components at boundary
!         rer     = density * total internal energy at boundary
!         pr      = pressure at boundary
!
! Notes: 
!   1. This condition is valid only for thermally and calorically perfect  
!      gas.
!   2. Important to avoid division by MakeNonZero(sl) when computing eta 
!      because that computation can become undefined for reservoir inflow
!      conditions, i.e., for a vanishing velocity vector. 
!
!******************************************************************************
!
! $Id: BcondInflowPerf.F90,v 1.3 2008/12/06 08:44:08 mtcampbe Exp $
!
! Copyright: (c) 2002-2004 by the University of Illinois
!
!******************************************************************************

SUBROUTINE bcondinflowperf(bcOptType,bcOptFixed,ptot,ttot,betah,betav,mach, & 
                           sxn,syn,szn,cpgas,mm,rl,rul,rvl,rwl,rr,rur,rvr, & 
                           rwr,rer,pr)

  USE moddatatypes
  USE modparameters
  USE modtools, ONLY     : makenonzero
  USE modinterfaces, ONLY: mixtperf_c_co2guvw, mixtperf_c_dgp, mixtperf_c_grt, &
        mixtperf_co2_cguvw, mixtperf_c2_grt,mixtperf_d_prt, & 
        mixtperf_eo_dgpuvw, mixtperf_eo_dgpvm, mixtperf_g_cpr, &
        mixtperf_p_gmapo, mixtperf_p_gpotto, mixtperf_po_gptto, & 
        mixtperf_po_cgpuvw, mixtperf_r_m, mixtperf_t_cgr, mixtperf_t_gmato, & 
        mixtperf_vm_c2co2g

  IMPLICIT NONE

! ... parameters
  INTEGER, INTENT(IN) :: bcoptfixed,bcopttype

  REAL(RFREAL), INTENT(IN) :: betah, betav, cpgas, mach, mm, sxn, syn, szn, &
                              ptot, rl, rul, rvl, rwl, ttot
  REAL(RFREAL), INTENT(OUT) :: rer, rr, rur, rvr, rwr, pr

! ... local variables
  REAL(RFREAL) :: al, ar, a02, cp, disc, eta, g, gm1, igm1, ql, rgas, rm, &
                  sl, sr , tr, ul, ur, vl, vr, wl, wr

!******************************************************************************
! gas properties

  rgas = mixtperf_r_m(mm)
  g    = mixtperf_g_cpr(cpgas,rgas)

! subsonic or mixed -----------------------------------------------------------   

  IF ( bcopttype == bcopt_subsonic .OR. bcopttype == bcopt_mixed ) THEN     
    gm1  = g - 1.0_rfreal
    igm1 = 1.0_rfreal/gm1

    ul = rul/rl
    vl = rvl/rl
    wl = rwl/rl

    a02 = mixtperf_c2_grt(g,rgas,ttot)

    al = mixtperf_c_co2guvw(a02,g,ul,vl,wl) ! make al consistent with a02
    ql = ul*sxn + vl*syn + wl*szn

! - subsonic

    IF ( abs(ql) < al ) THEN
      sl = sqrt(ul*ul + vl*vl + wl*wl)
      
      IF ( bcoptfixed == bcopt_fixed_no ) THEN 
        IF ( sl > 1.0e-6_rfreal ) THEN ! Avoid ill-defined angle computation
          eta = ql/sl        
        ELSE 
          eta = -1.0_rfreal
        END IF ! sl
      ELSE 
        eta = -1.0_rfreal
      END IF ! bcOptFixed

      rm   = al - 0.5_rfreal*gm1*ql
      disc = 0.5_rfreal*gm1*eta**2* & 
             (a02/(rm*rm)*(1.0_rfreal + 0.5_rfreal*gm1*eta**2) - 1.0_rfreal)     
        
      IF ( disc < 0.0_rfreal ) THEN ! discriminant cannot be negative
        ar = sqrt(a02)
        tr = ttot
        pr = ptot
        sr = 0.0_rfreal
      ELSE
        ar = rm/(1.0_rfreal + 0.5_rfreal*gm1*eta*eta)*(1.0_rfreal+sqrt(disc))                   
                      
        tr = mixtperf_t_cgr(ar,g,rgas)
        pr = mixtperf_p_gpotto(g,ptot,tr,ttot)
        sr = mixtperf_vm_c2co2g(ar*ar,a02,g)
      END IF ! disc    
                  
      rr = mixtperf_d_prt( pr,rgas,tr )
      ur = sr*cos(betah)*cos(betav)
      vr = sr*sin(betah)
      wr = sr*cos(betah)*sin(betav)

      rer = rr*mixtperf_eo_dgpvm(rr,g,pr,sr)
      rur = rr*ur
      rvr = rr*vr
      rwr = rr*wr

! - supersonic

    ELSE
      IF ( bcopttype == bcopt_mixed ) THEN
        pr = mixtperf_p_gmapo(g,mach,ptot)
        tr = mixtperf_t_gmato(g,mach,ttot)
        rr = mixtperf_d_prt(pr,rgas,tr)
        ar = mixtperf_c_grt(g,rgas,tr)

        ur = mach*ar*cos(betah)*cos(betav)
        vr = mach*ar*sin(betah)
        wr = mach*ar*cos(betah)*sin(betav)

        rer = rr*mixtperf_eo_dgpuvw(rr,g,pr,ur,vr,wr)
        rur = rr*ur
        rvr = rr*vr
        rwr = rr*wr
      END IF ! bcOptType
    END IF ! ql < al

! supersonic ------------------------------------------------------------------

  ELSE ! bcOptType == BCOPT_SUPERSONIC
    pr = mixtperf_p_gmapo(g,mach,ptot)
    tr = mixtperf_t_gmato(g,mach,ttot)
    rr = mixtperf_d_prt(pr,rgas,tr)
    ar = mixtperf_c_grt(g,rgas,tr)

    ur = mach*ar*cos(betah)*cos(betav)
    vr = mach*ar*sin(betah)
    wr = mach*ar*cos(betah)*sin(betav) 

    rer = rr*mixtperf_eo_dgpuvw(rr,g,pr,ur,vr,wr)
    rur = rr*ur
    rvr = rr*vr
    rwr = rr*wr
  END IF ! bcOptType

END SUBROUTINE bcondinflowperf

!******************************************************************************
!
! RCS Revision history:
!
! $Log: BcondInflowPerf.F90,v $
! Revision 1.3  2008/12/06 08:44:08  mtcampbe
! Updated license.
!
! Revision 1.2  2008/11/19 22:17:22  mtcampbe
! Added Illinois Open Source License/Copyright
!
! Revision 1.1  2004/12/01 16:47:56  haselbac
! Initial revision after changing case
!
! Revision 1.7  2004/01/29 22:52:36  haselbac
! Added bcOptFixed, fixed bug, clean-up
!
! Revision 1.6  2003/12/04 03:22:56  haselbac
! Fixed bug in formulation, added partial fix for eta
!
! Revision 1.5  2003/11/20 16:40:35  mdbrandy
! Backing out RocfluidMP changes from 11-17-03
!
! Revision 1.2  2002/06/22 00:49:50  jblazek
! Modified interfaces to BC routines.
!
! Revision 1.1  2002/06/10 21:19:34  haselbac
! Initial revision
!
!******************************************************************************