Specified Monin-Obukhov Length

amr-wind/equation_systems/icns/source_terms/DragForcing.H

@@ -48,6 +48,30 @@ private:
     int m_sponge_south{0};
     int m_sponge_north{1};
     bool m_is_laminar{false};
+    std::string m_wall_het_model{"none"};
+    amrex::Real m_mol_length{-1e30};
+    amrex::Real m_z0{0.1};
+    amrex::Real m_kappa{0.41};
+    amrex::Real m_gamma_m{5.0};
+    amrex::Real m_beta_m{16.0};
+
+    static amrex::Real stability(
+        const amrex::Real z,
+        const amrex::Real mol_length,
+        const amrex::Real gamma_m,
+        const amrex::Real beta_m)
+    {
+        const amrex::Real zeta = z / mol_length;
+        amrex::Real psi_m = 0;
+        if (zeta > 0) {
+            psi_m = -gamma_m * zeta;
+        } else {
+            amrex::Real x = std::sqrt(std::sqrt(1 - beta_m * zeta));
+            psi_m = 2.0 * std::log(0.5 * (1.0 + x)) + log(0.5 * (1 + x * x)) -
+                    2.0 * std::atan(x) + 2.0 * std::atan(1.0);
+        }
+        return psi_m;
+    }
 };
 
 } // namespace amr_wind::pde::icns

amr-wind/equation_systems/icns/source_terms/DragForcing.cpp

@@ -39,6 +39,13 @@ DragForcing::DragForcing(const CFDSim& sim)
         amrex::Gpu::copy(
             amrex::Gpu::hostToDevice, fa_velocity.line_average().begin(),
             fa_velocity.line_average().end(), m_device_vel_vals.begin());
+        amrex::ParmParse pp_abl("ABL");
+        pp_abl.query("wall_het_model", m_wall_het_model);
+        pp_abl.query("mol_length", m_mol_length);
+        pp_abl.query("surface_roughness_z0", m_z0);
+        pp_abl.query("kappa", m_kappa);
+        pp_abl.query("mo_gamma_m", m_gamma_m);
+        pp_abl.query("mo_beta_m", m_beta_m);
     } else {
         m_sponge_strength = 0.0;
     }
@@ -96,6 +103,14 @@ void DragForcing::operator()(
     const auto tiny = std::numeric_limits<amrex::Real>::epsilon();
     const amrex::Real kappa = 0.41;
     const amrex::Real cd_max = 1000.0;
+    amrex::Real non_neutral_neighbour = 0.0;
+    amrex::Real non_neutral_cell = 0.0;
+    if (m_wall_het_model == "mol") {
+        non_neutral_neighbour =
+            stability(1.5 * dx[2], m_mol_length, m_gamma_m, m_beta_m);
+        non_neutral_cell =
+            stability(0.5 * dx[2], m_mol_length, m_gamma_m, m_beta_m);
+    }
     amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
         const amrex::Real x = prob_lo[0] + (i + 0.5) * dx[0];
         const amrex::Real y = prob_lo[1] + (j + 0.5) * dx[1];
@@ -144,9 +159,11 @@ void DragForcing::operator()(
             const amrex::Real uy2 = vel(i, j, k + 1, 1);
             const amrex::Real m2 = std::sqrt(ux2 * ux2 + uy2 * uy2);
             const amrex::Real z0 = std::max(terrainz0(i, j, k), z0_min);
-            const amrex::Real ustar = m2 * kappa / std::log(1.5 * dx[2] / z0);
+            const amrex::Real ustar =
+                m2 * kappa /
+                (std::log(1.5 * dx[2] / z0) - non_neutral_neighbour);
             const amrex::Real uTarget =
-                ustar / kappa * std::log(0.5 * dx[2] / z0);
+                ustar / kappa * (std::log(0.5 * dx[2] / z0) - non_neutral_cell);
             const amrex::Real uxTarget =
                 uTarget * ux2 / (tiny + std::sqrt(ux2 * ux2 + uy2 * uy2));
             const amrex::Real uyTarget =

amr-wind/equation_systems/temperature/source_terms/DragTempForcing.H

@@ -4,7 +4,6 @@
 #include "amr-wind/equation_systems/temperature/TemperatureSource.H"
 #include "amr-wind/core/SimTime.H"
 #include "amr-wind/CFDSim.H"
-#include "amr-wind/transport_models/TransportModel.H"
 
 namespace amr_wind::pde::temperature {
 
@@ -25,17 +24,56 @@ public:
         const amrex::Array4<amrex::Real>& src_term) const override;
 
 private:
+    const SimTime& m_time;
     const CFDSim& m_sim;
     const amrex::AmrCore& m_mesh;
     const Field& m_velocity;
     const Field& m_temperature;
+    amrex::Real m_internal_temperature{300};
     amrex::Real m_drag_coefficient{1.0};
+    std::string m_wall_het_model{"none"};
+    amrex::Real m_mol_length{-1e30};
+    amrex::Real m_z0{0.1};
+    amrex::Real m_kappa{0.41};
+    amrex::Real m_gamma_m{5.0};
+    amrex::Real m_beta_m{16.0};
+    amrex::Real m_gamma_h{5.0};
+    amrex::Real m_beta_h{16.0};
 
-    //! Transport model
-    const transport::TransportModel& m_transport;
+    static amrex::Real stability(
+        const amrex::Real z,
+        const amrex::Real mol_length,
+        const amrex::Real gamma_m,
+        const amrex::Real beta_m)
+    {
+        const amrex::Real zeta = z / mol_length;
+        amrex::Real psi_m = 0;
+        if (zeta > 0) {
+            psi_m = -gamma_m * zeta;
+        } else {
+            amrex::Real x = std::sqrt(std::sqrt(1 - beta_m * zeta));
+            psi_m = 2.0 * std::log(0.5 * (1.0 + x)) + log(0.5 * (1 + x * x)) -
+                    2.0 * std::atan(x) + 2.0 * std::atan(1.0);
+        }
+        return psi_m;
+    }
 
-    //! Reference temperature
-    std::unique_ptr<ScratchField> m_ref_theta;
+    static amrex::Real thermal_stability(
+        const amrex::Real z,
+        const amrex::Real mol_length,
+        const amrex::Real gamma_h,
+        const amrex::Real beta_h)
+    {
+        const amrex::Real zeta = z / mol_length;
+        amrex::Real psi_h = 0;
+        if (zeta > 0) {
+            psi_h = -gamma_h * zeta;
+        } else {
+            amrex::Real x = std::sqrt(1 - beta_h * zeta);
+            psi_h = 2.0 * std::log(0.5 * (1 + x));
+        }
+        return psi_h;
+    }
 };
 
 } // namespace amr_wind::pde::temperature

amr-wind/equation_systems/temperature/source_terms/DragTempForcing.cpp

@@ -9,20 +9,24 @@
 namespace amr_wind::pde::temperature {
 
 DragTempForcing::DragTempForcing(const CFDSim& sim)
-    : m_sim(sim)
+    : m_time(sim.time())
+    , m_sim(sim)
     , m_mesh(sim.mesh())
     , m_velocity(sim.repo().get_field("velocity"))
     , m_temperature(sim.repo().get_field("temperature"))
-    , m_transport(sim.transport_model())
 {
     amrex::ParmParse pp("DragTempForcing");
     pp.query("drag_coefficient", m_drag_coefficient);
-    m_ref_theta = m_transport.ref_theta();
-    if (pp.contains("reference_temperature")) {
-        amrex::Abort(
-            "DragTempForcing.reference_temperature has been deprecated. Please "
-            "replace with transport.reference_temperature.");
-    }
+    pp.query("internal_temperature", m_internal_temperature);
+    amrex::ParmParse pp_abl("ABL");
+    pp_abl.query("wall_het_model", m_wall_het_model);
+    pp_abl.query("mol_length", m_mol_length);
+    pp_abl.query("surface_roughness_z0", m_z0);
+    pp_abl.query("kappa", m_kappa);
+    pp_abl.query("mo_gamma_m", m_gamma_m);
+    pp_abl.query("mo_beta_m", m_beta_m);
+    pp_abl.query("mo_gamma_m", m_gamma_h);
+    pp_abl.query("mo_beta_m", m_beta_h);
 }
 
 DragTempForcing::~DragTempForcing() = default;
@@ -47,22 +51,55 @@ void DragTempForcing::operator()(
     auto* const m_terrain_blank =
         &this->m_sim.repo().get_int_field("terrain_blank");
     const auto& blank = (*m_terrain_blank)(lev).const_array(mfi);
+    auto* const m_terrain_drag =
+        &this->m_sim.repo().get_int_field("terrain_drag");
+    const auto& drag = (*m_terrain_drag)(lev).const_array(mfi);
     const auto& geom = m_mesh.Geom(lev);
     const auto& dx = geom.CellSizeArray();
     const amrex::Real drag_coefficient = m_drag_coefficient / dx[2];
-    const auto& ref_theta = (*m_ref_theta)(lev).const_array(mfi);
+    const amrex::Real internal_temperature = m_internal_temperature;
+    const amrex::Real gravity_mod = 9.81;
+    amrex::Real psi_m = 0.0;
+    amrex::Real psi_h_neighbour = 0.0;
+    amrex::Real psi_h_cell = 0.0;
+    const amrex::Real kappa = m_kappa;
+    const amrex::Real z0 = m_z0;
+    const amrex::Real mol_length = m_mol_length;
+    const auto& dt = m_time.delta_t();
+    if (m_wall_het_model == "mol") {
+        psi_m = stability(1.5 * dx[2], m_mol_length, m_gamma_m, m_beta_m);
+        psi_h_neighbour =
+            thermal_stability(1.5 * dx[2], m_mol_length, m_gamma_h, m_beta_h);
+        psi_h_cell =
+            thermal_stability(0.5 * dx[2], m_mol_length, m_gamma_h, m_beta_h);
+    }
     const auto tiny = std::numeric_limits<amrex::Real>::epsilon();
     const amrex::Real cd_max = 10.0;
     amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
         const amrex::Real ux1 = vel(i, j, k, 0);
         const amrex::Real uy1 = vel(i, j, k, 1);
         const amrex::Real uz1 = vel(i, j, k, 2);
+        const amrex::Real theta = temperature(i, j, k, 0);
+        const amrex::Real theta2 = temperature(i, j, k + 1, 0);
+        const amrex::Real wspd = std::sqrt(ux1 * ux1 + uy1 * uy1);
+        const amrex::Real ustar =
+            wspd * kappa / (std::log(1.5 * dx[2] / z0) - psi_m);
+        //! We do not know the actual temperature so use cell above
+        const amrex::Real thetastar =
+            theta * ustar * ustar / (kappa * gravity_mod * mol_length);
+        const amrex::Real surf_temp =
+            theta2 -
+            thetastar / kappa * (std::log(1.5 * dx[2] / z0) - psi_h_neighbour);
+        const amrex::Real tTarget =
+            surf_temp +
+            thetastar / kappa * (std::log(0.5 * dx[2] / z0) - psi_h_cell);
+        const amrex::Real bc_forcing_t = -(tTarget - theta) / dt;
         const amrex::Real m = std::sqrt(ux1 * ux1 + uy1 * uy1 + uz1 * uz1);
         const amrex::Real Cd =
             std::min(drag_coefficient / (m + tiny), cd_max / dx[2]);
-        const amrex::Real T0 = ref_theta(i, j, k);
         src_term(i, j, k, 0) -=
-            (Cd * (temperature(i, j, k, 0) - T0) * blank(i, j, k, 0));
+            (Cd * (theta - internal_temperature) * blank(i, j, k, 0) +
+             bc_forcing_t * drag(i, j, k));
     });
 }
 

amr-wind/equation_systems/tke/source_terms/KransAxell.H

@@ -51,6 +51,28 @@ private:
 
     //! Reference temperature
     std::unique_ptr<ScratchField> m_ref_theta;
+    std::string m_wall_het_model{"none"};
+    amrex::Real m_mol_length{10000};
+    amrex::Real m_gamma_m{5.0};
+    amrex::Real m_beta_m{16.0};
+
+    static amrex::Real stability(
+        const amrex::Real z,
+        const amrex::Real mol_length,
+        const amrex::Real gamma_m,
+        const amrex::Real beta_m)
+    {
+        const amrex::Real zeta = z / mol_length;
+        amrex::Real psi_m = 0;
+        if (zeta > 0) {
+            psi_m = -gamma_m * zeta;
+        } else {
+            amrex::Real x = std::sqrt(std::sqrt(1 - beta_m * zeta));
+            psi_m = 2.0 * std::log(0.5 * (1.0 + x)) + log(0.5 * (1 + x * x)) -
+                    2.0 * std::atan(x) + 2.0 * std::atan(1.0);
+        }
+        return psi_m;
+    }
 };
 
 } // namespace amr_wind::pde::tke

amr-wind/equation_systems/tke/source_terms/KransAxell.cpp

@@ -26,6 +26,10 @@ KransAxell::KransAxell(const CFDSim& sim)
     pp.query("surface_roughness_z0", m_z0);
     pp.query("surface_temp_flux", m_heat_flux);
     pp.query("meso_sponge_start", m_sponge_start);
+    pp.query("wall_het_model", m_wall_het_model);
+    pp.query("mol_length", m_mol_length);
+    pp.query("mo_gamma_m", m_gamma_m);
+    pp.query("mo_beta_m", m_beta_m);
     {
         amrex::ParmParse pp_incflow("incflo");
         pp_incflow.queryarr("gravity", m_gravity);
@@ -64,14 +68,18 @@ void KransAxell::operator()(
     const amrex::Real z0 = m_z0;
     const amrex::GpuArray<amrex::Real, AMREX_SPACEDIM> gravity{
         m_gravity[0], m_gravity[1], m_gravity[2]};
+    amrex::Real psi_m = 0.0;
+    if (m_wall_het_model == "mol") {
+        psi_m = stability(0.5 * dx[2], m_mol_length, m_gamma_m, m_beta_m);
+    }
     amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
         amrex::Real bcforcing = 0;
         const amrex::Real ux = vel(i, j, k, 0);
         const amrex::Real uy = vel(i, j, k, 1);
         const amrex::Real z = problo[2] + (k + 0.5) * dx[2];
         if (k == 0) {
             const amrex::Real m = std::sqrt(ux * ux + uy * uy);
-            const amrex::Real ustar = m * kappa / std::log(z / z0);
+            const amrex::Real ustar = m * kappa / (std::log(z / z0) - psi_m);
             const amrex::Real T0 = ref_theta_arr(i, j, k);
             const amrex::Real hf = std::abs(gravity[2]) / T0 * heat_flux;
             const amrex::Real rans_b = std::pow(
@@ -107,7 +115,8 @@ void KransAxell::operator()(
                 const amrex::Real uy = vel(i, j, k, 1);
                 const amrex::Real z = 0.5 * dx[2];
                 amrex::Real m = std::sqrt(ux * ux + uy * uy);
-                const amrex::Real ustar = m * kappa / std::log(z / z0);
+                const amrex::Real ustar =
+                    m * kappa / (std::log(z / z0) - psi_m);
                 const amrex::Real T0 = ref_theta_arr(i, j, k);
                 const amrex::Real hf = std::abs(gravity[2]) / T0 * heat_flux;
                 const amrex::Real rans_b = std::pow(

amr-wind/wind_energy/ABLWallFunction.H

@@ -82,6 +82,12 @@ public:
 private:
     const ABLWallFunction& m_wall_func;
     std::string m_wall_shear_stress_type{"moeng"};
+    std::string m_wall_het_model{"none"};
+    amrex::Real m_mol_length{10000};
+    amrex::Real m_z0{0.1};
+    amrex::Real m_kappa{0.41};
+    amrex::Real m_gamma_m{5.0};
+    amrex::Real m_beta_m{16.0};
 };
 
 class ABLTempWallFunc : public FieldBCIface
@@ -98,6 +104,14 @@ public:
 private:
     const ABLWallFunction& m_wall_func;
     std::string m_wall_shear_stress_type{"moeng"};
+    std::string m_wall_het_model{"none"};
+    amrex::Real m_mol_length{10000};
+    amrex::Real m_z0{0.1};
+    amrex::Real m_kappa{0.41};
+    amrex::Real m_gamma_m{5.0};
+    amrex::Real m_gamma_h{5.0};
+    amrex::Real m_beta_m{16.0};
+    amrex::Real m_beta_h{16.0};
 };
 
 } // namespace amr_wind