absolute_relative.f File Reference

Go to the source code of this file.

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

Function/Subroutine Documentation

absolute_relative()

subroutine absolute_relative	(	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(*)	mi,
		real*8	ttime,
		real*8	time,
		integer	iaxial
	)

!     
!     orifice element
!
!     author: Yannick Muller
!     
      implicit none
!     
      logical identity
      character*8 lakon(*)
      character*81 set(*)
!     
      integer nelem,nactdog(0:3,*),node1,node2,nodem,numf,
     &     ielprop(*),nodef(*),idirf(*),index,iflag,
     &     ipkon(*),kon(*),nelemswirl,mi(*),iaxial
!     
      real*8 prop(*),v(0:mi(2),*),xflow,f,df(*),kappa,r,
     &     p1,p2,t1,t2,cp,physcon(*),km1,kp1,kdkm1,
     &     kdkp1,u,pi,qred_crit,pt1,pt2,tt1,tt2,ct,fact,
     &     cp_cor,ttime,time
!     
      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
         xflow=0.d0
!     
      elseif(iflag.eq.2) then
!     
         numf=4
         kappa=(cp/(cp-r))
         km1=kappa-1.d0
         kp1=kappa+1.d0
         kdkm1=kappa/km1
         kdkp1=kappa/kp1
!     
         index=ielprop(nelem)
!         
         u=prop(index+1)
         ct=prop(index+2)
         nelemswirl=nint(prop(index+3))
!
        if(nelemswirl.ne.0) then
!
!     previous element is a preswirl nozzle
!
            if(lakon(nelemswirl)(2:5).eq.'ORPN') then
               ct=prop(ielprop(nelemswirl)+5)
!
!     previous element is a forced vortex
!
            elseif(lakon(nelemswirl)(2:5).eq.'VOFO') then
               ct=prop(ielprop(nelemswirl)+7)
!
!     previous element is a free vortex
!
            elseif(lakon(nelemswirl)(2:5).eq.'VOFR') then
               ct=prop(ielprop(nelemswirl)+9)
            endif
         endif                
!     
         pt1=v(2,node1)
         pt2=v(2,node2)
!
         if(lakon(nelem)(2:4).eq.'ATR') then
!     
            xflow=v(1,nodem)*iaxial
            tt1=v(0,node1)-physcon(1)
            tt2=v(0,node2)-physcon(1)
!     
            nodef(1)=node1
            nodef(2)=node1
            nodef(3)=nodem
            nodef(4)=node2
!     
!     in the case of a negative flow direction
!     
            if(xflow.le.0d0) then
               write(*,*)''
               write(*,*)'*WARNING:'
               write(*,*)'in element',nelem
               write(*,*)'TYPE=ABSOLUTE TO RELATIVE'
               write(*,*)'mass flow negative!'
               write(*,*)'check results and element definition'
            endif
!       
         elseif(lakon(nelem)(2:4).eq.'RTA') then
!     
            xflow=v(1,nodem)*iaxial
            tt1=v(0,node1)-physcon(1)
            tt2=v(0,node2)-physcon(1)
!     
            nodef(1)=node1
            nodef(2)=node1
            nodef(3)=nodem
            nodef(4)=node2
!     
            if(xflow.le.0d0) then
               write(*,*)''
               write(*,*)'*WARNING:'
               write(*,*)'in element',nelem
               write(*,*)'TYPE=RELATIVE TO ABSOLUTE'
               write(*,*)'mass flow negative!'
               write(*,*)'check results and element definition'
            endif
         endif
!     
         idirf(1)=2
         idirf(2)=0
         idirf(3)=1
         idirf(4)=2
!     
!     computing temperature corrected Cp=Cp(T) coefficient 
         call cp_corrected(cp,tt1,tt2,cp_cor)
!     
            if(tt1.lt.273) then
               tt1= tt2
            endif
!
            if(cp_cor.eq.0.d0) then
               cp_cor=cp
            endif
!     
!     transformation from absolute system to relative system
!     
         if(lakon(nelem)(2:4).eq.'ATR') then
!     
            fact=1+(u**2-2*u*ct)/(2*cp_cor*tt1)
!     
            f=pt2-pt1*(fact)**kdkm1
!     
!     pressure node 1
!
            df(1)=-fact**kdkm1
!
!     temperature node1
!
            df(2)=-pt1*kdkm1*(-(u**2-2*u*ct)/(2*cp_cor*tt1**2))
     &           *fact**(kdkm1-1)
!
!     mass flow node m
!
            df(3)=0
!
!     pressure node 2
!     
            df(4)=1
!     
!     transformation from relative system to absolute system
!     
         elseif(lakon(nelem)(2:4).eq.'RTA') then
!     
            fact=1-(u**2-2*u*ct)/(2*cp*tt1)
!     
            f=pt2-pt1*(fact)**kdkm1
!     
            df(1)=-fact**kdkm1
!     
            df(2)=-pt1*kdkm1*((u**2-2*u*ct)/(2*cp*tt1**2))
     &           *fact**(kdkm1-1)
!     
            df(3)=0
!     
            df(4)=1
!     
         endif

      elseif(iflag.eq.3) then
            
         kappa=(cp/(cp-r))
         km1=kappa-1.d0
         kp1=kappa+1.d0
         kdkm1=kappa/km1
         kdkp1=kappa/kp1
!     
         index=ielprop(nelem)
!         
         u=prop(index+1)
         ct=prop(index+2)
         nelemswirl=nint(prop(index+3))
!
        if(nelemswirl.ne.0) then
!
!     previous element is a preswirl nozzle
!
            if(lakon(nelemswirl)(2:5).eq.'ORPN') then
               ct=prop(ielprop(nelemswirl)+5)
!
!     previous element is a forced vortex
!
            elseif(lakon(nelemswirl)(2:5).eq.'VOFO') then
               ct=prop(ielprop(nelemswirl)+7)
!
!     previous element is a free vortex
!
            elseif(lakon(nelemswirl)(2:5).eq.'VOFR') then
               ct=prop(ielprop(nelemswirl)+9)
            endif
         endif                
!     
         pt1=v(2,node1)
         pt2=v(2,node2)
!
         if(lakon(nelem)(2:4).eq.'ATR') then
!     
            xflow=v(1,nodem)*iaxial
            tt1=v(0,node1)-physcon(1)
            tt2=v(0,node2)-physcon(1)
!     
            nodef(1)=node1
            nodef(2)=node1
            nodef(3)=nodem
            nodef(4)=node2
!     
!     in the case of a negative flow direction
!     
            if(xflow.le.0d0) then
               write(*,*)''
               write(*,*)'*WARNING:'
               write(*,*)'in element',nelem
               write(*,*)'TYPE=ABSOLUTE TO RELATIVE'
               write(*,*)'mass flow negative!'
               write(*,*)'check results and element definition'
            endif
!       
         elseif(lakon(nelem)(2:4).eq.'RTA') then
!     
            xflow=v(1,nodem)*iaxial
            tt1=v(0,node1)-physcon(1)
            tt2=v(0,node2)-physcon(1)
!     
            nodef(1)=node1
            nodef(2)=node1
            nodef(3)=nodem
            nodef(4)=node2
!     
            if(xflow.le.0d0) then
               write(*,*)''
               write(*,*)'*WARNING:'
               write(*,*)'in element',nelem
               write(*,*)'TYPE=RELATIVE TO ABSOLUTE'
               write(*,*)'mass flow negative!'
               write(*,*)'check results and element definition'
            endif
         endif
!     
         idirf(1)=2
         idirf(2)=0
         idirf(3)=1
         idirf(4)=2
!     
!     computing temperature corrected Cp=Cp(T) coefficient 
         call cp_corrected(cp,tt1,tt2,cp_cor)
!     
         if(tt1.lt.273) then
            tt1= tt2
         endif
!     
         if(cp_cor.eq.0.d0) then
            cp_cor=cp
         endif

                  write(1,*) ''
         write(1,55) ' from node',node1,
     &' to node', node2,':   air massflow rate=',xflow,''
 55      FORMAT(1x,a,i6,a,i6,a,e11.4,a,a,e11.4,a)

            write(1,56)'       Inlet node ',node1,':     Tt1= ',tt1,
     &           ', Ts1= ',tt1,', Pt1= ',pt1
            write(1,*)'             Element ',nelem,lakon(nelem)
            write(1,57)'             u= ',u,' ,Ct= ',ct,''
            write(1,56)'      Outlet node ',node2,':     Tt2= ',tt2,
     &           ', Ts2= ',tt2,', Pt2= ',pt2
!     
 56      FORMAT(1x,a,i6,a,e11.4,a,e11.4,a,e11.4,a,e11.4)  
 57      FORMAT(1x,a,e11.4,a,e11.4,a)

      endif
!     
      xflow=xflow/iaxial
      df(3)=df(3)*iaxial
!     
      return