Rocstar-legacy/ZN__calc__burning__rate_8f90_source.html

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

 ! * 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.                          *

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

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

 !

 !   SUBROUTINE  : ZN_calc_burning_rate

 !

 !   purpose     : calculate burning rate

 !

 !   Authors          :  J. Weber, K.C. Tang, and M. Q. Brewster

 !

 !   Creation Date    :  Sep. 04, 2002

 !

 !   Modifications    :

 !

 !    No.     Date         Authors       Description

 !

 !

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

 !

 !   arguments   :

 !

 !

 !      G_ZN     : Global variables for Rocburn_1D_ZN

 !      P        : pressure (atm)

 !      qr       : radiant flux (cal/cm^2/s)

 !      To       : initial temperature (K)

 !      rhoc     : condensed phase density (g/cm^3)

 !      qr_old   : radiant flux at previous time step (cal/cm^2/s)

 !      fr_old   : fr at previous time step

 !      Toa      : apparent initial temperature (K)

 !      rb       : burning rate (cm/s)

 !      Ts       : surface temperature (K)

 !      fr       : fraction of radiation absorbed below surface reaction

 !                 zone

 !      Tn       : temperature profile at previous time step

 !      Tnp1     : temperature profile

 !

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

 !

   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


     IMPLICIT NONE


     include 'mpif.h'


     TYPE(g_burn_1d), POINTER :: g_zn

     REAL(DBL), INTENT(IN)    :: delt, p, qr, to, rhoc

     REAL(DBL), INTENT(IN)    :: qr_old, fr_old

     REAL(DBL), INTENT(INOUT) :: toa

     REAL(DBL), INTENT(INOUT) :: rb, ts, fr

     REAL(DBL), INTENT(IN)    :: tn(:)

     REAL(DBL), INTENT(OUT)   :: tnp1(:)


 !

 ! ============================================================================

 !

 !     delcare local variables

 !

     INTEGER   :: iter, itrmx, i, mpi_comm_rocburn

     REAL(DBL) :: rblast, tslast

     REAL(DBL) :: fs, rb_old

     REAL(DBL) :: first(1:g_zn%nxmax), second(1:g_zn%nxmax)

     REAL(DBL) :: da(1:g_zn%nxmax), db(1:g_zn%nxmax), dc(1:g_zn%nxmax)

     INTEGER   :: ierror


 !

 ! ============================================================================

 !


     IF(g_zn%Model_combustion == 3) THEN

 !

 !        Model_combustion = 3

 !

 !        rb=a*P**n

 !


          rb=g_zn%a_p*p**g_zn%n_p

          tnp1(1)=1000.0


          RETURN


     END IF


     rb_old = rb


     fs=rb*(ts-toa)/g_zn%alfac-fr_old*qr_old/g_zn%lamc


     DO iter=0,g_zn%itermax


 !

 !      get 1st and 2nd derivatives at last time step (n)

 !      with this iterations surface temperature B.C.


        CALL cmpct1(g_zn%nx,g_zn%delz,tn,first, da, db, dc)

        CALL cmpct2(g_zn%nx,g_zn%delz,tn,second, da, db, dc)


 !

 !      determine T profile at this time step (n+1) based on

 !      properties a last time step (n) (explicit differencing)

 !

        DO i=2,g_zn%nx-1

           tnp1(i)=tn(i)+delt*(                                 &

                  (g_zn%alfac*g_zn%zxx(i)-rb_old*g_zn%zx(i))*first(i) +        &

                   g_zn%alfac*g_zn%zx(i)**2*second(i) +                   &

                   fr_old*qr_old/rhoc/g_zn%C*g_zn%Ka*exp(g_zn%Ka*g_zn%x(i)))

        END DO ! i


        tnp1(g_zn%nx)=g_zn%To


 !        update suface B.C. for this iteration

 !        Tnp1(1)=(48.0*Tnp1(2)-36.0*Tnp1(3)+16.0*Tnp1(4)                 ! 4th

 !    +            -3.0*Tnp1(5)-12.0*G_ZN%delz*fs/G_ZN%zx(1))/25.0        ! 4th

          tnp1(1)=(18.0*tnp1(2)-9.0*tnp1(3)+2.0*tnp1(4)                   &

                  -6.0*g_zn%delz*fs/g_zn%zx(1))/11.0                      ! 3rd

 !        Tnp1(1)=(Tnp1(2)-G_ZN%delz*fs/G_ZN%zx(1))                       ! 1st


          ts=tnp1(1)

 !        Tn(1)=Tnp1(1)


 !

 !        calculate rb, Toa, and fr for given  P, qr, and Ts

 !

          CALL zn(p, qr, to, ts, toa, rb, fr)


 !

 !        calculate fs for this iteration of (n+1)

 !

          fs=rb*(ts-toa)/g_zn%alfac-fr*qr/g_zn%lamc


          IF(iter.ge.1) THEN

             IF(abs(ts-tslast).lt.g_zn%tol_Ts.and.abs(rb-rblast).lt.g_zn%tol_Ts)  THEN

                RETURN

             END IF

          END IF


          tslast=ts

          rblast=rb


     END DO ! iter


     print *,' ROCBURN_ZN: rank=',g_zn%rank

     print *,' ROCBURN_ZN: Error - burning rate: Maximum # of iterations '

     print *,' ROCBURN_ZN: itrmx=', g_zn%itermax, 'reached'

     print *,' ROCBURN_ZN: tol_Ts=',g_zn%tol_Ts

     print *,' ROCBURN_ZN: Tn=',tn

     print *,' ROCBURN_ZN: Tnp1=',tnp1

     stop


     CONTAINS


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

 !                             INTERNAL PROCEDURES

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

 !

 !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

 !                                                                             c

 !   subroutine  : ZN                                                          c

 !                                                                             c

 !   purpose     :  alculate rb,Toa given Ts,P,qr                              c

 !                                                                             c

 ! --------------------------------------------------------------------------- c

 !                                                                             c

 !                            written by                                       c

 !                      K.C. Tang and M.Q. Brewster                            c

 !                                                                             c

 !                       Last modified : 09/14/2000                            c

 ! --------------------------------------------------------------------------- c

 !                                                                             c

 !   input       :                                                             c

 !                                                                             c

 !      P        : pressure (atm)                                              c

 !      qr       : radiative heat flux (al/cm^2-s)                             c

 !      To       : propellant temperatuer at deep into the propellant (K)      c

 !      Ts       : surface temperature (K)                                     c

 !      Toa      : ZN temperature (K) for initial guess                        c

 !      rb       : burning rate (cm/s) for initial guess                       c

 !                                                                             c

 !                                                                             c

 !   output      :                                                             c

 !                                                                             c

 !      Toa      : ZN temperature (K)                                          c

 !      rb       : burning rate (cm/s)                                         c

 !      fr       : transmissivity of surface reaction zone                     c

 !                                                                             c

 !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

 !


 !

       SUBROUTINE zn(P, qr, To, Ts, Toa, rb, fr)


       IMPLICIT NONE


       REAL(DBL) :: p, qr, to

       REAL(DBL) :: ts, rb, toa, fr

       REAL(DBL) :: tf

       REAL(DBL) :: delta_t2, ag


       REAL(DBL), parameter :: tol=1.d-5


 ! ============================================================================

 !

 !     delcare local variables

 !


       REAL(DBL) :: j(2,2),xx(2),e(2)

       REAL(DBL) :: xreac

       REAL(DBL) :: term0, term1, term2, term3, term4, term5

       REAL(DBL) :: xcd, da, xg


       integer iter


       term1=g_zn%Ac*g_zn%R*(ts*rhoc)**2*g_zn%alfac*g_zn%C* &

             exp(-g_zn%Ec/g_zn%R/ts)/g_zn%Ec


 !

 !     initial guess

 !


       xx(1)=rb*rhoc

       xx(2)=toa


       IF(g_zn% Model_combustion == 1) THEN


 !        Model_combustion = 1

 !

 !        WSB homogeneouse propellant combustion model

 !


          ag=g_zn%Bg/(g_zn%R*g_zn%R)*g_zn%MW*g_zn%MW


          DO iter=0,g_zn%itermax


 !

 !           set up residual matrix and Jacobian derivative matrix for

 !           Ibiricu and Williams condensed phase expression


             xreac=2.0*rhoc*g_zn%alfac*g_zn%R*ts/xx(1)/g_zn%Ec

             term0=1.0-g_zn%Ka*xreac


             fr=dexp(-xreac*g_zn%Ka)

             term2=g_zn%C*(ts-xx(2))-g_zn%Qc/2.0-fr*qr/xx(1)

             term3=term1/term2


             e(1)=xx(1)**2-term3

             j(1,1)=2.0*xx(1)+term3/term2*fr*qr/xx(1)**2*term0

             j(1,2)=-term3*g_zn%C/term2


 !

 !           set up residual matrix and Jacobian derivative matrix for

 !           Ward-Son-Brewster gas phase expression


             xcd=g_zn%lamg/g_zn%C/xx(1)

             da=4.0*g_zn%Bg*g_zn%C/g_zn%lamg*(p*g_zn%MW*xcd/g_zn%R)**2

             term4=(1.0+da)**0.5

             xg=2.0*xcd/(term4-1.0)

             term5=g_zn%Qc-g_zn%C*(ts-xx(2))

             tf=(qr/xx(1)+g_zn%Qc+g_zn%Qg)/g_zn%C+xx(2)


             e(2)=(tf-ts)/xg+(qr+xx(1)*term5)/g_zn%lamg

             j(2,1)=-qr/xg/xx(1)**2/g_zn%C-(tf-ts)/xg/xx(1)/term4+term5/g_zn%lamg

             j(2,2)=1/xg+1/xcd


 !

 !           update guesses for m and Toa using Newton-Rhapson

 !           & guard against a negative (erroneus) rb solution

 !

             call newton(xx,j,e)

             if(xx(1).lt.0.0) then

                print *,' ROCBURN_ZN: rank=',g_zn%rank,   &

                        ' rb negative xx(1)=',xx(1),' Ts=',ts

                xx(1)=abs(xx(1))

             end if


             if(abs(e(1)).lt.tol .and. abs(e(2)).lt.tol) then


 !

 !              reset final converged values

 !

                rb=xx(1)/rhoc

                toa=xx(2)

                fr=dexp(-xreac*g_zn%Ka)

                return


             end if

          END DO ! iter


          WRITE(*,*) ' ROCBURN_ZN: in ZN rank=',g_zn%rank,  &

                     ' Error-WSB: Maximum # of iterations reached'

          stop


       end if


       if(g_zn%Model_combustion.eq.2) then


 !        Model_combustion = 2

 !

 !        Zeldovich-Novozhilov phenomenological model for composite

 !        propellant

 !


          DO iter=0,g_zn%itermax


 !

 !           set up residual matrix and Jacobian derivative matrix for

 !           Ibiricu and Williams condensed phase expression


             xreac=2.0*rhoc*g_zn%alfac*g_zn%R*ts/xx(1)/g_zn%Ec

             term0=1.0-g_zn%Ka*xreac


             fr=dexp(-xreac*g_zn%Ka)

             term2=g_zn%C*(ts-xx(2))-g_zn%Qc/2.0-fr*qr/xx(1)

             term3=term1/term2


             e(1)=xx(1)**2-term3

             j(1,1)=2.0*xx(1)+term3/term2*fr*qr/xx(1)**2*term0

             j(1,2)=-term3*g_zn%C/term2


 !

 !           set up residual matrix and Jacobian derivative matrix for

 !           Ward-Son-Brewster gas phase expression

 !


             delta_t2 = ts - 300.0

             e(2)=g_zn%a_T*((ts-xx(2)-qr/(xx(1)*g_zn%C))/delta_t2)**g_zn%n_T -    &

                  rhoc*g_zn%a_p*p**g_zn%n_p/xx(1)

             j(2,1)=rhoc*g_zn%a_p*p**g_zn%n_p/(xx(1)*xx(1))

             j(2,2)=-g_zn%n_T*g_zn%a_T*(ts-xx(2)-qr/(xx(1)*g_zn%C))**(g_zn%n_T-1.0)/   &

                     (delta_t2**g_zn%n_T)


 !

 !           update guesses for m and Toa using Newton-Rhapson

 !           & guard against a negative (erroneus) rb solution

 !


             call newton(xx,j,e)


             if(xx(1).lt.0.0) then

                print *,' ROCBURN_ZN: inside ZN rank=',g_zn%rank,  &

                        ' rb negative xx(1)=',xx(1),' Ts=',ts

                xx(1)=abs(xx(1))

             end if


             if(abs(e(1)).lt.tol .and. abs(e(2)).lt.tol) then


 !

 !              reset final converged values

 !

                rb=xx(1)/rhoc

                toa=xx(2)

                fr=dexp(-xreac*g_zn%Ka)

                RETURN


             END IF

          END DO ! iter


          WRITE(*,*) ' ROCBURN_ZN: in ZN rank=',g_zn%rank,    &

                     ' Error-Composite: Maximum # of iterations reached'

          stop


       end if


       IF(g_zn%Model_combustion == 3) THEN


 !

 !

 !        Model_combustion = 3

 !

 !        rb=a*P**n

 !


          rb=g_zn%a_p*p**g_zn%n_p

          tnp1(1)=1000.0


          RETURN


       END IF


       WRITE(*,*) ' ROCBURN_ZN: in ZN rank=',g_zn%rank,  &

                  ' Error-ZN: No appropriate combustion model'

       WRITE(*,*) ' ROCBURN_ZN: Combustion Model= ', g_zn%Model_combustion

       stop


       END SUBROUTINE zn


 !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

 !                                                                             c

 !   subroutine  : Newton                                                      c

 !                                                                             c

 !                                                                             c

 !   Author:                                                                   c

 !                                                                             c

 !      Paul Loner                                                             c

 !                                                                             c

 !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

 !                                                                             c

       SUBROUTINE newton(x,J,e)


       implicit none


       REAL(DBL) :: x(2),j(2,2),e(2),detj,delx(2),dum

       integer i


       detj=j(1,1)*j(2,2)-j(2,1)*j(1,2)


       IF(detj.eq.0) then

         print *,' ROCBURN_ZN: in Newton rank=',g_zn%rank,   &

                 ' Error: singular matrix encountered'

         print *,' ROCBURN_ZN: P=', p

         print *,' ROCBURN_ZN: J=', j

         stop

       ENDIF

 !

 !     invert the 2x2 Jacobian matrix

 !

       dum=j(1,1)

       j(1,1)=j(2,2)

       j(2,2)=dum

       j(1,1)=j(1,1)/detj

       j(1,2)=-j(1,2)/detj

       j(2,1)=-j(2,1)/detj

       j(2,2)=j(2,2)/detj

 !

 !     multiply J^-1()*e() and update guesses

 !

       do i=1,2

         delx(i)=j(i,1)*e(1)+j(i,2)*e(2)

         x(i)=x(i)-delx(i)

       end do ! i


       RETURN

       END SUBROUTINE newton


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

 !

       SUBROUTINE cmpct1(imax,h,u,f, a, b, cc)

 !

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

 !  Modifications:

 !

 !    No.     Date         Programmer    Description

 !    001  May  03, 2001   K. C. Tang    change memory allocation method for

 !                                       temporary arrarys inside a

 !                                       subroutine


       IMPLICIT NONE


       integer i,imax

       REAL(DBL) :: h,u(imax),f(imax)


       REAL(DBL) :: a(:),b(:),cc(:)


 !     set up equations for first derivatives

 !      print *,' in cmpct1 imax=',imax,' h=',h,' u=',u,' f=',f


 !     ALLOCATE(a(imax))

 !     ALLOCATE(b(imax))

 !     ALLOCATE(cc(imax))


       a(1)=0.0

       b(1)=1.0

       cc(1)=0.0

       f(1)=(-25.0*u(1)+48.0*u(2)-36.0*u(3)+16.0*u(4)-3.0*u(5))/(12.0*h) !4th


       do i=2,imax-1

         a(i)=1.0

         b(i)=4.0

         cc(i)=1.0

         f(i)=3.0*(u(i+1)-u(i-1))/h

       end do ! i


       a(imax)=0.0

       b(imax)=1.0

       cc(imax)=0.0

       f(imax)=0.0


 !     solve for first derivatives


       call tridg(imax,a,b,cc,f)


 !     DEALLOCATE(a)

 !     DEALLOCATE(b)

 !     DEALLOCATE(cc)

       RETURN

       END SUBROUTINE cmpct1


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

 !

       SUBROUTINE cmpct2(imax,h,u,s, a, b, cc )

 !

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

 !  Modifications:

 !

 !    No.     Date         Programmer    Description

 !    001  May  03, 2001   K. C. Tang    change memory allocation method for

 !                                       temporary arrarys inside a

 !                                       subroutine


       implicit none


       integer i,imax

       REAL(DBL) :: h,h2,u(imax),s(imax)

       REAL(DBL) :: a(:),b(:),cc(:)

 !     REAL(DBL), ALLOCATABLE :: a(:),b(:),cc(:)


       h2=h*h


 !     ALLOCATE(a(imax))

 !     ALLOCATE(b(imax))

 !     ALLOCATE(cc(imax))


 !     set up equations for second derivatives


       a(1)=0.0

       b(1)=1.0

       cc(1)=0.0

       s(1)=(45.0*u(1)-154.0*u(2)+214.0*u(3)-156.0*u(4)+61.0*u(5)-     &

            10.0*u(6))/12.0/h2                                        ! 4th


       DO i=2,imax-1

         a(i)=1.0

         b(i)=10.0

         cc(i)=1.0

         s(i)=12.0*(u(i+1)-2.0*u(i)+u(i-1))/h2

       END DO ! i


       a(imax)=0.0

       b(imax)=1.0

       cc(imax)=0.0

       s(imax)=0.0


 !     solve for second derivatives


       call tridg(imax,a,b,cc,s)


 !     DEALLOCATE(a)

 !     DEALLOCATE(b)

 !     DEALLOCATE(cc)


       RETURN

       END SUBROUTINE cmpct2


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

 !

       SUBROUTINE tridg (imax,a,b,c,d)

 !

 !     solves a tridiagonal matrix equation of the form:

 !

 !       a(i)*u(i-1) + b(i)*u(i) + c(i)*u(i+1) = d(i)

 !

 !     using the Thomas algorithm.

 !

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


       IMPLICIT NONE


       INTEGER :: imax

       REAL(DBL)  :: a(:), b(:), c(:), d(:)

 !     REAL(DBL), DIMENSION(imax) :: a, b, c, d


       INTEGER ::  i


       DO i=2,imax

         a(i)=a(i)/b(i-1)

         b(i)=b(i)-a(i)*c(i-1)

         d(i)=d(i)-a(i)*d(i-1)

       END DO ! i


       d(imax)=d(imax)/b(imax)


       DO     i=imax-1,1,-1

         d(i)=(d(i)-c(i)*d(i+1))/b(i)

       END DO ! i

       RETURN

       END SUBROUTINE tridg


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

 !                       END OF INTERNAL PROCEDURES

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


   END SUBROUTINE zn_calc_burning_rate


d
const NT & d
Definition: rational_rotation.h:73

tridg
subroutine tridg(imax, a, b, c, d)
Definition: ZN_calc_burning_rate.f90:584

s
double s
Definition: blastest.C:80

Color::b
unsigned char b() const
Definition: Color.h:70

Mesquite::StopWatch::stop
void stop()
Definition: Misc/MsqTimer.cpp:115

Line_2::c
RT c() const
Definition: Line_2.h:150

m_rocburn_zn_global_data
Definition: 02.burn_ZN_global_data.f90:23

i
blockLoc i
Definition: read.cpp:79

x
void int int REAL * x
Definition: read.cpp:74

zn
subroutine zn(P, qr, To, Ts, Toa, rb, fr)
Definition: ZN_calc_burning_rate.f90:216

m_rocburn_zn_global_data::g_burn_1d
Definition: 02.burn_ZN_global_data.f90:48

cmpct1
subroutine cmpct1(imax, h, u, f, a, b, cc)
Definition: ZN_calc_burning_rate.f90:471

newton
subroutine newton(x, J, e)
Definition: ZN_calc_burning_rate.f90:432

Aff_transformation_repS2::print
virtual std::ostream & print(std::ostream &os) const
Definition: Aff_transformationS2.h:277

j
j indices j
Definition: Indexing.h:6

cmpct2
subroutine cmpct2(imax, h, u, s, a, b, cc)
Definition: ZN_calc_burning_rate.f90:526

zn_calc_burning_rate
subroutine zn_calc_burning_rate(G_ZN, delt, P, qr, To, rhoc, qr_old, fr_old, Toa, rb, Ts, fr, Tn, Tnp1)
Definition: ZN_calc_burning_rate.f90:60

Line_2::a
RT a() const
Definition: Line_2.h:140