calc_external_forces.f

Go to the documentation of this file.

!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
!                                                                      C
!  Module name: calc_external_forces(IER)                              C
!  Purpose: Calculate the 3-components of external forces at cell      C
!           faces                                                      C
!                                                                      C
!  Author: S. Benyahia                              Date: 13-MAY-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:   JoiX  JoiY   JoiZ                             C
!                                                                      C
!     Local variables: all terms in mass flux                          C
!                                                                      C
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
!
      SUBROUTINE calc_external_forces()
!
!-----------------------------------------------
!     Modules
!-----------------------------------------------
      USE param
      USE param1
      USE geometry
      USE compar
      USE fldvar
      USE indices
      USE ghdtheory
      USE physprop
      USE run
      USE constant
      USE drag
      USE bc
      use scales
      use bodyforce
      USE fun_avg
      USE functions
      IMPLICIT NONE
!-----------------------------------------------
!     Local variables
!-----------------------------------------------
!                      Index
      INTEGER          IJK, I, J, K
!
!                      Index
      INTEGER          IJKE, IJKN, IJKT
!
!                      Solids phase
      INTEGER          M ,IM
!
!     number densities to compute del(Nj)
      DOUBLE PRECISION NjC, NjE, NjN, NjT
!
!     mass, volume of species
      DOUBLE PRECISION Mj, Vj
!
!     drag force on a particle
      DOUBLE PRECISION dragFc, dragFe, dragFn, dragFt
!
!     pressure terms in mass mobility
      DOUBLE PRECISION PGE, PGN, PGT, SDPx, SDPy, SDPz

!     off-diagonal terms for HYS drag relation
      DOUBLE PRECISION  avgDragx, avgDragy, avgDragz
!
!-----------------------------------------------
!     Function subroutines
!-----------------------------------------------

      DO 200 ijk = ijkstart3, ijkend3
          i = i_of(ijk)
          j = j_of(ijk)
          k = k_of(ijk)


          IF ( fluid_at(ijk) ) THEN

               ijke = east_of(ijk)
               ijkn = north_of(ijk)
               ijkt = top_of(ijk)

! pressure term (no need to have it inside smax loop)
               pge = p_g(ijke)
               pgn = p_g(ijkn)
               pgt = p_g(ijkt)
               IF (cyclic_x_pd) THEN
                 IF (cyclic_at_e(ijk)) pge = p_g(ijke) - delp_x
               ENDIF
               IF (cyclic_y_pd) THEN
                 IF (cyclic_at_n(ijk)) pgn = p_g(ijkn) - delp_y
               ENDIF
               IF (cyclic_z_pd) THEN
                 IF (cyclic_at_t(ijk)) pgt = p_g(ijkt) - delp_z
               ENDIF
               sdpx = -p_scale*(pge - p_g(ijk))  * odx_e(i)
               sdpy = -p_scale*(pgn - p_g(ijk))  * ody_n(j)
               sdpz = -p_scale*(pgt - p_g(ijk))  * (ox_e(i)*odz_t(k))

               DO m = 1, smax
                 vj = (pi/6.d0)*d_p(ijk,m)**3 ! particle volume
                 mj = vj * ro_s(ijk,m)            ! particle mass

                 njc = rop_s(ijk,m) / mj
                 nje = rop_s(ijke,m) / mj
                 njn = rop_s(ijkn,m) / mj
                 njt = rop_s(ijkt,m) / mj

! drag force on a particle in -x -y -z directions
                 dragfc = zero
                 dragfe = zero
                 dragfn = zero
                 dragft = zero
                 if(njc > zero) dragfc = f_gs(ijk ,m)/njc
                 if(nje > zero) dragfe = f_gs(ijke,m)/nje
                 if(njn > zero) dragfn = f_gs(ijkn,m)/njn
                 if(njt > zero) dragft = f_gs(ijkt,m)/njt

                 dragfxflux(ijk,m) = avg_x(dragfc,dragfe,i) * (u_g(ijk) - u_s(ijk,m))
                 dragfyflux(ijk,m) = avg_y(dragfc,dragfn,j) * (v_g(ijk) - v_s(ijk,m))
                 dragfzflux(ijk,m) = avg_z(dragfc,dragft,k) * (w_g(ijk) - w_s(ijk,m))

                 dragfx(ijk,m) = avg_x(dragfc,dragfe,i) * (u_g(ijk))
                 dragfy(ijk,m) = avg_y(dragfc,dragfn,j) * (v_g(ijk))
                 dragfz(ijk,m) = avg_z(dragfc,dragft,k) * (w_g(ijk))

                 beta_cell_x(ijk,m)=avg_x(dragfc,dragfe,i)
                 beta_cell_y(ijk,m)=avg_y(dragfc,dragfn,j)
                 beta_cell_z(ijk,m)=avg_z(dragfc,dragft,k)

                 IF(drag_type_enum == hys)THEN
                    DO im=1,smax
                       IF(im /= m)THEN

! HYS additional drag force on a particle in -x -y -z directions
                          dragfc = zero
                          dragfe = zero
                          dragfn = zero
                          dragft = zero
                          if(njc > zero) dragfc = beta_ij(ijk,m,im) /njc
                          if(nje > zero) dragfe = beta_ij(ijke,m,im)/nje
                          if(njn > zero) dragfn = beta_ij(ijkn,m,im)/njn
                          if(njt > zero) dragft = beta_ij(ijkt,m,im)/njt

                          avgdragx = avg_x(dragfc,dragfe,i)
                          dragfx(ijk,m) = dragfx(ijk,m) - avgdragx*(u_g(ijk))
                          dragfxflux(ijk,m) = dragfxflux(ijk,m) - avgdragx*(u_g(ijk) - u_s(ijk,im))
                          beta_ij_cell_x(ijk,m,im)= avg_x(dragfc,dragfe,i)

                          avgdragy = avg_y(dragfc,dragfn,j)
                          dragfy(ijk,m) = dragfy(ijk,m) - avgdragy*(v_g(ijk))
                          dragfyflux(ijk,m) = dragfyflux(ijk,m) - avgdragy*(v_g(ijk) - v_s(ijk,im))
                          beta_ij_cell_y(ijk,m,im)=avg_y(dragfc,dragfn,j)

                          avgdragz = avg_z(dragfc,dragft,k)
                          dragfz(ijk,m) = dragfz(ijk,m) - avgdragz*(w_g(ijk))
                          dragfzflux(ijk,m) = dragfzflux(ijk,m) - avgdragz*(w_g(ijk) - w_s(ijk,im))
                          beta_ij_cell_z(ijk,m,im)=avg_z(dragfc,dragft,k)

                       ENDIF
                    ENDDO
                 ENDIF


                 fixvel(ijk,m) =  (mj * bfx_s(ijk,m)+dragfx(ijk,m) +vj*sdpx)
                 fiyvel(ijk,m) =  (mj * bfy_s(ijk,m)+dragfy(ijk,m) +vj*sdpy)
                 fizvel(ijk,m) =  (mj * bfz_s(ijk,m)+dragfz(ijk,m) +vj*sdpz)

                 fix(ijk,m) =  (mj * bfx_s(ijk,m)+dragfxflux(ijk,m) +vj*sdpx)
                 fiy(ijk,m) =  (mj * bfy_s(ijk,m)+dragfyflux(ijk,m) +vj*sdpy)
                 fiz(ijk,m) =  (mj * bfz_s(ijk,m)+dragfzflux(ijk,m) +vj*sdpz)

                 fiminusdragx(ijk,m) =  (mj * bfx_s(ijk,m) + vj*sdpx)
                 fiminusdragy(ijk,m) =  (mj * bfy_s(ijk,m) + vj*sdpy)
                 fiminusdragz(ijk,m) =  (mj * bfz_s(ijk,m) + vj*sdpz)
               ENDDO
          ENDIF     ! Fluid_at
 200  CONTINUE     ! outer IJK loop

      RETURN
      END SUBROUTINE calc_external_forces