pt2zpt1_crit.f File Reference

Go to the source code of this file.

Functions/Subroutines
subroutine	pt2zpt1_crit (pt2, pt1, Tt1, lambda, kappa, r, l, d, pt2zpt1_c, Qred1_crit, crit, icase, M1)

Function/Subroutine Documentation

pt2zpt1_crit()

subroutine pt2zpt1_crit	(	real*8, intent(in)	pt2,
		real*8, intent(in)	pt1,
		real*8, intent(in)	Tt1,
		real*8, intent(in)	lambda,
		real*8, intent(in)	kappa,
		real*8, intent(in)	r,
		real*8, intent(in)	l,
		real*8, intent(in)	d,
		real*8, intent(inout)	pt2zpt1_c,
		real*8, intent(inout)	Qred1_crit,
		logical, intent(inout)	crit,
		integer, intent(in)	icase,
		real*8	M1
	)

!     
      implicit none
!
      logical crit
!
      integer icase,i
!     
      real*8 pt2,pt1,lambda,kappa,l,d,m1,pt2zpt1,pt2zpt1_c,
     &     km1,kp1,kp1zk,tt1,r,xflow_crit,qred1_crit,fmin,
     &     f,fmax,m1_min,m1_max,lld
!
      intent(in) pt2,pt1,tt1,lambda,kappa,r,l,d,
     &     icase
!
      intent(inout) qred1_crit,crit,pt2zpt1_c
!
      crit=.false.
!
!     useful variables and constants
! 
      km1=kappa-1.d0
      kp1=kappa+1.d0
      kp1zk=kp1/kappa
      lld=lambda*l/d
!
!     adiabatic case
!
      if(icase.eq.0) then
!     
!        computing M1 using dichotomy method (dividing the interval with the function
!        root iteratively by 2)
!     
         i=1
!
         m1_min=0.001d0
         m1_max=1
!
         fmin=(1.d0-m1_min**2)*(kappa*m1_min**2)**(-1)
     &        +0.5d0*kp1zk*log((0.5d0*kp1)*m1_min**2
     &        *(1+0.5d0*km1*m1_min**2)**(-1))-lld
!     
         fmax=(1.d0-m1_max**2)*(kappa*m1_max**2)**(-1)
     &        +0.5d0*kp1zk*log((0.5d0*kp1)*m1_max**2
     &        *(1+0.5d0*km1*m1_max**2)**(-1))-lld
         do
            i=i+1
            m1=(m1_min+m1_max)*0.5d0
!     
            f=(1.d0-m1**2)*(kappa*m1**2)**(-1)
     &           +0.5d0*kp1zk*log((0.5d0*kp1)*m1**2
     &           *(1+0.5d0*km1*m1**2)**(-1))-lld
!     
            if(abs(f).le.1e-6) then
               exit
            endif
            if(i.gt.50) then
               exit
            endif
!     
            if(fmin*f.le.0.d0) then
               m1_max=m1
               fmax=f
            else
               m1_min=m1
               fmin=f
            endif
         enddo
!     
         pt2zpt1_c=m1*(0.5d0*kp1)**(0.5*kp1/km1)
     &        *(1+0.5d0*km1*m1**2)**(-0.5d0*kp1/km1)
!     
!     isothermal case
!
      elseif (icase.eq.1) then
!     
!        computing M1 using dichotomy method for choked conditions M2=1/dsqrt(kappa)
!        (1.d0-kappa*M1**2)/(kappa*M1**2)+log(kappa*M1**2)-lambda*l/d=0
!     
         i=1
!     
         m1_min=0.001d0
         m1_max=1/dsqrt(kappa)
!     
         fmin=(1.d0-kappa*m1_min**2)/(kappa*m1_min**2)
     &        +log(kappa*m1_min**2)-lambda*l/d
!     
         fmax=(1.d0-kappa*m1_max**2)/(kappa*m1_max**2)
     &        +log(kappa*m1_max**2)-lambda*l/d
!     
         do
            i=i+1
            m1=(m1_min+m1_max)*0.5d0
!     
            f=(1.d0-kappa*m1**2)/(kappa*m1**2)
     &           +log(kappa*m1**2)-lambda*l/d
!     
            if((abs(f).le.1e-5).or.(i.ge.50)) then
               exit
            endif
!     
            if(fmin*f.le.0.d0) then
               m1_max=m1
               fmax=f
            else
               m1_min=m1
               fmin=f
            endif
         enddo
!     
!        computing the critical pressure ratio in the isothermal case
!        pt=A*dsqrt(kappa)/(xflow*dsqrt(kappa Tt))*
!           M*(1+0.5d0*(kappa-1)M**2)**(-0.5d0*(kappa+1)/(kappa-1))
!        and forming the pressure ratio between inlet and outlet(choked)
!     
c         Pt2zPt1_c=dsqrt(Tt2/Tt1)*M1*dsqrt(kappa)*((1+0.5d0*km1/kappa)
c     &        /(1+0.5d0*km1*M1**2))**(0.5d0*(kappa+1)/km1)
         pt2zpt1_c=m1*dsqrt(kappa)*((1+0.5d0*km1/kappa)
     &        /(1+0.5d0*km1*m1**2))**(0.5d0*(kappa+1)/km1+0.5d0)
!     
      endif
!     
      pt2zpt1=pt2/pt1
      if(pt2zpt1.le.pt2zpt1_c) then
         crit=.true.
      endif
!     
      qred1_crit=m1*dsqrt(kappa/r)
     &     *(1+0.5d0*km1*m1**2)**(-0.5d0*kp1/km1)
!
c      write(*,*) 'pt2zpt1_crit ',dsqrt(kappa/r)
c     &     *(1+0.5d0*km1)**(-0.5d0*kp1/km1)
!     
      return