solve_k_epsilon_eq.f

Go to the documentation of this file.

!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
!                                                                      C
!  Subroutine: SOLVE_K_Epsilon_EQ                                      C
!  Purpose: Solve K & Epsilon equations for a turbulent flow           C
!                                                                      C
!                                                                      C
!  Author: S. Benyahia                                Date: MAY-13-04  C
!                                                                      C
!                                                                      C
!  Literature/Document References:                                     C
!                                                                      C
!                                                                      C
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
      SUBROUTINE solve_k_epsilon_eq(IER)

! Modules
!---------------------------------------------------------------------//
      use ambm, only: a_m, b_m, lock_ambm, unlock_ambm
      use bc, only: bc_k_turb_g, bc_e_turb_g
      use compar, only: ijkstart3, ijkend3
      use constant, only: to_si
      use fldvar, only: rop_g, ep_g, u_g, v_g, w_g
      use fldvar, only: k_turb_g, k_turb_go
      use fldvar, only: e_turb_g, e_turb_go
      use functions, only: fluid_at, zmax
      use indices, only: i_of, j_of, k_of
      use geometry, only: vol, ijkmax2
      use leqsol, only: leq_it, leq_method, leq_sweep, leq_pc, leq_tol
      use mflux, only: flux_ge, flux_gn, flux_gt
      use mflux, only: flux_gse, flux_gsn, flux_gst
      use param, only: dimension_bc, dimension_3
      use param1, only: zero
      use residual, only: resid, max_resid, ijk_resid
      use residual, only: num_resid, den_resid
      use residual, only: resid_ke
      use run, only: discretize, k_epsilon, added_mass, odt
      use rxns, only: sum_r_g
      use toleranc, only: zero_ep_s
      use turb, only: dif_k_turb_g, k_turb_g_c, k_turb_g_p
      use turb, only: dif_e_turb_g, e_turb_g_c, e_turb_g_p
      use ur_facs, only: ur_fac
      use usr_src, only: call_usr_source, calc_usr_source
      use usr_src, only: k_epsilon_e, k_epsilon_k
      IMPLICIT NONE

! Dummy arguments
!---------------------------------------------------------------------//
! Error index
      INTEGER, INTENT(INOUT) :: IER

! Local variables
!---------------------------------------------------------------------//
! phase index
      INTEGER :: M
! indices
      INTEGER :: IJK, I, J, K
! loop counter
      INTEGER :: LC
!
      DOUBLE PRECISION :: apo
! linear equation solver method and iterations
      INTEGER :: LEQM, LEQI
! temporary variables in residual computation
      DOUBLE PRECISION :: res1, mres1, num_res, den_res
      INTEGER :: ires1
! A default zero flux will be defined for both K & Epsilon at walls
      DOUBLE PRECISION :: BC_hw_K_Turb_G (dimension_bc)
      DOUBLE PRECISION :: BC_hw_E_Turb_G (dimension_bc)
      DOUBLE PRECISION :: BC_K_Turb_GW (dimension_bc)
      DOUBLE PRECISION :: BC_E_Turb_GW (dimension_bc)
      DOUBLE PRECISION :: BC_C_K_Turb_G (dimension_bc)
      DOUBLE PRECISION :: BC_C_E_Turb_G (dimension_bc)

! small value of K or E, 1 cm2/s2 = 1e-4 m2/s2 = 1e-4 m2/s3
      DOUBLE PRECISION smallTheta

! source vector: coefficient of dependent variable
! becomes part of a_m matrix; must be positive
      DOUBLE PRECISION :: S_P(dimension_3)
! source vector: constant part becomes part of b_m vector
      DOUBLE PRECISION :: S_C(dimension_3)
!
      character(LEN=8) :: Vname
!---------------------------------------------------------------------//

      IF( .NOT. k_epsilon) RETURN

      call lock_ambm

      smalltheta = (to_si)**4 * zero_ep_s
      resid(resid_ke,0) = zero
      num_resid(resid_ke,0) = zero
      den_resid(resid_ke,0) = zero
      max_resid(resid_ke,0) = zero
      ijk_resid(resid_ke,0) = 0

! Setting default zero flux for K & Epsilon since we use wall functions.
! If an expert user want to use Low Re K-Epilon model and needs to set
! the turbulence quantities to zero at walls, then set the hw's to UNDEFINE will
! do it. All the variables below can be changed in the same way as in the
! MFIX data file in the boundary conditions section.
      DO lc = 1, dimension_bc
        bc_hw_k_turb_g(lc) = zero
        bc_k_turb_gw(lc) = zero
        bc_c_k_turb_g(lc) = zero
        bc_hw_e_turb_g(lc) = zero
        bc_e_turb_gw(lc) = zero
        bc_c_e_turb_g(lc) = zero
      ENDDO
! End of setting default zero flux for K & Epsilon wall boundary conditions

! Equations solved for gas phase, thus M = 0
      m = 0
      CALL init_ab_m (a_m, b_m, ijkmax2, m)

! Solve the K_Turb_G equation first
! ---------------------------------------------------------------->>>
      DO ijk = ijkstart3, ijkend3
         IF (fluid_at(ijk)) THEN
            apo = rop_g(ijk)*vol(ijk)*odt
            s_p(ijk) = apo + (zmax(sum_r_g(ijk)) + k_turb_g_p(ijk))*&
                       vol(ijk)
            s_c(ijk) = apo*k_turb_go(ijk) + &
                       k_turb_g(ijk)*zmax((-sum_r_g(ijk)))*vol(ijk) + &
                       k_turb_g_c(ijk)*vol(ijk)
         ELSE
            s_p(ijk) = zero
            s_c(ijk) = zero
         ENDIF
      ENDDO

      IF(.NOT.added_mass) THEN
         CALL conv_dif_phi (k_turb_g, dif_k_turb_g, discretize(9), &
                            u_g, v_g, w_g, flux_ge, flux_gn, flux_gt, &
                            m, a_m, b_m)
      ELSE
         CALL conv_dif_phi (k_turb_g, dif_k_turb_g, discretize(9), &
                            u_g, v_g, w_g, flux_gse, flux_gsn, flux_gst, &
                            m, a_m, b_m)
      ENDIF

      CALL bc_phi (k_turb_g, bc_k_turb_g, bc_k_turb_gw, bc_hw_k_turb_g, &
                   bc_c_k_turb_g, m, a_m, b_m)
      CALL source_phi (s_p, s_c, ep_g, k_turb_g, m, a_m, b_m)

! calculate any usr source terms
      IF (call_usr_source(9)) CALL calc_usr_source(k_epsilon_k,&
                            a_m, b_m, lm=0)

      CALL calc_resid_s (k_turb_g, a_m, b_m, m, num_res, den_res, res1, &
                         mres1, ires1, zero)

      resid(resid_ke,0) = resid(resid_ke,0)+res1
      num_resid(resid_ke,0) = num_resid(resid_ke,0)+num_res
      den_resid(resid_ke,0) = den_resid(resid_ke,0)+den_res
      if(mres1 .gt. max_resid(resid_ke,0)) then
         max_resid(resid_ke,0) = mres1
         ijk_resid(resid_ke,0) = ires1
      endif

      CALL under_relax_s (k_turb_g, a_m, b_m, m, ur_fac(9))

!      call check_ab_m(a_m, b_m, m, .false., ier)
!      call write_ab_m(a_m, b_m, ijkmax2, m, ier)
!      write(*,*) res1, mres1, ires1
!      call test_lin_eq(ijkmax2, ijmax2, imax2, a_m(1, -3, 0), 1, &
!                       DO_K, ier)

      CALL adjust_leq (resid(resid_ke,0), leq_it(9), leq_method(9), &
               leqi, leqm)

      write(vname, '(A,I2)')'K_Turb_G'
      CALL solve_lin_eq (vname, 9, k_turb_g, a_m, b_m, m, leqi, leqm, &
                         leq_sweep(9), leq_tol(9),  leq_pc(9), ier)

!      call out_array(K_Turb_G, Vname)

! remove small negative K values generated by linear solver
! same as adjust_theta.f
      DO ijk = ijkstart3, ijkend3
         IF (fluid_at(ijk)) THEN
            IF(k_turb_g(ijk) < smalltheta) k_turb_g(ijk) = smalltheta
         ENDIF
      ENDDO


! Now solve the E_Turb_G (dissipation) equation.
! ---------------------------------------------------------------->>>
! Initiate (again) the Am Bm matrix.
! This has to be done for every scalar equation.
      CALL init_ab_m (a_m, b_m, ijkmax2, m)

      DO ijk = ijkstart3, ijkend3
         i = i_of(ijk)
         j = j_of(ijk)
         k = k_of(ijk)
         IF (fluid_at(ijk)) THEN
            apo = rop_g(ijk)*vol(ijk)*odt
            s_p(ijk) = apo + (zmax(sum_r_g(ijk)) + e_turb_g_p(ijk))*&
                       vol(ijk)
            s_c(ijk) = apo*e_turb_go(ijk) + &
                       e_turb_g(ijk)*zmax((-sum_r_g(ijk)))*vol(ijk) + &
                       e_turb_g_c(ijk)*vol(ijk)
         ELSE
            s_p(ijk) = zero
            s_c(ijk) = zero
         ENDIF
      ENDDO

      IF(.NOT.added_mass) THEN
         CALL conv_dif_phi (e_turb_g, dif_e_turb_g, discretize(9), &
                            u_g, v_g, w_g, flux_ge, flux_gn, flux_gt, &
                            m, a_m, b_m)
      ELSE
         CALL conv_dif_phi (e_turb_g, dif_e_turb_g, discretize(9), &
                            u_g, v_g, w_g, flux_gse, flux_gsn, flux_gst, &
                            m, a_m, b_m)
      ENDIF

      CALL bc_phi (e_turb_g, bc_e_turb_g, bc_e_turb_gw, bc_hw_e_turb_g, &
                   bc_c_e_turb_g, m, a_m, b_m)

      CALL source_phi (s_p, s_c, ep_g, e_turb_g, m, a_m, b_m)

! calculate any usr source terms
      IF (call_usr_source(9)) CALL calc_usr_source(k_epsilon_e,&
                            a_m, b_m, lm=0)

! set epsilon in fluid cells adjacent to wall cells
      CALL source_k_epsilon_bc(a_m, b_m, m)

      CALL calc_resid_s (e_turb_g, a_m, b_m, m, num_res, den_res, res1, &
                         mres1, ires1, zero)

      resid(resid_ke,0) = resid(resid_ke,0)+res1
      num_resid(resid_ke,0) = num_resid(resid_ke,0)+num_res
      den_resid(resid_ke,0) = den_resid(resid_ke,0)+den_res
      if(mres1 .gt. max_resid(resid_ke,0))then
        max_resid(resid_ke,0) = mres1
        ijk_resid(resid_ke,0) = ires1
      endif

      CALL under_relax_s (e_turb_g, a_m, b_m, m, ur_fac(9))

!      call check_ab_m(a_m, b_m, m, .false., ier)
!      call write_ab_m(a_m, b_m, ijkmax2, m, ier)
!      write(*,*) res1, mres1, ires1
!      call test_lin_eq(ijkmax2, ijmax2, imax2, a_m(1, -3, 0), 1, &
!                       DO_K, ier)

       CALL adjust_leq (resid(resid_ke,0), leq_it(9), leq_method(9), &
                        leqi, leqm)

      write(vname, '(A,I2)')'E_Turb_G'
      CALL solve_lin_eq (vname, 9, e_turb_g, a_m, b_m, m, leqi, leqm, &
                         leq_sweep(9), leq_tol(9), leq_pc(9), ier)

!      call out_array(E_Turb_G, Vname)

! remove small negative Epsilon values generated by linear solver
! same as adjust_theta.f
      DO ijk = ijkstart3, ijkend3
         IF (fluid_at(ijk)) THEN
            IF(e_turb_g(ijk) < smalltheta) e_turb_g(ijk) = smalltheta
         ENDIF
      ENDDO

      call unlock_ambm

      RETURN
      END SUBROUTINE solve_k_epsilon_eq

!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
!                                                                      C
!  Subroutine: SOURCE_K_EPSILON_BC                                     C
!  Purpose: When implementing the wall functions, the epsilon          C
!  (dissipation) value at the fluid cell near the walls needs to be    C
!  set.                                                                C
!                                                                      C
!                                                                      C
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
      SUBROUTINE source_k_epsilon_bc(A_M, B_M, M)

! Modules
!---------------------------------------------------------------------//
      use compar, only: ijkstart3, ijkend3
      use cutcell, only: cut_cell_at, delh_scalar
      use fldvar, only: k_turb_g, e_turb_g
      use functions, only: fluid_at, wall_at
      use functions, only: im_of, jm_of, km_of
      use functions, only: ip_of, jp_of, kp_of
      use geometry, only: dy, odz, ox, dx, cylindrical
      use indices, only: i_of, j_of, k_of
      use param, only: dimension_3, dimension_m
      use param1, only: zero, one
      IMPLICIT NONE

! Dummy Arguments
!---------------------------------------------------------------------//
! Septadiagonal matrix A_m
      DOUBLE PRECISION, INTENT(INOUT) :: A_m(dimension_3, -3:3, 0:dimension_m)
! Vector b_m
      DOUBLE PRECISION, INTENT(INOUT) :: B_m(dimension_3, 0:dimension_m)
! phase index
      INTEGER, INTENT(IN) :: M
! Local variables
!---------------------------------------------------------------------//
      INTEGER :: IJK, I, J, K
!---------------------------------------------------------------------//

      DO ijk = ijkstart3, ijkend3
         i = i_of(ijk)
         j = j_of(ijk)
         k = k_of(ijk)
         IF (fluid_at(ijk)) THEN
            IF(wall_at(jp_of(ijk)).OR.wall_at(jm_of(ijk))) THEN
               a_m(ijk,1,m) = zero
               a_m(ijk,-1,m) = zero
               a_m(ijk,2,m) = zero
               a_m(ijk,-2,m) = zero
               a_m(ijk,3,m) = zero
               a_m(ijk,-3,m) = zero
               a_m(ijk,0,m) = -one
               b_m(ijk,m) =-((0.09d+0)**0.75*k_turb_g(ijk)**1.5)/dy(j) &
                           *2.0d+0/0.42d+0
            ELSEIF(wall_at(kp_of(ijk)).OR.wall_at(km_of(ijk))) THEN
               a_m(ijk,1,m) = zero
               a_m(ijk,-1,m) = zero
               a_m(ijk,2,m) = zero
               a_m(ijk,-2,m) = zero
               a_m(ijk,3,m) = zero
               a_m(ijk,-3,m) = zero
               a_m(ijk,0,m) = -one
               b_m(ijk,m) =-((0.09d+0)**0.75*k_turb_g(ijk)**1.5)* &
                            (odz(k)*ox(i)*2.0d+0)/0.42d+0
            ENDIF  !for identifying wall cells in J or K direction

            IF(cylindrical) THEN
               IF (wall_at(ip_of(ijk)))  THEN
                  a_m(ijk,1,m) = zero
                  a_m(ijk,-1,m) = zero
                  a_m(ijk,2,m) = zero
                  a_m(ijk,-2,m) = zero
                  a_m(ijk,3,m) = zero
                  a_m(ijk,-3,m) = zero
                  a_m(ijk,0,m) = -one
                  b_m(ijk,m) =-((0.09d+0)**0.75*k_turb_g(ijk)**1.5)/&
                               dx(i)*2.0d+0/0.42d+0

               ENDif! for wall cells in I direction
            ELSEIF (wall_at(ip_of(ijk)).OR.wall_at(im_of(ijk))) THEN
               a_m(ijk,1,m) = zero
               a_m(ijk,-1,m) = zero
               a_m(ijk,2,m) = zero
               a_m(ijk,-2,m) = zero
               a_m(ijk,3,m) = zero
               a_m(ijk,-3,m) = zero
               a_m(ijk,0,m) = -one
               b_m(ijk,m) =-((0.09d+0)**0.75*k_turb_g(ijk)**1.5)/dx(i) &
                             *2.0d+0/0.42d+0
            ENDIF ! for cylindrical

            IF(cut_cell_at(ijk)) THEN
               a_m(ijk,1,m) = zero
               a_m(ijk,-1,m) = zero
               a_m(ijk,2,m) = zero
               a_m(ijk,-2,m) = zero
               a_m(ijk,3,m) = zero
               a_m(ijk,-3,m) = zero
               a_m(ijk,0,m) = -one
               b_m(ijk,m) =-((0.09d+0)**0.75*k_turb_g(ijk)**1.5)/&
                            (0.42d+0*delh_scalar(ijk))
            ENDIF
         ENDIF  !for fluid at ijk
      ENDDO   ! enddo ijk
      RETURN
      END SUBROUTINE source_k_epsilon_bc