source_ghd_granular_energy.f

Go to the documentation of this file.

!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
!                                                                      C
!  Module name: SOURCE_GHD_GRANULAR_ENERGY(sourcelhs,sourcerhs,IJK,IER)
!  Purpose: Calculate the source terms in the granular energy equation C
!           for GHD theory                                             C
!                                                                      C
!  Author: S. Benyahia, J. Galvin, C. Hrenya        Date: 22-JAN-09    C
!  Reviewer:                                          Date:            C
!                                                                      C
!  Literature/Document References:                                     C
!     C. Hrenya handnotes and Garzo, Hrenya, Dufty papers (PRE, 2007)  C
!                                                                      C
!  Variables referenced:                                               C
!                                                                      C
!  Variables modified:                                                 C
!                                                                      C
!     Local variables: sourcelhs, sourcerhs                            C
!                                                                      C
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
      SUBROUTINE source_ghd_granular_energy(SOURCELHS, SOURCERHS, IJK)
!...Translated by Pacific-Sierra Research VAST-90 2.06G5  12:17:31  12/09/98
!...Switches: -xf
!
!     Include param.inc file to specify parameter values
!
!-----------------------------------------------
!   M o d u l e s
!-----------------------------------------------
      USE param
      USE param1
      USE parallel
      USE physprop
      USE run
      USE drag
      USE geometry
      USE fldvar
      USE visc_s
      USE ghdtheory
      USE trace
      USE indices
      USE constant
      USE toleranc
      USE compar        !//d
      USE fun_avg
      USE functions
      IMPLICIT NONE
!-----------------------------------------------
!   G l o b a l   P a r a m e t e r s
!-----------------------------------------------
!-----------------------------------------------
!   D u m m y   A r g u m e n t s
!-----------------------------------------------
!
!                      Indices
      INTEGER          IJK, I, J, K, IM, JM, KM, IMJK, IJMK, IJKM,&
                       IJKE, IJKW, IJKN, IJKS, IJKT, IJKB
!
!                      number densities
      DOUBLE PRECISION NiE, NiW, NiN, NiS, NiT,&
                       NiB, Nip, Ntotal
!
!                      mass of particles & non-zero gran. temp.
      DOUBLE PRECISION Mi, NonZeroTheta
!
!                      Thermal mobility-related source terms
      DOUBLE PRECISION ThermMobilityX, ThermMobilityY, ThermMobilityZ
!
!                      del.Joi and Fi.Joi terms
      DOUBLE PRECISION DelDotJoi, FiDotJoi, JoiXC, JoiYC, JoiZC
      DOUBLE PRECISION FiXC, FiYC, FiZC,ni(smax), SIGMAI(smax)
!
!                      phase index
      INTEGER          M, L
!
!                      Dufour-related source terms
      DOUBLE PRECISION DufourX, DufourY, DufourZ, DijQTerm, &
                       DijQTermE, DijQTermW, DijQTermN, DijQTermS, DijQTermT, DijQTermB,&
                       LijTermW,LijTermE,LijTermN,LijTermS,LijTermT,LijTermB

      DOUBLE PRECISION DijQTermE_H,DijQTermE_A,DijQTermW_H,DijQTermW_A,DijQTermN_H,DijQTermN_A,DijQTermS_H,&
                       DijQTermS_A,DijQTermT_H,DijQTermT_A,DijQTermB_H,DijQTermB_A,LijTermE_H,LijTermE_A, &
                       LijTermW_H,LijTermW_A,LijTermN_H,LijTermN_A,LijTermS_H,LijTermS_A,LijTermT_H, &
                       LijTermT_A,LijTermB_H,LijTermB_A
!
!                      Source terms to be kept on RHS
      DOUBLE PRECISION SOURCERHS, PressureRhs, ShearProduction, BulkViscRhs, DissDivURhs, phi_tot, &
                       SOURCE_FLUID,SINK_FLUID
!    DOUBLE PRECISION chi_ij
!
!                      Source terms to be kept on LHS
      DOUBLE PRECISION SOURCELHS, PressureLhs, CollDissipation, BulkViscLhs, DissDivULhs,VSLIP
      DOUBLE PRECISION UGC, USCM, VGC, VSCM, WGC, WSCM

      i = i_of(ijk)
      j = j_of(ijk)
      k = k_of(ijk)
      im = im1(i)
      jm = jm1(j)
      km = km1(k)
      imjk = im_of(ijk)
      ijmk = jm_of(ijk)
      ijkm = km_of(ijk)
      ijke = east_of(ijk)
      ijkw = west_of(ijk)
      ijkn = north_of(ijk)
      ijks = south_of(ijk)
      ijkt = top_of(ijk)
      ijkb = bottom_of(ijk)

      nonzerotheta = max(theta_m(ijk,mmax), small_number)

      ntotal = zero
      phi_tot = zero
      DO m = 1,smax
          ntotal = ntotal + rop_s(ijk,m)*6.d0/(pi*d_p(ijk,m)**3 * ro_s(ijk,m))
          ni(m) = rop_s(ijk,m)*6.d0/(pi*d_p(ijk,m)**3 * ro_s(ijk,m))
          phi_tot = phi_tot + rop_s(ijk,m)/ro_s(ijk,m)
          sigmai(m) = d_p(ijk,m)
      ENDDO

! Production by shear: (S:grad(v))
!     p_s*tr(D)
      pressurerhs = p_s_c(ijk,mmax)*zmax(( -trd_s_c(ijk,mmax) ))
      pressurelhs = p_s_c(ijk,mmax)*zmax((  trd_s_c(ijk,mmax) ))

!     mu_s*tr(D^2)
      shearproduction = 2d0*mu_s_c(ijk,mmax)*trd_s2(ijk,mmax)

!     lambda_s*tr(D)^2
      bulkviscrhs = zmax(  lambda_s_c(ijk,mmax) ) * trd_s_c(ijk,mmax)**2
      bulkvisclhs = zmax( -lambda_s_c(ijk,mmax) ) * trd_s_c(ijk,mmax)**2


! Energy dissipation by collisions: (3/2)*n*T*zeta
!     (3/2)*n*T*zeta0; zeroth order cooling rate term
      colldissipation = 1.5d0*ntotal*zeta0(ijk)
!     (3/2)*n*T*zetaU*div(U) :
      dissdivurhs = 1.5d0*ntotal*theta_m(ijk,mmax)* zmax( -zetau(ijk)*trd_s_c(ijk,mmax) )
      dissdivulhs = 1.5d0*ntotal* zmax(  zetau(ijk)*trd_s_c(ijk,mmax) )


      dufourx = zero
      dufoury = zero
      dufourz = zero
      thermmobilityx = zero
      thermmobilityy = zero
      thermmobilityz = zero
      deldotjoi = zero
      fidotjoi  = zero
      source_fluid = zero
      sink_fluid = zero
      DO m = 1,smax

! Part of heat flux: div (q)
!     Sum_ij [ div( T^2*DijQ/nj*grad(nj)) ] -> Dufour term

          DO l = 1,smax
              mi = (pi/6.d0)*d_p(ijk,l)**3 * ro_s(ijk,l)

              nip = rop_s(ijk,l) /mi
              nie = rop_s(ijke,l)/mi
              niw = rop_s(ijkw,l)/mi
              nin = rop_s(ijkn,l)/mi
              nis = rop_s(ijks,l)/mi

              dijqterm  = zero
              dijqterme = zero
              dijqtermw = zero
              dijqtermn = zero
              dijqterms = zero

              if(rop_s(ijk,l)/ro_s(ijk,l) > dil_ep_s) dijqterm = theta_m(ijk,mmax)**2*dijq(ijk,m,l) / nip
              if(rop_s(ijke,l)/ro_s(ijke,l) > dil_ep_s) dijqterme =theta_m(ijke,mmax)**2*dijq(ijke,m,l) / nie
              if(rop_s(ijkw,l)/ro_s(ijkw,l) > dil_ep_s) dijqtermw =theta_m(ijkw,mmax)**2*dijq(ijkw,m,l) / niw
              if(rop_s(ijkn,l)/ro_s(ijkn,l) > dil_ep_s) dijqtermn =theta_m(ijkn,mmax)**2*dijq(ijkn,m,l) / nin
              if(rop_s(ijks,l)/ro_s(ijks,l) > dil_ep_s) dijqterms =theta_m(ijks,mmax)**2*dijq(ijks,m,l) / nis

                dijqterme_h = avg_x_s(dijqterm, dijqterme, i)
                dijqterme_a = avg_x(dijqterm, dijqterme, i)
              IF(abs(dijqterme_h) .le. abs(dijqterme_a))THEN
                dijqterme = dijqterme_h
              ELSE
                dijqterme = dijqterme_a
              ENDIF

                dijqtermw_h = avg_x_s(dijqtermw, dijqterm, im)
                dijqtermw_a = avg_x(dijqtermw, dijqterm, im)
              IF(abs(dijqtermw_h) .le. abs(dijqtermw_a))THEN
                dijqtermw = dijqtermw_h
              ELSE
                dijqtermw = dijqtermw_a
              ENDIF

                dijqtermn_h = avg_y_s(dijqterm, dijqtermn, j)
                dijqtermn_a = avg_y(dijqterm, dijqtermn, j)
              IF(abs(dijqtermn_h) .le. abs(dijqtermn_a))THEN
                dijqtermn = dijqtermn_h
              ELSE
                dijqtermn = dijqtermn_a
              ENDIF

                dijqterms_h = avg_y_s(dijqterms, dijqterm, jm)
                dijqterms_a = avg_y(dijqterms, dijqterm, jm)
              IF(abs(dijqterms_h) .le. abs(dijqterms_a))THEN
                dijqterms = dijqterms_h
              ELSE
                dijqterms = dijqterms_a
              ENDIF

              dufourx = dufourx + ( dijqterme*(nie-nip)*odx_e(i) - &
                                    dijqtermw*(nip-niw)*odx_e(im) )*ayz(ijk)
              dufoury = dufoury + ( dijqtermn*(nin-nip)*ody_n(j) - &
                                    dijqterms*(nip-nis)*ody_n(jm) )*axz(ijk)

              IF(.NOT. no_k) THEN
                 nit = rop_s(ijkt,l)/mi
                 nib = rop_s(ijkb,l)/mi

                 dijqtermt = zero
                 dijqtermb = zero

                 if(rop_s(ijkt,l)/ro_s(ijkt,l) > dil_ep_s) dijqtermt = theta_m(ijk,mmax)**2*dijq(ijk,m,l) / nit
                 if(rop_s(ijkb,l)/ro_s(ijkb,l) > dil_ep_s) dijqtermb = theta_m(ijk,mmax)**2*dijq(ijk,m,l) / nib

                 dijqtermt_h = avg_z_s(dijqterm , dijqtermt, k)
                 dijqtermt_a = avg_z(dijqterm , dijqtermt, k)

                 IF(abs(dijqtermt_h) .le. abs(dijqtermt_a))THEN
                   dijqtermt=dijqtermt_h
                 ELSE
                   dijqtermt=dijqtermt_a
                 ENDIF

                 dijqtermb_h = avg_z_s(dijqtermb, dijqterm, km)
                 dijqtermb_a = avg_z(dijqtermb, dijqterm, km)

                 IF(abs(dijqtermb_h) .le. abs(dijqtermb_a))THEN
                   dijqtermb=dijqtermb_h
                 ELSE
                   dijqtermb=dijqtermb_a
                 ENDIF

                 dufourz = dufourz + ( dijqtermt*(nit-nip)*odz_t(k)*ox(i) - &
                                       dijqtermb*(nip-nib)*odz_t(km)*ox(i) )*axy(ijk)
              ENDIF

!     Sum_ij [ div( Lij*Fj) ]; thermal mobility term
!     Where Fj = Body Force

                lijtermw_h = avg_x_s(lij(ijkw,m,l),lij(ijk,m,l),im)
                lijtermw_a = avg_x(lij(ijkw,m,l),lij(ijk,m,l),im)
              IF(abs(lijtermw_h) .le. abs(lijtermw_a))THEN
                lijtermw = lijtermw_h
              ELSE
                lijtermw = lijtermw_a
              ENDIF

                lijterme_h = avg_x_s(lij(ijk,m,l),lij(ijke,m,l),i)
                lijterme_a = avg_x(lij(ijk,m,l),lij(ijke,m,l),i)

              IF(abs(lijterme_h) .le. abs(lijterme_a))THEN
                lijterme = lijterme_h
              ELSE
                lijterme = lijterme_a
              ENDIF

                lijtermn_h = avg_y_s(lij(ijk,m,l),lij(ijkn,m,l),j)
                lijtermn_a = avg_y(lij(ijk,m,l),lij(ijkn,m,l),j)
              IF(abs(lijtermn_h) .le. abs(lijtermn_a))THEN
                lijtermn = lijtermn_h
              ELSE
                lijtermn = lijtermn_a
              ENDIF

                lijterms_h = avg_y_s(lij(ijks,m,l),lij(ijk,m,l),jm)
                lijterms_a = avg_y(lij(ijks,m,l),lij(ijk,m,l),jm)
              IF(abs(lijterms_h) .le. abs(lijterms_a))THEN
                lijterms = lijterms_h
              ELSE
                lijterms = lijterms_a
              ENDIF

              thermmobilityx = thermmobilityx + ( &
                       fix(ijk,l) *lijterme - fix(imjk,l)*lijtermw )*ayz(ijk)

              thermmobilityy = thermmobilityy + ( &
                       fiy(ijk,l) *lijtermn - fiy(ijmk,l)*lijterms )*axz(ijk)

              IF(.NOT. no_k) THEN

              lijtermt_h = avg_z_s(lij(ijk,m,l),lij(ijkt,m,l),k)
              lijtermt_a = avg_z(lij(ijk,m,l),lij(ijkt,m,l),k)
              IF(abs(lijtermt_h) .le. abs(lijtermt_a))THEN
                lijtermt = lijtermt_h
              ELSE
                lijtermt = lijtermt_a
              ENDIF

              lijtermb_h = avg_z_s(lij(ijkb,m,l),lij(ijk,m,l),km)
              lijtermb_a = avg_z(lij(ijkb,m,l),lij(ijk,m,l),km)
              IF(abs(lijtermb_h) .le. abs(lijtermb_a))THEN
                lijtermb = lijtermb_h
              ELSE
                lijtermb = lijtermb_a
              ENDIF

              thermmobilityz = thermmobilityz + ( &
                       fiz(ijk,l) *lijtermt  - fiz(ijkm,l)*lijtermb )*axy(ijk)
              ENDIF
          ENDDO ! for L = 1, smax

! Additional term arising from subtraction of 3/2*T*continuity
!     + (3/2)*T* Sum_i [ div (joi/mi) ]

          mi = (pi/6.d0)*d_p(ijk,m)**3 * ro_s(ijk,m)

          deldotjoi = deldotjoi + 1.5d0*theta_m(ijk,mmax)/mi * ( &
                      (joix(ijk,m) - joix(imjk,m))*ayz(ijk) + &
                      (joiy(ijk,m) - joiy(ijmk,m))*axz(ijk) + &
                      (joiz(ijk,m) - joiz(ijkm,m))*axy(ijk) )


! Species force dot species mass flux
!     Sum_i [ Fi dot joi/mi ]

!     Calculate species mass flux components at cell center
          joixc = avg_x_e(joix(imjk,m),joix(ijk,m),i)
          joiyc = avg_y_n(joiy(ijmk,m),joiy(ijk,m))
          joizc = avg_z_t(joiz(ijkm,m),joiz(ijk,m))

! external forces evaluated at cell center
          fixc = avg_x_e(fix(imjk,m),fix(ijk,m),i)
          fiyc = avg_y_n(fiy(ijmk,m),fiy(ijk,m))
          fizc = avg_z_t(fiz(ijkm,m),fiz(ijk,m))

          fidotjoi  = fidotjoi  + ( joixc*fixc + joiyc*fiyc + joizc*fizc ) / mi


          IF(switch > zero .AND. abs(ro_g0) .gt. zero) THEN ! do nothing for gran. flow
            ugc = avg_x_e(u_g(imjk),u_g(ijk),i)
            vgc = avg_y_n(v_g(ijmk),v_g(ijk))
            wgc = avg_z_t(w_g(ijkm),w_g(ijk))
            uscm = avg_x_e(u_s(imjk,m),u_s(ijk,m),i)
            vscm = avg_y_n(v_s(ijmk,m),v_s(ijk,m))
            wscm = avg_z_t(w_s(ijkm,m),w_s(ijk,m))

            vslip = (uscm-ugc)**2 + (vscm-vgc)**2 + (wscm-wgc)**2
            vslip = dsqrt(vslip)

! Source/Sink due to fluid do not work well with fluid-solids cases that we run
! uncomment the lines of code below to use (W. Holloway and S. Benyahia).
!
!            call chi_ij_GHD(smax,M,M,SIGMAi,phi_tot,ni,chi_ij)

!            SOURCE_FLUID = SOURCE_FLUID + (81D0*EP_S(IJK,M)*(MU_G(IJK)*VSLIP)**2D0/ &
!                   (chi_ij*(D_P(IJK,M)**3D0*RO_S(IJK,M)*THETA_M(IJK,M))**0.5D0))*VOL(IJK)

!            SINK_FLUID = SINK_FLUID + 3.d0*F_GS(IJK,M)*THETA_M(IJK,M)/Mi
          ENDIF

      ENDDO ! for M = 1, smax

      sink_fluid = sink_fluid/nonzerotheta

      sourcerhs = (pressurerhs + shearproduction + bulkviscrhs + dissdivurhs)*vol(ijk) &
                 + zmax(dufourx)+zmax(dufoury)+zmax(dufourz) &
                 + zmax(thermmobilityx)+zmax(thermmobilityy)+zmax(thermmobilityz)      &
                 + zmax(deldotjoi) + zmax(fidotjoi)*vol(ijk)

      sourcerhs = sourcerhs + source_fluid

      sourcelhs = ( (pressurelhs + bulkvisclhs)/nonzerotheta   + &
                  (colldissipation + dissdivulhs + sink_fluid) ) * vol(ijk) + &
                  ( zmax(-dufourx)+zmax(-dufoury)+zmax(-dufourz) + &
                   zmax(-thermmobilityx)+zmax(-thermmobilityy)+zmax(-thermmobilityz) + &
                   zmax(-deldotjoi) + zmax(-fidotjoi)*vol(ijk) )/ nonzerotheta

      RETURN
      END SUBROUTINE source_ghd_granular_energy