Skip to content

Commit

Permalink
Some changes to actuator disk model
Browse files Browse the repository at this point in the history
  • Loading branch information
Mahesh Natarajan committed Oct 28, 2024
1 parent 820d4f1 commit 7068523
Show file tree
Hide file tree
Showing 3 changed files with 25 additions and 27 deletions.
20 changes: 11 additions & 9 deletions Source/WindFarmParametrization/ERF_InitWindFarm.cpp
View file
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
14 changes: 8 additions & 6 deletions Source/WindFarmParametrization/GeneralActuatorDisk/ERF_AdvanceGeneralAD.cpp
View file
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 @@ -468,17 +468,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
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

0 comments on commit 7068523

Please sign in to comment.