Skip to content

Commit c0d63f1

Browse files
hgopalanmarchdf
hgopalan
and
marchdf
authored
Improve Robustness of Immersed Forcing Wall Function and Method (#1770)
* First * Second * Fourth * Fifth * Sixth * Clang messed * Update * Bug Fix * Clean * refMOL * Wall Function --------- Co-authored-by: Marc T. Henry de Frahan <marc.henrydefrahan@nrel.gov>
1 parent 464ce08 commit c0d63f1

File tree

15 files changed

+92
-21
lines changed

15 files changed

+92
-21
lines changed

‎amr-wind/core/SimTime.H‎

Lines changed: 3 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -199,6 +199,9 @@ private:
199199
//! Max time for simulation
200200
amrex::Real m_stop_time{-1.0};
201201

202+
//! Delayed adaptive stepping
203+
amrex::Real m_delay_time{-1.0};
204+
202205
//! Maximum CFL constraint
203206
amrex::Real m_max_cfl{1.0};
204207

‎amr-wind/core/SimTime.cpp‎

Lines changed: 3 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -16,6 +16,7 @@ void SimTime::parse_parameters()
1616
amrex::ParmParse pp("time");
1717
pp.query("stop_time", m_stop_time);
1818
pp.query("max_step", m_stop_time_index);
19+
pp.query("delay_time", m_delay_time);
1920
pp.query("fixed_dt", m_fixed_dt);
2021
pp.query("initial_dt", m_initial_dt);
2122
pp.query("max_dt", m_max_dt);
@@ -248,12 +249,12 @@ void SimTime::set_current_cfl(
248249
dt_new = amrex::min(dt_new, m_initial_dt);
249250
}
250251

251-
m_dt[0] = dt_new;
252+
m_dt[0] = (m_cur_time < m_delay_time) ? m_initial_dt : dt_new;
252253

253254
} else {
254255
// Ensure that we use user-specified dt. Checkpoint restart might have
255256
// overridden this
256-
m_dt[0] = m_fixed_dt;
257+
m_dt[0] = (m_cur_time < m_delay_time) ? m_initial_dt : m_fixed_dt;
257258
}
258259

259260
m_current_cfl = 0.5 * cfl_unit_time * m_dt[0];

‎amr-wind/turbulence/LES/Kosovic.H‎

Lines changed: 6 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -58,11 +58,16 @@ private:
5858
amrex::Real m_surfaceRANSExp{2};
5959
amrex::Real m_LESTurnOff{1e15};
6060
bool m_writeTerms{false};
61-
amrex::Real m_refMOL{constants::LOW_NUM};
6261
const Field& m_vel;
6362
const Field& m_rho;
6463
Field& m_Nij;
6564
Field& m_divNij;
65+
std::string m_wall_het_model{"none"};
66+
amrex::Real m_monin_obukhov_length{constants::LARGE_NUM};
67+
amrex::Real m_kappa{0.41};
68+
amrex::Real m_gamma_m{5.0};
69+
amrex::Real m_beta_m{16.0};
70+
amrex::Real m_surface_roughness_z0{0.1};
6671
};
6772

6873
} // namespace amr_wind::turbulence

‎amr-wind/turbulence/LES/Kosovic.cpp‎

Lines changed: 69 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -7,6 +7,7 @@
77
#include "AMReX_REAL.H"
88
#include "AMReX_MultiFab.H"
99
#include "AMReX_ParmParse.H"
10+
#include "amr-wind/wind_energy/ABL.H"
1011

1112
namespace amr_wind {
1213
namespace turbulence {
@@ -25,7 +26,6 @@ Kosovic<Transport>::Kosovic(CFDSim& sim)
2526
m_Cs = std::sqrt(8 * (1 + m_Cb) / (27 * M_PI * M_PI));
2627
m_C1 = std::sqrt(960) * m_Cb / (7 * (1 + m_Cb) * m_Sk);
2728
m_C2 = m_C1;
28-
pp.query("refMOL", m_refMOL);
2929
pp.query("surfaceRANS", m_surfaceRANS);
3030
if (m_surfaceRANS) {
3131
m_surfaceFactor = 1;
@@ -40,6 +40,13 @@ Kosovic<Transport>::Kosovic(CFDSim& sim)
4040
this->m_sim.io_manager().register_io_var("divNij");
4141
}
4242
pp.query("LESOff", m_LESTurnOff);
43+
amrex::ParmParse pp_abl("ABL");
44+
pp_abl.query("wall_het_model", m_wall_het_model);
45+
pp_abl.query("monin_obukhov_length", m_monin_obukhov_length);
46+
pp_abl.query("kappa", m_kappa);
47+
pp_abl.query("mo_gamma_m", m_gamma_m);
48+
pp_abl.query("mo_beta_m", m_beta_m);
49+
pp_abl.query("surface_roughness_z0", m_surface_roughness_z0);
4350
}
4451
template <typename Transport>
4552
void Kosovic<Transport>::update_turbulent_viscosity(
@@ -59,6 +66,13 @@ void Kosovic<Transport>::update_turbulent_viscosity(
5966
const auto* m_terrain_blank =
6067
has_terrain ? &this->m_sim.repo().get_int_field("terrain_blank")
6168
: nullptr;
69+
const auto* m_terrain_drag =
70+
has_terrain ? &this->m_sim.repo().get_int_field("terrain_drag")
71+
: nullptr;
72+
const auto* m_terrain_height =
73+
has_terrain ? &this->m_sim.repo().get_field("terrain_height") : nullptr;
74+
const auto* m_terrain_z0 =
75+
has_terrain ? &this->m_sim.repo().get_field("terrainz0") : nullptr;
6276
// Populate strainrate into the turbulent viscosity arrays to avoid creating
6377
// a temporary buffer
6478
fvm::strainrate(mu_turb, vel);
@@ -75,29 +89,54 @@ void Kosovic<Transport>::update_turbulent_viscosity(
7589
const amrex::Real ds = std::cbrt(dx * dy * dz);
7690
const amrex::Real ds_sqr = ds * ds;
7791
const amrex::Real smag_factor = Cs_sqr * ds_sqr;
78-
const amrex::Real locMOL = m_refMOL;
7992
const amrex::Real locLESTurnOff = m_LESTurnOff;
8093
const amrex::Real locSwitchLoc = m_switchLoc;
8194
const amrex::Real locSurfaceRANSExp = m_surfaceRANSExp;
8295
const amrex::Real locSurfaceFactor = m_surfaceFactor;
8396
const amrex::Real locC1 = m_C1;
8497
const auto& mu_arrs = mu_turb(lev).arrays();
8598
const auto& rho_arrs = den(lev).const_arrays();
99+
const auto& vel_arrs = vel(lev).const_arrays();
86100
const auto& divNij_arrs = (this->m_divNij)(lev).arrays();
87101
const auto& blank_arrs = has_terrain
88102
? (*m_terrain_blank)(lev).const_arrays()
89103
: amrex::MultiArray4<const int>();
104+
const auto& drag_arrs = has_terrain
105+
? (*m_terrain_drag)(lev).const_arrays()
106+
: amrex::MultiArray4<const int>();
107+
const auto& height_arrs = has_terrain
108+
? (*m_terrain_height)(lev).const_arrays()
109+
: amrex::MultiArray4<const double>();
110+
const auto& z0_arrs = has_terrain ? (*m_terrain_z0)(lev).const_arrays()
111+
: amrex::MultiArray4<const double>();
112+
const amrex::Real monin_obukhov_length = m_monin_obukhov_length;
113+
const amrex::Real kappa = m_kappa;
114+
const amrex::Real surface_roughness_z0 = m_surface_roughness_z0;
115+
const amrex::Real non_neutral_neighbour =
116+
(m_wall_het_model == "mol")
117+
? MOData::calc_psi_m(
118+
1.5 * dz / monin_obukhov_length, m_beta_m, m_gamma_m)
119+
: 0.0;
90120
amrex::ParallelFor(
91121
mu_turb(lev),
92122
[=] AMREX_GPU_DEVICE(int nbx, int i, int j, int k) noexcept {
93123
const amrex::Real rho = rho_arrs[nbx](i, j, k);
94-
const amrex::Real x3 = problo[2] + (k + 0.5) * dz;
124+
amrex::Real x3 = problo[2] + (k + 0.5) * dz;
125+
x3 = (has_terrain)
126+
? std::max(x3 - height_arrs[nbx](i, j, k, 0), 0.5 * dz)
127+
: x3;
95128
const amrex::Real fmu = std::exp(-x3 / locSwitchLoc);
96129
const amrex::Real phiM =
97-
(locMOL < 0) ? std::pow(1 - 16 * x3 / locMOL, -0.25)
98-
: 1 + 5 * x3 / locMOL;
130+
(monin_obukhov_length < 0)
131+
? std::pow(1 - 16 * x3 / monin_obukhov_length, -0.25)
132+
: 1 + 5 * x3 / monin_obukhov_length;
133+
const amrex::Real wall_distance =
134+
(has_terrain)
135+
? std::max(
136+
(k + 1) * dz - height_arrs[nbx](i, j, k, 0), dz)
137+
: (k + 1) * dz;
99138
const amrex::Real ransL =
100-
std::pow(0.41 * (k + 1) * dz / phiM, 2);
139+
std::pow(0.41 * wall_distance / phiM, 2);
101140
amrex::Real turnOff = std::exp(-x3 / locLESTurnOff);
102141
amrex::Real viscosityScale =
103142
locSurfaceFactor *
@@ -108,6 +147,29 @@ void Kosovic<Transport>::update_turbulent_viscosity(
108147
(has_terrain) ? 1 - blank_arrs[nbx](i, j, k, 0) : 1.0;
109148
mu_arrs[nbx](i, j, k) *=
110149
rho * viscosityScale * turnOff * blankTerrain;
150+
// log-law
151+
const amrex::Real ux = vel_arrs[nbx](i, j, k + 1, 0);
152+
const amrex::Real uy = vel_arrs[nbx](i, j, k + 1, 1);
153+
const amrex::Real m = std::sqrt(ux * ux + uy * uy);
154+
const amrex::Real local_z0 =
155+
(has_terrain) ? std::max(z0_arrs[nbx](i, j, k, 0), 1e-4)
156+
: surface_roughness_z0;
157+
// ustar from neighbor cell above
158+
const amrex::Real ustar =
159+
m * kappa /
160+
(std::log(1.5 * dz / local_z0) - non_neutral_neighbour);
161+
const amrex::Real ux0 = vel_arrs[nbx](i, j, k, 0);
162+
const amrex::Real uy0 = vel_arrs[nbx](i, j, k, 1);
163+
const amrex::Real m0 = std::sqrt(ux0 * ux0 + uy0 * uy0);
164+
const amrex::Real uxm1 = vel_arrs[nbx](i, j, k - 1, 0);
165+
const amrex::Real uym1 = vel_arrs[nbx](i, j, k - 1, 1);
166+
const amrex::Real mm1 = std::sqrt(uxm1 * uxm1 + uym1 * uym1);
167+
const amrex::Real dMdz = std::max((m0 - mm1) / dz, 0.01);
168+
amrex::Real mut_loglaw = 2 * ustar * ustar * rho / dMdz;
169+
const amrex::Real drag =
170+
(has_terrain) ? drag_arrs[nbx](i, j, k, 0) : 0.0;
171+
mu_arrs[nbx](i, j, k) =
172+
mu_arrs[nbx](i, j, k) * (1 - drag) + drag * mut_loglaw;
111173
amrex::Real stressScale =
112174
locSurfaceFactor *
113175
(std::pow(1 - fmu, locSurfaceRANSExp) * smag_factor *
@@ -135,7 +197,7 @@ void Kosovic<Transport>::update_alphaeff(Field& alphaeff)
135197
auto lam_alpha = (this->m_transport).alpha();
136198
auto& mu_turb = this->m_mu_turb;
137199
auto& repo = mu_turb.repo();
138-
const amrex::Real muCoeff = (m_refMOL < 0) ? 3 : 1;
200+
const amrex::Real muCoeff = (m_monin_obukhov_length < 0) ? 3 : 1;
139201
const int nlevels = repo.num_active_levels();
140202
for (int lev = 0; lev < nlevels; ++lev) {
141203
const auto& muturb_arrs = mu_turb(lev).const_arrays();

‎test/test_files/abl_kosovic_neutral/abl_kosovic_neutral.inp‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -47,7 +47,7 @@ Kosovic.writeTerms = false
4747
Kosovic.switchLoc = 48
4848
# Reference Monin-Obukhov length for the RANS length scale at surface.
4949
# Can be replaced with the computed length scale in the code. However, adds code complications.
50-
Kosovic.refMOL = -1e30
50+
Kosovic.monin_obukhov_length = -1e30
5151
# Experimental feature to turn off the LES model and use the numerical dissipation as a sub-grid scale model
5252
# Initial results showed a value of 100 m improved the results by 3-5% for neutral stratification when dx<16 m
5353
Kosovic.LESOff = 100

‎test/test_files/abl_kosovic_neutral_ib/abl_kosovic_neutral_ib.inp‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -39,7 +39,7 @@ transport.thermal_expansion_coefficient = 0.003333
3939
turbulence.model = Kosovic
4040
Kosovic.surfaceRANS = false
4141
Kosovic.writeTerms = false
42-
Kosovic.refMOL = -1e30
42+
Kosovic.monin_obukhov_length = -1e30
4343
# Atmospheric boundary layer
4444
ABL.Uperiods = 12.0
4545
ABL.Vperiods = 12.0

‎test/test_files/abl_waves_terrain/abl_waves_terrain.inp‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -24,7 +24,7 @@ transport.viscosity = 1e-5
2424
transport.laminar_prandtl = 0.7
2525
transport.turbulent_prandtl = 0.333
2626
turbulence.model = Kosovic
27-
Kosovic.refMOL = -1e30
27+
Kosovic.monin_obukhov_length = -1e30
2828
incflo.physics = ABL OceanWaves TerrainDrag
2929
ICNS.source_terms = BoussinesqBuoyancy CoriolisForcing GeostrophicForcing NonLinearSGSTerm DragForcing
3030
BoussinesqBuoyancy.reference_temperature = 263.5

‎test/test_files/act_abl_joukowskydisk/act_abl_joukowskydisk.inp‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -35,7 +35,7 @@ transport.reference_temperature = 300
3535
transport.thermal_expansion_coefficient = 0.00333333
3636
# turbulence equation parameters
3737
turbulence.model = Kosovic
38-
Kosovic.refMOL = -1e30
38+
Kosovic.monin_obukhov_length = -1e30
3939
# Atmospheric boundary layer
4040
ABL.kappa = .41
4141
ABL.normal_direction = 2

‎test/test_files/act_abl_uniformctdisk/act_abl_uniformctdisk.inp‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -35,7 +35,7 @@ transport.reference_temperature = 300
3535
transport.thermal_expansion_coefficient = 0.00333333
3636
# turbulence equation parameters
3737
turbulence.model = Kosovic
38-
Kosovic.refMOL = -1e30
38+
Kosovic.monin_obukhov_length = -1e30
3939
# Atmospheric boundary layer
4040
ABL.kappa = .41
4141
ABL.normal_direction = 2

‎test/test_files/forest_drag/forest_drag.inp‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -47,7 +47,7 @@ Kosovic.writeTerms = false
4747
Kosovic.switchLoc = 48
4848
# Reference Monin-Obukhov length for the RANS length scale at surface.
4949
# Can be replaced with the computed length scale in the code. However, adds code complications.
50-
Kosovic.refMOL = -1e30
50+
Kosovic.monin_obukhov_length = -1e30
5151
# Experimental feature to turn off the LES model and use the numerical dissipation as a sub-grid scale model
5252
# Initial results showed a value of 100 m improved the results by 3-5% for neutral stratification when dx<16 m
5353
Kosovic.LESOff = 100

0 commit comments

Comments
 (0)