Kosovic Model - Multiphase Support

amr-wind/turbulence/LES/CMakeLists.txt

@@ -3,5 +3,6 @@ target_sources(${amr_wind_lib_name} PRIVATE
   OneEqKsgs.cpp
   AMD.cpp
   AMDNoTherm.cpp
-  Kosovic.cpp		
+  Kosovic.cpp
+  MultiPhaseKosovic.cpp
   )

amr-wind/turbulence/LES/MultiPhaseKosovic.H

@@ -0,0 +1,64 @@
+#ifndef MULTIPHASEKOSOVIC_H
+#define MULTIPHASEKOSOVIC_H
+
+#include <string>
+#include "amr-wind/turbulence/TurbModelBase.H"
+#include "amr-wind/core/FieldRepo.H"
+
+namespace amr_wind::turbulence {
+
+/** MultiPhaseKosovic LES Model
+ *  \ingroup turb_model
+ */
+template <typename Transport>
+class MultiPhaseKosovic : public TurbModelBase<Transport>
+{
+public:
+    static std::string identifier()
+    {
+        return "MultiPhaseKosovic-" + Transport::identifier();
+    }
+
+    explicit MultiPhaseKosovic(CFDSim& sim);
+
+    //! Model name for debugging purposes
+    std::string model_name() const override { return "MultiPhaseKosovic"; }
+
+    //! Update the turbulent viscosity field
+    void update_turbulent_viscosity(
+        const FieldState fstate, const DiffusionType /*unused*/) override;
+
+    //! No post advance work for this model
+    void post_advance_work() override {}
+
+    //! Update the effective thermal diffusivity field
+    void update_alphaeff(Field& alphaeff) override;
+
+    //! Parse turbulence model coefficients
+    void parse_model_coeffs() override;
+
+    //! Return model coefficients dictionary
+    TurbulenceModel::CoeffsDictType model_coeffs() const override;
+
+private:
+    //! MultiPhaseKosovic coefficient (default value set for ABL simulations)
+    // Ref: Mirocha et. al "Implementation of a Nonlinear Subfilter Turbulence
+    //  Stress Model for Large-Eddy Simulation in the Advanced Research WRF
+    //  Model"
+    // , MWR 2012.
+    amrex::Real m_Cb{0.36};
+    amrex::Real m_Cs{0.135};
+    amrex::Real m_C1{2.1};
+    amrex::Real m_C2{2.1};
+    amrex::Real m_Sk{0.5};
+    bool m_writeTerms{false};
+    amrex::Vector<amrex::Real> m_gravity{0.0, 0.0, -9.81};
+    const Field& m_vel;
+    const Field& m_rho;
+    Field& m_Nij;
+    Field& m_divNij;
+};
+
+} // namespace amr_wind::turbulence
+
+#endif /* MULTIPHASEKOSOVIC_H */

amr-wind/turbulence/LES/MultiPhaseKosovic.cpp

@@ -0,0 +1,158 @@
+#include <cmath>
+#include "amr-wind/turbulence/LES/MultiPhaseKosovic.H"
+#include "amr-wind/turbulence/TurbModelDefs.H"
+#include "amr-wind/fvm/nonLinearSum.H"
+#include "amr-wind/fvm/strainrate.H"
+#include "amr-wind/fvm/divergence.H"
+#include "amr-wind/fvm/gradient.H"
+#include "AMReX_REAL.H"
+#include "AMReX_MultiFab.H"
+#include "AMReX_ParmParse.H"
+
+namespace amr_wind {
+namespace turbulence {
+
+template <typename Transport>
+// cppcheck-suppress uninitMemberVar
+MultiPhaseKosovic<Transport>::MultiPhaseKosovic(CFDSim& sim)
+    : TurbModelBase<Transport>(sim)
+    , m_vel(sim.repo().get_field("velocity"))
+    , m_rho(sim.repo().get_field("density"))
+    , m_Nij(sim.repo().declare_field("Nij", 9, 1, 1))
+    , m_divNij(sim.repo().declare_field("divNij", 3))
+{
+    amrex::ParmParse pp("MultiPhaseKosovic");
+    pp.query("Cb", m_Cb);
+    m_Cs = std::sqrt(8 * (1 + m_Cb) / (27 * M_PI * M_PI));
+    m_C1 = std::sqrt(960) * m_Cb / (7 * (1 + m_Cb) * m_Sk);
+    m_C2 = m_C1;
+    pp.query("writeTerms", m_writeTerms);
+    if (m_writeTerms) {
+        this->m_sim.io_manager().register_io_var("Nij");
+        this->m_sim.io_manager().register_io_var("divNij");
+    }
+    this->m_sim.io_manager().register_io_var("density_prod");
+    {
+        amrex::ParmParse pp_incflow("incflo");
+        pp_incflow.queryarr("gravity", m_gravity);
+    }
+    if (!sim.repo().field_exists("vof")) {
+        amrex::Abort(
+            " Use the single-phase Kosovic model. This model is only for "
+            "multiphase");
+    }
+}
+template <typename Transport>
+void MultiPhaseKosovic<Transport>::update_turbulent_viscosity(
+    const FieldState fstate, const DiffusionType /*unused*/)
+{
+    BL_PROFILE(
+        "amr-wind::" + this->identifier() + "::update_turbulent_viscosity");
+
+    auto& mu_turb = this->mu_turb();
+    const auto& repo = mu_turb.repo();
+    const auto& vel = m_vel.state(fstate);
+    const auto& den = m_rho.state(fstate);
+    const auto& geom_vec = repo.mesh().Geom();
+    const amrex::Real Cs_sqr = this->m_Cs * this->m_Cs;
+    const auto* m_vof = &this->m_sim.repo().get_field("vof");
+    // Populate strainrate into the turbulent viscosity arrays to avoid creating
+    // a temporary buffer
+    fvm::strainrate(mu_turb, vel);
+    // Non-linear component Nij is computed here and goes into Body Forcing
+    fvm::nonlinearsum(m_Nij, vel);
+    fvm::divergence(m_divNij, m_Nij);
+    const int nlevels = repo.num_active_levels();
+    for (int lev = 0; lev < nlevels; ++lev) {
+        const auto& geom = geom_vec[lev];
+        const amrex::Real dx = geom.CellSize()[0];
+        const amrex::Real dy = geom.CellSize()[1];
+        const amrex::Real dz = geom.CellSize()[2];
+        const amrex::Real ds = std::cbrt(dx * dy * dz);
+        const amrex::Real ds_sqr = ds * ds;
+        const amrex::Real smag_factor = Cs_sqr * ds_sqr;
+        const amrex::Real locC1 = m_C1;
+        const auto& mu_arrs = mu_turb(lev).arrays();
+        const auto& rho_arrs = den(lev).const_arrays();
+        const auto& divNij_arrs = (this->m_divNij)(lev).arrays();
+        const auto& vof_arrs = (*m_vof)(lev).const_arrays();
+        amrex::ParallelFor(
+            mu_turb(lev),
+            [=] AMREX_GPU_DEVICE(int nbx, int i, int j, int k) noexcept {
+                const amrex::Real rho = rho_arrs[nbx](i, j, k);
+                //! Turn off turbulence in water as interface not handled with
+                //! temperature
+                const amrex::Real water_diff =
+                    (vof_arrs[nbx](i, j, k) == 0) ? 1 : 0;
+                mu_arrs[nbx](i, j, k) *= rho * smag_factor * water_diff;
+                //! Non-linear model is not tested in water and interface
+                //! effects are approximate Turning off non-linear terms in
+                //! water based on density
+                amrex::Real water_mixing_length =
+                    (vof_arrs[nbx](i, j, k) == 0) ? 1 : 0;
+                divNij_arrs[nbx](i, j, k, 0) *=
+                    rho * smag_factor * 0.25 * locC1 * water_mixing_length;
+                divNij_arrs[nbx](i, j, k, 1) *=
+                    rho * smag_factor * 0.25 * locC1 * water_mixing_length;
+                divNij_arrs[nbx](i, j, k, 2) *=
+                    rho * smag_factor * 0.25 * locC1 * water_mixing_length;
+            });
+    }
+    amrex::Gpu::streamSynchronize();
+
+    mu_turb.fillpatch(this->m_sim.time().current_time());
+}
+
+template <typename Transport>
+void MultiPhaseKosovic<Transport>::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();
+    const amrex::Real muCoeff = 1.0;
+    const auto* m_vof = &this->m_sim.repo().get_field("vof");
+    const int nlevels = repo.num_active_levels();
+    for (int lev = 0; lev < nlevels; ++lev) {
+        const auto& muturb_arrs = mu_turb(lev).const_arrays();
+        const auto& alphaeff_arrs = alphaeff(lev).arrays();
+        const auto& lam_diff_arrs = (*lam_alpha)(lev).const_arrays();
+        const auto& vof_arrs = (*m_vof)(lev).const_arrays();
+        amrex::ParallelFor(
+            mu_turb(lev),
+            [=] AMREX_GPU_DEVICE(int nbx, int i, int j, int k) noexcept {
+                //! Current air-water interface setup causes oscillation in
+                //! temperature with turbulent diffusion. Turning it off inside
+                //! water until a better way to handle it can be found.
+                const amrex::Real water_diff =
+                    (vof_arrs[nbx](i, j, k) == 0) ? 1 : 0;
+                alphaeff_arrs[nbx](i, j, k) =
+                    lam_diff_arrs[nbx](i, j, k) +
+                    muCoeff * muturb_arrs[nbx](i, j, k) * water_diff;
+            });
+    }
+    amrex::Gpu::streamSynchronize();
+
+    alphaeff.fillpatch(this->m_sim.time().current_time());
+}
+template <typename Transport>
+void MultiPhaseKosovic<Transport>::parse_model_coeffs()
+{
+    const std::string coeffs_dict = this->model_name() + "_coeffs";
+    amrex::ParmParse pp(coeffs_dict);
+    pp.query("Cs", this->m_Cs);
+}
+
+template <typename Transport>
+TurbulenceModel::CoeffsDictType
+MultiPhaseKosovic<Transport>::model_coeffs() const
+{
+    return TurbulenceModel::CoeffsDictType{{"Cb", this->m_Cb}};
+}
+
+} // namespace turbulence
+
+INSTANTIATE_TURBULENCE_MODEL(MultiPhaseKosovic);
+
+} // namespace amr_wind

docs/sphinx/theory/include/turbulence.rst

@@ -163,4 +163,7 @@ is the atmospheric stability function. Currently, the implementation for the sta
 The implementation of the non-linear model is split into two parts. The subgrid-scale viscosity term is directly used 
 within the AMR-Wind diffusion framework. The last two terms in :math:`M_{ij}` are added as source-terms in the momentum equation. 
 
+- **Multiphase Variant**
 
+The non-linear model can also be used for multiphase flows. However, the user has to change the model name in input file to 
+``turbulence.model = MultiPhaseKosovic`` instead of ``turbulence.model= Kosovic``. 

docs/sphinx/theory/theory.rst

@@ -20,3 +20,4 @@ Theory Manual
 .. include:: include/forest_model.rst
 
 .. include:: include/turbine_models.rst
+

test/CMakeLists.txt

@@ -207,6 +207,7 @@ add_test_re(abl_unstable_schumann_wall_model)
 add_test_re(abl_anelastic)
 add_test_re(abl_multiphase_laminar)
 add_test_re(abl_multiphase_laminar_2planes)
+add_test_re(abl_multiphase_kosovic)
 add_test_re(abl_waves_terrain)
 add_test_re(act_abl_joukowskydisk)
 add_test_re(act_abl_uniformctdisk)

test/test_files/abl_multiphase_kosovic/abl_multiphase_kosovic.inp

@@ -0,0 +1,126 @@
+#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨#
+#            SIMULATION CONTROL         #
+#.......................................#
+time.stop_time                           = 7200
+time.max_step                            = -1          # Max number of time steps
+time.fixed_dt                            = -1          # Use this constant dt if > 0
+time.initial_dt                          = 0.05
+time.cfl                                 = 0.9         # CFL factor
+
+time.plot_interval                       = 1000       # Steps between plot files
+time.checkpoint_interval                 = 25000       # Steps between checkpoint files
+io.output_default_variables              = false
+io.outputs                               = velocity terrain_blank  ow_velocity
+incflo.use_godunov                       = 1
+incflo.diffusion_type                    = 1
+incflo.godunov_type                      = weno_z
+incflo.mflux_type                        = upwind
+turbulence.model                         = MultiPhaseKosovic
+incflo.gravity                           = 0.  0. -9.81  # Gravitational force (3D)
+transport.model                          = TwoPhaseTransport
+transport.viscosity_fluid1               = 1e-3
+transport.viscosity_fluid2               = 1e-5
+transport.laminar_prandtl_fluid1         = 7.2
+transport.laminar_prandtl_fluid2         = 0.7
+transport.turbulent_prandtl              = 0.3333
+MultiPhase.density_fluid1                = 1020.
+MultiPhase.density_fluid2                = 1.225
+MultiPhase.water_level                   = 0.
+incflo.verbose                           =   2          # incflo_level
+
+#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨#
+#            GEOMETRY & BCs             #
+#.......................................#
+geometry.prob_lo                         = 0.       0.    -200.  # Lo corner coordinates
+geometry.prob_hi                         = 3000 3000 1600  # Hi corner coordinates; 600 / 128 = 4.6875 m resolution
+amr.n_cell                               = 384 384 904
+amr.max_level                            = 0
+geometry.is_periodic                     = 1   1   0  
+zlo.type                                 = slip_wall
+zlo.temperature_type                     = fixed_gradient
+zlo.temperature                          = 0.0										 
+zhi.type                                 = slip_wall
+zhi.temperature_type                     = fixed_gradient
+zhi.temperature                          = 0.003
+
+#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨#
+#               PHYSICS                 #
+#.......................................#
+incflo.physics                           = MultiPhase ABL OceanWaves
+ICNS.source_terms                        = BoussinesqBuoyancy CoriolisForcing GeostrophicForcing GravityForcing NonLinearSGSTerm
+ICNS.use_perturb_pressure                = true
+incflo.velocity                          = 10.0 0.0 0.0
+GeostrophicForcing.geostrophic_wind      = 10.0 0.0 0.0
+GeostrophicForcing.wind_forcing_off_height  = 12.0
+GeostrophicForcing.wind_forcing_ramp_height = 16.0
+CoriolisForcing.latitude                 = 90.0
+CoriolisForcing.north_vector             = 0.0 1.0 0.0
+CoriolisForcing.east_vector              = 1.0 0.0 0.0
+BoussinesqBuoyancy.reference_temperature = 290
+RayleighDamping.reference_velocity       = 10.0 0 0
+RayleighDamping.length_sloped_damping    = 400
+RayleighDamping.length_complete_damping  = 200
+RayleighDamping.time_scale               = 20.0
+RayleighDamping.force_coord_directions   = 0 0 1								
+ABL.reference_temperature                = 290
+ABL.temperature_heights                  = 0.0 300 1300 2300
+ABL.temperature_values                   = 290 290 293 296
+ABL.perturb_temperature                  = false
+ABL.cutoff_height                        = 50.0
+ABL.perturb_velocity                     = true
+ABL.perturb_ref_height                   = 50.0
+ABL.Uperiods                             = 4.0
+ABL.Vperiods                             = 4.0
+ABL.deltaU                               = 1e-5
+ABL.deltaV                               = 1e-5
+ABL.kappa                                = .41
+ABL.surface_roughness_z0                 = 0.0001
+ABL.surface_temp_flux                    = 0.0
+ABL.wall_shear_stress_type               = local
+
+OceanWaves.label                         = Wave1
+OceanWaves.Wave1.type                    = LinearWaves
+OceanWaves.Wave1.wave_height             = 6
+OceanWaves.Wave1.wave_length             = 400.0
+OceanWaves.Wave1.water_depth             = 200.0
+OceanWaves.Wave1.relax_zone_gen_length   = 400.0
+OceanWaves.Wave1.relax_zone_out_length   = 800.0
+OceanWaves.Wave1.initialize_wave_field   = true
+
+#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨#
+#          POST-Processing              #
+#.......................................#
+
+
+#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨#
+#              AVERAGING                #
+#.......................................#
+
+
+#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨#
+#            MESH REFINEMENT            #
+#.......................................#
+tagging.labels                           = static
+tagging.static.type                      = CartBoxRefinement
+tagging.static.static_refinement_def     = static_box.txt
+
+#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨#
+#               TURBINES                #
+#.......................................#
+mac_proj.num_pre_smooth                            = 8
+mac_proj.num_post_smooth                           = 8
+mac_proj.mg_rtol                                   = -1
+mac_proj.mg_atol                                   = 1e-4
+mac_proj.maxiter                                   = 25
+mac_proj.fmg_maxiter                               = 4
+nodal_proj.num_pre_smooth                          = 8
+nodal_proj.num_post_smooth                         = 8
+nodal_proj.mg_rtol                                 = -1
+nodal_proj.mg_atol                                 = 1e-4
+nodal_proj.maxiter                                 = 25
+nodal_proj.fmg_maxiter                             = 4
+diffusion.mg_rtol                                  = -1
+diffusion.mg_atol                                  = 1e-4
+temperature_diffusion.mg_rtol                      = -1
+temperature_diffusion.mg_atol                      = 1e-4
+