cross_split.f File Reference

Go to the source code of this file.

Functions/Subroutines
subroutine	cross_split (node1, node2, nodem, nelem, lakon, kon, ipkon, nactdog, identity, ielprop, prop, iflag, v, xflow, f, nodef, idirf, df, cp, r, physcon, numf, set, mi, ider, ttime, time, iaxial)

Function/Subroutine Documentation

cross_split()

subroutine cross_split	(	integer	node1,
		integer	node2,
		integer	nodem,
		integer	nelem,
		character*8, dimension(*)	lakon,
		integer, dimension(*)	kon,
		integer, dimension(*)	ipkon,
		integer, dimension(0:3,*)	nactdog,
		logical	identity,
		integer, dimension(*)	ielprop,
		real*8, dimension(*)	prop,
		integer	iflag,
		real*8, dimension(0:mi(2),*)	v,
		real*8	xflow,
		real*8	f,
		integer, dimension(*)	nodef,
		integer, dimension(*)	idirf,
		real*8, dimension(*)	df,
		real*8	cp,
		real*8	r,
		real*8, dimension(*)	physcon,
		integer	numf,
		character*81, dimension(*)	set,
		integer, dimension(2)	mi,
		integer	ider,
		real*8	ttime,
		real*8	time,
		integer	iaxial
	)

!
!     A cross split element(zeta calculation according to Idel'chik)
!     Written by Yavor Dobrev
!     For an explanation of the parameters see tee.f
!
      implicit none
!
      logical identity
!
      character*8 lakon(*)
      character*81 set(*)
!
      integer
     &nelem,nactdog(0:3,*),node1,node2,nodem,numf,kon(*),ipkon(*),
     &ielprop(*),nodef(*),idirf(*),iflag,inv,mi(2),nodem1,nelem1,
     &ichan_num,ider,icase,iaxial,index
!
      real*8 prop(*),v(0:mi(2),*),xflow,f,df(*),kappa,r,tt1,tt2,pt1,
     &pt2,cp,physcon(*),km1,kp1,kdkm1,pt2pt1,pt1pt2,pt2pt1_crit,
     &tdkp1,a,a1,a2,calc_residual_cross_split,dh1,dh2,alpha,xflow1,
     &xflow2,zeta_fac,a_s,ts0,pspt0,pspt2,m1,m2,ts2,ttime,time
!
      index=ielprop(nelem)
!
      if(iflag.eq.0) then
         identity=.true.
         if(nactdog(2,node1).ne.0)then
            identity=.false.
         elseif(nactdog(2,node2).ne.0)then
            identity=.false.
         elseif(nactdog(1,nodem).ne.0)then
            identity=.false.
         endif
!
      elseif(iflag.eq.1)then
!
         kappa=(cp/(cp-r))
         kp1=kappa+1d0
         km1=kappa-1d0
         kdkm1=kappa/km1
         tdkp1=2.d0/kp1
!
         if(nelem.eq.nint(prop(index+2))) then
            a=prop(index+5)
         elseif(nelem.eq.nint(prop(index+3)))then
            a=prop(index+6)
         elseif(nelem.eq.nint(prop(index+4)))then
            a=prop(index+6)
         endif
!
         pt1=v(2,node1)
         pt2=v(2,node2)
!
         if(pt1.ge.pt2) then
            inv=1
            tt1=v(0,node1)-physcon(1)
            tt2=v(0,node2)-physcon(1)
         else
            inv=-1
            pt1=v(2,node2)
            pt2=v(2,node1)
            tt1=v(0,node2)-physcon(1)
            tt2=v(0,node1)-physcon(1)
         endif
!     
         pt1pt2=pt1/pt2
         pt2pt1=1/pt1pt2
!
         pt2pt1_crit=tdkp1**kdkm1
!
         if(pt2pt1.gt.pt2pt1_crit) then
            xflow=inv*pt1*a*dsqrt(2.d0*kdkm1*pt2pt1**(2.d0/kappa)
     &           *(1.d0-pt2pt1**(1.d0/kdkm1))/r)/dsqrt(tt1)
         else
            xflow=inv*pt1*a*dsqrt(kappa/r)*tdkp1**(kp1/(2.d0*km1))/
     &           dsqrt(tt1)
         endif
!     
      elseif(iflag.eq.2)then
!     
         numf=6
!
         kappa=(cp/(cp-r))
!
!     setting icase (always adiabatic)
!     
         icase=0;
!
!        Inlet conditions are the same for both branches
!
!        Determining the previous element as
!        the incoming mass flow is defined there
         nelem1=nint(prop(index+1))
         nodem1=kon(ipkon(nelem1)+2)
!
!        Inlet conditions
         pt1=v(2,node1)
         tt1=v(0,node1)-physcon(1)
         xflow1=v(1,nodem1)*iaxial
         a1=prop(index+5)
         dh1=prop(index+9)
!
!        Set the node numbers for the degrees of freedom
         nodef(1)=node1
         nodef(2)=node1
         nodef(3)=nodem1
         nodef(4)=nodem
         nodef(5)=node2
         nodef(6)=node2
!
         idirf(1)=2
         idirf(2)=0
         idirf(3)=1
         idirf(4)=1
         idirf(5)=2
         idirf(6)=0
!
         if(nelem.eq.nint(prop(index+2))) then
            ichan_num=1
!
            tt2=v(0,node2)
            xflow2=v(1,nodem)*iaxial
            pt2=v(2,node2)
            a2=a1
            a_s=prop(index+6)
            dh2=dh1
            alpha=0
            zeta_fac=prop(index+11)
!
         elseif(nelem.eq.nint(prop(index+3))) then
            ichan_num=2
!
            tt2=v(0,node2)
            xflow2=v(1,nodem)*iaxial
            pt2=v(2,node2)
            a2=prop(index+6)
            dh2=prop(index+10)
            alpha=prop(index+7)
            zeta_fac=prop(index+12)
!
         elseif(nelem.eq.nint(prop(index+4))) then
            ichan_num=3
!
            tt2=v(0,node2)
            xflow2=v(1,nodem)*iaxial
            pt2=v(2,node2)
            a1=prop(index+5)
            a2=prop(index+6)
            dh2=prop(index+10)
            alpha=prop(index+8)
            zeta_fac=prop(index+12)
!
         endif
!
         if(ider.eq.0)then
!           Residual
            f=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)
         else
!           Derivatives
            call calc_ider_cross_split(df,pt1,tt1,xflow1,xflow2,pt2,
     &    tt2,ichan_num,a1,a2,a_s,dh1,dh2,alpha,zeta_fac,
     &    kappa,r,ider,iflag)
         endif
!
      elseif(iflag.eq.3) then
!
         kappa=(cp/(cp-r))
!
!         setting icase (always adiabatic)
!
         icase=0;
!
!        Inlet conditions are the same for both branches
!
!        Determining the previous element as
!        the incoming mass flow is defined there
!
         nelem1=nint(prop(index+1))
         nodem1=kon(ipkon(nelem1)+2)
!
!        Inlet conditions
!
         pt1=v(2,node1)
         tt1=v(0,node1)-physcon(1)
         xflow1=v(1,nodem1)*iaxial
         a1=prop(index+5)
         dh1=prop(index+9)      
!
         if(nelem.eq.nint(prop(index+2))) then
            ichan_num=1
!
            tt2=v(0,node2)
            xflow2=v(1,nodem)*iaxial
            pt2=v(2,node2)
            a2=a1
            a_s=prop(index+6)
            dh2=dh1
            alpha=0
            zeta_fac=prop(index+11)
!
         elseif(nelem.eq.nint(prop(index+3))) then
            ichan_num=2
!
            tt2=v(0,node2)
            xflow2=v(1,nodem)*iaxial
            pt2=v(2,node2)
            a2=prop(index+6)
            dh2=prop(index+10)
            alpha=prop(index+7)
            zeta_fac=prop(index+12)
!
         elseif(nelem.eq.nint(prop(index+4))) then
            ichan_num=3
!
            tt2=v(0,node2)
            xflow2=v(1,nodem)*iaxial
            pt2=v(2,node2)
            a1=prop(index+5)
            a2=prop(index+6)
            dh2=prop(index+10)
            alpha=prop(index+8)
            zeta_fac=prop(index+12)
!
         endif
!
!        Write the main information about the element
!
!        Flow velocity at inlet
         call ts_calc(xflow1,tt1,pt1,kappa,r,a1,ts0,icase)
         pspt0=(ts0/tt1)**(kappa/(kappa-1))
!        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,*) ''
         write(1,55) ' from node ',node1,
     &        ' to node ', node2,':   air massflow rate= ',xflow
!     
         write(1,56)'       Inlet node ',node1,':    Tt1= ',tt1,
     &        ', Ts1= ',ts0,', Pt1= ',pt1,
     &        ', M1= ',m1
         write(1,*)'             Element ',nelem,lakon(nelem)
     &        ,', Branch ',ichan_num
!     
 55      format(1x,a,i6,a,i6,a,e11.4,a)
 56      format(1x,a,i6,a,e11.4,a,e11.4,a,e11.4,a,e11.4)
!     
!     Set ider to calculate the residual
         ider=0
!     
!     Calculate the element one last time with enabled output
         f=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)
!     
         write(1,56)'      Outlet node ',node2,':   Tt2= ',tt2,
     &        ' , Ts2= ',ts2,' , Pt2= ',pt2,
     &        ', M2= ',m2
!     
      endif
!     
      xflow=xflow/iaxial
      df(3)=df(3)*iaxial
      df(4)=df(4)*iaxial
!     
      return