Donelan Heat Flux Calculation to Work in HEAT_FLUX Mode (#979)

* A fix to the Donelan wall flux model such that it works in both HEAT_FLUX and SURFACE_TEMPERATURE mode. Previously, it did not work properly in HEAT_FLUX mode. Co-authored-by: Matt Churchfield <[email protected]> Co-authored-by: gdeskos <[email protected]>
Exawind · Feb 27, 2024 · d4dd236 · d4dd236
1 parent 7291737
commit d4dd236
Show file tree

Hide file tree

Showing 5 changed files with 125 additions and 8 deletions.
diff --git a/amr-wind/wind_energy/ABLWallFunction.cpp b/amr-wind/wind_energy/ABLWallFunction.cpp
@@ -85,8 +85,8 @@ ABLWallFunction::ABLWallFunction(const CFDSim& sim)
         }
     }
 
-    m_mo.alg_type =
-        m_tempflux ? MOData::HEAT_FLUX : MOData::SURFACE_TEMPERATURE;
+    m_mo.alg_type = m_tempflux ? MOData::ThetaCalcType::HEAT_FLUX
+                               : MOData::ThetaCalcType::SURFACE_TEMPERATURE;
     m_mo.gravity = utils::vec_mag(m_gravity.data());
 }
 

diff --git a/amr-wind/wind_energy/MOData.H b/amr-wind/wind_energy/MOData.H
@@ -19,7 +19,7 @@ namespace amr_wind {
  */
 struct MOData
 {
-    enum ThetaCalcType {
+    enum class ThetaCalcType {
         HEAT_FLUX = 0,      ///< Heat-flux specified
         SURFACE_TEMPERATURE ///< Surface temperature specified
     };
@@ -46,7 +46,7 @@ struct MOData
     amrex::Real beta_m{16.0};
     amrex::Real beta_h{16.0};
 
-    ThetaCalcType alg_type{HEAT_FLUX};
+    ThetaCalcType alg_type{ThetaCalcType::HEAT_FLUX};
 
     amrex::Real phi_m() const
     {

diff --git a/amr-wind/wind_energy/MOData.cpp b/amr-wind/wind_energy/MOData.cpp
@@ -48,13 +48,13 @@ void MOData::update_fluxes(int max_iters)
     do {
         utau_iter = utau;
         switch (alg_type) {
-        case HEAT_FLUX:
+        case ThetaCalcType::HEAT_FLUX:
             surf_temp = surf_temp_flux * (std::log(zref / z0) - psi_h) /
                             (utau * kappa) +
                         theta_mean;
             break;
 
-        case SURFACE_TEMPERATURE:
+        case ThetaCalcType::SURFACE_TEMPERATURE:
             surf_temp_flux = -(theta_mean - surf_temp) * utau * kappa /
                              (std::log(zref / z0) - psi_h);
             break;

diff --git a/amr-wind/wind_energy/ShearStress.H b/amr-wind/wind_energy/ShearStress.H
@@ -14,6 +14,8 @@
  * \ingroup we_abl
  */
 
+namespace amr_wind {
+
 struct ShearStressConstant
 {
     explicit ShearStressConstant(const amr_wind::MOData& mo)
@@ -167,6 +169,8 @@ struct ShearStressDonelan
         : wspd_mean(mo.vmag_mean)
         , theta_mean(mo.theta_mean)
         , theta_surface(mo.surf_temp)
+        , temp_flux_surface(mo.surf_temp_flux)
+        , alg_type(mo.alg_type)
     {}
 
     AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real
@@ -204,12 +208,25 @@ struct ShearStressDonelan
     AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real
     calc_theta(amrex::Real /* wspd */, amrex::Real /* theta */) const
     {
-        return 0.0012 * wspd_mean * (theta_surface - theta_mean);
+        amrex::Real flux = 0.0;
+        switch (alg_type) {
+        case amr_wind::MOData::ThetaCalcType::HEAT_FLUX:
+            flux = temp_flux_surface;
+            break;
+        case amr_wind::MOData::ThetaCalcType::SURFACE_TEMPERATURE:
+            flux = 0.0012 * wspd_mean * (theta_surface - theta_mean);
+            break;
+        }
+        return flux;
     };
 
     amrex::Real wspd_mean;
     amrex::Real theta_mean;
     amrex::Real theta_surface;
+    amrex::Real temp_flux_surface;
+    amr_wind::MOData::ThetaCalcType alg_type;
 };
 
+} // namespace amr_wind
+
 #endif /* ShearStress_H */
diff --git a/unit_tests/wind_energy/test_abl_bc.cpp b/unit_tests/wind_energy/test_abl_bc.cpp
@@ -71,7 +71,7 @@ void init_velocity(amr_wind::Field& fld, amrex::Real vval, int dir)
 using ICNSFields =
     amr_wind::pde::FieldRegOp<amr_wind::pde::ICNS, amr_wind::fvm::Godunov>;
 
-TEST_F(ABLMeshTest, abl_wall_model)
+TEST_F(ABLMeshTest, abl_local_wall_model)
 {
     constexpr amrex::Real tol = 1.0e-12;
     constexpr amrex::Real mu = 0.01;
@@ -173,4 +173,104 @@ TEST_F(ABLMeshTest, abl_wall_model)
     EXPECT_NEAR(vexpct, vbase, tol);
 }
 
+TEST_F(ABLMeshTest, abl_donelan_wall_model)
+{
+    constexpr amrex::Real tol = 1.0e-12;
+    constexpr amrex::Real mu = 0.01;
+    constexpr amrex::Real vval = 35.0;
+    constexpr amrex::Real dt = 0.1;
+    constexpr amrex::Real kappa = 0.4;
+    constexpr amrex::Real z0 = 0.11;
+    int dir = 0;
+    populate_parameters();
+    {
+        amrex::ParmParse pp("geometry");
+        amrex::Vector<int> periodic{{1, 1, 0}};
+        pp.addarr("is_periodic", periodic);
+    }
+    {
+        amrex::ParmParse pp("zlo");
+        pp.add("type", (std::string) "wall_model");
+    }
+    {
+        amrex::ParmParse pp("zhi");
+        pp.add("type", (std::string) "slip_wall");
+    }
+    {
+        amrex::ParmParse pp("incflo");
+        pp.add("diffusion_type", 0);
+    }
+    {
+        amrex::ParmParse pp("transport");
+        pp.add("viscosity", mu);
+    }
+    {
+        amrex::ParmParse pp("time");
+        pp.add("fixed_dt", dt);
+    }
+    {
+        amrex::ParmParse pp("ABL");
+        pp.add("wall_shear_stress_type", (std::string) "donelan");
+        pp.add("kappa", kappa);
+        pp.add("surface_roughness_z0", z0);
+    }
+    initialize_mesh();
+
+    // Set up solver-related routines
+    auto& pde_mgr = sim().pde_manager();
+    pde_mgr.register_icns();
+    sim().create_turbulence_model();
+    sim().init_physics();
+
+    // Specify velocity as uniform in x direction
+    auto& velocity = sim().repo().get_field("velocity");
+    init_velocity(velocity, vval, dir);
+    // Specify density as unity
+    auto& density = sim().repo().get_field("density");
+    density.setVal(1.0);
+    // Perform post init for physics: turns on wall model
+    for (auto& pp : sim().physics()) {
+        pp->post_init_actions();
+    }
+
+    // Advance states to prepare for time step
+    pde_mgr.advance_states();
+
+    // Initialize icns pde
+    auto& icns_eq = pde_mgr.icns();
+    icns_eq.initialize();
+    // Initialize viscosity
+    sim().turbulence_model().update_turbulent_viscosity(
+        amr_wind::FieldState::Old, DiffusionType::Crank_Nicolson);
+    icns_eq.compute_mueff(amr_wind::FieldState::Old);
+
+    // Check test setup by verifying mu
+    const auto& viscosity = sim().repo().get_field("velocity_mueff");
+    EXPECT_NEAR(mu, utils::field_max(viscosity), tol);
+    EXPECT_NEAR(mu, utils::field_min(viscosity), tol);
+
+    // Zero source term and convection term to focus on diffusion
+    auto& src = icns_eq.fields().src_term;
+    auto& adv = icns_eq.fields().conv_term;
+    src.setVal(0.0);
+    adv.setVal(0.0);
+
+    // Calculate diffusion term
+    icns_eq.compute_diffusion_term(amr_wind::FieldState::Old);
+    // Setup mask_cell array to avoid errors in solve
+    auto& mask_cell = sim().repo().declare_int_field("mask_cell", 1, 1);
+    mask_cell.setVal(1);
+    // Compute result with just diffusion term
+    icns_eq.compute_predictor_rhs(DiffusionType::Explicit);
+
+    // Get resulting velocity in first cell
+    const amrex::Real vbase = get_val_at_kindex(velocity, dir, 0) / 8 / 8;
+
+    // Calculate expected velocity after one step
+    const amrex::Real dz = sim().mesh().Geom(0).CellSizeArray()[2];
+    const amrex::Real tau_wall = 0.0024 * std::pow(vval, 2);
+    const amrex::Real vexpct = vval + dt * (0.0 - tau_wall) / dz;
+    EXPECT_NEAR(vexpct, vbase, tol);
+}
+
 } // namespace amr_wind_tests