diff --git a/Docs/sphinx_doc/theory/WindFarmModels.rst b/Docs/sphinx_doc/theory/WindFarmModels.rst index 3fff9f6dd..d4ab0a183 100644 --- a/Docs/sphinx_doc/theory/WindFarmModels.rst +++ b/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: diff --git a/Source/DataStructs/ERF_DataStruct.H b/Source/DataStructs/ERF_DataStruct.H index e44795be1..b3e2992ae 100644 --- a/Source/DataStructs/ERF_DataStruct.H +++ b/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; diff --git a/Source/ERF_make_new_arrays.cpp b/Source/ERF_make_new_arrays.cpp index 68630c53a..a20f4517e 100644 --- a/Source/ERF_make_new_arrays.cpp +++ b/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 diff --git a/Source/IO/ERF_Plotfile.cpp b/Source/IO/ERF_Plotfile.cpp index e26ee2ec3..75557acd4 100644 --- a/Source/IO/ERF_Plotfile.cpp +++ b/Source/IO/ERF_Plotfile.cpp @@ -101,13 +101,9 @@ ERF::setPlotVariables (const std::string& pp_plot_var_names, Vector 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 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(); diff --git a/Source/Initialization/ERF_init_windfarm.cpp b/Source/Initialization/ERF_init_windfarm.cpp index bf7d3dcc2..bc41acf57 100644 --- a/Source/Initialization/ERF_init_windfarm.cpp +++ b/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); } } diff --git a/Source/WindFarmParametrization/ERF_InitWindFarm.cpp b/Source/WindFarmParametrization/ERF_InitWindFarm.cpp index 32fae15ce..99a6f2f20 100644 --- a/Source/WindFarmParametrization/ERF_InitWindFarm.cpp +++ b/Source/WindFarmParametrization/ERF_InitWindFarm.cpp @@ -3,6 +3,8 @@ */ #include +#include +#include // 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> 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 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(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 diff --git a/Source/WindFarmParametrization/ERF_WindFarm.H b/Source/WindFarmParametrization/ERF_WindFarm.H index ee7742de9..100a81f0a 100644 --- a/Source/WindFarmParametrization/ERF_WindFarm.H +++ b/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& a_bld_rad_loc, + const amrex::Vector& a_bld_twist, + const amrex::Vector& 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>& a_bld_airfoil_aoa, + const amrex::Vector>& a_bld_airfoil_Cl, + const amrex::Vector>& 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& a_velocity, + const amrex::Vector& a_C_P, + const amrex::Vector& a_C_T, + const amrex::Vector& a_rotor_RPM, + const amrex::Vector& 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 wind_speed, thrust_coeff, power; + amrex::Vector bld_rad_loc, bld_twist, bld_chord; + amrex::Vector> bld_airfoil_aoa, bld_airfoil_Cl, bld_airfoil_Cd; + int n_bld_sections; + amrex::Vector velocity, C_P, C_T, rotor_RPM, blade_pitch; /*! \brief Create and set the specified windfarm model */ template diff --git a/Source/WindFarmParametrization/GeneralActuatorDisk/ERF_AdvanceGeneralAD.cpp b/Source/WindFarmParametrization/GeneralActuatorDisk/ERF_AdvanceGeneralAD.cpp index 96243d29b..42f9d4be3 100644 --- a/Source/WindFarmParametrization/GeneralActuatorDisk/ERF_AdvanceGeneralAD.cpp +++ b/Source/WindFarmParametrization/GeneralActuatorDisk/ERF_AdvanceGeneralAD.cpp @@ -1,6 +1,7 @@ #include #include #include +#include using namespace amrex; @@ -18,14 +19,17 @@ GeneralAD::advance (const Geometry& geom, AMREX_ALWAYS_ASSERT(W_old.nComp() > 0); AMREX_ALWAYS_ASSERT(mf_Nturb.nComp() > 0); AMREX_ALWAYS_ASSERT(mf_vars_generalAD.nComp() > 0); + compute_freestream_velocity(cons_in, U_old, V_old, mf_SMark); source_terms_cellcentered(geom, cons_in, mf_SMark, mf_vars_generalAD); - update(dt_advance, cons_in, U_old, V_old, mf_vars_generalAD); + update(dt_advance, cons_in, U_old, V_old, W_old, mf_vars_generalAD); } void GeneralAD::update (const Real& dt_advance, MultiFab& cons_in, - MultiFab& U_old, MultiFab& V_old, + MultiFab& U_old, + MultiFab& V_old, + MultiFab& W_old, const MultiFab& mf_vars_generalAD) { @@ -33,12 +37,14 @@ GeneralAD::update (const Real& dt_advance, Box tbx = mfi.nodaltilebox(0); Box tby = mfi.nodaltilebox(1); + Box tbz = mfi.nodaltilebox(2); auto generalAD_array = mf_vars_generalAD.array(mfi); auto u_vel = U_old.array(mfi); auto v_vel = V_old.array(mfi); + auto w_vel = W_old.array(mfi); - ParallelFor(tbx, tby, + ParallelFor(tbx, tby, tbz, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { u_vel(i,j,k) = u_vel(i,j,k) + (generalAD_array(i-1,j,k,0) + generalAD_array(i,j,k,0))/2.0*dt_advance; @@ -46,10 +52,235 @@ GeneralAD::update (const Real& dt_advance, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { v_vel(i,j,k) = v_vel(i,j,k) + (generalAD_array(i,j-1,k,1) + generalAD_array(i,j,k,1))/2.0*dt_advance; + }, + [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept + { + w_vel(i,j,k) = w_vel(i,j,k) + (generalAD_array(i,j,k-1,2) + generalAD_array(i,j,k,2))/2.0*dt_advance; }); + } } +void GeneralAD::compute_freestream_velocity(const MultiFab& cons_in, + const MultiFab& U_old, + const MultiFab& V_old, + const MultiFab& mf_SMark) +{ + get_turb_loc(xloc, yloc); + freestream_velocity.clear(); + freestream_phi.clear(); + disk_cell_count.clear(); + freestream_velocity.resize(xloc.size(),0.0); + freestream_phi.resize(xloc.size(),0.0); + disk_cell_count.resize(xloc.size(),0.0); + + Gpu::DeviceVector d_freestream_velocity(xloc.size()); + Gpu::DeviceVector d_freestream_phi(yloc.size()); + Gpu::DeviceVector d_disk_cell_count(yloc.size()); + Gpu::copy(Gpu::hostToDevice, freestream_velocity.begin(), freestream_velocity.end(), d_freestream_velocity.begin()); + Gpu::copy(Gpu::hostToDevice, freestream_phi.begin(), freestream_phi.end(), d_freestream_phi.begin()); + Gpu::copy(Gpu::hostToDevice, disk_cell_count.begin(), disk_cell_count.end(), d_disk_cell_count.begin()); + + Real* d_freestream_velocity_ptr = d_freestream_velocity.data(); + Real* d_freestream_phi_ptr = d_freestream_phi.data(); + Real* d_disk_cell_count_ptr = d_disk_cell_count.data(); + + + for ( MFIter mfi(cons_in,TilingIfNotGPU()); mfi.isValid(); ++mfi) { + + auto SMark_array = mf_SMark.array(mfi); + auto u_vel = U_old.array(mfi); + auto v_vel = V_old.array(mfi); + Box tbx = mfi.nodaltilebox(0); + + ParallelFor(tbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + + if(SMark_array(i,j,k,0) != -1.0) { + int turb_index = static_cast(SMark_array(i,j,k,0)); + Real phi = std::atan2(v_vel(i,j,k),u_vel(i,j,k)); // Wind direction w.r.t the x-dreiction + Gpu::Atomic::Add(&d_freestream_velocity_ptr[turb_index],std::pow(u_vel(i,j,k)*u_vel(i,j,k) + v_vel(i,j,k)*v_vel(i,j,k),0.5)); + Gpu::Atomic::Add(&d_disk_cell_count_ptr[turb_index],1.0); + Gpu::Atomic::Add(&d_freestream_phi_ptr[turb_index],phi); + } + }); + } + + // Copy back to host + Gpu::copy(Gpu::deviceToHost, d_freestream_velocity.begin(), d_freestream_velocity.end(), freestream_velocity.begin()); + Gpu::copy(Gpu::deviceToHost, d_freestream_phi.begin(), d_freestream_phi.end(), freestream_phi.begin()); + Gpu::copy(Gpu::deviceToHost, d_disk_cell_count.begin(), d_disk_cell_count.end(), disk_cell_count.begin()); + + // Reduce the data on every processor + amrex::ParallelAllReduce::Sum(freestream_velocity.data(), + freestream_velocity.size(), + amrex::ParallelContext::CommunicatorAll()); + + amrex::ParallelAllReduce::Sum(freestream_phi.data(), + freestream_phi.size(), + amrex::ParallelContext::CommunicatorAll()); + + + amrex::ParallelAllReduce::Sum(disk_cell_count.data(), + disk_cell_count.size(), + amrex::ParallelContext::CommunicatorAll()); + + get_turb_loc(xloc, yloc); + + + if (ParallelDescriptor::IOProcessor()){ + for(int it=0; it rad) { + index = i; + break; + } + } + } + if(index == -1 and rad > bld_rad_loc[n_bld_sections-1]) { + index = n_bld_sections-1; + } + if(index == -1) { + //printf("The radial section is at %0.15g m\n",rad); + Abort("Could not find index of the radial section."); + } + + return index; +} + +AMREX_FORCE_INLINE +AMREX_GPU_DEVICE +std::array +compute_source_terms_Fn_Ft (const Real rad, + const Real avg_vel, + const Real* bld_rad_loc, + const Real* bld_twist, + const Real* bld_chord, + int n_bld_sections, + const Real* bld_airfoil_aoa, + const Real* bld_airfoil_Cl, + const Real* bld_airfoil_Cd, + const int n_pts_airfoil, + const Real* velocity, + const Real* rotor_RPM, + const Real* blade_pitch, + const int n_spec_extra) +{ + + Real rpm = interpolate_1d(velocity, rotor_RPM, avg_vel, n_spec_extra); + Real pitch = interpolate_1d(velocity, blade_pitch, avg_vel, n_spec_extra); + + Real Omega = rpm/60.0*2.0*PI; + Real rho = 1.226; + + Real B = 3.0; + Real rhub = 2.0; + Real rtip = 63.5; + + Real twist = interpolate_1d(bld_rad_loc, bld_twist, rad, n_bld_sections); + Real c = interpolate_1d(bld_rad_loc, bld_chord, rad, n_bld_sections); + + // Iteration procedure + + Real s = 0.5*c*B/(PI*rad); + + Real at, an, V1, Vt, Vr, psi, L, D, Cn, Ct; + Real ftip, fhub, F, Cl, Cd, at_new, an_new; + + at = 0.1; + an = 0.1; + + bool is_converged = false; + + for(int i=0;i<100;i++) { + V1 = avg_vel*(1-an); + Vt = Omega*(1.0+at)*rad; + Vr = std::pow(V1*V1+Vt*Vt,0.5); + + psi = std::atan2(V1,Vt); + + Real aoa = psi*180.0/PI - twist + pitch; + + Cl = interpolate_1d(bld_airfoil_aoa, bld_airfoil_Cl, aoa, n_pts_airfoil); + Cd = interpolate_1d(bld_airfoil_aoa, bld_airfoil_Cd, aoa, n_pts_airfoil); + + //Cl = 1.37; + //Cd = 0.014; + + //printf("rad, aoa, Cl, Cd = %0.15g %0.15g %0.15g %0.15g\n", rad, aoa, Cl, Cd); + + Cn = Cl*std::cos(psi) + Cd*std::sin(psi); + Ct = Cl*std::sin(psi) - Cd*std::cos(psi); + + ftip = B*(rtip-rad)/(2.0*rad*std::sin(psi)+1e-10); + fhub = B*(rad-rhub)/(2.0*rad*std::sin(psi)+1e-10); + + AMREX_ALWAYS_ASSERT(std::fabs(std::exp(-fhub))<=1.0); + AMREX_ALWAYS_ASSERT(std::fabs(std::exp(-ftip))<=1.0); + + F = 1.0;//2.0/PI*(std::acos(std::exp(-ftip)) + std::acos(std::exp(-fhub)) ); + + at_new = 1.0/ ( 4.0*F*std::sin(psi)*std::cos(psi)/(s*Ct+1e-10) - 1.0 ); + an_new = 1.0/ ( 1.0 + 4.0*F*std::pow(std::sin(psi),2)/(s*Cn + 1e-10) ); + at_new = std::max(0.0, at_new); + + if(std::fabs(at_new-at) < 1e-5 and std::fabs(an_new-an) < 1e-5) { + //printf("Converged at, an = %d %0.15g %0.15g %0.15g\n",i, at, an, psi); + at = at_new; + an = an_new; + is_converged = true; + break; + } + at = at_new; + an = an_new; + //printf("Iteration, at, an = %0.15g %0.15g %0.15g\n",at, an, psi); + } + + if(!is_converged) { + Abort("The iteration procedure for the generalized actuator disk did not converge. Exiting..."); + } + + // Iterations converged. Now compute Ft, Fn + + L = 0.5*rho*Vr*Vr*c*Cl; + D = 0.5*rho*Vr*Vr*c*Cd; + + Real Fn = L*std::cos(psi) + D*std::sin(psi); + Real Ft = L*std::sin(psi) - D*std::cos(psi); + + //printf("Fn and Ft %0.15g %0.15g %0.15g %0.15g\n", L, D, std::cos(psi), std::sin(psi)); + + std::array Fn_and_Ft; + Fn_and_Ft[0] = Fn; + Fn_and_Ft[1] = Ft; + + return Fn_and_Ft; + + //exit(0); +} + + void GeneralAD::source_terms_cellcentered (const Geometry& geom, const MultiFab& cons_in, @@ -58,9 +289,19 @@ GeneralAD::source_terms_cellcentered (const Geometry& geom, { get_turb_loc(xloc, yloc); + get_turb_spec(rotor_rad, hub_height, thrust_coeff_standing, wind_speed, thrust_coeff, power); + get_blade_spec(bld_rad_loc,bld_twist,bld_chord); + + get_blade_airfoil_spec(bld_airfoil_aoa, bld_airfoil_Cl, bld_airfoil_Cd); + + get_turb_spec_extra(velocity, C_P, C_T, rotor_RPM, blade_pitch); + + Real d_hub_height = hub_height; + Real d_rotor_rad = rotor_rad; + Gpu::DeviceVector d_xloc(xloc.size()); Gpu::DeviceVector d_yloc(yloc.size()); Gpu::copy(Gpu::hostToDevice, xloc.begin(), xloc.end(), d_xloc.begin()); @@ -70,6 +311,7 @@ GeneralAD::source_terms_cellcentered (const Geometry& geom, // Domain valid box const amrex::Box& domain = geom.Domain(); + auto ProbLoArr = geom.ProbLoArray(); int domlo_x = domain.smallEnd(0); int domhi_x = domain.bigEnd(0) + 1; int domlo_y = domain.smallEnd(1); @@ -80,23 +322,85 @@ GeneralAD::source_terms_cellcentered (const Geometry& geom, // The order of variables are - Vabs dVabsdt, dudt, dvdt, dTKEdt mf_vars_generalAD.setVal(0.0); - long unsigned int nturbs = xloc.size(); + long unsigned int nturbs = xloc.size(); + // This is the angle phi in Fig. 10 in Mirocha et. al. 2014 + // set_turb_disk angle in ERF_InitWindFarm.cpp sets this phi as + // the turb_disk_angle get_turb_disk_angle(turb_disk_angle); - Real nx = -std::cos(turb_disk_angle); - Real ny = -std::sin(turb_disk_angle); Real d_turb_disk_angle = turb_disk_angle; - Gpu::DeviceVector d_wind_speed(wind_speed.size()); - Gpu::DeviceVector d_thrust_coeff(thrust_coeff.size()); + Gpu::DeviceVector d_freestream_velocity(nturbs); + Gpu::DeviceVector d_disk_cell_count(nturbs); + Gpu::copy(Gpu::hostToDevice, freestream_velocity.begin(), freestream_velocity.end(), d_freestream_velocity.begin()); + Gpu::copy(Gpu::hostToDevice, disk_cell_count.begin(), disk_cell_count.end(), d_disk_cell_count.begin()); + + Real* d_xloc_ptr = d_xloc.data(); + Real* d_yloc_ptr = d_yloc.data(); + Real* d_freestream_velocity_ptr = d_freestream_velocity.data(); + Real* d_disk_cell_count_ptr = d_disk_cell_count.data(); + + int n_bld_sections = bld_rad_loc.size(); + + Gpu::DeviceVector d_bld_rad_loc(n_bld_sections); + Gpu::DeviceVector d_bld_twist(n_bld_sections); + Gpu::DeviceVector d_bld_chord(n_bld_sections); + + Gpu::copy(Gpu::hostToDevice, bld_rad_loc.begin(), bld_rad_loc.end(), d_bld_rad_loc.begin()); + Gpu::copy(Gpu::hostToDevice, bld_twist.begin(), bld_twist.end(), d_bld_twist.begin()); + Gpu::copy(Gpu::hostToDevice, bld_chord.begin(), bld_chord.end(), d_bld_chord.begin()); - // Copy data from host vectors to device vectors - Gpu::copy(Gpu::hostToDevice, wind_speed.begin(), wind_speed.end(), d_wind_speed.begin()); - Gpu::copy(Gpu::hostToDevice, thrust_coeff.begin(), thrust_coeff.end(), d_thrust_coeff.begin()); + Real* bld_rad_loc_ptr = d_bld_rad_loc.data(); + Real* bld_twist_ptr = d_bld_twist.data(); + Real* bld_chord_ptr = d_bld_chord.data(); - const Real* wind_speed_d = d_wind_speed.dataPtr(); - const Real* thrust_coeff_d = d_thrust_coeff.dataPtr(); - const int n_spec_table = d_wind_speed.size(); + Vector> d_bld_airfoil_aoa(n_bld_sections); + Vector> d_bld_airfoil_Cl(n_bld_sections); + Vector> d_bld_airfoil_Cd(n_bld_sections); + + int n_pts_airfoil = bld_airfoil_aoa[0].size(); + + for(int i=0;i hp_bld_airfoil_aoa, hp_bld_airfoil_Cl, hp_bld_airfoil_Cd; + for (auto & v :d_bld_airfoil_aoa) { + hp_bld_airfoil_aoa.push_back(v.data()); + } + for (auto & v :d_bld_airfoil_Cl) { + hp_bld_airfoil_Cl.push_back(v.data()); + } + for (auto & v :d_bld_airfoil_Cd) { + hp_bld_airfoil_Cd.push_back(v.data()); + } + + Gpu::AsyncArray aoa(hp_bld_airfoil_aoa.data(), n_bld_sections); + Gpu::AsyncArray Cl(hp_bld_airfoil_Cl.data(), n_bld_sections); + Gpu::AsyncArray Cd(hp_bld_airfoil_Cd.data(), n_bld_sections); + + auto d_bld_airfoil_aoa_ptr = aoa.data(); + auto d_bld_airfoil_Cl_ptr = Cl.data(); + auto d_bld_airfoil_Cd_ptr = Cd.data(); + + int n_spec_extra = velocity.size(); + + Gpu::DeviceVector d_velocity(n_spec_extra); + Gpu::DeviceVector d_rotor_RPM(n_spec_extra); + Gpu::DeviceVector d_blade_pitch(n_spec_extra); + + Gpu::copy(Gpu::hostToDevice, velocity.begin(), velocity.end(), d_velocity.begin()); + Gpu::copy(Gpu::hostToDevice, rotor_RPM.begin(), rotor_RPM.end(), d_rotor_RPM.begin()); + Gpu::copy(Gpu::hostToDevice, blade_pitch.begin(), blade_pitch.end(), d_blade_pitch.begin()); + + auto d_velocity_ptr = d_velocity.data(); + auto d_rotor_RPM_ptr = d_rotor_RPM.data(); + auto d_blade_pitch_ptr = d_blade_pitch.data(); for ( MFIter mfi(cons_in,TilingIfNotGPU()); mfi.isValid(); ++mfi) { @@ -109,34 +413,79 @@ GeneralAD::source_terms_cellcentered (const Geometry& geom, int jj = amrex::min(amrex::max(j, domlo_y), domhi_y); int kk = amrex::min(amrex::max(k, domlo_z), domhi_z); + Real x = ProbLoArr[0] + (ii+0.5)*dx[0]; + Real y = ProbLoArr[1] + (jj+0.5)*dx[1]; + Real z = ProbLoArr[2] + (kk+0.5)*dx[2]; + // ?? Density needed here + Real inv_dens_vol = 1.0/(1.0*dx[0]*dx[1]*dx[2]); + int check_int = 0; - Real source_x = 0.0; - Real source_y = 0.0; + Real source_x = 0.0, source_y = 0.0, source_z = 0.0; + std::array Fn_and_Ft; for(long unsigned int it=0;it(it)) { check_int++; - if(C_T <= 1) { - source_x = -2.0*std::pow(Uinfty_dot_nhat, 2.0)*a*(1.0-a)*dx[1]*dx[2]*std::cos(d_turb_disk_angle)/(dx[0]*dx[1]*dx[2])*std::cos(phi); - source_y = -2.0*std::pow(Uinfty_dot_nhat, 2.0)*a*(1.0-a)*dx[1]*dx[2]*std::cos(d_turb_disk_angle)/(dx[0]*dx[1]*dx[2])*std::sin(phi); - } - else { - source_x = -0.5*C_T*std::pow(Uinfty_dot_nhat, 2.0)*dx[1]*dx[2]*std::cos(d_turb_disk_angle)/(dx[0]*dx[1]*dx[2])*std::cos(phi); - source_y = -0.5*C_T*std::pow(Uinfty_dot_nhat, 2.0)*dx[1]*dx[2]*std::cos(d_turb_disk_angle)/(dx[0]*dx[1]*dx[2])*std::sin(phi); + + // Find radial distance of the point and the zeta angle + Real rad = std::pow( (x-d_xloc_ptr[it])*(x-d_xloc_ptr[it]) + + (y-d_yloc_ptr[it])*(y-d_yloc_ptr[it]) + + (z-d_hub_height)*(z-d_hub_height), 0.5 ); + + int index = find_rad_loc_index(rad, bld_rad_loc_ptr, n_bld_sections); + + // This if check makes sure it is a point with radial distance + // between the hub radius and the rotor radius. + // ?? hub radius needed here + if(rad >= 2.0 and rad <= d_rotor_rad) { + //AMREX_ASSERT( (z-d_hub_height) <= rad ); + // Consider the vector that joines the point and the turbine center. + // Dot it on to the vector that joins the turbine center and along + // the plane of the disk. See fig. 10 in Mirocha et. al. 2014. + + Real vec_proj = (x-d_xloc_ptr[it])*(std::sin(phi)) + + (y-d_yloc_ptr[it])*(-std::cos(phi)); + + + Real zeta = std::atan2(z-d_hub_height, vec_proj); + //printf("zeta val is %0.15g\n", zeta*180.0/PI); + Fn_and_Ft = compute_source_terms_Fn_Ft(rad, avg_vel, + bld_rad_loc_ptr, + bld_twist_ptr, + bld_chord_ptr, + n_bld_sections, + d_bld_airfoil_aoa_ptr[index], + d_bld_airfoil_Cl_ptr[index], + d_bld_airfoil_Cd_ptr[index], + n_pts_airfoil, + d_velocity_ptr, + d_rotor_RPM_ptr, + d_blade_pitch_ptr, + n_spec_extra); + + Real Fn = Fn_and_Ft[0]; + Real Ft = Fn_and_Ft[1]; + // Compute the source terms - pass in radial distance, free stream velocity + + Real Fx = Fn*std::cos(phi) + Ft*std::sin(zeta)*std::sin(phi); + Real Fy = Fn*std::sin(phi) - Ft*std::sin(zeta)*std::cos(phi); + Real Fz = -Ft*std::cos(zeta); + + source_x = -Fx*inv_dens_vol; + source_y = -Fy*inv_dens_vol; + source_z = -Fz*inv_dens_vol; + + + //printf("Val source_x, is %0.15g, %0.15g, %0.15g %0.15g %0.15g %0.15g\n", rad, Fn, Ft, source_x, source_y, source_z); } } } + if(check_int > 1){ amrex::Error("Actuator disks are overlapping. Visualize actuator_disks.vtk " "and check the windturbine locations input file. Exiting.."); @@ -144,6 +493,7 @@ GeneralAD::source_terms_cellcentered (const Geometry& geom, generalAD_array(i,j,k,0) = source_x; generalAD_array(i,j,k,1) = source_y; + generalAD_array(i,j,k,2) = source_z; }); } } diff --git a/Source/WindFarmParametrization/GeneralActuatorDisk/ERF_GeneralAD.H b/Source/WindFarmParametrization/GeneralActuatorDisk/ERF_GeneralAD.H index f3ab40227..593b16a48 100644 --- a/Source/WindFarmParametrization/GeneralActuatorDisk/ERF_GeneralAD.H +++ b/Source/WindFarmParametrization/GeneralActuatorDisk/ERF_GeneralAD.H @@ -23,6 +23,11 @@ public: const amrex::MultiFab& mf_Nturb, const amrex::MultiFab& mf_SMark) override; + void compute_freestream_velocity (const amrex::MultiFab& cons_in, + const amrex::MultiFab& U_old, + const amrex::MultiFab& V_old, + const amrex::MultiFab& mf_SMark); + void source_terms_cellcentered (const amrex::Geometry& geom, const amrex::MultiFab& cons_in, const amrex::MultiFab& mf_Smark, @@ -32,6 +37,7 @@ public: amrex::MultiFab& cons_in, amrex::MultiFab& U_old, amrex::MultiFab& V_old, + amrex::MultiFab& W_old, const amrex::MultiFab& mf_vars); protected: @@ -40,6 +46,9 @@ protected: amrex::Real hub_height, rotor_rad, thrust_coeff_standing, nominal_power; amrex::Vector wind_speed, thrust_coeff, power; amrex::Vector freestream_velocity, freestream_phi, disk_cell_count; + amrex::Vector bld_rad_loc, bld_twist, bld_chord; + amrex::Vector> bld_airfoil_aoa, bld_airfoil_Cl, bld_airfoil_Cd; + amrex::Vector velocity, C_P, C_T, rotor_RPM, blade_pitch; }; #endif diff --git a/Source/WindFarmParametrization/Null/ERF_NullWindFarm.H b/Source/WindFarmParametrization/Null/ERF_NullWindFarm.H index 934832b19..2daea0aa5 100644 --- a/Source/WindFarmParametrization/Null/ERF_NullWindFarm.H +++ b/Source/WindFarmParametrization/Null/ERF_NullWindFarm.H @@ -49,6 +49,37 @@ public: m_turb_disk_angle = turb_disk_angle; } + virtual void set_blade_spec(const amrex::Vector& bld_rad_loc, + const amrex::Vector& bld_twist, + const amrex::Vector& bld_chord) + { + m_bld_rad_loc = bld_rad_loc; + m_bld_twist = bld_twist; + m_bld_chord = bld_chord; + } + + virtual void set_blade_airfoil_spec(const amrex::Vector>& bld_airfoil_aoa, + const amrex::Vector>& bld_airfoil_Cl, + const amrex::Vector>& bld_airfoil_Cd) + { + m_bld_airfoil_aoa = bld_airfoil_aoa; + m_bld_airfoil_Cl = bld_airfoil_Cl; + m_bld_airfoil_Cd = bld_airfoil_Cd; + } + + virtual void set_turb_spec_extra(const amrex::Vector& velocity, + const amrex::Vector& C_P, + const amrex::Vector& C_T, + const amrex::Vector& rotor_RPM, + const amrex::Vector& blade_pitch) + { + m_velocity = velocity; + m_C_P = C_P; + m_C_T = C_T; + m_rotor_RPM = rotor_RPM; + m_blade_pitch = blade_pitch; + } + void get_turb_spec (amrex::Real& rotor_rad, amrex::Real& hub_height, amrex::Real& thrust_coeff_standing, amrex::Vector& wind_speed, amrex::Vector& thrust_coeff, amrex::Vector& power) @@ -73,6 +104,36 @@ public: turb_disk_angle = m_turb_disk_angle; } + void get_blade_spec(amrex::Vector& bld_rad_loc, + amrex::Vector& bld_twist, + amrex::Vector& bld_chord) + { + bld_rad_loc = m_bld_rad_loc; + bld_twist = m_bld_twist; + bld_chord = m_bld_chord; + } + + void get_blade_airfoil_spec(amrex::Vector>& bld_airfoil_aoa, + amrex::Vector>& bld_airfoil_Cl, + amrex::Vector>& bld_airfoil_Cd) + { + bld_airfoil_aoa = m_bld_airfoil_aoa; + bld_airfoil_Cl = m_bld_airfoil_Cl; + bld_airfoil_Cd = m_bld_airfoil_Cd; + } + + void get_turb_spec_extra(amrex::Vector& velocity, + amrex::Vector& C_P, + amrex::Vector& C_T, + amrex::Vector& rotor_RPM, + amrex::Vector& blade_pitch) + { + velocity = m_velocity; + C_P = m_C_P; + C_T = m_C_T; + rotor_RPM = m_rotor_RPM; + blade_pitch = m_blade_pitch; + } static AMREX_GPU_DEVICE bool find_if_marked(amrex::Real x1, amrex::Real x2, amrex::Real y1, amrex::Real y2, @@ -121,6 +182,9 @@ protected: amrex::Real m_turb_disk_angle; amrex::Real m_hub_height, m_rotor_rad, m_thrust_coeff_standing, m_nominal_power; amrex::Vector m_wind_speed, m_thrust_coeff, m_power; + amrex::Vector m_bld_rad_loc, m_bld_twist, m_bld_chord; + amrex::Vector> m_bld_airfoil_aoa, m_bld_airfoil_Cl, m_bld_airfoil_Cd; + amrex::Vector m_velocity, m_C_P, m_C_T, m_rotor_RPM, m_blade_pitch; };