moehring.f File Reference

Go to the source code of this file.

Functions/Subroutines
subroutine	moehring (node1, node2, nodem, nelem, lakon, kon, ipkon, nactdog, identity, ielprop, prop, iflag, v, xflow, f, nodef, idirf, df, cp, R, dvi, numf, set, mi, ttime, time, iaxial)
real *8 function	f_k (x, phi, lambda1, zk0, Pup, Tup, rurd, xflow, kup)
real *8 function	f_p (x, phi, lambda1, zk0, Pup, Tup, rurd, xflow, kup)
real *8 function	f_t (x, phi, lambda1, zk0, Pup, Tup, rurd, xflow, kup)
real *8 function	f_m (x, phi, lambda1, zk0, Pup, Tup, rurd, xflow, kup)
real *8 function	f_cm (x, phi, lambda1, zk0, Pup, Tup, rurd, xflow, kup)

Function/Subroutine Documentation

f_cm()

real*8 function f_cm	(	real*8	x,
		real*8	phi,
		real*8	lambda1,
		real*8	zk0,
		real*8	Pup,
		real*8	Tup,
		real*8	rurd,
		real*8	xflow,
		real*8	kup
	)

!     
      implicit none
      integer neq,idid,ipar,j,iwork(100),lrw,liw
      real*8 f_cm,x,rpar(8),rtol,atol,y(1),info(15),
     &     rurd,zk0,lambda1,kup,xflow,pup,tup,phi,t,rwork(160)
!     
      external dkdx
!  
!     storing the parameters
      rpar(1)=phi
      rpar(2)=lambda1
      rpar(3)=zk0
      rpar(4)=pup
      rpar(5)=tup
      rpar(6)=rurd
      rpar(7)=xflow
      rpar(8)=kup
!   
!     relative error
      rtol=1.d-7
!     absolute error
      atol=1.d-7
!     initial value
!
      if(xflow.lt.0d0) then
         t=rurd
      else
         t=1.d0
      endif
!
      neq=1
!     
!     initialisation info field
      do j=1,15
         info(j)=0
      enddo
!     initial condition f(0)
!     core swirl ratio at Rup
!     
      y(1)=kup     
!     
      lrw=160
      liw=60
!     
!     solving the differential equation Möhring 3.35
!     dK/dX=f(K(X))
!     
      call ddeabm(dkdx,neq,t,y,x,info,rtol,atol,idid,
     &     rwork,lrw,iwork,liw,rpar,ipar)
!     
      f_cm=dabs(1-y(1))/(1-y(1))*dabs(1-y(1))**(1.75d0)
     &     *x**(3.6d0)
!     
      return

f_k()

real*8 function f_k	(	real*8	x,
		real*8	phi,
		real*8	lambda1,
		real*8	zk0,
		real*8	Pup,
		real*8	Tup,
		real*8	rurd,
		real*8	xflow,
		real*8	kup
	)

!
      implicit none
      integer neq,idid,ipar,iwork(100),lrw,liw,j
      real*8 f_k,x,rpar(8),rtol,atol,y(1),info(15),
     &     rurd,zk0,lambda1,kup,xflow,pup,tup,phi,t,rwork(160)
!     
      external dkdx
!
!     storing the parameters
      rpar(1)=phi
      rpar(3)=zk0
!
!    relative error
         rtol=1.d-7
!     absolute error
         atol=1.d-7
!
!     initial value
         if(xflow.lt.0d0) then
            t=rurd
         else 
            t=1.d0
         endif

         neq=1
!
!     initialisation info field
         do j=1,15
            info(j)=0
         enddo
!     initial condition f(0)
!     core swirl ratio at Rup repectively Rdown depending 
!     on the type of element centrifugal or centripetal
!     
         y(1)= kup 
!
         lrw=160
         liw=60
!
!     solving the differential equation Möhring 3.35
!     dK/dX=f(K(X))
!
         if(dabs(xflow).gt.1e-6) then
            call ddeabm(dkdx,neq,t,y,x,info,rtol,atol,idid,
     &        rwork,lrw,iwork,liw,rpar,ipar)
         else
            y(1)=1/(zk0+1)
         endif
!
       f_k=y(1)**2*x
!
       return

f_m()

real*8 function f_m	(	real*8	x,
		real*8	phi,
		real*8	lambda1,
		real*8	zk0,
		real*8	Pup,
		real*8	Tup,
		real*8	rurd,
		real*8	xflow,
		real*8	kup
	)

!     
      implicit none
      integer neq,idid,ipar,j,iwork(100),lrw,liw
      real*8 f_m,x,rpar(8),rtol,atol,y(1),info(15),
     &     rurd,zk0,lambda1,kup,xflow,pup,tup,phi,t,rwork(160)
!     
      external dkdm
!
!     storing the parameters
      rpar(1)=phi
      rpar(2)=lambda1
      rpar(3)=zk0
      rpar(4)=pup
      rpar(5)=tup
      rpar(6)=rurd
      rpar(7)=xflow
      rpar(8)=kup
!     
!     relative error
      rtol=1.d-7
!     absolute error
      atol=1.d-7
! 
!     initial value
      if(xflow.lt.0d0) then
         t=rurd
      else 
         t=1.d0
      endif    
      neq=1
!     
!     initialisation info field
      do j=1,15
         info(j)=0
      enddo
!     initial condition f(0)
!     core swirl ratio at Rup
!     
      y(1)=0     
!     
      lrw=160
      liw=60
!     
!     solving the differential equation Möhring 3.35
!     dK/dX=f(K(X))
!     
      call ddeabm(dkdm,neq,t,y,x,info,rtol,atol,idid,
     &     rwork,lrw,iwork,liw,rpar,ipar)
!     
      f_m=2*y(1)*x
!     
      return

f_p()

real*8 function f_p	(	real*8	x,
		real*8	phi,
		real*8	lambda1,
		real*8	zk0,
		real*8	Pup,
		real*8	Tup,
		real*8	rurd,
		real*8	xflow,
		real*8	kup
	)

!
      implicit none
      integer neq,idid,ipar,iwork(100),lrw,liw,j
      real*8 f_p,x,rpar(8),rtol,atol,y(1),info(15),rurd,
     &     zk0,lambda1,kup,xflow,pup,tup,phi,t,rwork(160)
!     
      external dkdp
!     storing the parameters
      rpar(1)=phi
      rpar(2)=lambda1
      rpar(3)=zk0
      rpar(4)=pup
      rpar(5)=tup
      rpar(6)=rurd
      rpar(7)=xflow
      rpar(8)=kup
!

!    relative error
      rtol=1.d-7
!     absolute error
      atol=1.d-7
!
!     initial value
      if(xflow.lt.0d0) then
         t=rurd
      else
         t=1.d0
      endif
      neq=1
!     
!     initialisation info field
      do j=1,15
         info(j)=0
      enddo
!     initial condition f(0)
!     core swirl ratio at Rup
!     
      y(1)= 0d0    
!      
      lrw=160
      liw=60
!          
      call ddeabm(dkdp,neq,t,y,x,info,rtol,atol,idid,
     &     rwork,lrw,iwork,liw,rpar,ipar)
!     
      f_p=2*y(1)*x
!     
      return

f_t()

real*8 function f_t	(	real*8	x,
		real*8	phi,
		real*8	lambda1,
		real*8	zk0,
		real*8	Pup,
		real*8	Tup,
		real*8	rurd,
		real*8	xflow,
		real*8	kup
	)

!     
      implicit none
      integer neq,idid,ipar,iwork(100),lrw,liw,j
      real*8 f_t,x,rpar(8),rtol,atol,y(1),info(15),rurd,
     &     zk0,lambda1,kup,xflow, pup,tup,phi,t,rwork(160)
!     
      external dkdt
!
!     storing the parameters
      rpar(1)=phi
      rpar(2)=lambda1
      rpar(3)=zk0
      rpar(4)=pup
      rpar(5)=tup
      rpar(6)=rurd
      rpar(7)=xflow
      rpar(8)=kup
!             
!     relative error
      rtol=1.d-7
!     absolute error
      atol=1.d-7
!
!     initial value
      if(xflow.lt.0d0) then
         t=rurd
      else 
         t=1.d0
      endif
!      if(xflow.lt.0d0) then
!         t=rurd
!      else 
!         t=1.d0
!      endif     
      neq=1
!     
!     initialisation info field
      do j=1,15
         info(j)=0
      enddo
!     initial condition f(0)
!     core swirl ratio at Rup
!     
      y(1)= 0d0    
!     
       lrw=160
       liw=60
!     
      call ddeabm(dkdt,neq,t,y,x,info,rtol,atol,idid,
     &     rwork,lrw,iwork,liw,rpar,ipar)
!     
      f_t=2d0*y(1)*x
!     
      return

moehring()

subroutine moehring	(	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	dvi,
		integer	numf,
		character*81, dimension(*)	set,
		integer, dimension(*)	mi,
		real*8	ttime,
		real*8	time,
		integer	iaxial
	)

!     
!     moehring element
!     This subroutines computes the evolution of the core swirl ratio 
!     for a disc stator system with either centrifugal or centripetal
!     flow.
!     Theoretical explanations can be found in
!     "Untersuchung dfes radialen Druckverlaufes und des übertragenen
!     drehmomentes im Radseitenraum von Kreiselpumpen bei glatter,
!     ebene Radwand und bei Anvendung von Rückenschaufeln"
!     Uwe Klaus Möhring , Dissertation, 
!     An der Üniversität Carolo-Wilhelmina zu Braunschweig 1976
!
!     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,iaxial,
     &     ipkon(*),kon(*),kgas,key,neval,ier,limit,lenw,last,
     &     iwork2(400),node_up,node_down,mi(*)
!     
      real*8 prop(*),v(0:mi(2),*),xflow,f,df(*),r,dvi,pi,
     &     r_min, r_max,cr, r_shroud,rsrmax,gap,swirl_up,
     &     pup,pdown,tup,tdown,kappa,cp,ttime,time,
     &     rup,rdown,k0,kup,cq,re_phi,phi,lambda1, lambda2,
     &     a,b,epsabs,
     &     epsrel,result,abserr,work(1200),zk0,t1,t2,p1,p2,pr,
     &     qred_crit,omega,rurd,c_p,cm,mr,f_k,f_t,f_p,f_m,f_cm,
     &     pdiff,pdiff_min,xflow_0,cq_0,phi_0
!     
      external dkdx,dkdp,dkdt,dkdm,f_k,f_t,f_p,f_m,f_cm
!     
!      numf=4
!     
      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
!     
         pi=4.d0*datan(1.d0)
         kappa=(cp/(cp-r))
         index=ielprop(nelem)
         qred_crit=dsqrt(kappa/r)*
     &        (1+0.5d0*(kappa-1))**(-0.5*(kappa+1)/(kappa-1))
!     
!     Because there is no explicit expression relating massflow
!     to pressure loss for möhrings
!     initial mass flow is set to arbitrarily
!     with consideration to flow direction
!     
         node1=kon(ipkon(nelem)+1)
         node2=kon(ipkon(nelem)+3)
         p1=v(2,node1)
         p2=v(2,node2)
         t1=v(0,node1)
         t2=v(0,node2)
!     
!     fictious cross section
!         A=1E-5
         if(p1.gt.p2) then
            xflow=1/dsqrt(t1)*p1*qred_crit*0.5d0
         else
            xflow=-1/dsqrt(t1)*p1*qred_crit*0.5d0
         endif
!     
      elseif(iflag.eq.2)then
!     
         numf=4
         index=ielprop(nelem) 
         pi=4.d0*datan(1.d0)
!     
!     Cr
         cr=0.315
!     
!     minimal disc radius
         r_min=prop(index+1)
!     
!     maximal disc radius
         r_max=prop(index+2)
!     
!     R_min/R_max
         rurd=r_min/r_max
!     
!     disc/stator gap
         gap=prop(index+3)
!     
!     shroud radius
         r_shroud=prop(index+4)
!     
!     R_schroud/R_max
         rsrmax=r_shroud/r_max      
!     
!     defining flow parameters
!     
!     massflow
         xflow=v(1,nodem)*iaxial
!     
!     upstream node
         node_up=nint(prop(index+5))
!     
!     downstream node
         node_down=nint(prop(index+6))
!   
!     centripetal
         if(lakon(nelem)(2:5).eq.'MRGP') then
            if(xflow.lt.0d0) then
               xflow=-xflow
            endif
!     
            rup=r_max
            rdown=r_min
!     
            nodef(1)=node_up
            nodef(2)=node_up
            nodef(3)=nodem
            nodef(4)=node_down
!     
!     centrifugal
         elseif(lakon(nelem)(2:5).eq.'MRGF') then
            if(xflow.gt.0d0) then
               xflow=-xflow
            endif
!     
            rup=r_min
            rdown=r_max
! 
!            if(xflow.gt.0) then
            nodef(1)=node_up
            nodef(2)=node_up
            nodef(3)=nodem
            nodef(4)=node_down
         endif
!     
!     upstream pressure
         pup=v(2,node_up)
!
!     downstream pressure
         pdown=v(2,node_down)
!
!     Upstream temperature
         tup=v(0,node_up)
!
!     downstream temperature
         tdown=v(0,node_down)
!     
         idirf(1)=2
         idirf(2)=0
         idirf(3)=1
         idirf(4)=2
!     
!     rotation related parameters
!     
!     rotation
!
         omega=prop(index+7)
!     
!     swirl at R_upstream
         swirl_up=prop(index+8)
!     
!     core swirl ratio when xflow=0
         k0=1/(1+(rsrmax**3.6*
     &        (rsrmax+4.6*gap/r_max))**(4.d0/7.d0))
!     
!     core swirl ratio at R_inlet
         kup=swirl_up/(omega*rup)
!     
!     dynamic_viscosity
!         kgas=0
!         call dynamic_viscosity(kgas,Tup,dvi)
         if(dabs(dvi).lt.1e-30) then
            write(*,*) '*ERROR in moehring: '
            write(*,*) '       no dynamic viscosity defined'
            write(*,*) '       dvi= ',dvi
            call exit(201)
c            kgas=0
c            call dynamic_viscosity(kgas,Tup,dvi)
         endif 
!     
!     defining common coefficients
!     
         cq=xflow*r*tup/(pup*omega*(r_max)**3)
!     
         re_phi=(omega*r_max**2*pup)/(dvi*r*tup)
!     
         phi=cq*(re_phi)**0.2d0
!
         zk0=(1-k0)/k0
!     
!     lambda1
         lambda1=(r_max-r_min)/dabs(r_max-r_min)*pi*cr/4*
     &        (dvi*r/(omega*r_max**2)**0.2d0*(omega*r_max**3)/r)
!     
!     lambda2
         lambda2=2d0*r/(omega**2*r_max**2)
!     
!*************************************************************************
!     integration of K(X),dKdp(X),dKdT(X),dKdm(X)

         limit=201
         lenw=5*limit
!     
!         if(lakon(nelem)(2:5).eq.'MRGF') then
!xflow.lt.0d0) then
!     
!     lower integration boundary
            a=rurd
!     
!     upper integration boundary
            b=1d0
!     
!         elseif(lakon(nelem)(2:5).eq.'MRGP') then
!     
!     lower integration boundary
!            a=1d0
!     
!     upper integration boundary
!            b=rurd
!         endif         
!     
!     absolute error
         epsabs=1e-7
!     
!     relative error
         epsrel=1e-7
!     
!     choice for local integration rule
         key=1
!
!     determining minimal pressure difference for xflow<<1
!
         if(lakon(nelem)(2:5).eq.'MRGF')then
            xflow_0=-0.00000003e-3
         elseif(lakon(nelem)(2:5).eq.'MRGP') then
            xflow_0=0.00000003e-3
         endif
!
         cq_0=xflow_0*r*tup/(pup*omega*(r_max)**3)
!     
         phi_0=cq_0*(re_phi)**0.2d0
!
         call dqag(f_k,a,b,epsabs,epsrel,key,result,abserr,neval,ier,
     &        limit,lenw,last,iwork2,work,phi_0,lambda1,zk0,pup,tup,
     &        rurd,xflow_0,kup)
!
!     pressure coefficient
         c_p=2*result
         pdiff_min=c_p*pup/(4*r*tup)*omega**2*r_max**2
         pdiff=dabs(pdown-pup)
!
!     K(x)
!     
         call dqag(f_k,a,b,epsabs,epsrel,key,result,abserr,neval,ier,
     &        limit,lenw,last,iwork2,work,phi,lambda1,zk0,pup,tup,rurd,
     &        xflow,kup)
!        
!     residual
!
         if(lakon(nelem)(2:5).eq.'MRGF') then
            f=lambda2*(pdown-pup)/(pdown+pup)*tup-result
         elseif(lakon(nelem)(2:5).eq.'MRGP') then
            f=lambda2*(pup-pdown)/(pup+pdown)*tup-result
         endif
!     
!     pressure coefficient
         c_p=2*result
!     
!     moment coefficient
         call dqag(f_cm,a,b,epsabs,epsrel,key,result,abserr,neval,ier,
     &        limit,lenw,last,iwork2,work,phi,lambda1,zk0,pup,tup,rurd,
     &        xflow,kup)
!         
         cm=0.5d0*pi*cr*re_phi**(-0.2d0)*result
         mr=0.5d0*cm*(pup/1000d0/(r*tup))*omega**2*r_max**5
         pr=mr*omega
!     
!     pressure
!     
         call dqag(f_p,a,b,epsabs,epsrel,key,result,abserr,neval,ier,
     &        limit,lenw,last,iwork2,work,phi,lambda1,zk0,pup,tup,rurd,
     &        xflow,kup)
!     
!     partial derivative (upstream pressure)
! 
         if(lakon(nelem)(2:5).eq.'MRGF') then
            df(1)=-2*lambda2**pdown/(pdown+pup)**2*tup-result
         elseif(lakon(nelem)(2:5).eq.'MRGP') then
            df(1)=2*lambda2**pdown/(pdown+pup)**2*tup-result
         endif
!     
!     temperature
!     
         call dqag(f_t,a,b,epsabs,epsrel,key,result,abserr,neval,ier,
     &        limit,lenw,last,iwork2,work,phi,lambda1,zk0,pup,tup,rurd,
     &        xflow,kup)
!     
!     partial derivative (upstream temperature)       
         if(lakon(nelem)(2:5).eq.'MRGF') then
            df(2)=lambda2**(pdown-pup)/(pdown+pup)-result
         elseif(lakon(nelem)(2:5).eq.'MRGP') then  
            df(2)=lambda2**(pup-pdown)/(pdown+pup)-result
         endif
!     
!     mass flow
!     
         call dqag(f_m,a,b,epsabs,epsrel,key,result,abserr,neval,ier,
     &        limit,lenw,last,iwork2,work,phi,lambda1,zk0,pup,tup,rurd,
     &        xflow,kup)
!     
!     partial derivative (mass flow)             
         df(3)=-result
!     
!     pressure 
!     partial derivative (downstream pressure)
         if(lakon(nelem)(2:5).eq.'MRGF') then
            df(4)=2*lambda2**pup/(pdown+pup)**2*tup
         elseif(lakon(nelem)(2:5).eq.'MRGP') then
            df(4)=-2*lambda2**pup/(pdown+pup)**2*tup
         endif
!   
      endif
!     
      xflow=xflow/iaxial
      df(3)=df(3)*iaxial
!
      return