wye.f File Reference

Go to the source code of this file.

Functions/Subroutines
subroutine	wye (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

wye()

subroutine wye	(	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 wye split element(zeta calculation according to Idel'chik)
!     Written by Yavor Dobrev
!     For an explanation of the parameters see tee.f
!
!     author: Yannick Muller
!
      implicit none
!
      logical identity
!
      character*81 set(*)
      character*8 lakon(*)
!
      integer nelem,nactdog(0:3,*),node1,node2,nodem,nodem1,numf,kon(*),
     &ipkon(*),ielprop(*),nodef(*),idirf(*),index,iflag,inv,mi(2),
     &nelem1,ichan_num,ider,icase,i,iaxial
!
      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,a_s,calc_residual_wye,dh1,dh2,alpha,xflow1,xflow2,
     &pi,zeta_fac,ts0,pspt0,pspt2,m1,m2,ts2,ttime,time
!
      index=ielprop(nelem)
!
      if (iflag.eq.0.d0) 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+4)
         elseif(nelem.eq.nint(prop(index+3)))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
         xflow=xflow/50
!     
      elseif (iflag.eq.2)then
!     
         numf=6
!
         kappa=(cp/(cp-r))
         pi=4.d0*datan(1.d0)
!
!        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+4)
!
!        Outlet conditions
!
         tt2=v(0,node2)
         xflow2=v(1,nodem)*iaxial
         pt2=v(2,node2)
!
         if(nelem.eq.nint(prop(index+2))) then
!
            a2 = a1
            a_s = prop(index+6)
!
            dh1 = prop(index+9)
            if(dh1.eq.0.d0) then
               dh1 = dsqrt(4*a1/pi)
            endif
!
            dh2 = dh1
            alpha = 0
            ichan_num = 1
            zeta_fac = prop(index+13)
!
         elseif(nelem.eq.nint(prop(index+3))) then
!
            ichan_num = 2
            a2 = prop(index+6)
!
            dh1 = prop(index+9)
            if(dh1.eq.0.d0) then
               dh1 = dsqrt(4*a1/pi)
            endif
!
            dh2 = prop(index+10)
            if(dh2.eq.0.d0) then
               dh2 = dsqrt(4*a2/pi)
            endif
!
            alpha = prop(index+8)
            zeta_fac = prop(index+14)
!
         endif
!
!        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
!
!        Sets the types of the degrees of freedom
         idirf(1)=2
         idirf(2)=0
         idirf(3)=1
         idirf(4)=1
         idirf(5)=2
         idirf(6)=0
!
         if(ider.eq.0.d0) then
!           Residual
            f=calc_residual_wye(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_wye(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;
!
         pi=4.d0*datan(1.d0)
!
!        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+4)
!
!        Outlet conditions
         tt2=v(0,node2)
         xflow2=v(1,nodem)*iaxial
         pt2=v(2,node2)
!
         if(nelem.eq.nint(prop(index+2))) then
!
            a2 = a1
            a_s = prop(index+6)
!
            dh1 = prop(index+9)
            if(dh1.eq.0.d0) then
               dh1 = dsqrt(4*a1/pi)
            endif
!
            dh2 = dh1
            alpha = 0
            ichan_num = 1
            zeta_fac = prop(index+13)
!
         elseif(nelem.eq.nint(prop(index+3))) then
!
            ichan_num = 2
            a2 = prop(index+6)
!
            dh1 = prop(index+9)
            if(dh1.eq.0.d0) then
               dh1 = dsqrt(4*a1/pi)
            endif
!
            dh2 = prop(index+10)
            if(dh2.eq.0.d0) then
               dh2 = dsqrt(4*a2/pi)
            endif
!
            alpha = prop(index+8)
            zeta_fac = prop(index+14)
!
         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_wye(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