liquidpump.f File Reference

Go to the source code of this file.

Functions/Subroutines
subroutine	liquidpump (node1, node2, nodem, nelem, nactdog, identity, ielprop, prop, iflag, v, xflow, f, nodef, idirf, df, rho, g, co, numf, mi, ttime, time, iaxial)

Function/Subroutine Documentation

liquidpump()

subroutine liquidpump	(	integer	node1,
		integer	node2,
		integer	nodem,
		integer	nelem,
		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	rho,
		real*8, dimension(3)	g,
		real*8, dimension(3,*)	co,
		integer	numf,
		integer, dimension(*)	mi,
		real*8	ttime,
		real*8	time,
		integer	iaxial
	)

!
!     pump for incompressible media
!     
      implicit none
!     
      logical identity
!      
      integer nelem,nactdog(0:3,*),node1,node2,nodem,
     &     ielprop(*),nodef(*),idirf(*),index,iflag,
     &     inv,id,numf,npu,i,mi(*),iaxial
!      
      real*8 prop(*),v(0:mi(2),*),xflow,f,df(*),ttime,time,
     &     p1,p2,rho,g(3),dg,z1,z2,co(3,*),
     &     xpu(10),ypu(10),xxpu(10),yypu(10),dh
!
      numf=3
!
      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).or.(iflag.eq.2))then
!     
         index=ielprop(nelem)
!
         npu=nint(prop(index+1))
         do i=1,npu
            xpu(i)=prop(index+2*i)
            ypu(i)=prop(index+2*i+1)
         enddo
!     
         p1=v(2,node1)
         p2=v(2,node2)
!     
         z1=-g(1)*co(1,node1)-g(2)*co(2,node1)-g(3)*co(3,node1)
         z2=-g(1)*co(1,node2)-g(2)*co(2,node2)-g(3)*co(3,node2)
!     
         if(iflag.eq.2) then
            xflow=v(1,nodem)*iaxial
            if(xflow.ge.0.d0) then
               inv=1
            else
               inv=-1
            endif
            nodef(1)=node1
            nodef(2)=nodem
            nodef(3)=node2
            idirf(1)=2
            idirf(2)=1
            idirf(3)=2
         endif
!     
         dg=dsqrt(g(1)*g(1)+g(2)*g(2)+g(3)*g(3))
!     
         if(iflag.eq.1) then
            dh=(z2-z1+(p2-p1)/rho)/dg
!
!           reverting the order in xpu and ypu and storing the
!           result in xxpu and yypu
!
            do i=1,npu
               xxpu(i)=xpu(npu+1-i)
               yypu(i)=ypu(npu+1-i)
            enddo
            call ident(yypu,dh,npu,id)
            if(id.eq.0) then
               xflow=xxpu(1)
            elseif(id.eq.npu) then
               xflow=0.d0
            else
               xflow=xxpu(id)+(xxpu(id+1)-xxpu(id))*(dh-yypu(id))/
     &               (yypu(id+1)-yypu(id))
            endif
         else
            df(1)=1.d0/rho
            df(3)=-df(1)
            xflow=xflow/rho
            call ident(xpu,xflow,npu,id)
            if(id.eq.0) then
               if(xflow.ge.0.d0) then
                  f=z1-z2+(p1-p2)/rho+dg*ypu(1)
                  df(2)=0.d0
               else
                  df(2)=-1.d10
                  f=z1-z2+(p1-p2)/rho+dg*(ypu(1)+xflow*df(2))
                  df(2)=df(2)*dg/rho
               endif
            elseif(id.eq.npu) then
               df(2)=-1.d10
               f=z1-z2+(p1-p2)/rho+dg*(ypu(npu)+df(2)*(xflow-xpu(npu)))
               df(2)=df(2)*dg/rho
            else
               df(2)=(ypu(id+1)-ypu(id))/(xpu(id+1)-xpu(id))
               f=z1-z2+(p1-p2)/rho+dg*(ypu(id)+(xflow-xpu(id))*df(2))
               df(2)=df(2)*dg/rho
            endif
         endif
!     
      endif
!     
      xflow=xflow/iaxial
      df(2)=df(2)*iaxial
!     
      return