calc_residual_cross_split.f File Reference

Go to the source code of this file.

Functions/Subroutines
real *8 function	calc_residual_cross_split (pt1, Tt1, xflow1, xflow2, pt2, Tt2, ichan_num, A1, A2, A_s, dh1, dh2, alpha, zeta_fac, kappa, R, ider, iflag)

Function/Subroutine Documentation

calc_residual_cross_split()

real*8 function calc_residual_cross_split	(	real*8	pt1,
		real*8	Tt1,
		real*8	xflow1,
		real*8	xflow2,
		real*8	pt2,
		real*8	Tt2,
		integer	ichan_num,
		real*8	A1,
		real*8	A2,
		real*8	A_s,
		real*8	dh1,
		real*8	dh2,
		real*8	alpha,
		real*8	zeta_fac,
		real*8	kappa,
		real*8	R,
		integer	ider,
		integer	iflag
	)

!
      implicit none
!
      integer icase,ichan_num,ider,icrit1,icrit2,iflag,ier
!
      real*8 
!     In- and Output
     &f,r,
!
!     Kappa stuff
     &kappa,km1,kp1,
!
     &pt1,pt2,tt1,tt2,xflow1,xflow2,
!
     &pt2_lim,
!
     &zeta,
!
     &a1,a2,
!
     &ts0,ts1,ts2,dh1,dh2,alpha,q_crit,pspt_crit,q0,q1,q2,pspt0,
     &pspt1,pspt2,w1,w2,w1w2,w2w1,pi,z2d390,z1p090,z60,z90,hq,m1,m2,
     &zeta_fac,xflow_s,q_s,ts_s,pspt_s,w_s,wsw1,a_s,asa1,vsv1
!
      real*8 table_zeta(2,10)
!
      pi=4.d0*datan(1.d0)
!
      icrit1=0
      icrit2=0
!
!     setting icase (always adiabatic)
!     
      icase=0;
!
      km1=kappa-1.d0
      kp1=kappa+1.d0
      q_crit=dsqrt(kappa/r)*
     &   (1+0.5d0*(kappa-1))**(-0.5d0*(kappa+1)/(kappa-1))
      pspt_crit=(2./(kappa+1.))**(kappa/(kappa-1.))
!
      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 Cross Split:'
         write(*,*)'Critical conditions at 1'
      endif
      q2=xflow2*dsqrt(tt1)/pt2/a2
      if(q2.ge.q_crit) then
         q2=q_crit
         icrit2=1
         write(*,*)'*WARNING in Cross Split:'
         write(*,*)'Critical conditions at 2'
      endif
!
!     Flow velocity at inlet
      call ts_calc(xflow1,tt1,pt1,kappa,r,a1,ts0,icase)
      pspt0=(ts0/tt1)**(kappa/(kappa-1))
      call wpi(w1, pspt0, q0, 
     &      dsqrt(tt1),kappa,r)
!
!     Flow velocity at outlet
      call ts_calc(xflow2,tt1,pt1,kappa,r,a2,ts1,icase)
      pspt1=(ts1/tt1)**(kappa/(kappa-1))
      call wpi(w2, pspt1, q1, 
     &      dsqrt(tt2),kappa,r) 
!
      w2w1=w2/w1
      w1w2=w1/w2
!
!     Main branch
      if(ichan_num.eq.1) then
!
!          Zeta as in Calculix
         zeta=0.4d0*(1-w2w1)**2
!
         zeta=zeta*(w1w2)**2
!
!     First branch
      elseif((ichan_num.eq.2).or.(ichan_num.eq.3)) then            
         hq=dh2/dh1
         if(alpha.le.60.or.hq.le.2.d0/3.d0) then
            zeta=0.95d0*((w2w1-2d0*dcos(alpha*pi/180))
     &                 *w2w1+1.d0)
            zeta=zeta*(w1w2)**2
         else
            z2d390=0.95d0*((w2w1-2d0*dcos(90.d0*pi/180))
     &                 *w2w1+1.d0)
            z1p090=0.95*(0.34d0+w2w1**2)
            z90=z2d390+(3*hq-2.d0)*(z1p090-z2d390)
            z60=0.95d0*((w2w1-2d0*dcos(60.d0*pi/180))
     &                 *w2w1+1.d0)
            zeta=z60+(alpha/30.d0-2.d0)*(z90-z60)
            zeta=zeta*(w1w2)**2
         endif
!
      endif
!
!     zeta_fac for side branches are all =1
!     main branch can be set by the user in ACC Designer
      zeta=zeta*zeta_fac
!
      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_cross_split=f
!
      return