Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Correction to generalized actuator disk model #1912

Merged
merged 4 commits into from
Oct 28, 2024
Merged

Correction to generalized actuator disk model #1912

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
6e191b0
Changing actuator disk model names
Oct 22, 2024
820d4f1
Merge remote-tracking branch 'upstream/development' into AWAKEN_Refactor
Oct 28, 2024
7068523
Some changes to actuator disk model
Oct 28, 2024
c4d4c2e
Revmoing unused variable
Oct 28, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
20 changes: 11 additions & 9 deletions Source/WindFarmParametrization/ERF_InitWindFarm.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -125,14 +125,6 @@ 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
Expand Down Expand Up @@ -456,6 +448,8 @@ WindFarm::fill_SMark_multifab(const Geometry& geom,

Real z = ProbLoArr[2] + (kk+0.5) * dx[2];

int turb_indices_overlap[2];
int check_int = 0;
for(int it=0; it<num_turb; it++){
Real x0 = d_xloc_ptr[it] + d_sampling_distance*nx;
Real y0 = d_yloc_ptr[it] + d_sampling_distance*ny;
Expand All @@ -467,13 +461,21 @@ WindFarm::fill_SMark_multifab(const Geometry& geom,
}
x0 = d_xloc_ptr[it];
y0 = d_yloc_ptr[it];
//printf("Values are %d, %0.15g, %0.15g\n", it, x0, y0);

is_cell_marked = find_if_marked(x1, x2, y1, y2, x0, y0,
nx, ny, d_hub_height, d_rotor_rad, z);
if(is_cell_marked) {
SMark_array(i,j,k,1) = it;
turb_indices_overlap[check_int] = it;
check_int++;
if(check_int > 1){
printf("Actuator disks with indices %d and %d are overlapping\n",
turb_indices_overlap[0],turb_indices_overlap[1]);
amrex::Error("Actuator disks are overlapping. Visualize actuator_disks.vtk "
" and check the windturbine locations input file. Exiting..");
}
}

}
});
}
Expand Down
15 changes: 8 additions & 7 deletions Source/WindFarmParametrization/GeneralActuatorDisk/ERF_AdvanceGeneralAD.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -239,7 +239,7 @@ compute_source_terms_Fn_Ft (const Real rad,
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)) );
F = 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) );
Expand Down Expand Up @@ -417,7 +417,6 @@ GeneralAD::source_terms_cellcentered (const Geometry& geom,
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;

Expand Down Expand Up @@ -468,17 +467,19 @@ GeneralAD::source_terms_cellcentered (const Geometry& geom,
d_blade_pitch_ptr,
n_spec_extra);

Real Fn = Fn_and_Ft[0];
Real Ft = Fn_and_Ft[1];
Real Fn = 3.0*Fn_and_Ft[0];
Real Ft = 3.0*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;
//Real dn = (

source_x = -Fx/(2.0*PI*rad*dx[0])*1.0/std::pow(2.0*PI,0.5);
source_y = -Fy/(2.0*PI*rad*dx[0])*1.0/std::pow(2.0*PI,0.5);
source_z = -Fz/(2.0*PI*rad*dx[0])*1.0/std::pow(2.0*PI,0.5);


//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);
Expand Down
18 changes: 6 additions & 12 deletions Source/WindFarmParametrization/SimpleActuatorDisk/ERF_AdvanceSimpleAD.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -114,15 +114,15 @@ void SimpleAD::compute_freestream_velocity(const MultiFab& cons_in,
amrex::ParallelContext::CommunicatorAll());

get_turb_loc(xloc, yloc);
if (ParallelDescriptor::IOProcessor()){
/*if (ParallelDescriptor::IOProcessor()){
for(int it=0; it<xloc.size(); it++){

std::cout << "turbine index, freestream velocity is " << it << " " << freestream_velocity[it] << " " <<
disk_cell_count[it] << " " <<
freestream_velocity[it]/(disk_cell_count[it] + 1e-10) << " " <<
freestream_phi[it]/(disk_cell_count[it] + 1e-10) << "\n";
}
}
}*/
}

void
Expand Down Expand Up @@ -195,12 +195,13 @@ SimpleAD::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);

int check_int = 0;

Real source_x = 0.0;
Real source_y = 0.0;

for(long unsigned int it=0;it<nturbs;it++) {
int it = static_cast<int>(SMark_array(ii,jj,kk,1));

if(it != -1) {
Real avg_vel = d_freestream_velocity_ptr[it]/(d_disk_cell_count_ptr[it] + 1e-10);
Real phi = d_freestream_phi_ptr[it]/(d_disk_cell_count_ptr[it] + 1e-10);

Expand All @@ -210,8 +211,6 @@ SimpleAD::source_terms_cellcentered (const Geometry& geom,
a = 0.5 - 0.5*std::pow(1.0-C_T,0.5);
}
Real Uinfty_dot_nhat = avg_vel*(std::cos(phi)*nx + std::sin(phi)*ny);
if(SMark_array(ii,jj,kk,1) == static_cast<double>(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);
Expand All @@ -220,12 +219,7 @@ SimpleAD::source_terms_cellcentered (const Geometry& geom,
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);
}
}
}
if(check_int > 1){
amrex::Error("Actuator disks are overlapping. Visualize actuator_disks.vtk "
"and check the windturbine locations input file. Exiting..");
}
}

simpleAD_array(i,j,k,0) = source_x;
simpleAD_array(i,j,k,1) = source_y;
Expand Down
Loading