Joukowsky: Fix bug and add Cp (Exawind#633)

* Remove second theta division * Add cp computation * Use dimensional moment calculation * Fix error * Fix D vs R bug and add inline function to help in the future
stephan-rohr · Jun 16, 2022 · 9d851fc · 9d851fc
1 parent 88b7512
commit 9d851fc
Show file tree

Hide file tree

Showing 4 changed files with 28 additions and 7 deletions.
diff --git a/amr-wind/wind_energy/actuator/disk/ActuatorDisk.H b/amr-wind/wind_energy/actuator/disk/ActuatorDisk.H
@@ -29,6 +29,10 @@ struct DiskBaseData
     RealList thrust_coeff;
     RealList table_velocity;
     std::string spreading_type{"LinearBasis"};
+    AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real radius() const
+    {
+        return 0.5 * diameter;
+    }
 };
 
 } // namespace actuator

diff --git a/amr-wind/wind_energy/actuator/disk/Joukowsky.H b/amr-wind/wind_energy/actuator/disk/Joukowsky.H
@@ -15,6 +15,7 @@ struct JoukowskyData : public DiskBaseData
     RealList root_correction;
     amrex::Real current_angular_velocity{0.0};
     amrex::Real vortex_core_size;
+    amrex::Real current_cp;
     // --- Sorenson 2020 equation 10 constants ----
     amrex::Real root_correction_coefficient{1.256};
     amrex::Real root_correction_exponent{2.0};

diff --git a/amr-wind/wind_energy/actuator/disk/Joukowsky_ops.H b/amr-wind/wind_energy/actuator/disk/Joukowsky_ops.H
@@ -68,7 +68,7 @@ struct InitDataOp<Joukowsky, ActSrcDisk>
             const auto b = meta.root_correction_exponent;
 
             const auto delta = std::max(
-                meta.vortex_core_size * 2.0 / meta.diameter,
+                meta.vortex_core_size * meta.radius(),
                 vs::DTraits<amrex::Real>::eps());
 
             const auto factor = dx / delta;
@@ -122,14 +122,14 @@ struct ComputeForceOp<Joukowsky, ActSrcDisk>
         joukowsky::set_current_angular_velocity(ddata, uInfSqr);
 
         const amrex::Real Ct = ddata.current_ct;
-        const amrex::Real dx = ddata.dr * 2.0 / ddata.diameter;
+        const amrex::Real dx = ddata.dr / ddata.radius();
         const amrex::Real dTheta =
             amr_wind::utils::two_pi() / ddata.num_vel_pts_t;
 
         const amrex::Real geometry_factor = ddata.dr * ddata.dr * dTheta * 0.5;
 
         const amrex::Real lambda =
-            0.5 * ddata.diameter * ddata.current_angular_velocity / U_ref;
+            ddata.radius() * ddata.current_angular_velocity / U_ref;
         const amrex::Real lambda_2 = lambda * lambda;
 
         // step 0: compute tip correction based on current wind speed (eq 9)
@@ -172,11 +172,15 @@ struct ComputeForceOp<Joukowsky, ActSrcDisk>
         const amrex::Real q0 = numerator / denominator;
 
         // step 3: compute normal force (fz) and azimuthal force (f_theta) (eq
-        // 13)
+        // 13) and cp (eq 20)
 
         VecSlice disk_velocity = ::amr_wind::utils::slice(
             grid.vel, ddata.num_force_pts, ddata.num_force_pts);
+
+        amrex::Real moment = 0.0;
+
         int ip = 0;
+
         for (int i = 0; i < ddata.num_vel_pts_r; i++) {
             const amrex::Real& F = ddata.tip_correction[i];
             const amrex::Real& g = ddata.root_correction[i];
@@ -204,15 +208,25 @@ struct ComputeForceOp<Joukowsky, ActSrcDisk>
                 const amrex::Real dArea =
                     geometry_factor * (std::pow(i + 1.0, 2) - std::pow(i, 2));
 
+                // eq 18
+                moment += x * ddata.radius() * f_theta * dArea;
+
                 grid.force[ip] =
                     (f_z * ddata.normal_vec + f_theta * theta_vec) *
-                    ddata.density * uInfSqr * dArea / ddata.num_vel_pts_t;
+                    ddata.density * uInfSqr * dArea;
+
                 ddata.disk_force = ddata.disk_force + grid.force[ip];
             }
         }
 
-        // step 5: compute moment and power (eqs 18 and 19)
-        // step 6: compute cp
+        // equation 20
+        const amrex::Real eq_20_denominator = std::max(
+            0.5 * M_PI * ddata.density * std::pow(ddata.radius(), 2) *
+                std::pow(U_ref, 3),
+            machine_eps);
+
+        ddata.current_cp = ddata.density * uInfSqr * moment *
+                           ddata.current_angular_velocity / eq_20_denominator;
     }
 };
 

diff --git a/amr-wind/wind_energy/actuator/disk/Joukowsky_ops.cpp b/amr-wind/wind_energy/actuator/disk/Joukowsky_ops.cpp
@@ -85,6 +85,7 @@ void prepare_netcdf_file(
     grp.def_var("vref", NC_DOUBLE, {nt_name, "ndim"});
     grp.def_var("vdisk", NC_DOUBLE, {nt_name, "ndim"});
     grp.def_var("ct", NC_DOUBLE, {nt_name});
+    grp.def_var("cp", NC_DOUBLE, {nt_name});
     grp.def_var("density", NC_DOUBLE, {nt_name});
     grp.def_var("total_disk_force", NC_DOUBLE, {nt_name, "ndim"});
     grp.def_var("angular_velocity", NC_DOUBLE, {nt_name});
@@ -130,6 +131,7 @@ void write_netcdf(
     grp.var("vdisk").put(
         &(data.mean_disk_velocity[0]), {nt, 0}, {1, AMREX_SPACEDIM});
     grp.var("ct").put(&data.current_ct, {nt}, {1});
+    grp.var("cp").put(&data.current_cp, {nt}, {1});
     grp.var("density").put(&data.density, {nt}, {1});
     grp.var("total_disk_force")
         .put(&data.disk_force[0], {nt}, {1, AMREX_SPACEDIM});