Rocburn_1D_ZN.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.                          *
! *********************************************************************
  MODULE m_rocburn_1d_zn
! ----------------------------------------------------------------
!
! This module provides the following subroutines for ROCBURN_1D_ZN  
! to support ROCBURN_2D:
!
!     INITIALIZE_0D
!     INITIALIZE_1D
!     GET_FILM_COEFF_1D
!     GET_TIME_STEP_1D
!     GET_BURNING_RATE_1D
!
!
! The combustion model used in this module is Zeldovich-Novozhilov 
! model implemented by J. Weber, K.C. Tang, and M. Q. Brewster.
!
!   Creation Date    :  Sep. 2, 2002
!
!   Modifications    :
!
!    No.     Date         Authors       Description
!
! ----------------------------------------------------------------
!
!

  USE m_rocburn_zn_global_data

  IMPLICIT NONE

  include 'mpif.h'


! ----------------------------------------------------------------
! local variables

  CONTAINS

!
!     ------------------------------------------------------------------------
!                             INTERNAL PROCEDURES
!     ------------------------------------------------------------------------

  SUBROUTINE check_alloc( ierr)
    INTEGER, INTENT(IN) :: ierr
    IF(ierr /= 0) THEN
       print *, "ROCBRUN_ZN ERROR: unable to allocate memory"
       CALL mpi_abort( mpi_comm_world, -1)
    ENDIF
  END SUBROUTINE check_alloc

  SUBROUTINE initialize_0d(G_ZN, comm, Indir, nxmax, To_read)
      TYPE(g_burn_1d), POINTER :: g_zn
      INTEGER, INTENT(IN)      :: comm
      CHARACTER(*), INTENT(IN) :: indir
      INTEGER, INTENT(OUT)     :: nxmax
      REAL(DBL), INTENT (OUT)  :: to_read

! ---------------------------------------------------------------------------
!     local variables
!


      INTEGER :: ierror

      INTERFACE 
        SUBROUTINE zn_input_0d(G_ZN, Indir)
          USE m_rocburn_zn_global_data
          TYPE(g_burn_1d), POINTER :: g_zn
          CHARACTER(*), INTENT(IN)  :: indir
        END SUBROUTINE zn_input_0d

        SUBROUTINE zn_gen_grid(G_ZN, gridtype, numx, x, z, zx, zxx)
          USE m_rocburn_zn_global_data
          TYPE(g_burn_1d), POINTER :: g_zn
          INTEGER, INTENT(IN) :: gridtype, numx
          REAL(DBL), INTENT (OUT) :: x(:), z(:), zx(:), zxx(:)
        END SUBROUTINE zn_gen_grid
      END INTERFACE

      ALLOCATE(g_zn)

!
!     read Rocburn_ZN 0d data
!

      g_zn%MPI_COMM_ROCBURN     = comm
      CALL mpi_comm_rank(comm, g_zn%rank, ierror)

      CALL zn_input_0d(g_zn, trim(indir))

      nxmax   = g_zn%nxmax
      to_read = g_zn%To

!
!     allocate memory for variables related to grid generation
!
      ALLOCATE(g_zn%x(1:g_zn%nxmax), stat=ierror); CALL check_alloc(ierror)
      ALLOCATE(g_zn%z(1:g_zn%nxmax), stat=ierror); CALL check_alloc(ierror)
      ALLOCATE(g_zn%zx(1:g_zn%nxmax), stat=ierror); CALL check_alloc(ierror)
      ALLOCATE(g_zn%zxx(1:g_zn%nxmax), stat=ierror); CALL check_alloc(ierror)

!
!     grid generation
!
      CALL zn_gen_grid(g_zn, g_zn%igrid, g_zn%nx,                  &
                       g_zn%x, g_zn%z, g_zn%zx, g_zn%zxx)

      RETURN

    END SUBROUTINE initialize_0d
! ============================================================================






    SUBROUTINE initialize_1d( G_ZN, bflag, P_mks, To, rhoc_mks, p_coor, rb_mks, Toa, fr, Tn, Tflame)


      IMPLICIT NONE

      TYPE(g_burn_1d), POINTER :: g_zn
      INTEGER, INTENT(INOUT)  :: bflag
      REAL(DBL), INTENT (IN)  :: p_mks, to, rhoc_mks, p_coor(3)
      REAL(DBL), INTENT (OUT) :: rb_mks, toa, fr
      REAL(DBL), INTENT (OUT) :: tn(:)
      REAL(DBL), INTENT (OUT) :: tflame
     
      INTERFACE
        SUBROUTINE zn_sswsb(G_ZN, P, qr, To, rhoc, rb, Ts, Tf, fr, Tn)
          USE m_rocburn_zn_global_data
          TYPE(g_burn_1d), POINTER :: g_zn
          REAL(DBL), INTENT(IN) :: p, qr, to, rhoc
          REAL(DBL), INTENT(OUT) :: rb, ts, tf, fr
          REAL(DBL), INTENT(OUT) :: tn(:)
        END SUBROUTINE zn_sswsb
      END INTERFACE

! ---------------------------------------------------------------------------
!     local variables
      REAL(DBL) :: p, rhoc, rb, ts, tf
      REAL(DBL) :: qr

! ---------------------------------------------------------------------------

      IF( (bflag == 0).OR.(g_zn%ign_flag == 1) ) THEN
!
!        not burning or ignition simulation required
!        set condensed phase temperature profile to initial temperature
!
         tn = g_zn%To

      ELSE
!
!        burning from onset and no ignition simulation required, 
!        calculate steady state solution and use as initial condition
!        
         p    = p_mks *9.869232667e-6           ! Pa to atm
         rhoc = rhoc_mks*1.0e-3                 ! Kg/m^3 to g/cm^3
         qr   = 0.0

         CALL zn_sswsb(g_zn, p, qr  , g_zn%To, rhoc, rb, ts, tf, fr, tn)

         rb_mks = rb*0.01                       ! cm/s to m/s

      END IF

      tflame = g_zn%Tf_adiabatic
      toa    = g_zn%To

      RETURN

    END SUBROUTINE initialize_1d
!! ============================================================================






    SUBROUTINE get_film_coeff_1d(G_ZN, p_coor, Ts, T_euler, P, qc, qcPrime)
       TYPE (g_burn_1d), POINTER :: g_zn

       REAL(DBL), INTENT (IN)      :: p_coor(3), ts, t_euler, p
       REAL(DBL), INTENT (OUT)     :: qc, qcprime  ! do not change qcPrime
!
!      place holder for calculating convective heat flux qc for
!      Rocburn_1D_ZN
!
!      currently not available; set qc to 0
!
       qc = 0.0
       qcprime = 0.0

       RETURN

    END SUBROUTINE get_film_coeff_1d
! ============================================================================






    SUBROUTINE get_time_step_1d(G_ZN, rb, Toa, dt_max)
      TYPE (g_burn_1d), POINTER :: g_zn

       REAL(DBL), INTENT (IN)      :: rb, toa
       REAL(DBL), INTENT (OUT)     :: dt_max

! ---------------------------------------------------------------------------
!     local variables

      REAL(DBL) :: dt_c
! ---------------------------------------------------------------------------

      IF( abs(toa - g_zn%To) >= 0.9*g_zn%To) THEN
         dt_c  = 0.1*g_zn%alfac/(rb*rb*1.0e4)  ! G_ZN%alfac in cm^2/sec, rb in m/s
         dt_max= min(g_zn%delt_max, dt_c)      ! dt_max in sec
      ELSE
         dt_max= g_zn%delt_max                 ! dt_max in sec
      END IF

      RETURN

    END SUBROUTINE get_time_step_1d
! ============================================================================







    SUBROUTINE get_burning_rate_1d( G_ZN, delt, P_mks, To, Tn,   &
              qc_mks, qc_old_mks, qr_mks, qr_old_mks, rhoc_mks, &
              toa, rb_mks, fr, bflag, tnp1, tflame)
      TYPE (g_burn_1d), POINTER :: g_zn

      REAL(DBL), INTENT (IN)      :: delt, p_mks, to
      REAL(DBL), INTENT (IN)      :: tn(:)
      REAL(DBL), INTENT (IN)      :: qc_mks, qc_old_mks, qr_mks, qr_old_mks
      REAL(DBL), INTENT (IN)      :: rhoc_mks
      REAL(DBL), INTENT (INOUT)   :: toa, rb_mks, fr
      INTEGER,   INTENT (INOUT)   :: bflag
      REAL(DBL), INTENT (OUT)     :: tnp1(:)
      REAL(DBL), INTENT (OUT)     :: tflame

      INTERFACE
        SUBROUTINE zn_sswsb(G_ZN, P, qr, To, rhoc, rb, Ts, Tf, fr, Tn)
          USE m_rocburn_zn_global_data
          TYPE(g_burn_1d), POINTER :: g_zn
          REAL(DBL), INTENT(IN) :: p, qr, to, rhoc
          REAL(DBL), INTENT(OUT) :: rb, ts, tf, fr
          REAL(DBL), INTENT(OUT) :: tn(:)
        END SUBROUTINE zn_sswsb

        SUBROUTINE zn_calc_burning_rate(G_ZN, delt, P, qr, To, rhoc, qr_old, fr_old,  &
                                        toa, rb, ts, fr, tn, tnp1)
          USE m_rocburn_zn_global_data
          TYPE(g_burn_1d), POINTER :: g_zn
          REAL(DBL), INTENT(IN)  :: delt, p, qr, to, rhoc
          REAL(DBL), INTENT(IN)  :: qr_old, fr_old, toa
          REAL(DBL), INTENT(OUT) :: rb, ts, fr
          REAL(DBL), INTENT(IN)  :: tn(:)
          REAL(DBL), INTENT(OUT) :: tnp1(:)
        END SUBROUTINE zn_calc_burning_rate
      END INTERFACE


! ---------------------------------------------------------------------------
!     local variables
      REAL(DBL) :: p, qc, qc_old, qr, qr_old, rhoc, rb
      REAL(DBL) :: ts, tf
! ---------------------------------------------------------------------------


      tflame = g_zn%Tf_adiabatic

      IF(bflag/=0) THEN

        IF(tn(1) > g_zn%To) THEN
!
!         propellant burning already, calculate burning rate using ZN_cal_burning_rate
!
          p      = p_mks *9.869232667e-6           ! Pa to atm
          rb     = rb_mks * 100.0                  ! m/s to cm/s
          qr     = qr_mks * 0.2388459e-4           ! W/m^2 to cal/cm^2/s  1 J = 0.2388459 cal
          qc     = qc_mks *  0.2388459e-4          ! W/m^2 to cal/cm^2/s
          qr_old = qr_old_mks * 0.2388459e-4       ! W/m^2 to cal/cm^2/s
          qc_old = qc_old_mks *  0.2388459e-4      ! W/m^2 to cal/cm^2/s
          rhoc   = rhoc_mks*1.0e-3                 ! Kg/m^3 to g/cm^3
          ts     = tn(1)

          CALL zn_calc_burning_rate(g_zn, delt, p, qr, to, rhoc, qr_old, fr,     &
                               toa, rb, ts, fr, tn, tnp1)


          rb_mks = rb*0.01                         ! cm/s to m/s

        ELSE 
!
!         propellant burning for the first time , set initial condition using ZN_ssWSB
!
          p    = p_mks *9.869232667e-6             ! Pa to atm
          rhoc = rhoc_mks*1.0e-3                   ! Kg/m^3 to g/cm^3
          qr     = qr_mks * 0.2388459e-4           ! W/m^2 to cal/cm^2/s  1 J = 0.2388459 cal

          CALL zn_sswsb(g_zn, p, qr, g_zn%To, rhoc, rb, ts, tf, fr, tnp1)

          rb_mks = rb*0.01                         ! cm/s to m/s
        END IF


      ELSE
!
!       propellant not burning yet, check for ignition simulation requirement
!
        IF(g_zn%ign_flag == 1) THEN
!
!         ignition simulation required
!         
!         place holder for ignition model
!
          WRITE(*,*) 'ROCBURN_ZN: rank=',g_zn%rank
          WRITE(*,*) '  Error: igntion model not ready'
          WRITE(*,*) '  job aborted'
          CALL mpi_abort( mpi_comm_world, -1)
          stop

        ELSE
!
!         ignition simulation not required
!
          bflag = 1
          RETURN

        END IF

      END IF

      RETURN

    END SUBROUTINE get_burning_rate_1d

!***************************************************************************
    SUBROUTINE finalize_0d(G_ZN)
 
      TYPE (g_burn_1d), POINTER :: g_zn

      DEALLOCATE( g_zn)

    END SUBROUTINE finalize_0d


  END MODULE m_rocburn_1d_zn