One Equation RANS Model with Wall Function Support (Exawind#1299)

---------

Co-authored-by: Marc T. Henry de Frahan <[email protected]>

amr-wind/equation_systems/temperature/source_terms/CMakeLists.txt

@@ -2,4 +2,5 @@ target_sources(${amr_wind_lib_name} PRIVATE
 ABLMesoForcingTemp.cpp BodyForce.cpp
 HurricaneTempForcing.cpp
 DragTempForcing.cpp
+TempSpongeForcing.cpp
   )

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

@@ -0,0 +1,37 @@
+#ifndef TEMPSPONGEFORCING_H
+#define TEMPSPONGEFORCING_H
+
+#include "amr-wind/equation_systems/temperature/TemperatureSource.H"
+#include "amr-wind/core/SimTime.H"
+#include "amr-wind/CFDSim.H"
+
+namespace amr_wind::pde::temperature {
+
+class TempSpongeForcing : public TemperatureSource::Register<TempSpongeForcing>
+{
+public:
+    static std::string identifier() { return "TempSpongeForcing"; }
+
+    explicit TempSpongeForcing(const CFDSim& sim);
+
+    ~TempSpongeForcing() override;
+
+    void operator()(
+        const int lev,
+        const amrex::MFIter& mfi,
+        const amrex::Box& bx,
+        const FieldState /*fstate*/,
+        const amrex::Array4<amrex::Real>& src_term) const override;
+
+private:
+    const amrex::AmrCore& m_mesh;
+    const Field& m_temperature;
+    amrex::Vector<amrex::Real> m_theta_heights;
+    amrex::Vector<amrex::Real> m_theta_values;
+    amrex::Gpu::DeviceVector<amrex::Real> m_theta_heights_d;
+    amrex::Gpu::DeviceVector<amrex::Real> m_theta_values_d;
+    amrex::Real m_sponge_start{600};
+};
+
+} // namespace amr_wind::pde::temperature
+#endif

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

@@ -0,0 +1,68 @@
+
+#include "amr-wind/equation_systems/temperature/source_terms/TempSpongeForcing.H"
+#include "amr-wind/utilities/IOManager.H"
+#include "amr-wind/utilities/linear_interpolation.H"
+
+#include "AMReX_ParmParse.H"
+#include "AMReX_Gpu.H"
+#include "AMReX_Random.H"
+
+namespace amr_wind::pde::temperature {
+
+TempSpongeForcing::TempSpongeForcing(const CFDSim& sim)
+    : m_mesh(sim.mesh()), m_temperature(sim.repo().get_field("temperature"))
+{
+    amrex::ParmParse pp_abl("ABL");
+    //! Temperature variation as a function of height
+    pp_abl.query("meso_sponge_start", m_sponge_start);
+    pp_abl.getarr("temperature_heights", m_theta_heights);
+    pp_abl.getarr("temperature_values", m_theta_values);
+    AMREX_ALWAYS_ASSERT(m_theta_heights.size() == m_theta_values.size());
+    const int num_theta_values = static_cast<int>(m_theta_heights.size());
+    m_theta_heights_d.resize(num_theta_values);
+    m_theta_values_d.resize(num_theta_values);
+    amrex::Gpu::copy(
+        amrex::Gpu::hostToDevice, m_theta_heights.begin(),
+        m_theta_heights.end(), m_theta_heights_d.begin());
+    amrex::Gpu::copy(
+        amrex::Gpu::hostToDevice, m_theta_values.begin(), m_theta_values.end(),
+        m_theta_values_d.begin());
+}
+
+TempSpongeForcing::~TempSpongeForcing() = default;
+
+void TempSpongeForcing::operator()(
+    const int lev,
+    const amrex::MFIter& mfi,
+    const amrex::Box& bx,
+    const FieldState fstate,
+    const amrex::Array4<amrex::Real>& src_term) const
+{
+    const auto& geom = m_mesh.Geom(lev);
+    const auto& dx = geom.CellSizeArray();
+    const auto& prob_lo = geom.ProbLoArray();
+    const auto& prob_hi = geom.ProbHiArray();
+    const auto& temperature =
+        m_temperature.state(field_impl::dof_state(fstate))(lev).const_array(
+            mfi);
+    const amrex::Real sponge_start = m_sponge_start;
+    const auto vsize = m_theta_heights_d.size();
+    const auto* theta_heights_d = m_theta_heights_d.data();
+    const auto* theta_values_d = m_theta_values_d.data();
+    amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
+        const amrex::Real z = prob_lo[2] + (k + 0.5) * dx[2];
+        const amrex::Real zi =
+            std::max((z - sponge_start) / (prob_hi[2] - sponge_start), 0.0);
+        amrex::Real ref_temp = temperature(i, j, k);
+        if (zi > 0) {
+            ref_temp = (vsize > 0)
+                           ? interp::linear(
+                                 theta_heights_d, theta_heights_d + vsize,
+                                 theta_values_d, z)
+                           : temperature(i, j, k);
+        }
+        src_term(i, j, k, 0) -= zi * zi * (temperature(i, j, k) - ref_temp);
+    });
+}
+
+} // namespace amr_wind::pde::temperature

amr-wind/equation_systems/tke/source_terms/CMakeLists.txt

@@ -1,4 +1,5 @@
 target_sources(${amr_wind_lib_name} PRIVATE
   KsgsM84Src.cpp
   KwSSTSrc.cpp
+  KransAxell.cpp
   )

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

@@ -0,0 +1,52 @@
+#ifndef KRANSAXELL_H
+#define KRANSAXELL_H
+
+#include "amr-wind/equation_systems/tke/TKESource.H"
+
+namespace amr_wind::pde::tke {
+
+/** TKE source term based on Axell 2011 paper
+ * Axell, L. B., & Liungman, O. (2001). A one-equation turbulence model for
+ * geophysical applications: comparison with data and the k− ε model.
+ * Environmental Fluid Mechanics, 1, 71-106.
+ *  \ingroup tke_src turb_model we_abl
+ */
+class KransAxell : public TKESource::Register<KransAxell>
+{
+public:
+    static std::string identifier() { return "KransAxell"; }
+
+    explicit KransAxell(const CFDSim& /*sim*/);
+
+    ~KransAxell() override;
+
+    void operator()(
+        const int lev,
+        const amrex::MFIter& mfi,
+        const amrex::Box& bx,
+        const FieldState fstate,
+        const amrex::Array4<amrex::Real>& src_term) const override;
+
+private:
+    Field& m_turb_lscale;
+    Field& m_shear_prod;
+    Field& m_buoy_prod;
+    Field& m_dissip;
+    Field& m_tke;
+    amrex::Real m_Cmu{0.556};
+    amrex::Real m_heat_flux{0.0};
+    amrex::Real m_ref_temp{300.0};
+    amrex::Real m_z0{0.1};
+    amrex::Real m_kappa{0.41};
+    amrex::Real m_sponge_start{600};
+    amrex::Real m_ref_tke{1e-10};
+    amrex::Vector<amrex::Real> m_gravity{0.0, 0.0, -9.81};
+    const SimTime& m_time;
+    const CFDSim& m_sim;
+    const amrex::AmrCore& m_mesh;
+    const Field& m_velocity;
+};
+
+} // namespace amr_wind::pde::tke
+
+#endif /* KRANSAXELL_H */

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

@@ -0,0 +1,126 @@
+#include <AMReX_Orientation.H>
+
+#include "amr-wind/equation_systems/tke/source_terms/KransAxell.H"
+#include "amr-wind/CFDSim.H"
+#include "amr-wind/turbulence/TurbulenceModel.H"
+#include "amr-wind/utilities/linear_interpolation.H"
+namespace amr_wind::pde::tke {
+
+KransAxell::KransAxell(const CFDSim& sim)
+    : m_turb_lscale(sim.repo().get_field("turb_lscale"))
+    , m_shear_prod(sim.repo().get_field("shear_prod"))
+    , m_buoy_prod(sim.repo().get_field("buoy_prod"))
+    , m_dissip(sim.repo().get_field("dissipation"))
+    , m_tke(sim.repo().get_field("tke"))
+    , m_time(sim.time())
+    , m_sim(sim)
+    , m_mesh(sim.mesh())
+    , m_velocity(sim.repo().get_field("velocity"))
+{
+    AMREX_ALWAYS_ASSERT(sim.turbulence_model().model_name() == "KLAxell");
+    auto coeffs = sim.turbulence_model().model_coeffs();
+    amrex::ParmParse pp("ABL");
+    pp.query("Cmu", m_Cmu);
+    pp.query("kappa", m_kappa);
+    pp.query("surface_roughness_z0", m_z0);
+    pp.query("reference_temperature", m_ref_temp);
+    pp.query("surface_temp_flux", m_heat_flux);
+    pp.query("meso_sponge_start", m_sponge_start);
+    {
+        amrex::ParmParse pp_incflow("incflo");
+        pp_incflow.queryarr("gravity", m_gravity);
+    }
+}
+
+KransAxell::~KransAxell() = default;
+
+void KransAxell::operator()(
+    const int lev,
+    const amrex::MFIter& mfi,
+    const amrex::Box& bx,
+    const FieldState fstate,
+    const amrex::Array4<amrex::Real>& src_term) const
+{
+    const auto& vel =
+        m_velocity.state(field_impl::dof_state(fstate))(lev).const_array(mfi);
+    const auto& tlscale_arr = (this->m_turb_lscale)(lev).array(mfi);
+    const auto& shear_prod_arr = (this->m_shear_prod)(lev).array(mfi);
+    const auto& buoy_prod_arr = (this->m_buoy_prod)(lev).array(mfi);
+    const auto& dissip_arr = (this->m_dissip)(lev).array(mfi);
+    const auto& tke_arr = m_tke(lev).array(mfi);
+    const auto& geom = m_mesh.Geom(lev);
+    const auto& problo = m_mesh.Geom(lev).ProbLoArray();
+    const auto& probhi = m_mesh.Geom(lev).ProbHiArray();
+    const auto& dx = geom.CellSizeArray();
+    const auto& dt = m_time.delta_t();
+    const amrex::Real ref_temp = m_ref_temp;
+    const amrex::Real heat_flux =
+        std::abs(m_gravity[2]) / ref_temp * m_heat_flux;
+    const amrex::Real Cmu = m_Cmu;
+    const amrex::Real sponge_start = m_sponge_start;
+    const amrex::Real ref_tke = m_ref_tke;
+    const auto tiny = std::numeric_limits<amrex::Real>::epsilon();
+    const amrex::Real kappa = m_kappa;
+    const amrex::Real z0 = m_z0;
+    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 rans_b = std::pow(
+                std::max(heat_flux, 0.0) * kappa * z / std::pow(Cmu, 3),
+                (2.0 / 3.0));
+            bcforcing =
+                (ustar * ustar / (Cmu * Cmu) + rans_b - tke_arr(i, j, k)) / dt;
+        }
+        const amrex::Real zi =
+            std::max((z - sponge_start) / (probhi[2] - sponge_start), 0.0);
+        const amrex::Real sponge_forcing =
+            zi * zi * (tke_arr(i, j, k) - ref_tke);
+        dissip_arr(i, j, k) = std::pow(Cmu, 3) *
+                              std::pow(tke_arr(i, j, k), 1.5) /
+                              (tlscale_arr(i, j, k) + tiny);
+        src_term(i, j, k) += shear_prod_arr(i, j, k) + buoy_prod_arr(i, j, k) -
+                             dissip_arr(i, j, k) - sponge_forcing + bcforcing;
+    });
+    // Add terrain components
+    const bool has_terrain =
+        this->m_sim.repo().int_field_exists("terrain_blank");
+    if (has_terrain) {
+        const auto* const m_terrain_blank =
+            &this->m_sim.repo().get_int_field("terrain_blank");
+        const auto* const m_terrain_drag =
+            &this->m_sim.repo().get_int_field("terrain_drag");
+        const auto& blank_arr = (*m_terrain_blank)(lev).const_array(mfi);
+        const auto& drag_arr = (*m_terrain_drag)(lev).const_array(mfi);
+        amrex::ParallelFor(
+            bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
+                amrex::Real terrainforcing = 0;
+                amrex::Real dragforcing = 0;
+                const amrex::Real ux = vel(i, j, k, 0);
+                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 rans_b = std::pow(
+                    std::max(heat_flux, 0.0) * kappa * z / std::pow(Cmu, 3),
+                    (2.0 / 3.0));
+                terrainforcing =
+                    (ustar * ustar / (Cmu * Cmu) + rans_b - tke_arr(i, j, k)) /
+                    dt;
+                const amrex::Real uz = vel(i, j, k, 2);
+                m = std::sqrt(ux * ux + uy * uy + uz * uz);
+                const amrex::Real Cd =
+                    std::min(10 / (dx[2] * m + tiny), 100 / dx[2]);
+                dragforcing = -Cd * m * tke_arr(i, j, k, 0);
+
+                src_term(i, j, k) += drag_arr(i, j, k) * terrainforcing +
+                                     blank_arr(i, j, k) * dragforcing;
+            });
+    }
+}
+
+} // namespace amr_wind::pde::tke

amr-wind/turbulence/RANS/CMakeLists.txt

@@ -1,4 +1,5 @@
 target_sources(${amr_wind_lib_name} PRIVATE
   KOmegaSST.cpp
   KOmegaSSTIDDES.cpp
+  KLAxell.cpp
   )

amr-wind/turbulence/RANS/KLAxell.H

@@ -0,0 +1,69 @@
+#ifndef KLAXELL_H
+#define KLAXELL_H
+
+#include <string>
+#include "amr-wind/turbulence/TurbModelBase.H"
+
+namespace amr_wind::turbulence {
+
+/** Base class for 1-Equation RANS TKE turbulence model
+ *  \ingroup turb_model
+ */
+template <typename Transport>
+class KLAxell : public TurbModelBase<Transport>
+{
+public:
+    static std::string identifier()
+    {
+        return "KLAxell-" + Transport::identifier();
+    }
+
+    explicit KLAxell(CFDSim& sim);
+
+    std::string model_name() const override { return "KLAxell"; }
+
+    //! Update the turbulent viscosity field
+    void update_turbulent_viscosity(
+        const FieldState fstate, const DiffusionType /*unused*/) override;
+
+    //! Do any post advance work
+    void post_advance_work() override;
+
+    //! Update the effective thermal diffusivity field
+    void update_alphaeff(Field& alphaeff) override;
+
+    //! Update the effective scalar diffusivity field
+    void update_scalar_diff(Field& deff, const std::string& name) override;
+
+    //! Parse turbulence model coefficients
+    void parse_model_coeffs() override;
+
+    //! Return turbulence model coefficients
+    TurbulenceModel::CoeffsDictType model_coeffs() const override;
+
+private:
+    Field& m_vel;
+    Field& m_turb_lscale;
+    Field& m_shear_prod;
+    Field& m_buoy_prod;
+    Field& m_dissip;
+    Field& m_rho;
+    Field* m_tke{nullptr};
+    //! Turbulence constant
+    amrex::Real m_Cmu{0.556};
+    amrex::Real m_Cmu_prime{0.556};
+    amrex::Real m_Cb_stable{0.25};
+    amrex::Real m_Cb_unstable{0.35};
+    amrex::Real m_prandtl{1.0};
+    Field& m_temperature;
+    //! Gravity vector (m/s^2)
+    amrex::Vector<amrex::Real> m_gravity{0.0, 0.0, -9.81};
+    //! Reference temperature (Kelvin)
+    amrex::Real m_ref_theta{300.0};
+    amrex::Real m_surf_flux{0};
+    amrex::Real m_lengthscale_switch{800};
+};
+
+} // namespace amr_wind::turbulence
+
+#endif /* KLAXELL_H */