Cleanup Moeng shear-stress model for clarity (#1426)

* Rename num1, num2 to match Moeng 1984 * Clarify denominator in Moeng shear-stress model
erf-model · Feb 5, 2024 · e336d57 · e336d57
1 parent f9bc46d
commit e336d57
Showing 1 changed file with 20 additions and 12 deletions.
diff --git a/Source/BoundaryConditions/MOSTStress.H b/Source/BoundaryConditions/MOSTStress.H
@@ -648,8 +648,8 @@ struct moeng_flux
         amrex::Real theta_surf  = t_surf_arr(ic,jc,zlo);
 
         amrex::Real wsp     = sqrt(velx*velx+vely*vely);
-        amrex::Real num1    = (theta-theta_mean)*wsp_mean;
-        amrex::Real num2    = (theta_mean-theta_surf)*wsp;
+        amrex::Real num1    = wsp * (theta_mean-theta_surf);
+        amrex::Real num2    = wsp_mean * (theta-theta_mean);
         amrex::Real moflux  = (std::abs(tstar) > eps) ?
                               tstar*ustar*(num1+num2)/((theta_mean-theta_surf)*wsp_mean) : 0.0;
         amrex::Real deltaz  = dz * (zlo - k);
@@ -709,17 +709,21 @@ struct moeng_flux
         amrex::Real wsp_mean = 0.5 * ( umm_arr(ic-1,jc,zlo) + umm_arr(ic,jc,zlo) );
         amrex::Real ustar    = 0.5 * ( u_star_arr(ic-1,jc,zlo) + u_star_arr(ic,jc,zlo) );
 
+        // Note: The surface mean shear stress is decomposed into tau_xz by
+        //       multiplying the modeled shear stress (rho*ustar^2) with
+        //       a factor of umean/wsp_mean for directionality; this factor
+        //       modifies the demoninator from what is in Moeng 1984.
         amrex::Real wsp     = sqrt(velx*velx+vely*vely);
-        amrex::Real num1    = (velx-umean)*wsp_mean;
-        amrex::Real num2    = umean*wsp;
-        amrex::Real stressx = ustar*ustar*(num1+num2)/(wsp_mean*wsp_mean);
+        amrex::Real num1    = wsp * umean;
+        amrex::Real num2    = wsp_mean * (velx-umean);
+        amrex::Real stressx = rho*ustar*ustar * (num1+num2)/(wsp_mean*wsp_mean);
         amrex::Real deltaz  = dz * (zlo - k);
 
         if (var_idx == Vars::xmom) {
-            dest_arr(i,j,k,icomp) = dest_arr(i,j,zlo,icomp) - stressx*rho*rho/eta*deltaz;
+            dest_arr(i,j,k,icomp) = dest_arr(i,j,zlo,icomp) - rho*stressx/eta*deltaz;
         } else {
             AMREX_ALWAYS_ASSERT(var_idx == Vars::xvel);
-            dest_arr(i,j,k,icomp) = dest_arr(i,j,zlo,icomp) - stressx*rho/eta*deltaz;
+            dest_arr(i,j,k,icomp) = dest_arr(i,j,zlo,icomp) - stressx/eta*deltaz;
         }
     }
 
@@ -773,17 +777,21 @@ struct moeng_flux
         amrex::Real wsp_mean  = 0.5 * ( umm_arr(ic,jc-1,zlo) + umm_arr(ic,jc,zlo) );
         amrex::Real ustar     = 0.5 * ( u_star_arr(ic,jc-1,zlo) + u_star_arr(ic,jc,zlo) );
 
+        // Note: The surface mean shear stress is decomposed into tau_yz by
+        //       multiplying the modeled shear stress (rho*ustar^2) with
+        //       a factor of vmean/wsp_mean for directionality; this factor
+        //       modifies the demoninator from what is in Moeng 1984.
         amrex::Real wsp     = sqrt(velx*velx+vely*vely);
-        amrex::Real num1    = (vely-vmean)*wsp_mean;
-        amrex::Real num2    = vmean*wsp;
-        amrex::Real stressy = ustar*ustar*(num1+num2)/(wsp_mean*wsp_mean);
+        amrex::Real num1    = wsp * vmean;
+        amrex::Real num2    = wsp_mean * (vely-vmean);
+        amrex::Real stressy = rho*ustar*ustar * (num1+num2)/(wsp_mean*wsp_mean);
         amrex::Real deltaz  = dz * (zlo - k);
 
         if (var_idx == Vars::ymom) {
-            dest_arr(i,j,k,icomp) = dest_arr(i,j,zlo,icomp) - stressy*rho*rho/eta*deltaz;
+            dest_arr(i,j,k,icomp) = dest_arr(i,j,zlo,icomp) - rho*stressy/eta*deltaz;
         } else {
             AMREX_ALWAYS_ASSERT(var_idx == Vars::yvel);
-            dest_arr(i,j,k,icomp) = dest_arr(i,j,zlo,icomp) - stressy*rho/eta*deltaz;
+            dest_arr(i,j,k,icomp) = dest_arr(i,j,zlo,icomp) - stressy/eta*deltaz;
         }
     }