k_epsilon_prop.f

Go to the documentation of this file.

!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
!                                                                      C
!  Subroutine: K_Epsilon_PROP                                          C
!  Purpose: Calculate diffusion coefficeint and sources for K &        C
!           Epsilon equations                                          C
!                                                                      C
!  Author:                                       Date:                 C
!  Modified: S. Benyahia                         Date:May-13-04        C
!                                                                      C
!                                                                      C
!  Literature/Document References:                                     C
!  Wilcox, D.C., Turbulence Modeling for CFD. DCW Industries, Inc.     C
!     La Canada, Ca. 1994.                                             C
!                                                                      C
!                                                                      C
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
      SUBROUTINE k_epsilon_prop()

! Modules
!---------------------------------------------------------------------//
      USE param
      USE param1
      USE parallel
      USE physprop
      USE drag
      USE run
      USE output
      USE geometry
      USE fldvar
      USE visc_g
      USE visc_s
      USE trace
      USE indices
      USE constant
      Use vshear
      USE turb
      USE toleranc
      USE compar
      USE tau_g
      USE sendrecv

      USE cutcell
      USE fun_avg
      USE functions

      IMPLICIT NONE

! Dummy arguments
!---------------------------------------------------------------------//

! Local parameters
!---------------------------------------------------------------------//
      DOUBLE PRECISION, PARAMETER :: F2O3 = 2.d0/3.d0

! Local variables
!---------------------------------------------------------------------//
! Indices
      INTEGER :: I, J, K, IJK
      INTEGER :: IMJK, IPJK, IJMK, IJPK, IJKM, IJKP
! Loop indices
      INTEGER :: I1, J1, K1
! Solids phase index
      INTEGER :: M
! U_g at the east (i+1/2, j, k) and west face (i-1/2, j, k)
      DOUBLE PRECISION :: UgE, UgW
! U_g at the north (i, j+1/2, k) and south face (i, j-1/2, k)
      DOUBLE PRECISION :: UgN, UgS
! U_g at the top (i, j, k+1/2) and bottom face (i, j, k-1/2)
      DOUBLE PRECISION :: UgT, UgB
! U_g at the center of a scalar cell (i, j, k)
! Calculated for Cylindrical coordinates only.
      DOUBLE PRECISION :: UgcC
! Cell center value of U_g
      DOUBLE PRECISION UGC
! V_g at the east (i+1/2, j, k) and west face (i-1/2, j, k)
      DOUBLE PRECISION :: VgE, VgW
! V_g at the north (i, j+1/2, k) and south face (i, j-1/2, k)
      DOUBLE PRECISION :: VgN, VgS
! V_g at the top (i, j, k+1/2) and bottom face (i, j, k-1/2)
      DOUBLE PRECISION :: VgT, VgB
! Cell center value of V_g
      DOUBLE PRECISION VGC
! W_g at the east (i+1/2, j, k) and west face (i-1/2, j, k)
      DOUBLE PRECISION :: WgE, WgW
! W_g at the north (i, j+1/2, k) and south face (i, j-1/2, k)
      DOUBLE PRECISION :: WgN, WgS
! W_g at the top (i, j, k+1/2) and bottom face (i, j, k-1/2)
      DOUBLE PRECISION :: WgT, WgB
! W_g at the center of a scalar cell (i, j, k).
! Calculated for Cylindrical coordinates only.
      DOUBLE PRECISION :: WgcC
! Cell center value of W_g
      DOUBLE PRECISION WGC
! trace_g and eddy viscosity
      DOUBLE PRECISION :: Trace_G, Mu_gas_t
      DOUBLE PRECISION :: C_Eps_Pope, Xsi_Pope
! gradient of velocity
      DOUBLE PRECISION :: DelV_G(3,3)
! rate of strain tensor
      DOUBLE PRECISION :: D_g(3,3)

!  Production of Turb. Due to shear, Turb Visc, and Constants.
!  See Wilcox PP. 89
      DOUBLE PRECISION Tauij_gDUi_gODxj, C_MU, Sigma_k, Sigma_E, Kappa
      DOUBLE PRECISION Pos_Tauij_gDUi_gODxj, Neg_Tauij_gDUi_gODxj
      DOUBLE PRECISION Ceps_1, Ceps_2, C_Eps_3, Check_Log
      DOUBLE PRECISION Pos_PI_kq_2, Neg_PI_kq_2
! Modif. for Sof Local Var.
      INTEGER :: P,Q
!---------------------------------------------------------------------//

      IF( .NOT. k_epsilon) RETURN

! M should be forced to be equal to one to get some info.
! from solids phase.
      m = 1

! Add constants. Most of these constants have the same names and values
! as the ones defined in Wilcox book (turbulence modeling for CFD).
! Some are necessary only for Simonin turbulence model

      c_mu = 9.0d-02
      kappa = 0.42d+0
      sigma_k = 1.0d0
      sigma_e = 1.3d0
      ceps_1 = 1.44d0 !should be 1.6 for axisymmetric cases with no Pope Correction.
      ceps_2 = 1.92d0
      c_eps_3 = 1.2d0 ! for Simonin model
      c_eps_pope = 0.79d0 ! for Pope Correction.

!!!$omp  parallel do private(ijk, L)
      DO ijk = ijkstart3, ijkend3
         IF (fluid_at(ijk)) THEN
            i = i_of(ijk)
            j = j_of(ijk)
            k = k_of(ijk)

! Velocity components at the scalar faces
            CALL get_face_vel_gas(ijk, uge, ugw, ugn, ugs, ugt, ugb, ugcc, &
                                       vge, vgw, vgn, vgs, vgt, vgb, &
                                       wge, wgw, wgn, wgs, wgt, wgb, wgcc)

! Velocity derivatives (gradient and rate of strain tensor)
            CALL calc_deriv_vel_gas(ijk, delv_g, d_g)

! Value of trace is also available via trd_g but the latter is
! calculated only once per time step
            trace_g = d_g(1,1) + d_g(2,2) + d_g(3,3)

! Pope's correction in 2-D to the Epsilon Equation for a round-Jet
! from Wilcox book, Page 103.
            xsi_pope = zero
            DO i1 = 1,3
               DO j1 = 1,3
                  DO k1 = 1,3
                     xsi_pope = xsi_pope + (delv_g(i1,j1) - delv_g(j1,i1))* &
                                (delv_g(j1,k1) - delv_g(k1,j1))*&
                                (delv_g(k1,i1) + delv_g(i1,k1))
                  ENDDO
               ENDDO
            ENDDO
            xsi_pope = xsi_pope/6. * k_turb_g(ijk)**2/(e_turb_g(ijk)+small_number)

! This IF statment is to ensure that we are using the turbulent viscosity
! and NOT the effective viscosity.
            IF (mu_gt(ijk) .GE. mu_g(ijk)) THEN
               mu_gas_t = mu_gt(ijk) - mu_g(ijk)
            ELSE
               mu_gas_t = zero
            ENDIF

! Calculate Tau(i,j)*dUi/dXj (production term in the K Equation
            IF(.NOT.cut_cell_at(ijk)) THEN
               tauij_gdui_godxj = 2d0*mu_gas_t*(             &
                  d_g(1,1) * d_g(1,1) + &
                  d_g(1,2) * (ugn-ugs)*ody(j) + &
                  d_g(1,3) * ((ugt-ugb)*(ox(i)*odz(k))-wgcc*ox(i)) + &
                  d_g(2,1) * (vge-vgw)*odx(i) + &
                  d_g(2,2) * d_g(2,2) + &
                  d_g(2,3) * (vgt-vgb)*(ox(i)*odz(k)) + &
                  d_g(3,1) * (wge-wgw)*odx(i) + &
                  d_g(3,2) * (wgn-wgs)*ody(j) + &
                  d_g(3,3) * d_g(3,3)) - &
                  f2o3 * ro_g(ijk) * k_turb_g(ijk)*trace_g - &
                  f2o3 * mu_gas_t * trace_g**2

            ELSE  ! CUT_CELL
! This is actually not used because of wall functions in cut cells
!               Tauij_gDUi_gODxj = 2D0*Mu_gas_t*(                             &
!                  D_G(1,1) * UG(1,1)  +                      &
!                  D_G(1,2) * UG(1,2)  +                      &
!                  D_G(1,3) * UG(1,3)  +                      &
!                  D_G(2,1) * UG(2,1)  +                      &
!                  D_G(2,2) * UG(2,2)  +                      &
!                  D_G(2,3) * UG(2,3)  +                      &
!                  D_G(3,1) * UG(3,1)  +                      &
!                  D_G(3,2) * UG(3,2)  +                      &
!                  D_G(3,3) * UG(3,3)) -                      &
!                  F2O3 * RO_G(IJK) * K_Turb_G(IJK)*Trace_g-  &
!                  F2O3 * Mu_gas_t * Trace_g**2
            ENDIF


! To avoid very small negative numbers
            IF (tauij_gdui_godxj .GE. zero) THEN
               pos_tauij_gdui_godxj = tauij_gdui_godxj
               neg_tauij_gdui_godxj = zero
            ELSE
               pos_tauij_gdui_godxj = zero
               neg_tauij_gdui_godxj = tauij_gdui_godxj
            ENDIF


! Interaction terms in the K-Epsilon equations FOR USE Simonin and Ahmadi models
            IF(simonin) THEN
               pos_pi_kq_2 = f_gs(ijk,1)*k_12(ijk)
               neg_pi_kq_2 = f_gs(ijk,1)*2.0d0* k_turb_g(ijk)
            ELSE IF(ahmadi) THEN
               pos_pi_kq_2 = f_gs(ijk,1)*3.0d0*theta_m(ijk,m)
               neg_pi_kq_2 = f_gs(ijk,1)*2.0d0* k_turb_g(ijk)
               c_eps_3 = zero ! no extra terms in epsilon equation for Ahmadi model !
            ELSE
               pos_pi_kq_2 = zero
               neg_pi_kq_2 = zero
            ENDIF


            IF(k_turb_g(ijk) > small_number .AND. &
               e_turb_g(ijk) > small_number) THEN

! Start Adding source terms to the K equation
               k_turb_g_c(ijk) = (ep_g(ijk)*pos_tauij_gdui_godxj +  &
                                   pos_pi_kq_2 )
               k_turb_g_p(ijk) =(-ep_g(ijk)*neg_tauij_gdui_godxj +  &
                  ep_g(ijk)*ro_g(ijk)*e_turb_g(ijk)+neg_pi_kq_2)/ &
                  k_turb_g(ijk)


! Calculate velocity components at i, j, k to be used in wall functions
            imjk = im_of(ijk)
            ipjk = ip_of(ijk)
            ijmk = jm_of(ijk)
            ijpk = jp_of(ijk)
            ijkm = km_of(ijk)
            ijkp = kp_of(ijk)
            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))

! Implementing wall functions targeted to fluid cells next to walls...
! Setting Source and sink terms in the K equation since the production
! in the K eq. due to shear should include the LOG law of the wall.
! North/South wall
               IF(wall_at(ijpk).OR.wall_at(ijmk)) THEN
                  check_log = log(9.81*c_mu**0.25*                     &
                     ro_g(ijk)*sqrt(k_turb_g(ijk))*dy(j)/2.0d+0/       &
                     mu_g(ijk))
                  IF(check_log .LE. zero) THEN
                     k_turb_g_c(ijk) = zero
                     k_turb_g_p(ijk) = zero
                  ELSE
                     k_turb_g_c(ijk) = sqrt(c_mu) * 2.d+0/dy(j)*      &
                        max(abs(ugc),abs(wgc))*ep_g(ijk)*ro_g(ijk)*    &
                        k_turb_g(ijk)/check_log
                     k_turb_g_p(ijk) = ((ep_g(ijk)*ro_g(ijk)*         &
                        c_mu**0.75*k_turb_g(ijk)**1.5/dy(j)*           &
                        2.0d+0/kappa))/k_turb_g(ijk)
                  ENDIF !for check_log less than zero
! Top/Bottom wall
               ELSEIF(wall_at(ijkp).OR.wall_at(ijkm)) THEN
                  check_log = log(9.81*c_mu**0.25*                     &
                     ro_g(ijk)*sqrt(k_turb_g(ijk))/                    &
                     (odz(k)*ox(i)*2.0d+0)/mu_g(ijk))
                  IF(check_log .LE. zero) THEN
                     k_turb_g_c(ijk) = zero
                     k_turb_g_p(ijk) = zero
                  ELSE
                     k_turb_g_c(ijk) = sqrt(c_mu)*(odz(k)*ox(i)*      &
                        2.0d+0)* max(abs(ugc),abs(vgc))*ep_g(ijk)*     &
                        ro_g(ijk)*k_turb_g(ijk)/check_log
                     k_turb_g_p(ijk) = ((ep_g(ijk)*ro_g(ijk)*         &
                        c_mu**0.75*k_turb_g(ijk)**1.5/kappa*           &
                        (odz(k)*ox(i)*2.0d+0)))/k_turb_g(ijk)
                  ENDIF !for check_log less than zero
               ENDIF !For walls

! For Cylindrical cases, wall_at (IP) is a wall cell, but wall_at (IM) is
! the axis of symmetry where wall functions obviously don't apply.
               IF(cylindrical) THEN
                  IF (wall_at(ipjk))  THEN
                     check_log = log(9.81*c_mu**0.25* ro_g(ijk)*       &
                     sqrt(k_turb_g(ijk))*dx(i)/2.0d+0/mu_g(ijk))
                     IF(check_log .LE. zero) THEN
                        k_turb_g_c(ijk) = zero
                        k_turb_g_p(ijk) = zero
                     ELSE
                        k_turb_g_c(ijk) = sqrt(c_mu)*2.d+0/dx(i)*     &
                           max(abs(vgc),abs(wgc)) *ep_g(ijk)*ro_g(ijk)*&
                           k_turb_g(ijk)/check_log
                        k_turb_g_p(ijk) = ((ep_g(ijk)*ro_g(ijk)*      &
                           c_mu**0.75*k_turb_g(ijk)**1.5/dx(i)*        &
                           2.0d+0/kappa))/ k_turb_g(ijk)
                     ENDIF !for check_log less than zero
                  ENDIF  ! for wall cells in I direction
! East/West wall
               ELSEIF (wall_at(ipjk).OR.wall_at(imjk)) THEN
                  check_log = log(9.81*c_mu**0.25*ro_g(ijk)*           &
                     sqrt(k_turb_g(ijk))*dx(i)/2.0d+0/mu_g(ijk))
                  IF(check_log .LE. zero) THEN
                     k_turb_g_c(ijk) = zero
                     k_turb_g_p(ijk) = zero
                  ELSE
                     k_turb_g_c(ijk) = sqrt(c_mu)*2.d+0/dx(i)*        &
                        max(abs(vgc),abs(wgc))*ep_g(ijk)*ro_g(ijk)*    &
                        k_turb_g(ijk)/check_log
                     k_turb_g_p(ijk) = ((ep_g(ijk)*ro_g(ijk)*         &
                        c_mu**0.75*k_turb_g(ijk)**1.5/dx(i)*           &
                        2.0d+0/kappa))/k_turb_g(ijk)
                  ENDIF !for check_log less than zero
               ENDIF ! for cylindrical

               IF(cut_cell_at(ijk)) THEN
                  check_log = log(9.81*c_mu**0.25* ro_g(ijk)*          &
                  sqrt(k_turb_g(ijk))*delh_scalar(ijk)/mu_g(ijk))
                  IF(check_log .LE. zero) THEN
                     k_turb_g_c(ijk) = zero
                     k_turb_g_p(ijk) = zero
                  ELSE
!                     K_Turb_G_c (IJK) = SQRT(C_mu)/DELH_Scalar(IJK)*   &
!                        MAX(ABS(UGC),ABS(VGC),ABS(WGC))*               &
!                        EP_g(IJK)*RO_G(IJK)*K_Turb_G(IJK) /Check_Log
                     k_turb_g_c(ijk) = sqrt(c_mu)/delh_scalar(ijk)*   &
                        dsqrt(ugc**2 + vgc**2 + wgc**2) *              &
                        ep_g(ijk)*ro_g(ijk)*k_turb_g(ijk)/check_log

                     k_turb_g_p(ijk) = (ep_g(ijk)*ro_g(ijk)*          &
                        c_mu**0.75*k_turb_g(ijk)**1.5)/                &
                       (delh_scalar(ijk)*kappa)/k_turb_g(ijk)
                  ENDIF !for check_log less than zero
               ENDIF


! Diffusion coefficient for turbulent kinetic energy (K)
               dif_k_turb_g(ijk) = ep_g(ijk)* &
                  (mu_g(ijk) + mu_gas_t /sigma_k)

! Add here Dissipation of Turbulence source terms
               e_turb_g_c(ijk) = (ceps_1 *&
                  ep_g(ijk)*pos_tauij_gdui_godxj*                      &
                  e_turb_g(ijk)/k_turb_g(ijk) +                        &
                  c_eps_3*pos_pi_kq_2*e_turb_g(ijk)/k_turb_g(ijk))

! Pope Correction in Xsi_Pope, Add it to E_Turb_G_c to use this option.
                  ! + C_Eps_Pope*RO_G(IJK)*EP_g(IJK)*&
                   !  ZMAX(Xsi_Pope)

               e_turb_g_p(ijk) = -ceps_1 * &
                  ep_g(ijk)*neg_tauij_gdui_godxj /k_turb_g(ijk) +      &
                  ceps_2 * ep_g(ijk) *ro_g(ijk)*                       &
                  e_turb_g(ijk)/k_turb_g(ijk) +                        &
                  c_eps_3*(neg_pi_kq_2)/k_turb_g(ijk)

! Diffusion coefficient for Dissipation of turbulent energy (Epsilon)
               dif_e_turb_g(ijk) =ep_g(ijk)* &
                  (mu_g(ijk) + mu_gas_t /sigma_e)

            ELSE
               k_turb_g_c(ijk) = zero
               k_turb_g_p(ijk) = zero
               e_turb_g_c(ijk) = zero
               e_turb_g_p(ijk) = zero
               dif_k_turb_g(ijk) = zero
               dif_e_turb_g(ijk) = zero
            ENDIF !for K_turb_G and E_Turb_G having very small numbers
         ENDIF
      ENDDO

      RETURN
      END SUBROUTINE k_epsilon_prop