WIP: Generalized actuator disk (#1875)

* First version of generalized actuator disk

* Reading in blade specification table

* Reading in files for generalized actuator disk

* Finalizing generalized actuator disk

* Correcting Windows error

* Correcting SYCL error

* Reading in table for rpm and blade pitch

* Correcting SYCL error

* Correcting unused variable error

* Correcting spelling error

* Correcting docs errors

* Correcting indexing error

---------

Co-authored-by: Mahesh Natarajan <[email protected]>
Co-authored-by: Mahesh Natarajan <[email protected]>
Co-authored-by: Mahesh Natarajan <[email protected]>
Co-authored-by: Mahesh Natarajan <[email protected]>
Co-authored-by: Ann Almgren <[email protected]>
Co-authored-by: Mahesh Natarajan <[email protected]>

Docs/sphinx_doc/theory/WindFarmModels.rst

@@ -127,7 +127,7 @@ The EWP model does not have a concept of intersected area by the turbine rotor l
 .. _actuator_disk_model_simplified:
 
 Simplified actuator disk model
-=================================
+-----------------------------------
 
 A simplified actuator disk model based on one-dimensional momentum theory is implemented (See Section 3.2 in `Wind Energy Handbook 2nd edition`_). A schematic of the actuator disk is shown in  Fig. :numref:`fig:ActuatorDisk_Schematic`.
 The model is implemented as source terms in the equations for the horizontal velocity components (ie. `x` and `y` directions). The thrust force from the one-dimensional momentum theory is given by
@@ -185,7 +185,7 @@ where :math:`dA` is the area of the actuator disk in the mesh cell (see Fig. :nu
 .. _generalized_actuator_disk_model:
 
 Generalized actuator disk model
-===============================
+------------------------------------
 
 The generalized actuator model (GAD) based on blade element theory (`Mirocha et. al. 2014`_, see Chapter 3 of `Small Wind Turbines`_) is implemented. Similar to the simplified actuator disk model, GAD also models the wind turbine as a disk, but takes into account the details of the blade geometry (see Fig. :numref:`fig:GAD_Schematic`). The forces on the blades in the x, y, z directions are computed, and that contributes to the source terms for the fluid momentum equations. The source terms in a mesh cell inside the actuator disk are given as:
 

Source/DataStructs/ERF_DataStruct.H

@@ -353,7 +353,8 @@ struct SolverChoice {
         pp.query("windfarm_loc_table",  windfarm_loc_table);
         pp.query("windfarm_spec_table", windfarm_spec_table);
         pp.query("windfarm_blade_table", windfarm_blade_table);
-        pp.query("windfarm_airofil_tables", windfarm_airfoil_tables);
+        pp.query("windfarm_airfoil_tables", windfarm_airfoil_tables);
+        pp.query("windfarm_spec_table_extra", windfarm_spec_table_extra);
 
         // Sampling distance upstream of the turbine to find the
         // incoming free stream velocity as a factor of the diameter of the
@@ -647,7 +648,7 @@ struct SolverChoice {
     // if     SAM, then it will be set to RhoQ4
     int RhoQr_comp {-1};
 
-    std::string windfarm_loc_table, windfarm_spec_table;
+    std::string windfarm_loc_table, windfarm_spec_table, windfarm_spec_table_extra;
     std::string windfarm_blade_table, windfarm_airfoil_tables;
     amrex::Real sampling_distance_by_D = -1.0;
     amrex::Real turb_disk_angle = -1.0;

Source/ERF_make_new_arrays.cpp

@@ -239,7 +239,7 @@ ERF::init_stuff (int lev, const BoxArray& ba, const DistributionMapping& dm,
         vars_windfarm[lev].define(ba, dm, 2, ngrow_state);// dudt, dvdt
     }
     if (solverChoice.windfarm_type == WindFarmType::GeneralAD) {
-        vars_windfarm[lev].define(ba, dm, 2, ngrow_state);// dudt, dvdt
+        vars_windfarm[lev].define(ba, dm, 3, ngrow_state);// dudt, dvdt, dwdt
     }
         Nturb[lev].define(ba, dm, 1, ngrow_state); // Number of turbines in a cell
         SMark[lev].define(ba, dm, 2, ngrow_state); // Free stream velocity/source term

Source/IO/ERF_Plotfile.cpp

@@ -101,13 +101,9 @@ ERF::setPlotVariables (const std::string& pp_plot_var_names, Vector<std::string>
                     tmp_plot_names.push_back(derived_names[i]);
                 }
             }
-            if(solverChoice.windfarm_type == WindFarmType::SimpleAD) {
-                if(derived_names[i] == "num_turb" or derived_names[i] == "SMark0" or derived_names[i] == "Smark1") {
-                    tmp_plot_names.push_back(derived_names[i]);
-                }
-            }
-            if(solverChoice.windfarm_type == WindFarmType::GeneralAD) {
-                if(derived_names[i] == "num_turb" or derived_names[i] == "SMark1") {
+            if( solverChoice.windfarm_type == WindFarmType::SimpleAD or
+                solverChoice.windfarm_type == WindFarmType::GeneralAD ) {
+                if(derived_names[i] == "num_turb" or derived_names[i] == "SMark0" or derived_names[i] == "SMark1") {
                     tmp_plot_names.push_back(derived_names[i]);
                 }
             }
@@ -486,8 +482,8 @@ ERF::WritePlotFile (int which, Vector<std::string> plot_var_names)
             mf_comp ++;
         }
 
-        if(containerHasElement(plot_var_names, "SMark0") and
-           solverChoice.windfarm_type == WindFarmType::SimpleAD) {
+        if( containerHasElement(plot_var_names, "SMark0") and
+           (solverChoice.windfarm_type == WindFarmType::SimpleAD or solverChoice.windfarm_type == WindFarmType::GeneralAD) ) {
              for ( MFIter mfi(mf[lev],TilingIfNotGPU()); mfi.isValid(); ++mfi)
             {
                 const Box& bx = mfi.tilebox();

Source/Initialization/ERF_init_windfarm.cpp

@@ -43,8 +43,10 @@ ERF::init_windfarm (int lev)
     }
 
     if(solverChoice.windfarm_type == WindFarmType::GeneralAD) {
-        //windfarm->read_windfarm_blade_table(solverChoice.windfarm_blade_table);
-        //windfarm->read_airfoil_tables
+        windfarm->read_windfarm_blade_table(solverChoice.windfarm_blade_table);
+        windfarm->read_windfarm_airfoil_tables(solverChoice.windfarm_airfoil_tables,
+                                               solverChoice.windfarm_blade_table);
+        windfarm->read_windfarm_spec_table_extra(solverChoice.windfarm_spec_table_extra);
     }
 }
 

Source/WindFarmParametrization/ERF_InitWindFarm.cpp

@@ -3,6 +3,8 @@
  */
 
 #include <ERF_WindFarm.H>
+#include <filesystem>
+#include <dirent.h>   // For POSIX directory handling
 
 using namespace amrex;
 
@@ -123,6 +125,14 @@ WindFarm::init_windfarm_lat_lon (const std::string windfarm_loc_table,
         xloc[it] = xloc[it] - xloc_min + windfarm_x_shift;
         yloc[it] = yloc[it] - yloc_min + windfarm_y_shift;
     }
+
+    FILE* file_xy_loc;
+    file_xy_loc = fopen("file_xy_loc_KingPlains.txt","w");
+
+    for(int it = 0;it<xloc.size(); it++){
+        fprintf(file_xy_loc,"%0.15g %0.15g %0.15g\n", xloc[it], yloc[it], 89.0);
+    }
+    fclose(file_xy_loc);
 }
 
 void
@@ -199,14 +209,152 @@ void
 WindFarm::read_windfarm_blade_table(const std::string windfarm_blade_table)
 {
     std::ifstream filename(windfarm_blade_table);
+    std::string line;
+    Real temp, var1, var2, var3;
     if (!filename.is_open()) {
         Error("You are using a generalized actuator disk model based on blade element theory. This needs info of blades."
                       " An entry erf.windfarm_blade_table is needed. Either the entry is missing or the file specified"
                       " in the entry - " + windfarm_blade_table + " is missing.");
     }
     else {
         Print() << "Reading in wind farm blade table: " << windfarm_blade_table << "\n";
+
+        // First 6 lines are comments
+
+        for (int i = 0; i < 6; ++i) {
+            if (std::getline(filename, line)) {  // Read one line into the array
+            }
+        }
+
+        while(filename >> var1 >> temp >> temp >> temp >> var2 >> var3 >> temp) {
+            bld_rad_loc.push_back(var1);
+            bld_twist.push_back(var2);
+            bld_chord.push_back(var3);
+            //int idx = bld_rad_loc.size()-1;
+            //printf("Values are = %0.15g %0.15g %0.15g\n", bld_rad_loc[idx], bld_twist[idx], bld_chord[idx]);
+        }
+        set_blade_spec(bld_rad_loc, bld_twist, bld_chord);
+        n_bld_sections = bld_rad_loc.size();
+    }
+}
+
+void
+WindFarm::read_windfarm_spec_table_extra(const std::string windfarm_spec_table_extra)
+{
+    // Open the file
+    std::ifstream file(windfarm_spec_table_extra);
+
+    // Check if file opened successfully
+    if (!file.is_open()) {
+        Abort("Error: You are using generalized wind farms option. This requires an input file erf.windfarm_spec_table_extra."
+              " Either this entry is missing in the inputs or the file specified -" + windfarm_spec_table_extra + " does"
+              " not exist. Exiting...");
+    } else {
+        printf("Reading in windfarm_spec_table_extra %s", windfarm_spec_table_extra.c_str());
+    }
+
+    // Ignore the first line (header)
+    std::string header;
+    std::getline(file, header);
+
+    // Variables to hold each row's values
+    double V, Cp, Ct, rpm, pitch, temp;
+
+    // Read the file row by row
+    while (file >> V) {
+        char comma;  // To ignore the commas
+        file >> comma >> Cp >> comma >> Ct >> comma >> temp >> comma >> temp >> comma
+             >> temp >> comma >> rpm >> comma >> pitch >> comma >> temp;
+
+        velocity.push_back(V);
+        C_P.push_back(Cp);
+        C_T.push_back(Ct);
+        rotor_RPM.push_back(rpm);
+        blade_pitch.push_back(pitch);
+    }
+
+    set_turb_spec_extra(velocity, C_P, C_T, rotor_RPM, blade_pitch);
+}
+
+
+void
+WindFarm::read_windfarm_airfoil_tables(const std::string windfarm_airfoil_tables,
+                                       const std::string windfarm_blade_table)
+{
+    DIR* dir;
+    struct dirent* entry;
+    std::vector<std::string> files;
+
+    // Check if directory exists
+    if ((dir = opendir(windfarm_airfoil_tables.c_str())) == nullptr) {
+       Abort("You are using a generalized actuator disk model based on blade element theory. This needs info of airfoil"
+             " cross sections over the span of the blade. There needs to be an entry erf.airfoil_tables which is the directory that"
+             " contains the angle of attack, Cl, Cd data for each airfoil cross-section. Either the entry is missing or the directory specified"
+             " in the entry - " + windfarm_airfoil_tables + " is missing. Exiting...");
+    }
+
+    // Loop through directory entries and collect filenames
+    while ((entry = readdir(dir)) != nullptr) {
+        // Skip special directory entries "." and ".."
+        if (std::string(entry->d_name) == "." || std::string(entry->d_name) == "..") {
+            continue;
+        }
+        files.emplace_back(windfarm_airfoil_tables + "/" + entry->d_name);  // Add file path to vector
     }
+
+    // Close the directory
+    closedir(dir);
+
+    if (files.empty()) {
+        Abort("It seems the directory containing the info of airfoil cross sections of the blades - " + windfarm_airfoil_tables +
+              " is empty. Exiting...");
+    }
+
+    if(files.size() != static_cast<long double>(n_bld_sections)) {
+        printf("There are %d airfoil sections in the last column of %s. But the number"
+               " of files in %s is only %ld.\n", n_bld_sections, windfarm_blade_table.c_str(),
+                windfarm_airfoil_tables.c_str(), files.size());
+        Abort("The number of blade sections from " + windfarm_blade_table + " should match the number of"
+              " files in " + windfarm_airfoil_tables + ". Exiting...");
+    }
+
+    // Sort filenames in lexicographical (alphabetical) order
+    std::sort(files.begin(), files.end());
+
+    // Process each file
+    int count = 0;
+    bld_airfoil_aoa.resize(n_bld_sections);
+    bld_airfoil_Cl.resize(n_bld_sections);
+    bld_airfoil_Cd.resize(n_bld_sections);
+    for (const auto& filePath : files) {
+        std::ifstream filename(filePath.c_str());
+
+        if (!filename.is_open()) {
+            std::cerr << "Failed to open file: " << filePath << std::endl;
+            continue;  // Move on to the next file
+        }
+
+           std::cout << "Reading file: " << filePath << std::endl;
+
+        std::string line;
+        for (int i = 0; i < 54; ++i) {
+            if (std::getline(filename, line)) {  // Read one line into the array
+            }
+        }
+
+        Real var1, var2, var3, temp;
+
+        while(filename >> var1 >> var2 >> var3 >> temp) {
+            bld_airfoil_aoa[count].push_back(var1);
+            bld_airfoil_Cl[count].push_back(var2);
+            bld_airfoil_Cd[count].push_back(var3);
+            //int idx = bld_airfoil_aoa.size()-1;
+            //printf("Values are = %0.15g %0.15g %0.15g\n", bld_airfoil_aoa[idx], bld_airfoil_Cl[idx], bld_airfoil_Cd[idx]);
+        }
+        count++;
+    }
+
+    set_blade_airfoil_spec(bld_airfoil_aoa, bld_airfoil_Cl, bld_airfoil_Cd);
 }
 
 void

Source/WindFarmParametrization/ERF_WindFarm.H

@@ -65,7 +65,10 @@ public:
 
     void read_windfarm_blade_table(const std::string windfarm_blade_table);
 
-    void read_windfarm_airofil_tables(const std::string windfarm_airfoils_tables);
+    void read_windfarm_airfoil_tables(const std::string windfarm_airfoil_tables,
+                                      const std::string windfarm_blade_table);
+
+    void read_windfarm_spec_table_extra(const std::string windfarm_spec_table_extra);
 
     void fill_Nturb_multifab(const amrex::Geometry& geom,
                              amrex::MultiFab& mf_Nturb);
@@ -108,9 +111,32 @@ public:
         m_windfarm_model[0]->set_turb_loc(a_xloc, a_yloc);
     }
 
-    void set_turb_disk_angle (const amrex::Real& turb_disk_angle) override
+    void set_turb_disk_angle (const amrex::Real& a_turb_disk_angle) override
+    {
+        m_windfarm_model[0]->set_turb_disk_angle(a_turb_disk_angle);
+    }
+
+    void set_blade_spec (const amrex::Vector<amrex::Real>& a_bld_rad_loc,
+                         const amrex::Vector<amrex::Real>& a_bld_twist,
+                         const amrex::Vector<amrex::Real>& a_bld_chord) override
+    {
+        m_windfarm_model[0]->set_blade_spec(a_bld_rad_loc, a_bld_twist, a_bld_chord);
+    }
+
+    void set_blade_airfoil_spec (const amrex::Vector<amrex::Vector<amrex::Real>>& a_bld_airfoil_aoa,
+                                 const amrex::Vector<amrex::Vector<amrex::Real>>& a_bld_airfoil_Cl,
+                                 const amrex::Vector<amrex::Vector<amrex::Real>>& a_bld_airfoil_Cd) override
+    {
+        m_windfarm_model[0]->set_blade_airfoil_spec(a_bld_airfoil_aoa, a_bld_airfoil_Cl, a_bld_airfoil_Cd);
+    }
+
+    void set_turb_spec_extra (const amrex::Vector<amrex::Real>& a_velocity,
+                              const amrex::Vector<amrex::Real>& a_C_P,
+                              const amrex::Vector<amrex::Real>& a_C_T,
+                              const amrex::Vector<amrex::Real>& a_rotor_RPM,
+                              const amrex::Vector<amrex::Real>& a_blade_pitch) override
     {
-        m_windfarm_model[0]->set_turb_disk_angle(turb_disk_angle);
+        m_windfarm_model[0]->set_turb_spec_extra(a_velocity, a_C_P, a_C_T, a_rotor_RPM, a_blade_pitch);
     }
 
 protected:
@@ -119,6 +145,10 @@ protected:
     amrex::Real my_turb_disk_angle;
     amrex::Real hub_height, rotor_rad, thrust_coeff_standing, nominal_power;
     amrex::Vector<amrex::Real> wind_speed, thrust_coeff, power;
+    amrex::Vector<amrex::Real> bld_rad_loc, bld_twist, bld_chord;
+    amrex::Vector<amrex::Vector<amrex::Real>> bld_airfoil_aoa, bld_airfoil_Cl, bld_airfoil_Cd;
+    int n_bld_sections;
+    amrex::Vector<amrex::Real> velocity, C_P, C_T, rotor_RPM, blade_pitch;
 
    /*! \brief Create and set the specified windfarm model */
     template<class NewWindFarmModel>