Mohamed particles

Sources/Process/Backup_Mod/Load.f90

@@ -184,12 +184,11 @@ subroutine Backup_Mod_Load(fld, time_step, time_step_stat, backup)
     ! Turbulence quantities connected with heat transfer
   end if
 
-  if(turbulence_model .eq. K_EPS_ZETA_F .and. heat_transfer) then
+  if( (turbulence_model .eq. K_EPS_ZETA_F .and. heat_transfer) .or. &
+      (turbulence_model .eq. HYBRID_LES_RANS .and. heat_transfer) ) then
     call Backup_Mod_Read_Variable(fh, d, vc, 't2',       t2)
     call Backup_Mod_Read_Cell_Bnd(fh, d, vc, 'p_t2',     p_t2    (-nb_s:nc_s))
     call Backup_Mod_Read_Cell_Bnd(fh, d, vc, 'con_wall', con_wall(-nb_s:nc_s))
-  else if (turbulence_model .eq. HYBRID_LES_RANS .and. heat_transfer) then
-    call Backup_Mod_Read_Cell_Bnd(fh, d, vc, 'con_wall', con_wall(-nb_s:nc_s))
   end if
 
   !----------------------------!

Sources/Process/Backup_Mod/Save.f90

@@ -163,12 +163,11 @@ subroutine Backup_Mod_Save(fld, time_step, time_step_stat, name_save)
     call Backup_Mod_Write_Cell_Bnd(fh, d, vc, 'l_scale',  l_scale(-nb_s:nc_s))
   end if
 
-  if(turbulence_model .eq. K_EPS_ZETA_F .and. heat_transfer) then
+  if( (turbulence_model .eq. K_EPS_ZETA_F .and. heat_transfer) .or. &
+      (turbulence_model .eq. HYBRID_LES_RANS .and. heat_transfer) ) then
     call Backup_Mod_Write_Variable(fh, d, vc, 't2',       t2)
     call Backup_Mod_Write_Cell_Bnd(fh, d, vc, 'p_t2',     p_t2    (-nb_s:nc_s))
     call Backup_Mod_Write_Cell_Bnd(fh, d, vc, 'con_wall', con_wall(-nb_s:nc_s))
-  else if (heat_transfer .and. turbulence_model .eq. HYBRID_LES_RANS) then
-    call Backup_Mod_Write_Cell_Bnd(fh, d, vc, 'con_wall', con_wall(-nb_s:nc_s))
   end if
 
 

Sources/Process/Initialize_Variables.f90

@@ -158,7 +158,8 @@ subroutine Initialize_Variables(flow)
             i=Key_Ind('F22', keys,nks);prof(k,0)=f22_def; f22 %n(c)=prof(k,i)
           end if
 
-          if(turbulence_model .eq. K_EPS_ZETA_F .and. heat_transfer) then
+          if( (turbulence_model .eq. K_EPS_ZETA_F .and. heat_transfer) .or. &
+              (turbulence_model .eq. HYBRID_LES_RANS .and. heat_transfer) ) then
             i=Key_Ind('T2', keys,nks);prof(k,0)=t2_def; t2 %n(c)=prof(k,i)
           end if
 
@@ -309,8 +310,8 @@ subroutine Initialize_Variables(flow)
           y_plus(c) = 0.001
         end if
 
-        if(turbulence_model .eq. K_EPS_ZETA_F .and. &
-           heat_transfer) then
+          if( (turbulence_model .eq. K_EPS_ZETA_F .and. heat_transfer) .or. &
+              (turbulence_model .eq. HYBRID_LES_RANS .and. heat_transfer) ) then
           vals(0) = t2_def;  t2 % n(c) = vals(Key_Ind('T2',  keys, nks))
           t2 % o(c)  = t2 % n(c)
           t2 % oo(c) = t2 % n(c)

Sources/Process/Load_Boundary_Conditions.f90

@@ -281,7 +281,10 @@ subroutine Load_Boundary_Conditions(flow, backup)
               i = Key_Ind('F22',  keys, nks); if(i > 0) f22  % n(c) = vals(i)
             end if
 
-            if(turbulence_model .eq. K_EPS_ZETA_F .and. heat_transfer) then
+            if( (turbulence_model .eq. K_EPS_ZETA_F .and.    &
+                 heat_transfer) .or.                         &
+                (turbulence_model .eq. HYBRID_LES_RANS .and. &
+                  heat_transfer) )                           then
               i = Key_Ind('T2',  keys, nks); if(i > 0) t2  % n(c) = vals(i)
             end if
 
@@ -420,7 +423,10 @@ subroutine Load_Boundary_Conditions(flow, backup)
                 i = Key_Ind('F22',  keys, nks); if(i>0) f22  % n(c) = prof(k,i)
               end if
 
-              if(turbulence_model .eq. K_EPS_ZETA_F .and. heat_transfer) then
+              if( (turbulence_model .eq. K_EPS_ZETA_F .and.    &
+                   heat_transfer) .or.                         &
+                  (turbulence_model .eq. HYBRID_LES_RANS .and. &
+                    heat_transfer) )                           then
                 i = Key_Ind('T2',  keys, nks); if(i>0) t2  % n(c) = prof(k,i)
               end if
 
@@ -640,7 +646,10 @@ subroutine Load_Boundary_Conditions(flow, backup)
                     if(i > 0) f22 % n(c) = wi*prof(m,i) + (1.-wi)*prof(m+1,i)
                   end if
 
-                  if(turbulence_model .eq. K_EPS_ZETA_F .and. heat_transfer) then
+                  if( (turbulence_model .eq. K_EPS_ZETA_F .and.    &
+                       heat_transfer) .or.                         &
+                      (turbulence_model .eq. HYBRID_LES_RANS .and. &
+                        heat_transfer) )                           then
                     i = Key_Ind('T2',keys,nks)
                     if(i > 0) t2 % n(c) = wi*prof(m,i) + (1.-wi)*prof(m+1,i)
                   end if

Sources/Process/Save_Cgns_Results.f90

@@ -186,7 +186,8 @@ subroutine Save_Results(flow, name_save)
                               f22 % n(1),  "TurbulentQuantityF22")
   end if
 
-  if(turbulence_model .eq. K_EPS_ZETA_F .and. heat_transfer) then
+  if( (turbulence_model .eq. K_EPS_ZETA_F .and. heat_transfer) .or. &
+      (turbulence_model .eq. HYBRID_LES_RANS .and. heat_transfer) ) then
     call Cgns_Mod_Write_Field(base, block, solution, field, grid, &
                               t2 % n(1),  "TurbulentQuantityT2")
     call Cgns_Mod_Write_Field(base, block, solution, field, grid, &

Sources/Process/Save_Vtu_Results.f90

@@ -263,10 +263,10 @@ subroutine Save_Results(flow, name_save)
     call Save_Vtu_Scalar(grid, IN_4, IN_5, "TurbulentQuantityF22",  f22  % n(1))
   end if
 
-  if (turbulence_model .eq. K_EPS_ZETA_F .and. heat_transfer) then
+  if( (turbulence_model .eq. K_EPS_ZETA_F .and. heat_transfer) .or. &
+      (turbulence_model .eq. HYBRID_LES_RANS .and. heat_transfer) ) then
     call Save_Vtu_Scalar(grid, IN_4, IN_5, "TurbulentQuantityT2", t2 % n(1))
-    call Save_Vtu_Scalar(grid, IN_4, IN_5, "TurbulentT2Production", &
-                                           p_t2(1))
+    call Save_Vtu_Scalar(grid, IN_4, IN_5, "TurbulentT2Production", p_t2(1))
   end if
   ! Save vis and vis_t
   if(turbulence_model .eq. DES_SPALART .or.  &

Sources/Process/Swarm_Mod.f90

@@ -9,7 +9,7 @@ module Swarm_Mod
 !------------------------------------------------------------------------------!
   implicit none
 !==============================================================================!
-
+! Flag for branch "Mohamed_Particles"
   !-------------------!
   !   Particle type   !
   !-------------------!

Sources/Process/Turbulence/Allocate.f90

@@ -123,7 +123,7 @@ subroutine Turbulence_Allocate(flow)
 
     ! Hydraulic roughness given by formula
     if(rough_walls) then
-      allocate(z_o_f(-grid % n_bnd_cells:grid % n_cells));  z_o_f   = 0.
+      allocate(z_o_f(-grid % n_bnd_cells:grid % n_cells));  z_o_f   = -1.
     end if
 
     if(heat_transfer) then
@@ -228,6 +228,13 @@ subroutine Turbulence_Allocate(flow)
 
     if(turbulence_statistics) then
 
+      ! First moments
+      call Var_Mod_Allocate_Statistics(u)
+      call Var_Mod_Allocate_Statistics(v)
+      call Var_Mod_Allocate_Statistics(w)
+      call Var_Mod_Allocate_Statistics(p)
+
+      ! Second moments
       call Var_Mod_Allocate_Statistics(uu)
       call Var_Mod_Allocate_Statistics(vv)
       call Var_Mod_Allocate_Statistics(ww)
@@ -565,10 +572,12 @@ subroutine Turbulence_Allocate(flow)
 
     if(heat_transfer) then
 
+      call Var_Mod_Allocate_Solution('T2', '', t2, grid)
       call Var_Mod_Allocate_New_Only('UT', ut, grid)
       call Var_Mod_Allocate_New_Only('VT', vt, grid)
       call Var_Mod_Allocate_New_Only('WT', wt, grid)
       allocate(con_wall(-grid % n_bnd_cells:grid % n_cells)); con_wall = 0.
+      allocate(p_t2    (-grid % n_bnd_cells:grid % n_cells)); p_t2     = 0.
 
     end if ! heat_transfer
 
@@ -603,6 +612,7 @@ subroutine Turbulence_Allocate(flow)
         call Var_Mod_Allocate_Statistics(wt)  ! new value allocated above
         call Var_Mod_Allocate_New_Only('TT', tt,  grid)
         call Var_Mod_Allocate_Statistics(tt)
+        call Var_Mod_Allocate_Statistics(t2)
 
       end if ! heat_transfer
 

Sources/Process/Turbulence/Calculate_Vis_T_K_Eps.f90

@@ -130,8 +130,7 @@ subroutine Calculate_Vis_T_K_Eps(flow)
                                           grid % wall_dist(c1),  &
                                           kin_vis)
           u_plus     = U_Plus_Rough_Walls(grid % wall_dist(c1))
-          vis_wall(c1) = y_plus(c1) * viscosity * kappa  &
-                       / log((grid % wall_dist(c1)+z_o)/z_o)  ! is this U+?
+          vis_wall(c1) = y_plus(c1) * viscosity / u_plus
         end if
 
         if(heat_transfer) then

Sources/Process/Turbulence/Calculate_Vis_T_K_Eps_Zeta_F.f90

@@ -119,12 +119,12 @@ subroutine Calculate_Vis_T_K_Eps_Zeta_F(flow)
 
         if(rough_walls) then
           z_o = Roughness_Coefficient(grid, z_o_f(c1), c1)      
+!          z_o = max(grid % wall_dist(c1)/(e_log*y_plus(c1)),z_o) 
           y_plus(c1) = Y_Plus_Rough_Walls(u_tau(c1),             &
                                           grid % wall_dist(c1),  &
                                           kin_vis)
           u_plus     = U_Plus_Rough_Walls(grid % wall_dist(c1))
-          vis_wall(c1) = y_plus(c1) * viscosity * kappa  &
-                       / log((grid % wall_dist(c1)+z_o)/z_o)  ! is this U+?
+          vis_wall(c1) = y_plus(c1) * viscosity / u_plus
         end if  
 
         if(heat_transfer) then

Sources/Process/Turbulence/Compute_Turbulent.f90

@@ -212,7 +212,8 @@ subroutine Compute_Turbulent(flow, sol, dt, ini, phi, n_step)
     if(phi % name .eq. 'ZETA') call Source_Zeta_K_Eps_Zeta_F(flow, sol, n_step)
   end if
 
-  if(turbulence_model .eq. K_EPS_ZETA_F .and. heat_transfer) then
+  if( (turbulence_model .eq. K_EPS_ZETA_F .and. heat_transfer) .or. &
+      (turbulence_model .eq. HYBRID_LES_RANS .and. heat_transfer) ) then
     if(phi % name .eq. 'T2')  call Source_T2(flow, sol)
   end if
 
@@ -242,7 +243,7 @@ subroutine Compute_Turbulent(flow, sol, dt, ini, phi, n_step)
 
   do c = 1, grid % n_cells
     if( phi % n(c) < 0.0 ) phi % n(c) = phi % o(c)
-    if(phi % name .eq. 'ZETA')  phi % n(c) = min(phi % n(c), 1.8) 
+    if(phi % name .eq. 'ZETA')  phi % n(c) = min(phi % n(c), 1.8)
   end do
 
   ! Print info on the screen
@@ -257,7 +258,8 @@ subroutine Compute_Turbulent(flow, sol, dt, ini, phi, n_step)
       call Info_Mod_Iter_Fill_At(3, 3, phi % name, exec_iter, phi % res)
   end if
 
-  if(turbulence_model .eq. K_EPS_ZETA_F .and. heat_transfer) then
+  if( (turbulence_model .eq. K_EPS_ZETA_F .and. heat_transfer) .or. &
+      (turbulence_model .eq. HYBRID_LES_RANS .and. heat_transfer) ) then
     if(phi % name .eq. 'T2')  &
       call Info_Mod_Iter_Fill_At(3, 5, phi % name, exec_iter, phi % res)
   end if

Sources/Process/Turbulence/Roughness_Coefficient.f90

@@ -17,7 +17,7 @@ real function Roughness_Coefficient(grid, z_o_function, c)
 
   Roughness_Coefficient = z_o
 
-  if(z_o_function > tiny) then
+  if(z_o_function > -tiny) then
     Roughness_Coefficient = z_o_function
   end if
 

Sources/Process/Turbulence/Source_Eps_K_Eps.f90

@@ -35,7 +35,7 @@ subroutine Source_Eps_K_Eps(flow, sol)
   integer                    :: s, c, c1, c2, j
   real                       :: u_tan, u_nor_sq, u_nor, u_tot_sq
   real                       :: re_t, f_mu, u_tau_new, fa, kin_vis
-  real                       :: eps_wf, eps_int, ebf, y_star
+  real                       :: eps_wf, eps_int, y_star
 !==============================================================================!
 !   Dimensions:                                                                !
 !                                                                              !
@@ -116,7 +116,7 @@ subroutine Source_Eps_K_Eps(flow, sol)
         end if
 
         if(rough_walls) then 
-          z_o = Roughness_Coefficient(grid, z_o_f(c1), c1)       
+          z_o = Roughness_Coefficient(grid, z_o_f(c1), c1)     
           eps % n(c1) = c_mu75 * kin % n(c1)**1.5  &
                       / ((grid % wall_dist(c1) + z_o) * kappa)
 
@@ -136,7 +136,6 @@ subroutine Source_Eps_K_Eps(flow, sol)
 
           u_tau_new = sqrt(tau_wall(c1)/density)
           y_plus(c1) = Y_Plus_Low_Re(u_tau_new, grid % wall_dist(c1), kin_vis)
-          ebf = 0.01 * y_plus(c1)**4 / (1.0 + 5.0*y_plus(c1))
 
           eps_int = 2.0*viscosity/density * kin % n(c1)    &
                   / grid % wall_dist(c1)**2

Sources/Process/Turbulence/Source_Eps_K_Eps_Zeta_F.f90

@@ -27,7 +27,7 @@ subroutine Source_Eps_K_Eps_Zeta_F(flow, sol)
   real,              pointer :: b(:)
   integer                    :: c, s, c1, c2, j
   real                       :: u_tan, u_nor_sq, u_nor, u_tot_sq
-  real                       :: e_sor, c_11e, ebf
+  real                       :: e_sor, c_11e
   real                       :: eps_wf, eps_int
   real                       :: fa, u_tau_new, kin_vis
 !==============================================================================!
@@ -135,8 +135,6 @@ subroutine Source_Eps_K_Eps_Zeta_F(flow, sol)
           u_tau_new = sqrt(tau_wall(c1)/density)
           y_plus(c1) = Y_Plus_Low_Re(u_tau_new, grid % wall_dist(c1), kin_vis)
 
-          ebf = 0.001 * y_plus(c1)**4 / (1.0 + y_plus(c1))
-
           eps_int = 2.0* kin_vis * kin % n(c1)  &
                   / grid % wall_dist(c1)**2
           eps_wf  = c_mu75 * kin % n(c1)**1.5            &
@@ -148,7 +146,6 @@ subroutine Source_Eps_K_Eps_Zeta_F(flow, sol)
                      / (kappa*grid % wall_dist(c1) * p_kin(c1)),  &
                      1.0)
             eps % n(c1) = (1.0 - fa) * eps_int + fa * eps_wf
-
             ! Adjusting coefficient to fix eps value in near wall calls
             do j = a % row(c1), a % row(c1 + 1) - 1
               a % val(j) = 0.0

Sources/Process/Turbulence/Source_Kin_K_Eps.f90

@@ -112,11 +112,11 @@ subroutine Source_Kin_K_Eps(flow, sol)
                      / kin_vis
 
           tau_wall(c1) = density*kappa*u_tau(c1)*u_tan  &
-                       / log(max(((grid % wall_dist(c1)+z_o)/z_o), 1.05))
+                       / log(((grid % wall_dist(c1)+z_o) / z_o))
 
-          p_kin(c1) = density * tau_wall(c1) * c_mu25 * sqrt(kin % n(c1))   &
+          p_kin(c1) = tau_wall(c1) * c_mu25 * sqrt(kin % n(c1))   &
                     / (kappa*(grid % wall_dist(c1) + z_o))
-          b(c1)     = b(c1) + (p_kin(c1) - p_kin(c1))   &
+          b(c1)     = b(c1) + (p_kin(c1) - vis_t(c1) * shear(c1)**2)   &
                     * grid % vol(c1)
         else
           u_tau(c1) = c_mu25 * sqrt(kin % n(c1))
@@ -127,7 +127,7 @@ subroutine Source_Kin_K_Eps(flow, sol)
 
           ebf = 0.01 * y_plus(c1)**4 / (1.0 + 5.0*y_plus(c1))
 
-          p_kin_wf  = tau_wall(c1) * 0.07**0.25 * sqrt(kin % n(c1))  &
+          p_kin_wf  = tau_wall(c1) * c_mu25 * sqrt(kin % n(c1))  &
                     / (grid % wall_dist(c1) * kappa)
 
           p_kin_int = vis_t(c1) * shear(c1)**2

Sources/Process/Turbulence/Source_Kin_K_Eps_Zeta_F.f90

@@ -132,7 +132,7 @@ subroutine Source_Kin_K_Eps_Zeta_F(flow, sol)
         end if
 
         if(rough_walls) then
-          z_o = Roughness_Coefficient(grid, z_o_f(c1), c1)       
+          z_o = Roughness_Coefficient(grid, z_o_f(c1), c1)    
           u_tau(c1)  = c_mu25 * sqrt(kin % n(c1))
           y_plus(c1) = Y_Plus_Rough_Walls(u_tau(c1), &
                        grid % wall_dist(c1), kin_vis) 
@@ -153,7 +153,7 @@ subroutine Source_Kin_K_Eps_Zeta_F(flow, sol)
 
           ebf = max(0.01 * y_plus(c1)**4 / (1.0 + 5.0*y_plus(c1)),tiny)
 
-          p_kin_wf  = tau_wall(c1) * 0.07**0.25 * sqrt(kin % n(c1))  &
+          p_kin_wf  = tau_wall(c1) * c_mu25 * sqrt(kin % n(c1))  &
                     / (grid % wall_dist(c1) * kappa)
 
           p_kin_int = vis_t(c1) * shear(c1)**2

Sources/Process/Turbulence/Source_Zeta_K_Eps_Zeta_F.f90

@@ -15,10 +15,10 @@ subroutine Source_Zeta_K_Eps_Zeta_F(grid, sol, n_step)
 !------------------------------------------------------------------------------!
   implicit none
 !--------------------------------[Arguments]-----------------------------------!
-  type(Grid_Type)           :: grid
   type(Solver_Type), target :: sol
   integer                   :: n_step
 !----------------------------------[Locals]------------------------------------!
+  type(Grid_Type),   pointer :: grid
   type(Matrix_Type), pointer :: a
   real,              pointer :: b(:)
   integer                    :: c
@@ -50,7 +50,7 @@ subroutine Source_Zeta_K_Eps_Zeta_F(grid, sol, n_step)
     if(n_step > 500) then
       b(c) = b(c) + f22 % n(c) * grid % vol(c) * density
     else
-      b(c) = b(c) + max(0.0, f22 % n(c)*grid % vol(c)) * density
+      b(c) = b(c) + max(0.0, f22 % n(c) * grid % vol(c)) * density
       a % val(a % dia(c)) = a % val(a % dia(c))                  &
                           + max(0.0, -f22 % n(c) * grid % vol(c) &
                           / (zeta % n(c) + TINY)) * density

Sources/Process/Turbulence/Time_And_Length_Scale.f90

@@ -45,7 +45,7 @@ subroutine Time_And_Length_Scale(grid)
 
       t_1(c) = kin % n(c)/eps_l(c)
       t_2(c) = c_t*sqrt(kin_vis/eps_l(c))
-      t_3(c) = 0.6/(sqrt(3.0)*c_mu_d * zeta % n(c) * shear(c) + TINY)
+      t_3(c) = 0.8/(sqrt(3.0)*c_mu_d * zeta % n(c) * shear(c) + TINY)
 
       l_1(c) = kin % n(c)**1.5/eps_l(c)
       l_2(c) = c_nu * (kin_vis**3 / eps_l(c))**0.25

Sources/Process/Update_Boundary_Values.f90

@@ -113,7 +113,8 @@ subroutine Update_Boundary_Values(flow)
         end if
       end if
 
-      if(turbulence_model .eq. K_EPS_ZETA_F .and. heat_transfer) then
+      if( (turbulence_model .eq. K_EPS_ZETA_F .and. heat_transfer) .or. &
+          (turbulence_model .eq. HYBRID_LES_RANS .and. heat_transfer) ) then
         if(Grid_Mod_Bnd_Cond_Type(grid,c2) .eq. OUTFLOW  .or.   &
            Grid_Mod_Bnd_Cond_Type(grid,c2) .eq. CONVECT  .or.   &
            Grid_Mod_Bnd_Cond_Type(grid,c2) .eq. PRESSURE .or.   &
@@ -236,8 +237,8 @@ subroutine Update_Boundary_Values(flow)
           f22  % n(c2) = f22  % n(grid % bnd_cond % copy_c(c2))
         end if
 
-        if(turbulence_model .eq. K_EPS_ZETA_F .and. &
-           heat_transfer) then 
+        if( (turbulence_model .eq. K_EPS_ZETA_F .and. heat_transfer) .or. &
+            (turbulence_model .eq. HYBRID_LES_RANS .and. heat_transfer) ) then
           t2 % n(c2) = t2 % n(grid % bnd_cond % copy_c(c2))
         end if
 

Tests/Hybrid_Les_Rans/Channel_Re_Tau_2000/Stretched_Mesh/Xmgrace/k_plus_Re_tau_2000.dat

@@ -0,0 +1 @@
+../../Uniform_Mesh/Xmgrace/k_plus_Re_tau_2000.dat