Adding a dynamic stress (wave model) to wall_models

amr-wind/boundary_conditions/wall_models/MOSD.H

@@ -0,0 +1,57 @@
+#ifndef MOSD_H
+#define MOSD_H
+
+#include "AMReX_AmrCore.H"
+#include <AMReX.H>
+namespace amr_wind {
+struct MOSD
+{
+    /*
+     * A dynamic wall model that calculates the stress from wave to wind
+     * based on geometric information of the wave.
+     */
+
+    amrex::Real amplitude{0.05};
+    amrex::Real wavenumber{4};
+    amrex::Real omega{0.8};
+    amrex::Real time;
+
+    AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real get_dyn_tau(
+        const amrex::Real u_dx,
+        const amrex::Real v_dx,
+        const amrex::Real xc,
+        const amrex::Real unit_nor) const
+    {
+
+        // Building the wave surface, gradients, wave velocities and unit normal
+        const amrex::Real dx_eta_wave =
+            -amplitude * wavenumber * std::sin(wavenumber * xc - omega * time);
+        const amrex::Real dy_eta_wave = 0;
+        const amrex::Real dt_eta_wave =
+            amplitude * omega * std::sin(wavenumber * xc - omega * time);
+        const amrex::Real grad_eta_wave =
+            std::sqrt(dx_eta_wave * dx_eta_wave + dy_eta_wave * dy_eta_wave);
+        const amrex::Real Cx_wave =
+            -dt_eta_wave * dx_eta_wave * (1 / (grad_eta_wave * grad_eta_wave));
+        const amrex::Real Cy_wave =
+            -dt_eta_wave * dy_eta_wave * (1 / (grad_eta_wave * grad_eta_wave));
+        const amrex::Real n_x = dx_eta_wave / grad_eta_wave;
+        const amrex::Real n_y = dy_eta_wave / grad_eta_wave;
+
+        // Calculating the relative velocity, heaviside function
+        const amrex::Real u_r = u_dx - Cx_wave;
+        const amrex::Real v_r = v_dx - Cy_wave;
+        const amrex::Real ur_mag =
+            std::sqrt(u_r * u_r * n_x * n_x + v_r * v_r * n_y * n_y);
+        const amrex::Real Heavi_arg = (u_r * dx_eta_wave + v_r * dy_eta_wave);
+        const amrex::Real Heavi =
+            (Heavi_arg + std::abs(Heavi_arg)) / (2 * Heavi_arg);
+
+        // Calculating the magnitude of the stress
+        return (1 / M_PI) * ur_mag * ur_mag * grad_eta_wave * grad_eta_wave *
+               Heavi * (unit_nor == 0 ? n_x : n_y);
+    }
+};
+
+} // namespace amr_wind
+#endif /* MOSD_H */

amr-wind/boundary_conditions/wall_models/ShearStressSimple.H

@@ -3,6 +3,7 @@
 
 #include "amr-wind/boundary_conditions/wall_models/LogLaw.H"
 #include "amr-wind/wind_energy/ShearStress.H"
+#include "amr-wind/boundary_conditions/wall_models/MOSD.H"
 
 namespace amr_wind {
 
@@ -12,8 +13,13 @@ struct SimpleShearSchumann
         : utau2(ll.utau_mean * ll.utau_mean), wspd_mean(ll.wspd_mean), m_ll(ll)
     {}
 
-    AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real
-    get_shear(amrex::Real u, amrex::Real /* wspd */) const
+    AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real get_shear(
+        amrex::Real u,
+        amrex::Real /* wspd */,
+        amrex::Real /*u_dx*/,
+        amrex::Real /*v_dx*/,
+        amrex::Real /*x_c*/,
+        amrex::Real /*unit_nor*/) const
     {
         return u / wspd_mean * utau2;
     };
@@ -26,15 +32,48 @@ struct SimpleShearLogLaw
 {
     explicit SimpleShearLogLaw(const amr_wind::LogLaw& ll) : m_ll(ll) {}
 
-    AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real
-    get_shear(amrex::Real u, amrex::Real wspd) const
+    AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real get_shear(
+        amrex::Real u,
+        amrex::Real wspd,
+        amrex::Real /*u_dx*/,
+        amrex::Real /*v_dx*/,
+        amrex::Real /*x_c*/,
+        amrex::Real /*unit_nor*/) const
     {
         amrex::Real utau = m_ll.get_utau(wspd);
         return utau * utau * u / wspd;
     };
 
     const amr_wind::LogLaw m_ll;
 };
+
+struct SimpleShearMOSD
+{
+    explicit SimpleShearMOSD(
+        const amr_wind::LogLaw& ll, const amr_wind::MOSD& md)
+        : m_ll(ll), m_md(md)
+    {}
+
+    AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real get_shear(
+        amrex::Real u,
+        amrex::Real wspd,
+        amrex::Real u_dx,
+        amrex::Real v_dx,
+        amrex::Real x_c,
+        amrex::Real unit_nor) const
+    {
+        amrex::Real utau = m_ll.get_utau(wspd);
+        amrex::Real tau_vis = utau * utau * u / wspd;
+
+        amrex::Real tau_wave = m_md.get_dyn_tau(u_dx, v_dx, x_c, unit_nor);
+
+        return tau_vis + tau_wave;
+    };
+
+    const amr_wind::LogLaw m_ll;
+    const amr_wind::MOSD m_md;
+};
+
 } // namespace amr_wind
 
 #endif /* SHEARSTRESSSIMPLE_H */

amr-wind/boundary_conditions/wall_models/WallFunction.H

@@ -5,6 +5,7 @@
 #include "amr-wind/core/FieldBCOps.H"
 #include "amr-wind/utilities/FieldPlaneAveragingFine.H"
 #include "amr-wind/boundary_conditions/wall_models/LogLaw.H"
+#include "amr-wind/boundary_conditions/wall_models/MOSD.H"
 
 namespace amr_wind {
 
@@ -22,10 +23,12 @@ public:
 
     amrex::Real utau() const { return m_log_law.utau_mean; }
     LogLaw log_law() const { return m_log_law; }
+    MOSD mosd() const { return m_mosd; }
 
     //! Update the mean velocity at a given timestep
     void update_umean();
     void update_utau_mean();
+    void update_time();
 
     ~WallFunction() = default;
 
@@ -38,6 +41,8 @@ private:
     LogLaw m_log_law;
     int m_direction{2}; ///< Direction normal to wall, hardcoded to z
     VelPlaneAveragingFine m_pa_vel;
+
+    MOSD m_mosd;
 };
 
 /** Applies a shear-stress value at the domain boundary

amr-wind/boundary_conditions/wall_models/WallFunction.cpp

@@ -17,6 +17,7 @@ WallFunction::WallFunction(CFDSim& sim)
 {
     amrex::Real mu;
     amrex::Real rho{1.0};
+
     {
         amrex::ParmParse pp("BodyForce");
         amrex::Vector<amrex::Real> body_force{0.0, 0.0, 0.0};
@@ -44,6 +45,13 @@ WallFunction::WallFunction(CFDSim& sim)
             (geom.ProbLo(m_direction) +
              (m_log_law.ref_index + 0.5) * geom.CellSize(m_direction));
     }
+    {
+        amrex::ParmParse pp(
+            "wave_mosd"); // "wave_mosd" is the prefix used in the input file
+        pp.query("amplitude", m_mosd.amplitude);
+        pp.query("wavenumber", m_mosd.wavenumber);
+        pp.query("frequency", m_mosd.omega);
+    }
 }
 
 VelWallFunc::VelWallFunc(Field& /*unused*/, WallFunction& wall_func)
@@ -55,7 +63,8 @@ VelWallFunc::VelWallFunc(Field& /*unused*/, WallFunction& wall_func)
 
     if (m_wall_shear_stress_type == "constant" ||
         m_wall_shear_stress_type == "log_law" ||
-        m_wall_shear_stress_type == "schumann") {
+        m_wall_shear_stress_type == "schumann" ||
+        m_wall_shear_stress_type == "mosd") {
         amrex::Print() << "Shear Stress model: " << m_wall_shear_stress_type
                        << std::endl;
     } else {
@@ -92,6 +101,9 @@ void VelWallFunc::wall_model(
         const auto& vold_lev = velocity.state(FieldState::Old)(lev);
         const auto& eta_lev = viscosity(lev);
 
+        const auto& problo = geom.ProbLoArray();
+        const auto& dx = geom.CellSizeArray();
+
         if (amrex::Gpu::notInLaunchRegion()) {
             mfi_info.SetDynamic(true);
         }
@@ -116,14 +128,103 @@ void VelWallFunc::wall_model(
                         const amrex::Real vv =
                             vold_arr(i, j, k + idx_offset, 1);
                         const amrex::Real wspd = std::sqrt(uu * uu + vv * vv);
-                        // Dirichlet BC
-                        varr(i, j, k - 1, 2) = 0.0;
+                        const amrex::Real xc = problo[0] + (i + 0.5) * dx[0];
 
-                        // Shear stress BC
-                        varr(i, j, k - 1, 0) =
-                            tau.get_shear(uu, wspd) / mu * den(i, j, k);
-                        varr(i, j, k - 1, 1) =
-                            tau.get_shear(vv, wspd) / mu * den(i, j, k);
+                        // Check if dx[0] > dx[2]
+                        if (dx[0] <= dx[2]) {
+                            // Handle the case where dx[0] <= dx[2]
+                            // Use the nearest cell value without interpolation
+                            const int k_nearest =
+                                static_cast<int>(std::round(dx[0] / dx[2]));
+                            const amrex::Real u_dx =
+                                vold_arr(i, j, k_nearest, 0);
+                            const amrex::Real v_dx =
+                                vold_arr(i, j, k_nearest, 1);
+
+                            varr(i, j, k - 1, 0) =
+                                tau.get_shear(uu, wspd, u_dx, v_dx, xc, 0) /
+                                mu * den(i, j, k);
+                            varr(i, j, k - 1, 1) =
+                                tau.get_shear(vv, wspd, u_dx, v_dx, xc, 1) /
+                                mu * den(i, j, k);
+                        } else {
+                            // Case where dx[0] > dx[2]
+                            const int k_lo = static_cast<int>(
+                                std::floor(dx[0] / dx[2] - 0.5));
+                            const int k_hi = k_lo + 1;
+
+                            // Ensure the interpolation points are within the
+                            // valid box
+                            if (bx.contains(amrex::IntVect(i, j, k_lo)) &&
+                                bx.contains(amrex::IntVect(i, j, k_hi))) {
+                                const amrex::Real z_lo = (k_lo + 0.5) * dx[2];
+                                const amrex::Real z_hi = (k_hi + 0.5) * dx[2];
+                                const amrex::Real z_interp = dx[0];
+
+                                // Check if the interpolation condition is
+                                // satisfied
+                                if (z_lo < z_interp && z_interp < z_hi) {
+                                    const amrex::Real weight =
+                                        (z_interp - z_lo) / (z_hi - z_lo);
+                                    const amrex::Real u_low =
+                                        vold_arr(i, j, k_lo, 0);
+                                    const amrex::Real u_up =
+                                        vold_arr(i, j, k_hi, 0);
+                                    const amrex::Real v_low =
+                                        vold_arr(i, j, k_lo, 1);
+                                    const amrex::Real v_up =
+                                        vold_arr(i, j, k_hi, 1);
+
+                                    const amrex::Real u_dx =
+                                        u_low + (u_up - u_low) * weight;
+                                    const amrex::Real v_dx =
+                                        v_low + (v_up - v_low) * weight;
+
+                                    varr(i, j, k - 1, 0) =
+                                        tau.get_shear(
+                                            uu, wspd, u_dx, v_dx, xc, 0) /
+                                        mu * den(i, j, k);
+                                    varr(i, j, k - 1, 1) =
+                                        tau.get_shear(
+                                            vv, wspd, u_dx, v_dx, xc, 1) /
+                                        mu * den(i, j, k);
+                                } else {
+                                    // Fallback: use nearest point if
+                                    // interpolation condition is not met
+                                    const int k_nearest =
+                                        (z_interp - z_lo < z_hi - z_interp)
+                                            ? k_lo
+                                            : k_hi;
+                                    const amrex::Real u_dx =
+                                        vold_arr(i, j, k_nearest, 0);
+                                    const amrex::Real v_dx =
+                                        vold_arr(i, j, k_nearest, 1);
+
+                                    varr(i, j, k - 1, 0) =
+                                        tau.get_shear(
+                                            uu, wspd, u_dx, v_dx, xc, 0) /
+                                        mu * den(i, j, k);
+                                    varr(i, j, k - 1, 1) =
+                                        tau.get_shear(
+                                            vv, wspd, u_dx, v_dx, xc, 1) /
+                                        mu * den(i, j, k);
+                                }
+                            } else {
+                                // Fallback: use the original cell value if
+                                // interpolation points are out of bounds
+                                const amrex::Real u_dx = vold_arr(i, j, k, 0);
+                                const amrex::Real v_dx = vold_arr(i, j, k, 1);
+
+                                varr(i, j, k - 1, 0) =
+                                    tau.get_shear(uu, wspd, u_dx, v_dx, xc, 0) /
+                                    mu * den(i, j, k);
+                                varr(i, j, k - 1, 1) =
+                                    tau.get_shear(vv, wspd, u_dx, v_dx, xc, 1) /
+                                    mu * den(i, j, k);
+                            }
+                        }
+
+                        varr(i, j, k - 1, 2) = 0.0;
                     });
             }
 
@@ -143,9 +244,11 @@ void VelWallFunc::wall_model(
 
                         // Shear stress BC
                         varr(i, j, k, 0) =
-                            -tau.get_shear(uu, wspd) / mu * den(i, j, k);
+                            -tau.get_shear(uu, wspd, 0, 0, 0, 0) / mu *
+                            den(i, j, k);
                         varr(i, j, k, 1) =
-                            -tau.get_shear(vv, wspd) / mu * den(i, j, k);
+                            -tau.get_shear(vv, wspd, 0, 0, 0, 1) / mu *
+                            den(i, j, k);
                     });
             }
         }
@@ -238,6 +341,12 @@ void VelWallFunc::operator()(Field& velocity, const FieldState rho_state)
         m_wall_func.update_umean();
         auto tau = SimpleShearSchumann(m_wall_func.log_law());
         wall_model(velocity, rho_state, tau);
+    } else if (m_wall_shear_stress_type == "mosd") {
+        m_wall_func.update_umean();
+        m_wall_func.update_utau_mean();
+        m_wall_func.update_time();
+        auto tau = SimpleShearMOSD(m_wall_func.log_law(), m_wall_func.mosd());
+        wall_model(velocity, rho_state, tau);
     }
 }
 
@@ -249,4 +358,6 @@ void WallFunction::update_umean()
 }
 
 void WallFunction::update_utau_mean() { m_log_law.update_utau_mean(); }
+
+void WallFunction::update_time() { m_mosd.time = m_sim.time().current_time(); }
 } // namespace amr_wind

docs/sphinx/theory/theory.rst

@@ -427,6 +427,23 @@ z direction (example: *half-channel* simulations) at the centerline.
        w &= 0
    \end{aligned}
 
+Dynamic wall model (Wave model)
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+This wall model is used to calculate the stress due to moving surfaces,
+like ocean waves. It aims to introduce wave phase-resolving physics 
+at a cost similar to using the Log-law wall model, without the need of using 
+wave adapting computational grids. The model developed by `Ayala et al (2024)
+<https://doi.org/10.1007/s10546-024-00884-8>`_ has two components:
+
+.. math:: \tau_{i3} = \frac{1}{\pi}|(\boldsymbol{u-C}) \cdot \boldsymbol{\hat{n}}|^2|\boldsymbol{\nabla} \eta|^2 \, \hat{n}_i  \, \text{H} \Bigl[ (u_j-C_j)\frac{\partial \eta}{\partial x_j} \Bigr] \, + \, \tau^{visc}_{i3}, \quad i = 1,2.
+
+The first component gives the form drag due to ocean waves, where, :math::`\boldsymbol{C}` 
+is the wave velocity vector, :math:: `\eta` is the surface height distribution and 
+:math:: `\hat{\boldsymbol n} =  \boldsymbol \nabla \eta /|\boldsymbol{\nabla} \eta|`. The 
+second component (:math:: `\tau^{visc}_{i3}`) is the stress due to unresolved effects, 
+like viscous effects. For this component, the ``Log-law wall model`` is used. 
+
 Navigating source code
 ------------------------
 

docs/sphinx/user/features.rst

@@ -61,7 +61,7 @@ Methods and models
 
    * Large Eddy Simulation: constant Smagorinsky,  AMD, one equation :math:`k_{sgs}`, Kosovic [:ref:`doc <turbulence>`, :ref:`inp <inputs_turbulence>`]
 
-   * Wall models: log-law, constant stress, Schumann [:ref:`doc <wall_models>`, :ref:`inp <inputs_abl>`]
+   * Wall models: log-law, constant stress, Schumann [:ref:`doc <wall_models>`, :ref:`inp <inputs_abl>`], dynamic (wave model) [:ref:`doc <wall_models>`, :ref:`inp <inputs_boundary_conditions>` 
 
    * Reynolds-Average Navier-Stokes: :math:`k`-:math:`\omega` SST (and IDDES variant) [:ref:`doc <turbulence>`, :ref:`inp <inputs_turbulence>`]