diff --git a/amr-wind/equation_systems/tke/source_terms/CMakeLists.txt b/amr-wind/equation_systems/tke/source_terms/CMakeLists.txt index beed97be50..507817e12d 100644 --- a/amr-wind/equation_systems/tke/source_terms/CMakeLists.txt +++ b/amr-wind/equation_systems/tke/source_terms/CMakeLists.txt @@ -1,4 +1,5 @@ target_sources(${amr_wind_lib_name} PRIVATE KsgsM84Src.cpp KwSSTSrc.cpp + Krans.cpp ) diff --git a/amr-wind/equation_systems/tke/source_terms/Krans.H b/amr-wind/equation_systems/tke/source_terms/Krans.H new file mode 100644 index 0000000000..bfbe151001 --- /dev/null +++ b/amr-wind/equation_systems/tke/source_terms/Krans.H @@ -0,0 +1,44 @@ +#ifndef Krans_H +#define Krans_H + +#include "amr-wind/equation_systems/tke/TKESource.H" + +namespace amr_wind::pde::tke { + +/** TKE source term based on Axell 2011 paper + * \ingroup tke_src turb_model we_abl + */ +class Krans : public TKESource::Register +{ +public: + static std::string identifier() { return "Krans"; } + + explicit Krans(const CFDSim& /*sim*/); + + ~Krans() override; + + void operator()( + const int lev, + const amrex::MFIter& mfi, + const amrex::Box& bx, + const FieldState fstate, + const amrex::Array4& 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_heat_flux{0.0}; + amrex::Real m_ref_temp{300.0}; + const SimTime& m_time; + const CFDSim& m_sim; + const amrex::AmrCore& m_mesh; + const Field& m_velocity; +}; + + +} // namespace amr_wind::pde::tke + +#endif /* Krans_H */ diff --git a/amr-wind/equation_systems/tke/source_terms/Krans.cpp b/amr-wind/equation_systems/tke/source_terms/Krans.cpp new file mode 100644 index 0000000000..e5349378d0 --- /dev/null +++ b/amr-wind/equation_systems/tke/source_terms/Krans.cpp @@ -0,0 +1,90 @@ +#include + +#include "amr-wind/equation_systems/tke/source_terms/Krans.H" +#include "amr-wind/CFDSim.H" +#include "amr-wind/turbulence/TurbulenceModel.H" + +namespace amr_wind::pde::tke { + +Krans::Krans(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() == "OneEqRANS"); + auto coeffs = sim.turbulence_model().model_coeffs(); + amrex::ParmParse pp("ABL"); + pp.query("reference_temperature",m_ref_temp); + pp.query("surface_temp_flux",m_heat_flux); +} + +Krans::~Krans() = default; + +void Krans::operator()( + const int lev, + const amrex::MFIter& mfi, + const amrex::Box& bx, + const FieldState fstate, + const amrex::Array4& src_term) const +{ + const auto& vel = + m_velocity.state(field_impl::dof_state(fstate))(lev).const_array(mfi); + const bool has_terrain = + this->m_sim.repo().int_field_exists("terrain_blank"); + const auto* const m_terrain_blank = has_terrain + ? &this->m_sim.repo().get_int_field("terrain_blank") + : nullptr; + const auto* const m_terrain_drag = has_terrain + ? &this->m_sim.repo().get_int_field("terrain_drag") + : nullptr; + const auto& blank_arr = has_terrain? (*m_terrain_blank)(lev).const_array(mfi) : amrex::Array4(); + const auto& drag_arr = has_terrain? (*m_terrain_drag)(lev).const_array(mfi) : amrex::Array4(); + 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& dx = geom.CellSizeArray(); + const auto& dt = m_time.delta_t(); + const amrex::Real ref_temp=m_ref_temp; + const amrex::Real heat_flux=9.81/ref_temp*m_heat_flux; + const amrex::Real Cmu=0.556; + amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + amrex::Real bcforcing=0; + amrex::Real dragforcing=0; + amrex::Real ux = vel(i, j, k, 0); + amrex::Real uy = vel(i, j, k, 1); + const int cell_drag = (has_terrain)? drag_arr(i,j,k):0; + const int cell_blank = (has_terrain)? blank_arr(i,j,k):0; + amrex::Real ustar=0.4; + if(k==0 && cell_drag==1){ + amrex::Real m = std::sqrt(ux * ux + uy * uy); + const amrex::Real kappa = 0.41; + const amrex::Real z0 = 0.1; + const amrex::Real x3 = (k==0)? problo[2] + (k + 0.5) * dx[2]:0.5*dx[2]; + ustar = m * kappa /std::log(x3/z0); + const amrex::Real rans_b=std::pow(std::max(heat_flux,0.0)*0.41*x3/std::pow(0.556,3),(2.0/3.0)); + bcforcing=(ustar*ustar/(0.556*0.556)+rans_b-tke_arr(i,j,k))/dt; + } + if(cell_blank==1){ + amrex::Real uz = vel(i, j, k, 2); + const amrex::Real m = std::sqrt(ux * ux + uy * uy + uz*uz); + const amrex::Real Cd = std::min(10/(dx[2]*m+1e-5),100/dx[2]); + dragforcing= - Cd * m * tke_arr(i,j,k,0);} + dissip_arr(i,j,k)=std::pow(Cmu,3)*std::pow(tke_arr(i,j,k), 1.5)/(tlscale_arr(i,j,k)+1e-3); + src_term(i, j, k) += shear_prod_arr(i, j, k) + buoy_prod_arr(i, j, k) - + dissip_arr(i, j, k)+bcforcing+dragforcing; + }); + + +} + +} // namespace amr_wind::pde::tke diff --git a/amr-wind/turbulence/RANS/CMakeLists.txt b/amr-wind/turbulence/RANS/CMakeLists.txt index 137db454c3..7061a582e6 100644 --- a/amr-wind/turbulence/RANS/CMakeLists.txt +++ b/amr-wind/turbulence/RANS/CMakeLists.txt @@ -1,4 +1,5 @@ target_sources(${amr_wind_lib_name} PRIVATE KOmegaSST.cpp KOmegaSSTIDDES.cpp + OneEqRANS.cpp ) diff --git a/amr-wind/turbulence/RANS/OneEqRANS.H b/amr-wind/turbulence/RANS/OneEqRANS.H new file mode 100644 index 0000000000..7fdd02cef5 --- /dev/null +++ b/amr-wind/turbulence/RANS/OneEqRANS.H @@ -0,0 +1,69 @@ +#ifndef ONEEQRANS_H +#define ONEEQRANS_H + +#include +#include "amr-wind/turbulence/TurbModelBase.H" + +namespace amr_wind::turbulence { + +/** Base class for 1-Equation RANS TKE turbulence model + * \ingroup turb_model + */ +template +class OneEqRANS : public TurbModelBase +{ +public: + + static std::string identifier() + { + return "OneEqRANS-" + Transport::identifier(); + } + + explicit OneEqRANS(CFDSim& sim); + + //~OneEqRANS() override; + + std::string model_name() const override { return "OneEqRANS"; } + + //! 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_prandtl{1.0}; + Field& m_temperature; + //! Gravity vector (m/s^2) + amrex::Vector m_gravity{0.0, 0.0, -9.81}; + //! Reference temperature (Kelvin) + amrex::Real m_ref_theta{300.0}; +}; + +} // namespace amr_wind::turbulence + +#endif /* ONEEQRANS_H */ diff --git a/amr-wind/turbulence/RANS/OneEqRANS.cpp b/amr-wind/turbulence/RANS/OneEqRANS.cpp new file mode 100644 index 0000000000..0ef1a4b0c1 --- /dev/null +++ b/amr-wind/turbulence/RANS/OneEqRANS.cpp @@ -0,0 +1,236 @@ +#include "amr-wind/turbulence/RANS/OneEqRANS.H" +#include "amr-wind/equation_systems/PDEBase.H" +#include "amr-wind/turbulence/TurbModelDefs.H" +#include "amr-wind/fvm/gradient.H" +#include "amr-wind/fvm/strainrate.H" +#include "amr-wind/turbulence/turb_utils.H" +#include "amr-wind/equation_systems/tke/TKE.H" + +#include "AMReX_ParmParse.H" + +namespace amr_wind { +namespace turbulence { + +template +OneEqRANS::OneEqRANS(CFDSim& sim) + : TurbModelBase(sim) + , m_vel(sim.repo().get_field("velocity")) + , m_turb_lscale(sim.repo().declare_field("turb_lscale", 1)) + , m_shear_prod(sim.repo().declare_field("shear_prod", 1)) + , m_buoy_prod(sim.repo().declare_field("buoy_prod", 1)) + , m_dissip(sim.repo().declare_field("dissipation", 1)) + , m_rho(sim.repo().get_field("density")) + , m_temperature(sim.repo().get_field("temperature")) +{ + auto& tke_eqn = + sim.pde_manager().register_transport_pde(pde::TKE::pde_name()); + m_tke = &(tke_eqn.fields().field); + auto& phy_mgr = this->m_sim.physics_manager(); + if (!phy_mgr.contains("ABL")) { + amrex::Abort("OneEqRANS model only works with ABL physics"); + } + { + amrex::ParmParse pp("ABL"); + pp.get("reference_temperature", m_ref_theta); + } + + { + amrex::ParmParse pp("incflo"); + pp.queryarr("gravity", m_gravity); + } + + // TKE source term to be added to PDE + turb_utils::inject_turbulence_src_terms( + pde::TKE::pde_name(), {"Krans"}); +} + + + +template +void OneEqRANS::parse_model_coeffs() +{ + const std::string coeffs_dict = this->model_name() + "_coeffs"; + amrex::ParmParse pp(coeffs_dict); + pp.query("Cmu", this->m_Cmu); + pp.query("Cmu_prime", this->m_Cmu_prime); + pp.query("prandtl",this->m_prandtl); +} + +template +TurbulenceModel::CoeffsDictType OneEqRANS::model_coeffs() const +{ + return TurbulenceModel::CoeffsDictType{ + {"Cmu", this->m_Cmu}, {"Cmu_prime", this->m_Cmu_prime},{"prandtl",this->m_prandtl}}; +} + +template +void OneEqRANS::update_turbulent_viscosity( + const FieldState fstate, const DiffusionType /*unused*/) +{ + BL_PROFILE( + "amr-wind::" + this->identifier() + "::update_turbulent_viscosity"); + + auto gradT = (this->m_sim.repo()).create_scratch_field(3, 0); + fvm::gradient(*gradT, m_temperature.state(fstate)); + + const auto& vel = this->m_vel.state(fstate); + fvm::strainrate(this->m_shear_prod, vel); + + const amrex::GpuArray gravity{ + m_gravity[0], m_gravity[1], m_gravity[2]}; + const amrex::Real beta = 1.0 / m_ref_theta; + const amrex::Real Cmu=m_Cmu; + const amrex::Real Cb_stable=0.25; + const amrex::Real Cb_unstable=0.35; + auto& mu_turb = this->mu_turb(); + const auto& den = this->m_rho.state(fstate); + const auto& repo = mu_turb.repo(); + const auto& geom_vec = repo.mesh().Geom(); + const int nlevels = repo.num_active_levels(); + const bool has_terrain_height = this->m_sim.repo().field_exists("terrain_height"); + const auto* m_terrain_height = has_terrain_height + ? &this->m_sim.repo().get_field("terrain_height") + : nullptr; + const auto* m_terrain_blank = has_terrain_height + ? &this->m_sim.repo().get_int_field("terrain_blank") + : nullptr; + for (int lev = 0; lev < nlevels; ++lev) { + const auto& geom = geom_vec[lev]; + const auto& problo = repo.mesh().Geom(lev).ProbLoArray(); + const auto& probhi = repo.mesh().Geom(lev).ProbHiArray(); + const amrex::Real dz = geom.CellSize()[2]; + for (amrex::MFIter mfi(mu_turb(lev)); mfi.isValid(); ++mfi) { + const auto& bx = mfi.tilebox(); + const auto& mu_arr = mu_turb(lev).array(mfi); + const auto& rho_arr = den(lev).const_array(mfi); + const auto& gradT_arr = (*gradT)(lev).array(mfi); + const auto& tlscale_arr = (this->m_turb_lscale)(lev).array(mfi); + const auto& tke_arr = (*this->m_tke)(lev).array(mfi); + const auto& buoy_prod_arr = (this->m_buoy_prod)(lev).array(mfi); + const auto& shear_prod_arr = (this->m_shear_prod)(lev).array(mfi); + const auto& ht_arr = has_terrain_height? (*m_terrain_height)(lev).const_array(mfi) : amrex::Array4(); + const auto& blank_arr = has_terrain_height? (*m_terrain_blank)(lev).const_array(mfi) : amrex::Array4(); + const amrex::Real Rtc=-1.0; + const amrex::Real Rtmin=-3.0; + const amrex::Real lambda=30.0; + const amrex::Real kappa=0.41; + amrex::ParallelFor( + bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + amrex::Real stratification = + -(gradT_arr(i, j, k, 0) * gravity[0] + + gradT_arr(i, j, k, 1) * gravity[1] + + gradT_arr(i, j, k, 2) * gravity[2]) * + beta; + amrex::Real x3 = problo[2] + (k + 0.5) * dz; + x3=(has_terrain_height)? std::max(x3-ht_arr(i,j,k),0.5*dz):x3; + const amrex::Real lscale_s=(lambda*kappa*x3)/(lambda+kappa*x3); + const amrex::Real lscale_b=Cb_stable*std::sqrt(tke_arr(i,j,k)/std::max(stratification,1e-10)); + amrex::Real epsilon=std::pow(Cmu,3)*std::pow(tke_arr(i,j,k), 1.5)/(tlscale_arr(i,j,k)+1e-3); + amrex::Real Rt= std::pow(tke_arr(i,j,k)/epsilon,2)*stratification; + Rt=(Rt>Rtc)? Rt:std::max(Rt,Rt - std::pow(Rt-Rtc,2)/(Rt+Rtmin-2*Rtc)); + tlscale_arr(i,j,k)= (stratification>0)? + std::sqrt(std::pow(lscale_s*lscale_b,2)/(std::pow(lscale_s,2) + std::pow(lscale_b,2))): + lscale_s*std::sqrt(1.0 - std::pow(Cmu, 6.0)*std::pow(Cb_unstable, -2.0)*Rt); + tlscale_arr(i,j,k)= (stratification>0)? std::min(tlscale_arr(i,j,k),std::sqrt(0.56*tke_arr(i,j,k)/stratification)):tlscale_arr(i,j,k); + const amrex::Real Cmu_Rt=(0.556+0.108*Rt)/(1+0.308*Rt+0.00837*std::pow(Rt,2)); + const int blank_terrain= (has_terrain_height)? + 1-blank_arr(i,j,k):1; + mu_arr(i, j, k) = (std::abs(x3-probhi[2])<=200)? 1e-10: rho_arr(i, j, k) * Cmu_Rt * + tlscale_arr(i, j, k) * + std::sqrt(tke_arr(i, j, k))*blank_terrain; + const amrex::Real Cmu_prime_Rt=0.556/(1+0.277*Rt); + const amrex::Real muPrime= (std::abs(x3-probhi[2])<=200)? 1e-10: rho_arr(i, j, k) * Cmu_prime_Rt * + tlscale_arr(i, j, k) * + std::sqrt(tke_arr(i, j, k))*blank_terrain; + buoy_prod_arr(i, j, k) = -muPrime *stratification; + shear_prod_arr(i, j, k) *= + shear_prod_arr(i, j, k) * mu_arr(i, j, k); + }); + } + } + + mu_turb.fillpatch(this->m_sim.time().current_time()); +} + +template +void OneEqRANS::update_alphaeff(Field& alphaeff) +{ + + BL_PROFILE("amr-wind::" + this->identifier() + "::update_alphaeff"); + auto lam_alpha = (this->m_transport).alpha(); + auto& mu_turb = this->m_mu_turb; + auto& repo = mu_turb.repo(); + auto gradT = (this->m_sim.repo()).create_scratch_field(3, 0); + fvm::gradient(*gradT, m_temperature); + const amrex::GpuArray gravity{ + m_gravity[0], m_gravity[1], m_gravity[2]}; + const amrex::Real beta = 1.0 / m_ref_theta; + const amrex::Real Cmu=m_Cmu; + const int nlevels = repo.num_active_levels(); + for (int lev = 0; lev < nlevels; ++lev) { + for (amrex::MFIter mfi(mu_turb(lev)); mfi.isValid(); ++mfi) { + const auto& bx = mfi.tilebox(); + const auto& muturb_arr = mu_turb(lev).array(mfi); + const auto& alphaeff_arr = alphaeff(lev).array(mfi); + const auto& lam_diff_arr = (*lam_alpha)(lev).array(mfi); + const auto& tke_arr = (*this->m_tke)(lev).array(mfi); + const auto& gradT_arr = (*gradT)(lev).array(mfi); + const auto& tlscale_arr = (this->m_turb_lscale)(lev).array(mfi); + amrex::ParallelFor( + bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + amrex::Real stratification = + -(gradT_arr(i, j, k, 0) * gravity[0] + + gradT_arr(i, j, k, 1) * gravity[1] + + gradT_arr(i, j, k, 2) * gravity[2]) * + beta; + amrex::Real epsilon=std::pow(Cmu,3)*std::pow(tke_arr(i,j,k), 1.5)/tlscale_arr(i,j,k); + amrex::Real Rt= std::pow(tke_arr(i,j,k)/epsilon,2)*stratification; + Rt=(Rt<-1)? std::max(Rt,Rt-pow(1+Rt,2)/(Rt-1)):Rt; + const amrex::Real prandtlRt=(1+0.193*Rt)/(1+0.0302*Rt); + alphaeff_arr(i, j, k) = + lam_diff_arr(i, j, k) + + muturb_arr(i, j, k)/prandtlRt; + }); + } + } + + alphaeff.fillpatch(this->m_sim.time().current_time()); +} + +template +void OneEqRANS::update_scalar_diff( + Field& deff, const std::string& name) +{ + + BL_PROFILE("amr-wind::" + this->identifier() + "::update_scalar_diff"); + + if (name == pde::TKE::var_name()) { + auto& mu_turb = this->mu_turb(); + deff.setVal(0.0); + field_ops::saxpy( + deff, 2.0, mu_turb, 0, 0, deff.num_comp(), deff.num_grow()); + } else { + amrex::Abort( + "OneEqRANS:update_scalar_diff not implemented for field " + + name); + } +} + +template +void OneEqRANS::post_advance_work() +{ + + + BL_PROFILE("amr-wind::" + this->identifier() + "::post_advance_work"); + +} + + + + + +} // namespace turbulence + +INSTANTIATE_TURBULENCE_MODEL(OneEqRANS); + +} // namespace amr_wind