calc_residual_tee.f File Reference

Go to the source code of this file.

Functions/Subroutines
real *8 function	calc_residual_tee (pt1, Tt1, xflow1, xflow2, pt2, Tt2, A1, A2, zeta_fac, kappa, R, ider, iflag)

Function/Subroutine Documentation

calc_residual_tee()

real*8 function calc_residual_tee	(	real*8	pt1,
		real*8	Tt1,
		real*8	xflow1,
		real*8	xflow2,
		real*8	pt2,
		real*8	Tt2,
		real*8	A1,
		real*8	A2,
		real*8	zeta_fac,
		real*8	kappa,
		real*8	R,
		integer	ider,
		integer	iflag
	)

!
!     Function for calculating the residual of both branches of a tee
!
!     author: Yannick Muller
!
      implicit none
!
      integer 
!     Isothermal/adiabatic
     &icase,
!     Residual/Derivatives
     &ider,
!     Critical conditions at inlet or outlet
     &icrit1,
     &icrit2,
!     Where we are in the program
     &iflag
!
      real*8 
!     Residual
     &f,
!     Kappa stuff
     &kappa,
     &r,
!     State variables
     &pt1,
     &pt2,
     &tt1,
     &tt2,
     &ts0,
     &ts1,
     &ts2,
     &xflow1,
     &xflow2,
!
     &pt2_lim,
!
     &zeta,
     &zeta_fac,
!
!     Areas
     &a1,
     &a2,
!     Reduced mass flows
     &q0,
     &q1,
     &q2,
     &q_crit,
!     Pressure ratios
     &pspt_crit,
     &pspt0,
     &pspt1,
     &pspt2,
!     Flow velocities
     &w1,
     &w2,
     &w1w2,
     &w2w1,
!    Mach numbers,
     &m1,
     &m2
!
      icrit1 = 0
      icrit2 = 0
!     setting icase (always adiabatic)
      icase=0;
!
!     Critical values
      q_crit = dsqrt(kappa/r)*
     &   (1+0.5d0*(kappa-1))**(-0.5d0*(kappa+1)/(kappa-1))
      pspt_crit = (2/(kappa+1)) ** (kappa/(kappa-1))
!
!     These reduced mass flows are equivalent
!     (reduced mass flow at inlet)
      q0 = xflow1*dsqrt(tt1)/pt1/a1
      q1 = xflow2*dsqrt(tt1)/pt1/a2
      if(q1.ge.q_crit) then
         q1 = q_crit
         icrit1 = 1
         write(*,*)'*WARNING in Tee:'
         write(*,*)'Critical conditions at 1'
      endif
!     Reduced mass flow at outlet
      q2 = xflow2*dsqrt(tt1)/pt2/a2
      if(q2.ge.q_crit) then
         q2 = q_crit
         icrit2 = 1
         write(*,*)'*WARNING in Tee:'
         write(*,*)'Critical conditions at 2'
      endif
!
!     Flow velocity at inlet
!     Static temperature
      ts0=tt1
      call ts_calc(xflow1,tt1,pt1,kappa,r,a1,ts0,icase)
!     Pressure ratio
      pspt0 = (ts0/tt1)**(kappa/(kappa-1))
!     Velocity
      call wpi(w1, pspt0, q0, 
     &      dsqrt(tt1),kappa,r) 
!
!     Flow velocity at outlet
!     Static temperature
      call ts_calc(xflow2,tt1,pt1,kappa,r,a2,ts1,icase)
!     Pressure ratio
      pspt1 = (ts1/tt1)**(kappa/(kappa-1))
!     Velocity      
      call wpi(w2, pspt1, q1, 
     &      dsqrt(tt1),kappa,r) 
!
!     Velocity ratio
      w2w1=w2/w1
      w1w2=w1/w2
!
!     Zeta calculation
      zeta=1.d0+0.3d0*w2w1**2
      zeta=zeta*(w1w2)**2
!
      zeta = zeta_fac*zeta
!
!     Residual calculation
      if(icrit2.ne.1) then
         if(icrit1.ne.1) then
            f = pt2 - pt1*pspt1**zeta
         else
            f = xflow2*dsqrt(tt1)/pt1/a2-q_crit
         endif
      else
         f = xflow2*dsqrt(tt1)/pt2/a2-q_crit
      endif
!
      if(iflag.eq.3) then
!
         write(1,57)'             zeta= ',zeta
 57      format(1x,a,f9.4)
!
      else if (iflag.eq.4) then

!        Calculate Mach numbers
         call machpi(m1,pspt0,kappa,r)
         call ts_calc(xflow2,tt2,pt2,kappa,r,a2,ts2,icase)
!        Pressure ratio
         pspt2 = (ts2/tt2)**(kappa/(kappa-1))
         call machpi(m2,pspt2,kappa,r)
      
         write(1,80)'Inlet: Tt1= ',tt1,
     &              ', pt1= ',pt1,', M1= ',m1
     
         write(1,77)'mass flow = ',xflow2,', kappa = ',kappa,
     &              ', zeta= ',zeta

         write(1,80)'Outlet: Tt2= ',tt2,
     &              ', pt2= ',pt2,', M2= ',m2
     
 80   format(3x,a,f10.6,a,f10.2,a,f10.6)
 77   format(3x,a,f10.6,a,f10.2,a,f10.6)
      endif
!
      calc_residual_tee=f
!
      return