diff --git a/amr-wind/equation_systems/icns/source_terms/ABLForcing.H b/amr-wind/equation_systems/icns/source_terms/ABLForcing.H index 53a9f0fbe6..739983e64b 100644 --- a/amr-wind/equation_systems/icns/source_terms/ABLForcing.H +++ b/amr-wind/equation_systems/icns/source_terms/ABLForcing.H @@ -5,6 +5,7 @@ #include "amr-wind/core/SimTime.H" #include "amr-wind/utilities/trig_ops.H" #include "amr-wind/utilities/linear_interpolation.H" +#include "amr-wind/utilities/FieldPlaneAveraging.H" namespace amr_wind::pde::icns { @@ -29,12 +30,18 @@ public: const FieldState fstate, const amrex::Array4& src_term) const override; + void mean_velocity_init(const VelPlaneAveraging& /*vavg*/); + + void mean_velocity_update(const VelPlaneAveraging& /*vavg*/); + inline void set_target_velocities(amrex::Real ux, amrex::Real uy) { m_target_vel[0] = ux; m_target_vel[1] = uy; } + inline void set_boundary_layer_height(amrex::Real zi) { m_zi = zi; } + inline void set_mean_velocities(amrex::Real ux, amrex::Real uy) { m_mean_vel[0] = ux; @@ -71,7 +78,9 @@ public: outfile.open(m_force_timetable, std::ios::out | std::ios::app); outfile << std::setprecision(17) << nph_time << "\t" << m_abl_forcing[0] << "\t" << m_abl_forcing[1] << "\t" - << 0.0 << std::endl; + << 0.0 << "\t" << m_abl_forcing[1] / m_coriolis_factor + << "\t" << -m_abl_forcing[0] / m_coriolis_factor << "\t" + << m_zi << std::endl; } } @@ -83,9 +92,36 @@ private: const SimTime& m_time; const amrex::AmrCore& m_mesh; + amrex::Gpu::DeviceVector m_vel_ht; + amrex::Gpu::DeviceVector m_vel_vals; + + //! Axis over which averages are computed. + int m_axis{2}; + //! Activated when water is present in domain bool m_use_phase_ramp{false}; + //! Activate free-atmosphere damping + bool m_fa_damping{false}; + + //! Detect free-atmosphere base height (set to zi from ABLStats). + bool m_fa_detect_height{false}; + + //! Free-atmosphere base height. + amrex::Real m_fa_height{750.0}; + + //! Start time for free-atmosphere damping. + amrex::Real m_fa_time_start{0.0}; + + //! End time for free-atmosphere damping. + amrex::Real m_fa_time_end{1.0E9}; + + //! Time scale for free-atmosphere damping. + amrex::Real m_fa_tau{100.0}; + + //! Coriolis factor used for free-atmosphere damping. + amrex::Real m_coriolis_factor; + //! Number of cells in band to prevent forcing near liquid int m_n_band{2}; @@ -97,10 +133,13 @@ private: //! Bool for writing forcing time table bool m_write_force_timetable{false}; + //! File name for forcing time table output std::string m_force_timetable; + //! Output frequency for forces int m_force_outfreq{1}; + //! Output start time for force amrex::Real m_force_outstart{0.0}; @@ -130,6 +169,9 @@ private: //! Local storage of interface location amrex::Real m_water_level; + //! boundary layer height. + amrex::Real m_zi; + //! VOF field, to avoid forcing on liquid above force-off height const Field* m_vof; }; diff --git a/amr-wind/equation_systems/icns/source_terms/ABLForcing.cpp b/amr-wind/equation_systems/icns/source_terms/ABLForcing.cpp index 79ad1fbbd6..bdc00fba77 100644 --- a/amr-wind/equation_systems/icns/source_terms/ABLForcing.cpp +++ b/amr-wind/equation_systems/icns/source_terms/ABLForcing.cpp @@ -47,16 +47,43 @@ ABLForcing::ABLForcing(const CFDSim& sim) m_write_force_timetable = pp_abl.contains("forcing_timetable_output_file"); if (m_write_force_timetable) { + amrex::Print() << "Here!!!!" << std::endl; pp_abl.get("forcing_timetable_output_file", m_force_timetable); pp_abl.query("forcing_timetable_frequency", m_force_outfreq); pp_abl.query("forcing_timetable_start_time", m_force_outstart); if (amrex::ParallelDescriptor::IOProcessor()) { + amrex::Print() << m_force_timetable << std::endl; std::ofstream outfile; outfile.open(m_force_timetable, std::ios::out); - outfile << "time\tfx\tfy\tfz\n"; + outfile << "time\tfx\tfy\tfz\tUgx\tUgy\tzi\n"; } } + // If free-atmosphere damping is used, read these inputs. + pp_abl.query("free_atmosphere_damping", m_fa_damping); + if (m_fa_damping) { + pp_abl.query("detect_free_atmosphere_height", m_fa_detect_height); + if (!m_fa_detect_height) { + pp_abl.get("free_atmosphere_height", m_fa_height); + } + pp_abl.query("free_atmosphere_damping_start_time", m_fa_time_start); + pp_abl.query("free_atmosphere_damping_end_time", m_fa_time_end); + pp_abl.query("free_atmosphere_damping_time_scale", m_fa_tau); + + amrex::ParmParse pp_coriolis("CoriolisForcing"); + amrex::Real rot_time_period = 86164.091; + pp_coriolis.query("rotational_time_period", rot_time_period); + + amrex::Real omega = utils::two_pi() / rot_time_period; + + amrex::Real latitude = 90.0; + pp_coriolis.query("latitude", latitude); + latitude = utils::radians(latitude); + amrex::Real sinphi = std::sin(latitude); + + m_coriolis_factor = 2.0 * omega * sinphi; + } + for (int i = 0; i < AMREX_SPACEDIM; ++i) { m_mean_vel[i] = m_target_vel[i]; } @@ -80,6 +107,8 @@ ABLForcing::ABLForcing(const CFDSim& sim) // Point to something, will not be used m_vof = &sim.repo().get_field("velocity"); } + + mean_velocity_init(abl.abl_statistics().vel_profile_coarse()); } ABLForcing::~ABLForcing() = default; @@ -95,6 +124,19 @@ void ABLForcing::operator()( const amrex::Real dudt = m_abl_forcing[0]; const amrex::Real dvdt = m_abl_forcing[1]; + const bool fa_damping = m_fa_damping; + const bool fa_detect_height = m_fa_detect_height; + const amrex::Real fa_height = fa_detect_height ? m_zi : m_fa_height; + const amrex::Real fa_time_start = m_fa_time_start; + const amrex::Real fa_time_end = m_fa_time_end; + const amrex::Real fa_tau = m_fa_tau; + const amrex::Real fa_u = dvdt / m_coriolis_factor; + const amrex::Real fa_v = -dudt / m_coriolis_factor; + + const auto& current_time = m_time.current_time(); + const auto& new_time = m_time.new_time(); + const auto& nph_time = 0.5 * (current_time + new_time); + const bool ph_ramp = m_use_phase_ramp; const int n_band = m_n_band; const amrex::Real wlev = m_water_level; @@ -103,11 +145,26 @@ void ABLForcing::operator()( const auto& problo = m_mesh.Geom(lev).ProbLoArray(); const auto& dx = m_mesh.Geom(lev).CellSizeArray(); + const int idir = m_axis; + const amrex::Real* heights = m_vel_ht.data(); + const amrex::Real* heights_end = m_vel_ht.end(); + const amrex::Real* vals = m_vel_vals.data(); + const auto& vof = (*m_vof)(lev).const_array(mfi); + amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + amrex::IntVect iv(i, j, k); + const amrex::Real z = problo[idir] + (iv[idir] + 0.5) * dx[idir]; + + const amrex::Real umean = + amr_wind::interp::linear(heights, heights_end, vals, z, 3, 0); + const amrex::Real vmean = + amr_wind::interp::linear(heights, heights_end, vals, z, 3, 1); + + // This part deals with air and water phase and only masks ABL forcing + // on the water phase. amrex::Real fac = 1.0; if (ph_ramp) { - const amrex::Real z = problo[2] + (k + 0.5) * dx[2]; if (z - wlev < wrht0 + wrht1) { if (z - wlev < wrht0) { // Apply no forcing within first interval @@ -127,11 +184,56 @@ void ABLForcing::operator()( fac = 0.0; } } - src_term(i, j, k, 0) += fac * dudt; - src_term(i, j, k, 1) += fac * dvdt; - // No forcing in z-direction + // This part applies the free atmosphere damping, which nudges + // the solution in the free atmosphere toward geostrophic as + // computed by the wind speed controller. (i.e., the wind speed + // controller computes a forcing term that is really a driving + // pressure gradient, so there is a corresponding geostrophic + // wind. Nudge the free atmosphere toward that geostrophic wind.) + amrex::Real src_fa_damping_x = 0.0; + amrex::Real src_fa_damping_y = 0.0; + if (fa_damping && + ((nph_time >= fa_time_start) && (nph_time <= fa_time_end)) && + (z >= fa_height)) { + src_fa_damping_x = (1.0 / fa_tau) * (fa_u - umean); + src_fa_damping_y = (1.0 / fa_tau) * (fa_v - vmean); + } + + // Sum up the source term + src_term(i, j, k, 0) += (fac * dudt) + src_fa_damping_x; + src_term(i, j, k, 1) += (fac * dvdt) + src_fa_damping_y; + src_term(i, j, k, 2) += 0.0; }); } +void ABLForcing::mean_velocity_init(const VelPlaneAveraging& vavg) +{ + m_axis = vavg.axis(); + + // The implementation depends the assumption that the ABL statistics class + // computes statistics at the cell-centeres only on level 0. If this + // assumption changes in future, the implementation will break... so put in + // a check here to catch this. + AMREX_ALWAYS_ASSERT( + m_mesh.Geom(0).Domain().length(m_axis) == + static_cast(vavg.line_centroids().size())); + + m_vel_ht.resize(vavg.line_centroids().size()); + m_vel_vals.resize(vavg.line_average().size()); + + amrex::Gpu::copy( + amrex::Gpu::hostToDevice, vavg.line_centroids().begin(), + vavg.line_centroids().end(), m_vel_ht.begin()); + + mean_velocity_update(vavg); +} + +void ABLForcing::mean_velocity_update(const VelPlaneAveraging& vavg) +{ + amrex::Gpu::copy( + amrex::Gpu::hostToDevice, vavg.line_average().begin(), + vavg.line_average().end(), m_vel_vals.begin()); +} + } // namespace amr_wind::pde::icns diff --git a/amr-wind/wind_energy/ABL.cpp b/amr-wind/wind_energy/ABL.cpp index 8ca419de57..7d1aa3b751 100644 --- a/amr-wind/wind_energy/ABL.cpp +++ b/amr-wind/wind_energy/ABL.cpp @@ -173,6 +173,7 @@ void ABL::post_regrid_actions() { m_abl_anelastic->post_regrid_actions(); } */ void ABL::pre_advance_work() { + const auto& zi = m_stats->zi(); const auto& vel_pa = m_stats->vel_profile(); m_abl_wall_func.update_umean( m_stats->vel_profile(), m_stats->theta_profile_fine()); @@ -201,6 +202,9 @@ void ABL::pre_advance_work() #endif m_abl_forcing->set_mean_velocities(vx, vy); + m_abl_forcing->set_boundary_layer_height(zi); + + m_abl_forcing->mean_velocity_update(m_stats->vel_profile_coarse()); } if (m_abl_mean_bous != nullptr) { diff --git a/amr-wind/wind_energy/ABLStats.H b/amr-wind/wind_energy/ABLStats.H index 1ba8269e3e..3bcc6490b5 100644 --- a/amr-wind/wind_energy/ABLStats.H +++ b/amr-wind/wind_energy/ABLStats.H @@ -58,6 +58,8 @@ public: //! Compute height of capping inversion void compute_zi(); + amrex::Real zi() const override { return m_zi; }; + //! Return vel plane averaging instance const VelPlaneAveragingFine& vel_profile() const override { diff --git a/amr-wind/wind_energy/ABLStats.cpp b/amr-wind/wind_energy/ABLStats.cpp index e855cdd78b..b19b0fbde7 100644 --- a/amr-wind/wind_energy/ABLStats.cpp +++ b/amr-wind/wind_energy/ABLStats.cpp @@ -230,13 +230,29 @@ void ABLStats::post_advance_work() void ABLStats::compute_zi() { - + // This finds the location of the capping inversion (z_i) using the + // method of Sullivan et al. in "Structure of the entrainment zone + // capping the convective atmospheric boundary layer," JAS, Vol. 55, + // 1998. + // + // In this method, for every x,y location in the computational domain, + // the maximum d(\theta)/dz is found. You then have an x,y array + // of x,y local z_i. z_i is then averaged over the x,y array to give + // . + + // Compute the potential temperature gradient. auto gradT = (this->m_sim.repo()) .create_scratch_field(3, m_temperature.num_grow()[0]); fvm::gradient(*gradT, m_temperature); // Only compute zi using coarsest level const int lev = 0; + + // Using the AMReX ReduceToPlane function, find the maximum + // d(\theta)/dz over the mesh along lines in the normal direction, + // which typically would be z (if z is up). We end up with this + // 2D "planar" array of max d(\theta)/dz over the horizontal + // extent of the domain. const int dir = m_normal_dir; const auto& geom = (this->m_sim.repo()).mesh().Geom(lev); auto const& domain_box = geom.Domain(); @@ -260,11 +276,21 @@ void ABLStats::compute_zi() auto& pinned_tg_fab = device_tg_fab; #endif + // Because of the parallel decomposition, there may be multiple + // of these 2D "planar" arrays local to different processes handling + // boxes of mesh at different heights. So, we need to parallel + // reduce these by picking the max value at each x,y location to get + // the final resultant planar array of max temperature gradient values. + // This final array is collected only to the IO processor. amrex::ParallelReduce::Max( pinned_tg_fab.dataPtr(), static_cast(pinned_tg_fab.size()), amrex::ParallelDescriptor::IOProcessorNumber(), amrex::ParallelDescriptor::Communicator()); + // To get the planar average over the non-normal directions (typically + // x and y), we sum up all the max temperature gradient values. Again, + // this only happens on the IO processor. Once we have the sum, divide + // by the number of coarse grid points in x and y. if (amrex::ParallelDescriptor::IOProcessor()) { const auto dnval = m_dn; auto* p = pinned_tg_fab.dataPtr(); @@ -276,6 +302,11 @@ void ABLStats::compute_zi() 0.0); m_zi /= static_cast(pinned_tg_fab.size()); } + + // It is necessary for other processors to know , so broadcast + // it back out to all processors. + amrex::ParallelDescriptor::Bcast( + &m_zi, 1, amrex::ParallelDescriptor::IOProcessorNumber()); } void ABLStats::process_output() diff --git a/amr-wind/wind_energy/ABLStatsBase.H b/amr-wind/wind_energy/ABLStatsBase.H index 7a10f5ab48..2fb08cd3ab 100644 --- a/amr-wind/wind_energy/ABLStatsBase.H +++ b/amr-wind/wind_energy/ABLStatsBase.H @@ -32,6 +32,8 @@ public: //! Flag indicating ABL simulation mode virtual ABLStatsMode abl_mode() const = 0; + virtual amrex::Real zi() const = 0; + //! Interpolating object for vertical velocity profile virtual const VelPlaneAveraging& vel_profile_coarse() const = 0; virtual const VelPlaneAveragingFine& vel_profile() const = 0; diff --git a/amr-wind/wind_energy/ABLWallFunction.H b/amr-wind/wind_energy/ABLWallFunction.H index 9d8e4b1e1c..79c5a0b163 100644 --- a/amr-wind/wind_energy/ABLWallFunction.H +++ b/amr-wind/wind_energy/ABLWallFunction.H @@ -53,13 +53,32 @@ private: amrex::Vector m_gravity{0.0, 0.0, -9.81}; + //! Ability to read in a table of surface temperature flux versus time. + std::string m_surf_temp_flux_timetable; + amrex::Vector m_surf_temp_flux_time; + amrex::Vector m_surf_temp_flux_value; + //! Ability to read in a table of surface temperature versus time. std::string m_surf_temp_timetable; amrex::Vector m_surf_temp_time; amrex::Vector m_surf_temp_value; - bool m_tempflux{true}; - bool m_temp_table{false}; + //! Ability to read in a table of near-surface temperature verus time. + std::string m_nearsurf_temp_timetable; + amrex::Vector m_nearsurf_temp_time; + amrex::Vector m_nearsurf_temp_value; + + bool m_surf_temp_flux{true}; + bool m_surf_temp{false}; + bool m_nearsurf_temp{false}; + + bool m_surf_temp_flux_use_table{false}; + + bool m_surf_temp_use_table{false}; + bool m_surf_temp_use_rate{false}; + + bool m_nearsurf_temp_use_table{false}; + amrex::Real m_surf_temp_rate{0.0}; amrex::Real m_surf_temp_rate_tstart{0.0}; amrex::Real m_surf_temp_init{300.0}; diff --git a/amr-wind/wind_energy/ABLWallFunction.cpp b/amr-wind/wind_energy/ABLWallFunction.cpp index 4438573751..70a8736508 100644 --- a/amr-wind/wind_energy/ABLWallFunction.cpp +++ b/amr-wind/wind_energy/ABLWallFunction.cpp @@ -26,6 +26,7 @@ ABLWallFunction::ABLWallFunction(const CFDSim& sim) amrex::ParmParse pp("ABL"); pp.query("kappa", m_mo.kappa); + pp.query("mo_alpha_h", m_mo.alpha_h); pp.query("mo_gamma_m", m_mo.gamma_m); pp.query("mo_gamma_h", m_mo.gamma_h); pp.query("mo_beta_m", m_mo.beta_m); @@ -33,6 +34,7 @@ ABLWallFunction::ABLWallFunction(const CFDSim& sim) const char* z0_same = "surface_roughness_z0"; const char* z0_aero = "aerodynamic_roughness_length"; const char* z0_therm = "thermal_roughness_length"; + pp.query(z0_same, m_mo.z0); if (!pp.contains(z0_same)) { pp.query(z0_aero, m_mo.z0); @@ -60,81 +62,194 @@ ABLWallFunction::ABLWallFunction(const CFDSim& sim) << std::endl; } - if (pp.contains("surface_temp_flux")) { - pp.query("surface_temp_flux", m_mo.surf_temp_flux); + // Check for type of surface heating and read in data. + // - Specified surface temperature flux mode: + if ((pp.contains("surface_temp_flux")) || + (pp.contains("surface_temp_flux_timetable"))) { + m_surf_temp_flux = true; + m_surf_temp = false; + m_nearsurf_temp = false; + + // - Fixed, time-invariant surface temperature flux mode: + if (pp.contains("surface_temp_flux")) { + m_surf_temp_flux_use_table = false; + + pp.query("surface_temp_flux", m_mo.surf_temp_flux); + + amrex::Print() + << "ABLWallFunction: Fixed, time-invariant surface temperature " + "flux mode is selected." + << std::endl; + } + // - surface temperature flux from time-table mode: + else if (pp.contains("surface_temp_flux_timetable")) { + m_surf_temp_flux_use_table = true; + + pp.query("surface_temp_flux_timetable", m_surf_temp_flux_timetable); + + amrex::Print() + << "ABLWallFunction: Surface temperature flux time table " + "mode is selected." + << std::endl; + if (!m_surf_temp_flux_timetable.empty()) { + std::ifstream ifh(m_surf_temp_flux_timetable, std::ios::in); + if (!ifh.good()) { + amrex::Abort( + "Cannot find surface_temp_flux_timetable file: " + + m_surf_temp_flux_timetable); + } + amrex::Real data_time, data_value; + while (ifh >> data_time >> data_value) { + m_surf_temp_flux_time.push_back(data_time); + m_surf_temp_flux_value.push_back(data_value); + ifh.ignore( + std::numeric_limits::max(), '\n'); + } + } + } + + // - Specified near-surface temperature mode: + } else if (pp.contains("near_surface_temp_timetable")) { + m_surf_temp_flux = false; + m_surf_temp = false; + m_nearsurf_temp = true; + + m_nearsurf_temp_use_table = true; + + if (pp.contains("near_surface_height")) { + pp.get("near_surface_height", m_mo.zNearSurf); + } else { + amrex::Abort( + "ABLWallFunction: near_surface_height must be specified when " + "specifying near-surface temperature."); + } + + pp.query("near_surface_temp_timetable", m_nearsurf_temp_timetable); + amrex::Print() - << "ABLWallFunction: Surface temperature flux mode is selected." + << "ABLWallFunction: Near-surface temperature time table " + "mode is selected." << std::endl; - } else if (pp.contains("surface_temp_timetable")) { - pp.query("surface_temp_timetable", m_surf_temp_timetable); - m_tempflux = false; - m_temp_table = true; - amrex::Print() << "ABLWallFunction: Surface temperature time table " - "mode is selected." - << std::endl; - if (!m_surf_temp_timetable.empty()) { - std::ifstream ifh(m_surf_temp_timetable, std::ios::in); + + if (!m_nearsurf_temp_timetable.empty()) { + std::ifstream ifh(m_nearsurf_temp_timetable, std::ios::in); if (!ifh.good()) { amrex::Abort( - "Cannot find surface_temp_timetable file: " + - m_surf_temp_timetable); + "Cannot find near_surface_temp_timetable file: " + + m_nearsurf_temp_timetable); } - amrex::Real data_time; - amrex::Real data_value; - ifh.ignore(std::numeric_limits::max(), '\n'); - while (ifh >> data_time) { - ifh >> data_value; - m_surf_temp_time.push_back(data_time); - m_surf_temp_value.push_back(data_value); + amrex::Real data_time, data_value; + while (ifh >> data_time >> data_value) { + m_nearsurf_temp_time.push_back(data_time); + m_nearsurf_temp_value.push_back(data_value); + ifh.ignore(std::numeric_limits::max(), '\n'); } } - } else if (pp.contains("surface_temp_rate")) { - m_tempflux = false; - pp.get("surface_temp_rate", m_surf_temp_rate); - amrex::Print() - << "ABLWallFunction: Surface temperature rate mode is selected." - << std::endl; - if (pp.contains("surface_temp_init")) { - pp.get("surface_temp_init", m_surf_temp_init); - } else { + + // - Specified surface temperature mode: + } else if ( + (pp.contains("surface_temp_timetable")) || + (pp.contains("surface_temp_rate"))) { + m_surf_temp_flux = false; + m_surf_temp = true; + m_nearsurf_temp = false; + + if (pp.contains("surface_temp_timetable")) { + m_surf_temp_use_table = true; + m_surf_temp_use_rate = false; + + pp.query("surface_temp_timetable", m_surf_temp_timetable); + amrex::Print() << "ABLWallFunction: Surface temperature time table " + "mode is selected." + << std::endl; + if (!m_surf_temp_timetable.empty()) { + std::ifstream ifh(m_surf_temp_timetable, std::ios::in); + if (!ifh.good()) { + amrex::Abort( + "Cannot find surface_temp_timetable file: " + + m_surf_temp_timetable); + } + amrex::Real data_time, data_value; + while (ifh >> data_time >> data_value) { + m_surf_temp_time.push_back(data_time); + m_surf_temp_value.push_back(data_value); + ifh.ignore( + std::numeric_limits::max(), '\n'); + } + } + // - Specified surface temperature rate mode: + } else if (pp.contains("surface_temp_rate")) { + m_surf_temp_use_table = false; + m_surf_temp_use_rate = true; + + pp.get("surface_temp_rate", m_surf_temp_rate); + amrex::Print() - << "ABLWallFunction: Initial surface temperature not found for " - "ABL. Assuming to be equal to the reference temperature" + << "ABLWallFunction: Surface temperature rate mode is selected." << std::endl; - m_surf_temp_init = sim.transport_model().reference_temperature(); - } - if (pp.contains("surface_temp_rate_tstart")) { - pp.get("surface_temp_rate_tstart", m_surf_temp_rate_tstart); - } else { - amrex::Print() - << "ABLWallFunction: Surface temperature heating/cooling start " - "time (surface_temp_rate_tstart) not found for ABL. " - "Assuming zero." - << m_surf_temp_rate_tstart << std::endl; + + if (pp.contains("surface_temp_init")) { + pp.get("surface_temp_init", m_surf_temp_init); + } else { + amrex::Print() + << "ABLWallFunction: Initial surface temperature not found " + "for " + "ABL. Assuming to be equal to the reference temperature" + << std::endl; + m_surf_temp_init = + sim.transport_model().reference_temperature(); + } + + if (pp.contains("surface_temp_rate_tstart")) { + pp.get("surface_temp_rate_tstart", m_surf_temp_rate_tstart); + } else { + amrex::Print() + << "ABLWallFunction: Surface temperature heating/cooling " + "start " + "time (surface_temp_rate_tstart) not found for ABL. " + "Assuming zero." + << m_surf_temp_rate_tstart << std::endl; + } } + // - If no surface heating mode is specified, default to no surface + // heating. } else { - amrex::Print() << "ABLWallFunction: Neither surface_temp_flux nor " - "surface_temp_rate specified for ABL physics. " - "Assuming Neutral Stratification" + amrex::Print() << "ABLWallFunction: None of 'surface_temp_flux,' " + "'surface_temp_rate', 'surface_temp_timetable,' " + "or 'near_surface_temp_timetable' are specified " + "for ABL physics. Assuming no heat flux at surface." << std::endl; } + // Check to see if the simulation is non-periodic. If so, recommend + // using the wall model locally instead of planar averaging. if (pp.contains("inflow_outflow_mode")) { pp.query("inflow_outflow_mode", m_inflow_outflow); if (m_inflow_outflow) { pp.query("wf_velocity", m_wf_vel); pp.query("wf_vmag", m_wf_vmag); pp.query("wf_theta", m_wf_theta); - amrex::Print() << "ABLWallFunction: Inflow/Outflow mode is turned " - "on. Please make sure wall shear stress type is " - "set to local." - << std::endl; + amrex::Print() + << "ABLWallFunction: Inflow/Outflow mode is turned " + "on. It is recommended that the wall shear stress " + "type is set to local IF the flow is heterogeneous." + << std::endl; } } - m_mo.alg_type = m_tempflux ? MOData::ThetaCalcType::HEAT_FLUX - : MOData::ThetaCalcType::SURFACE_TEMPERATURE; + // Set the Monin-Obukhov algorithm type. + if (m_surf_temp_flux) { + m_mo.alg_type = MOData::ThetaCalcType::HEAT_FLUX; + } else if (m_surf_temp) { + m_mo.alg_type = MOData::ThetaCalcType::SURFACE_TEMPERATURE; + } else if (m_nearsurf_temp) { + m_mo.alg_type = MOData::ThetaCalcType::NEAR_SURFACE_TEMPERATURE; + } else { + m_mo.alg_type = MOData::ThetaCalcType::HEAT_FLUX; + } + + // Set gravity strength in Monin-Obukhov. m_mo.gravity = utils::vec_mag(m_gravity.data()); } @@ -151,20 +266,35 @@ void ABLWallFunction::update_umean( { const auto& time = m_sim.time(); - if (!m_tempflux) { - if (!m_temp_table) { + // If not using a constant heat flux... + if (m_surf_temp_flux && m_surf_temp_flux_use_table) { + m_mo.surf_temp_flux = amr_wind::interp::linear( + m_surf_temp_flux_time, m_surf_temp_flux_value, time.current_time()); + amrex::Print() << "Current surface temperature flux: " + << m_mo.surf_temp_flux << std::endl; + } + // - if specifying surface temperature + if (m_surf_temp) { + if (m_surf_temp_use_rate) { m_mo.surf_temp = m_surf_temp_init + m_surf_temp_rate * amrex::max( time.current_time() - m_surf_temp_rate_tstart, 0.0) / 3600.0; - } else { + } else if (m_surf_temp_use_table) { m_mo.surf_temp = amr_wind::interp::linear( m_surf_temp_time, m_surf_temp_value, time.current_time()); } amrex::Print() << "Current surface temperature: " << m_mo.surf_temp << std::endl; + // - if specifying near-surface temperature + } else if (m_nearsurf_temp) { + if (m_nearsurf_temp_use_table) { + m_mo.near_surf_temp = amr_wind::interp::linear( + m_nearsurf_temp_time, m_nearsurf_temp_value, + time.current_time()); + } } if (m_inflow_outflow) { diff --git a/amr-wind/wind_energy/MOData.H b/amr-wind/wind_energy/MOData.H index 101a1c627a..b2a6342cd7 100644 --- a/amr-wind/wind_energy/MOData.H +++ b/amr-wind/wind_energy/MOData.H @@ -20,16 +20,19 @@ namespace amr_wind { struct MOData { enum class ThetaCalcType { - HEAT_FLUX = 0, ///< Heat-flux specified - SURFACE_TEMPERATURE ///< Surface temperature specified + HEAT_FLUX = 0, ///< Heat-flux specified + SURFACE_TEMPERATURE = 1, ///< Surface temperature specified + NEAR_SURFACE_TEMPERATURE = 2, ///< Temperature specified above but + // near the surface }; - amrex::Real zref{0.0}; ///< Reference height (m) - amrex::Real z0{0.1}; ///< Aerodynamic roughness height (m) - amrex::Real z0t{z0}; ///< Thermal roughness height (m) - amrex::Real utau; ///< Friction velocity (m/s) - amrex::Real kappa{0.41}; ///< von Karman constant - amrex::Real gravity{9.81}; ///< Acceleration due to gravity (m/s^2) + amrex::Real zref{0.0}; ///< Reference height (m) + amrex::Real zNearSurf{2.0}; ///< Near-surface temperature driving height (m) + amrex::Real z0{0.1}; ///< Aerodynamic roughness height (m) + amrex::Real z0t{z0}; ///< Thermal roughness height (m) + amrex::Real utau; ///< Friction velocity (m/s) + amrex::Real kappa{0.41}; ///< von Karman constant + amrex::Real gravity{9.81}; ///< Acceleration due to gravity (m/s^2) amrex::Real obukhov_len{1.0e16}; ///< Non-dimensional Obukhov length amrex::RealArray vel_mean; ///< Mean velocity (at zref) @@ -40,7 +43,10 @@ struct MOData amrex::Real surf_temp_flux{0.0}; ///< Heat flux amrex::Real surf_temp; ///< Instantaneous surface temperature + amrex::Real near_surf_temp; ///< Instantaneous near-surface temperature (at + ///< zNearSurf). + amrex::Real alpha_h{1.0}; amrex::Real gamma_m{5.0}; amrex::Real gamma_h{5.0}; amrex::Real beta_m{16.0}; diff --git a/amr-wind/wind_energy/MOData.cpp b/amr-wind/wind_energy/MOData.cpp index 75241cbc1f..bc31d46e3f 100644 --- a/amr-wind/wind_energy/MOData.cpp +++ b/amr-wind/wind_energy/MOData.cpp @@ -52,8 +52,8 @@ void MOData::update_fluxes(int max_iters) amrex::Real utau_iter = 0.0; // Initialize variables - amrex::Real psi_m = 0.0; - amrex::Real psi_h = 0.0; + amrex::Real psi_m_zref = 0.0; + amrex::Real psi_h_zref = 0.0; utau = kappa * vmag_mean / (std::log(zref / z0)); int iter = 0; @@ -61,14 +61,39 @@ void MOData::update_fluxes(int max_iters) utau_iter = utau; switch (alg_type) { case ThetaCalcType::HEAT_FLUX: - surf_temp = surf_temp_flux * (std::log(zref / z0t) - psi_h) / + surf_temp = alpha_h * surf_temp_flux * + (std::log(zref / z0t) - psi_h_zref) / (utau * kappa) + theta_mean; + near_surf_temp = alpha_h * surf_temp_flux * + (std::log(2.0 / z0t) - psi_h_zref) / + (utau * kappa) + + theta_mean; break; case ThetaCalcType::SURFACE_TEMPERATURE: surf_temp_flux = -(theta_mean - surf_temp) * utau * kappa / - (std::log(zref / z0t) - psi_h); + (alpha_h * (std::log(zref / z0t) - psi_h_zref)); + near_surf_temp = alpha_h * surf_temp_flux * + (std::log(2.0 / z0t) - psi_h_zref) / + (utau * kappa) + + theta_mean; + break; + + case ThetaCalcType::NEAR_SURFACE_TEMPERATURE: + amrex::Real psi_h_zNearSurf = calc_psi_h(zNearSurf / obukhov_len); + amrex::Real a11 = 1.0; + amrex::Real a12 = -(alpha_h / (kappa * utau)) * + (std::log(zref / z0t) - psi_h_zref); + amrex::Real a21 = 1.0; + amrex::Real a22 = -(alpha_h / (kappa * utau)) * + (std::log(zNearSurf / z0t) - psi_h_zNearSurf); + + amrex::Real coeff = 1.0 / (a11 * a22 - a12 * a21); + + surf_temp = coeff * (a22 * theta_mean - a12 * near_surf_temp); + surf_temp_flux = coeff * (-a21 * theta_mean + a11 * near_surf_temp); + break; } @@ -82,18 +107,27 @@ void MOData::update_fluxes(int max_iters) obukhov_len = std::numeric_limits::max(); zeta = 0.0; } - psi_m = calc_psi_m(zeta); - psi_h = calc_psi_h(zeta); - utau = kappa * vmag_mean / (std::log(zref / z0) - psi_m); + psi_m_zref = calc_psi_m(zeta); + psi_h_zref = calc_psi_h(zeta); + utau = kappa * vmag_mean / (std::log(zref / z0) - psi_m_zref); ++iter; + + amrex::Print() << "iteration: " << iter << ", T(" << zNearSurf + << " m) = " << near_surf_temp << " K, T(" << zref + << " m) = " << theta_mean << "K , L = " << obukhov_len + << " m, uStar = " << utau + << "m/s, Qs = " << surf_temp_flux + << " K-m/s, theta_0 = " << surf_temp << " K" + << std::endl; } while ((std::abs(utau_iter - utau) > 1e-5) && iter <= max_iters); if (iter >= max_iters) { amrex::Print() << "MOData::update_fluxes: Convergence criteria not met after " << max_iters << " iterations\nObuhov length = " << obukhov_len - << " zeta = " << zeta << "\npsi_m = " << psi_m - << " psi_h = " << psi_h << "\nutau = " << utau + << " zeta = (z_ref/L) = " << zeta + << "\npsi_m(z_ref/L) = " << psi_m_zref + << " psi_h(z_ref/L) = " << psi_h_zref << "\nutau = " << utau << " Tsurf = " << surf_temp << " q = " << surf_temp_flux << std::endl; } diff --git a/amr-wind/wind_energy/ShearStress.H b/amr-wind/wind_energy/ShearStress.H index 34d20f97cb..76aff3e26c 100644 --- a/amr-wind/wind_energy/ShearStress.H +++ b/amr-wind/wind_energy/ShearStress.H @@ -216,6 +216,9 @@ struct ShearStressDonelan case amr_wind::MOData::ThetaCalcType::SURFACE_TEMPERATURE: flux = 0.0012 * wspd_mean * (theta_surface - theta_mean); break; + case amr_wind::MOData::ThetaCalcType::NEAR_SURFACE_TEMPERATURE: + amrex::Print() << "To be implemented..." << std::endl; + break; } return flux; };