xref: /phasta/phSolver/incompressible/getdiff.f (revision 595995161822a203c8467e0e4a253d7bd7d6df32)

      subroutine getDiff( dwl,yl, shape, xmudmi, xl, rmu,  rho)
c-----------------------------------------------------------------------
c  compute and add the contribution of the turbulent
c  eddy viscosity to the molecular viscosity.
c-----------------------------------------------------------------------
      use     turbSA
      include "common.h"

      real*8  yl(npro,nshl,ndof), rmu(npro), xmudmi(npro,ngauss),
     &        shape(npro,nshl),   rho(npro),
     &        dwl(npro,nshl),     sclr(npro),
     &        xl(npro,nenl,nsd)
      integer n, e

      real*8  epsilon_ls, kay, epsilon
     &        h_param, prop_blend(npro),test_it(npro)
c
c
c.... get the material properties (2 fluid models will need to determine
c     the "interpolated in phase space" properties....constant for now.
c     two options exist in the interpolation 1) smooth (recommended)
c     interpolation of nodal data, 2) discontinuous "sampling" phase at
c     quadrature points.
c
CAD
CAD    prop_blend is a smoothing function to avoid possible large density
CAD   gradients, e.g., water and air flows where density ratios can approach
CAD   1000.
CAD
CAD    epsilon_ls is an adjustment to control the width of the band over which
CAD    the properties are blended.



      if (iLSet .eq. 0)then

         rho  = datmat(1,1,1)	! single fluid model, i.e., only 1 density
         rmu = datmat(1,2,1)

      else     !  two fluid properties used in this model

!        Smooth the tranistion of properties for a "distance" of epsilon_ls
!        around the interface.  Here "distance" is define as the value of the
!        levelset function.  If the levelset function is properly defined,
!        this is the true distance normal from the front.  Of course, the
!        distance is in a driection normal to the front.

         Sclr = zero
         isc=abs(iRANS)+6
         do n = 1, nshl
            Sclr = Sclr + shape(:,n) * yl(:,n,isc)
         enddo
         do i= 1, npro
            if (sclr(i) .lt. - epsilon_ls)then
               prop_blend(i) = zero
            elseif  (abs(sclr(i)) .le. epsilon_ls)then
               prop_blend(i) = 0.5*(one + Sclr(i)/epsilon_ls +
     &              (sin(pi*Sclr(i)/epsilon_ls))/pi )
            elseif (sclr(i) .gt. epsilon_ls) then
               prop_blend(i) = one
            endif
         enddo
c
        rho = datmat(1,1,2) + (datmat(1,1,1)-datmat(1,1,2))*prop_blend
        rmu = datmat(1,2,2) + (datmat(1,2,1)-datmat(1,2,2))*prop_blend

      endif

CAD	At this point we have a rho that is bounded by the two values for
CAD 	density 1, datmat(1,1,1), the fluid,  and density 2, datmat(1,1,2)
CAD     the gas

c
c  The above approach evaluates all intermediate quantities at the
c  quadrature point, then combines them to form the needed quantities there.
c  1 alternative is calculating all quanties (only rho here) at the nodes and
c  then interpolating the result to the quadrature points.  If this is done,
c  do not forget to do the same thing for rou in e3b!!!
c  ^^^^^^^^^^
c  ||||||||||
c  WARNING
c
c.... dynamic model
c
      if (iLES .gt. 0 .and. iRANS.eq.0) then   ! simple LES
         rmu = rmu + xmudmi(:,intp)
      else if (iRANS.lt.0) then
         if (iRANS .eq. -1) then ! RANS (Spalart-Allmaras)
            call AddEddyViscSA(yl, shape, rmu)
         else if(iRANS.eq.-2) then ! RANS-KE
            sigmaInv=1.0        ! use full eddy viscosity for flow equations
            call AddEddyViscKE(yl, dwl, shape, rho, sigmaInv, rmu)
         endif
         if (iLES.gt.0) then    ! this is DES so we have to blend in
                                ! xmudmi based on max edge length of
                                ! element
            call EviscDESIC (xl,rmu,xmudmi)
         endif
      endif                     ! check for LES or RANS
c
      return
      end

      subroutine EviscDESIC(xl,xmut,xmudmi)

      include "common.h"
      real*8 xmut(npro),xl(npro,nenl,nsd),xmudmi(npro,ngauss)


      do i=1,npro
         dx=maxval(xl(i,:,1))-minval(xl(i,:,1))
         dy=maxval(xl(i,:,2))-minval(xl(i,:,2))
         dz=maxval(xl(i,:,3))-minval(xl(i,:,3))
         emax=max(dx,max(dy,dz))
         if(emax.lt.eles) then  ! pure les
            xmut(i)=xmudmi(i,intp)
         else if(emax.lt.two*eles) then ! blend
            xi=(emax-eles)/(eles)
            xmut(i)=xi*xmut(i)+(one-xi)*(xmudmi(1,intp)+datmat(1,2,2))
         endif                  ! leave at RANS value as edge is twice pure les
      enddo
 !this was made messy by the addEddyVisc routines  Must clean up later.



      return
      end

      subroutine getdiffsclr(shape, dwl, yl, diffus)

      use turbSA
      use turbKE ! access to KE model constants
      include "common.h"

      real*8   diffus(npro), rho(npro)
      real*8   yl(npro,nshl,ndof), dwl(npro,nshl), shape(npro,nshl)
      integer n, e
      rho(:)  = datmat(1,1,1)	! single fluid model, i.e., only 1 density
      if(isclr.eq.0) then  ! solving the temperature equation
         diffus(:) = datmat(1,4,1)
      else if(iRANS.eq.-1) then ! solving SA model
         diffus(:) = datmat(1,2,1)
         call AddSAVar(yl, shape, diffus)
      else if(iRANS.eq.-2)then ! solving KE model
         diffus(:) = datmat(1,2,1)
         if(isclr.eq.2) then
            sigmaInv=1.0/ke_sigma ! different eddy viscosity for epsilon
         else
            sigmaInv=1.0 ! full eddy viscosity for solving kay equation
         endif
         call AddEddyViscKE(yl, dwl, shape, rho, sigmaInv, diffus)
      else                      ! solving scalar advection diffusion equations
         diffus = scdiff(isclr)
      endif
c
      return
      end

      function ev2sa(xmut,rm,cv1)
      implicit none
      real*8 err,ev2sa,rm,cv1,f,dfds,rat,efac
      real*8 pt5,kappa,B,xmut,chi3,denom,cv1_3
      integer iter
      pt5=0.5
      err=1.0d-6
      ev2sa=rm*cv1*1.2599       ! inflection point chi=cv1*cuberoot(2)
      kappa=0.4
c$$$        B=5.5
      efac=0.1108               ! exp(-kappa*B)
      do iter=1,50
         chi3=ev2sa/rm
         chi3=chi3*chi3*chi3
         cv1_3=cv1**3
         denom=chi3+cv1_3

         f=ev2sa*chi3/denom - xmut
         dfds=chi3*(chi3+4.0*cv1_3)/(denom**2)
         rat=-f/dfds
         ev2sa=ev2sa+rat
         if(abs(rat).le.err) goto 20
      enddo
      write(*,*)'ev2sa failed to converge'
      write(*,*) 'dfds,        rat,        ev2sa,        mu'
      write(*,*) dfds,rat,ev2sa,rm
 20   continue
      return
      end
c


      subroutine AddEddyViscSA(yl,shape,rmu)
      use turbSA
      include "common.h"
c     INPUTS
      double precision, intent(in), dimension(npro,nshl,ndof) ::
     &     yl
      double precision, intent(in), dimension(npro,nshl) ::
     &     shape
c     INPUT-OUTPUTS
      double precision, intent(inout), dimension(npro) ::
     &     rmu
c     LOCALS
      logical, dimension(nshl) ::
     &     wallnode
      integer e, n
      double precision xki, xki3, fv1, evisc
c
c     Loop over elements in this block
      do e = 1, npro
c        assume no wall nodes on this element
         wallnode(:) = .false.
         if(itwmod.eq.-2) then  ! effective viscosity
c           mark the wall nodes for this element, if there are any
            do n = 1, nshl
               u1=yl(e,n,2)
               u2=yl(e,n,3)
               u3=yl(e,n,4)
               if((u1.eq.zero).and.(u2.eq.zero).and.(u3.eq.zero))
     &              then
                  wallnode(n)=.true.
               endif
            enddo
         endif
c
         if( any(wallnode(:)) ) then
c if there are wall nodes for this elt, then we are using effective-
c viscosity near-wall modeling, and eddy viscosity has been stored
c at the wall nodes in place of the spalart-allmaras variable; the
c eddy viscosity for the whole element is taken to be the avg of the
c wall values
            evisc = zero
            nwnode=0
            do n = 1, nshl
               if(wallnode(n)) then
                  evisc = evisc + yl(e,n,6)
                  nwnode = nwnode + 1
               endif
            enddo
            evisc = evisc/nwnode
            rmu(e) = rmu(e) + abs(evisc)
c this is what we would use instead of the above if we were allowing
c the eddy viscosity to vary through the element based on non-wall nodes
c$$$               evisc = zero
c$$$               Turb = zero
c$$$               do n = 1, nshl
c$$$                  if(wallmask(n).eq.1) then
c$$$                     evisc = evisc + shape(e,n) * yl(e,n,6)
c$$$                  else
c$$$                     Turb = Turb + shape(e,n) * yl(e,n,6)
c$$$                  endif
c$$$               enddo
c$$$               xki    = abs(Turb)/rmu(e)
c$$$               xki3   = xki * xki * xki
c$$$               fv1    = xki3 / (xki3 + saCv1P3)
c$$$               rmu(e) = rmu(e) + fv1*abs(Turb)
c$$$               rmu(e) = rmu(e) + abs(evisc)
         else
c else one of the following is the case:
c   using effective-viscosity, but no wall nodes on this elt
c   using slip-velocity
c   using no model; walls are resolved
c in all of these cases, eddy viscosity is calculated normally
            Turb = zero
            do n = 1, nshl
               Turb = Turb + shape(e,n) * yl(e,n,6)
            enddo
            xki    = abs(Turb)/rmu(e)
            xki3   = xki * xki * xki
            fv1    = xki3 / (xki3 + saCv1P3)
            rmu(e) = rmu(e) + fv1*abs(Turb)
         endif
      enddo                     ! end loop over elts
      return
      end subroutine AddEddyViscSA



      subroutine AddSAVar(yl,shape,rmu)
      use turbSA
      include "common.h"
c     INPUTS
      double precision, intent(in), dimension(npro,nshl,ndof) ::
     &     yl
      double precision, intent(in), dimension(npro,nshl) ::
     &     shape
c     INPUT-OUTPUTS
      double precision, intent(inout), dimension(npro) ::
     &     rmu
c     LOCALS
      logical, dimension(nshl) ::
     &     wallnode
      integer e, n
      double precision savar, savarw
c     Loop over elements in this block
      do e = 1, npro
c        assume no wall nodes on this element
         wallnode(:) = .false.
         if(itwmod.eq.-2) then  ! effective viscosity
c           mark the wall nodes for this element, if there are any
            do n = 1, nshl
               u1=yl(e,n,2)
               u2=yl(e,n,3)
               u3=yl(e,n,4)
               if((u1.eq.zero).and.(u2.eq.zero).and.(u3.eq.zero))
     &              then
                  wallnode(n)=.true.
               endif
            enddo
         endif
c
         savar=zero
         do n = 1, nshl
            if( wallnode(n) ) then
c if wallmask was set, we're using effective-viscosity wall-model and
c this node is on a wall.  Eddy viscosity has been stored at the wall
c nodes in place of the S-A variable, so we must convert it
               savarw = ev2sa(yl(e,n,6),datmat(1,2,1),saCv1)
               savar  = savar + shape(e,n) * savarw
            else
c if wallmask wasn't set, then one of the following is the case:
c   using effective-viscosity, but this isn't a wall node
c   using slip-velocity
c   using no wall model; wall is resolved
c in all these cases, the S-A variable is calculated normally
               savar  = savar + shape(e,n) * yl(e,n,6)
            endif
         enddo
         rmu(e)=datmat(1,2,1)
         rmu(e) = (rmu(e) + abs(savar)) * saSigmaInv
      enddo                     ! end loop over elts
      return
      end subroutine AddSAVar



      subroutine AddEddyViscKE(yl, dwl, shape, rho, sigmaInv, rmu)
      use turbKE ! access to KE model constants
      include "common.h"
c     INPUTS
      double precision, intent(in), dimension(npro,nshl,ndof) ::
     &     yl
      double precision, intent(in), dimension(npro,nshl) ::
     &     shape, dwl
      double precision, intent(in), dimension(npro) ::
     &     rho
      double precision sigmaInv
c     INPUT-OUTPUTS
      double precision, intent(inout), dimension(npro) ::
     &     rmu
c     LOCALS
      double precision eviscKE, kay, epsilon, dw, CmuKE
      double precision epsInv, Rey, Ret, RetInv, tmp1, fmuKE
      integer e,n
c
      do e = 1, npro
         kay = 0.0
         epsilon = 0.0
         dw = 0.0
         do n = 1, nshl
            kay = kay + shape(e,n)*yl(e,n,6)
            epsilon = epsilon + shape(e,n)*yl(e,n,7)
            dw = dw + shape(e,n)*dwl(e,n)
         enddo
         kay = abs(kay)
         if(kay.lt.1.0e-32) kay=0.0
         epsInv	    = 0
         if ( abs(epsilon) .gt.1.e-32) then
            epsInv        = 1. / abs(epsilon)
         endif

         Rey                 = sqrt(kay) *    dw * rho(e) / rmu(e)
         Ret                 = kay*kay   * epsInv * rho(e) / rmu(e)
         RetInv              = 0
         if(Ret.lt.1.d100.AND.Ret.gt.zero) RetInv=1./Ret
         tmp1     = exp(-0.0165*Rey)
         fmuKE    = (1. -tmp1) ** 2 * (1.+20.5*RetInv) ! fmu of Lam-Bremhorst

         eviscKE=rho(e)*ke_C_mu*fmuKE*kay*kay*epsInv

         rmu(e) = rmu(e) + eviscKE*sigmaInv
      enddo
      return
      end subroutine AddEddyViscKE