diff --git a/windse/ProblemManager.py b/windse/ProblemManager.py index 0fb20eca..35953efb 100644 --- a/windse/ProblemManager.py +++ b/windse/ProblemManager.py @@ -594,9 +594,10 @@ def ComputeFunctional(self,inflow_angle): # Define fluid properties # FIXME: These should probably be set in params.yaml input filt # nu = 1/10000 - rho = 1 + rho = self.params["problem"]["density"] # User-defined change nu_c = Constant(self.viscosity, name="viscosity") rho_c = Constant(rho, name="rho") + MISSING_RHO_C = rho_c # Define time step size (this value is used only for step 1 if adaptive timestepping is used) # FIXME: change variable name to avoid confusion within dolfin adjoint @@ -623,7 +624,7 @@ def ComputeFunctional(self,inflow_angle): self.u_k1 = Function(self.fs.V, name="u_k1") self.u_k2 = Function(self.fs.V, name="u_k2") - # Seed previous velocity fields with the chosen initial condition + # Seed previous velocity fields with the chosen initial condition #DEBUGING: Following lines are commented so that velocity isn't assigned to all nodes self.u_k.assign(self.bd.bc_velocity) self.u_k1.assign(self.bd.bc_velocity) self.u_k2.assign(self.bd.bc_velocity) @@ -694,12 +695,19 @@ def ComputeFunctional(self,inflow_angle): # + (nu_c+self.nu_T)*inner(grad(U_CN), grad(v))*dx \ # + dot(nabla_grad(self.p_k1), v)*dx \ # - dot(self.tf, v)*dx + # Update: density included with pressure term of NS + #F1 = (1.0/self.dt_c)*inner(u - self.u_k1, v)*dx \ + # + inner(dot(U_AB, nabla_grad(U_CN)), v)*dx \ + # + (nu_c+self.nu_T)*inner(grad(U_CN), grad(v))*dx \ + # + inner(grad(self.p_k1), v)*dx \ + # - dot(self.tf, v)*dx F1 = (1.0/self.dt_c)*inner(u - self.u_k1, v)*dx \ + inner(dot(U_AB, nabla_grad(U_CN)), v)*dx \ + (nu_c+self.nu_T)*inner(grad(U_CN), grad(v))*dx \ - + inner(grad(self.p_k1), v)*dx \ - - dot(self.tf, v)*dx + + (1.0/MISSING_RHO_C)*inner(grad(self.p_k1), v)*dx \ + - (1.0/MISSING_RHO_C)*dot(self.tf, v)*dx + self.a1 = lhs(F1) self.L1 = rhs(F1) @@ -708,7 +716,8 @@ def ComputeFunctional(self,inflow_angle): # self.a2 = dot(nabla_grad(p), nabla_grad(q))*dx # self.L2 = dot(nabla_grad(self.p_k1), nabla_grad(q))*dx - (1.0/self.dt_c)*div(self.u_k)*q*dx self.a2 = inner(grad(p), grad(q))*dx - self.L2 = inner(grad(self.p_k1), grad(q))*dx - (1.0/self.dt_c)*div(self.u_k)*q*dx + #self.L2 = inner(grad(self.p_k1), grad(q))*dx - (1/self.dt_c)*div(self.u_k)*q*dx + self.L2 = inner(grad(self.p_k1), grad(q))*dx - (MISSING_RHO_C/self.dt_c)*div(self.u_k)*q*dx # phi = p - self.p_k # F2 = inner(grad(q), grad(phi))*dx - (1.0/self.dt_c)*div(u_k)*q*dx @@ -719,7 +728,8 @@ def ComputeFunctional(self,inflow_angle): # self.a3 = dot(u, v)*dx # self.L3 = dot(self.u_k, v)*dx - self.dt_c*dot(nabla_grad(self.p_k - self.p_k1), v)*dx self.a3 = inner(u, v)*dx - self.L3 = inner(self.u_k, v)*dx - self.dt_c*inner(grad(self.p_k - self.p_k1), v)*dx + #self.L3 = inner(self.u_k, v)*dx - self.dt_c*inner(grad(self.p_k - self.p_k1), v)*dx + self.L3 = inner(self.u_k, v)*dx - (self.dt_c/MISSING_RHO_C)*inner(grad(self.p_k - self.p_k1), v)*dx # F3 = inner(u, v)*dx - inner(self.u_k, v)*dx + self.dt_c*inner(phi, v)*dx # self.a3 = lhs(F3) @@ -744,4 +754,3 @@ def ComputeFunctional(self,inflow_angle): # self.mcd[k] = Constant(cd[k]) # # ================================================================ - diff --git a/windse/SolverManager.py b/windse/SolverManager.py index 3a9a24c9..cb5cd9ef 100644 --- a/windse/SolverManager.py +++ b/windse/SolverManager.py @@ -1467,13 +1467,13 @@ def rot_z(theta): #================================================================ # Set the density - rho = 1.0 + rho = self.params["problem"]["density"] # Set the hub height hub_height = self.problem.farm.HH[0] # For a SWIFT turbine # Get the hub-height velocity - u_inf = 8.0 + u_inf = 8.0 # DEBUGGING: This variable is set to a constant number. Is this right? # Set the rotational speed of the turbine RPM = 15.0 @@ -1482,7 +1482,7 @@ def rot_z(theta): yaw = self.problem.farm.yaw[0] # Set the number of blades in the turbine - num_blades = 3 + num_blades = self.params["turbines"]["num_blades"] # Blade length (turbine radius) L = self.problem.farm.radius[0] # For a SWIFT turbine 27 m in diameter @@ -1502,7 +1502,7 @@ def rot_z(theta): # print(self.problem.farm.radius) # Discretize each blade into separate nodes - num_actuator_nodes = 10 + num_actuator_nodes = 10 # self.params["turbines"]["blade_segments"] # DEBUGGING: This variable is updated to be a user-defined parameter #================================================================ # Set Derived Constants @@ -1666,7 +1666,7 @@ def UpdateActuatorLineForceOld(self, simTime): coords = np.copy(coords[0::self.problem.dom.dim, :]) # Set up the turbine geometry - num_blades = 3 + num_blades = self.params["turbines"]["num_blades"] theta_vec = np.linspace(0, 2.0*np.pi, num_blades+1) theta_vec = theta_vec[0:num_blades] @@ -1676,7 +1676,7 @@ def UpdateActuatorLineForceOld(self, simTime): cl = 1.0 cd = 2.0 - rho = 1 + rho = self.params["problem"]["density"] u_inf = 8 @@ -1687,7 +1687,7 @@ def UpdateActuatorLineForceOld(self, simTime): c = L/20 # Number of blade evaluation sections - num_actuator_nodes = 50 + num_actuator_nodes = 50 # self.params["turbines"]["blade_segments"] # DEBUGGING: Hard-coded parameter changed to user-defined parameter rdim = np.linspace(0, L, num_actuator_nodes) zdim = 0.0 + np.zeros(num_actuator_nodes) xblade = np.vstack((np.zeros(num_actuator_nodes), rdim, zdim)) diff --git a/windse/input_schema.yaml b/windse/input_schema.yaml index 9615376b..2157fc7f 100644 --- a/windse/input_schema.yaml +++ b/windse/input_schema.yaml @@ -563,6 +563,11 @@ properties: type: "number" description: "The radius of the hub. If non-zero, actuator nodes will still be placed in the full range [0, rotor_radius], but the lift/drag properties in the range [0, hub_rad] will be modified to reflect the blade root." units: "meter" + num_blades: + default: 3 + type: "integer" + description: "Integer number of blades for each turbine" + units: "dimensionless" chord_perturb: default: 0.0 type: "number" @@ -881,6 +886,11 @@ properties: type: "number" description: "Kinematic viscosity." units: "meter^2/second" + density: + default: 1.0 + type: "number" + description: "Fluid density" + units: "kg/meter^3" lmax: default: 15 type: "number" diff --git a/windse/turbine_types/ActuatorDisk.py b/windse/turbine_types/ActuatorDisk.py index d397119b..ca9ff21d 100644 --- a/windse/turbine_types/ActuatorDisk.py +++ b/windse/turbine_types/ActuatorDisk.py @@ -202,7 +202,7 @@ def _setup_chord_force(self): if self.mchord is not None: if self.chord is not None: if self.blade_segments == "computed": - self.num_actuator_nodes = 10 ##### FIX THIS #### + self.num_actuator_nodes = 10 ##### FIX THIS #### # DEBUGGING: self.num_actuator_nodes = int(2.0*self.radius/self.gaussian_width) # FromLine 150 of ActuatorLine.py self.blade_segments = self.num_actuator_nodes else: self.num_actuator_nodes = self.blade_segments diff --git a/windse/turbine_types/ActuatorLine.py b/windse/turbine_types/ActuatorLine.py index 2c48a1df..3b04fb34 100644 --- a/windse/turbine_types/ActuatorLine.py +++ b/windse/turbine_types/ActuatorLine.py @@ -133,7 +133,7 @@ def init_constant_alm_terms(self, fs): #================================================================ # Set the number of blades in the turbine - self.num_blades = 3 + self.num_blades = self.params["turbines"]["num_blades"] # Set the spacing pf each blade self.theta_0_vec = np.linspace(0.0, 2.0*np.pi, self.num_blades, endpoint = False) @@ -548,7 +548,7 @@ def get_angle_between_vectors(a, b, n): tip_loss_fac[k] = 1.0 else: - loss_exponent = 3.0/2.0*(self.radius-self.rdim[k])/(self.rdim[k]*np.sin(aoa)) + loss_exponent = self.num_blades/2.0*(self.radius-self.rdim[k])/(self.rdim[k]*np.sin(aoa)) #DEBUGGING: 3.0 is changed to self.num_blades for adjusted tip-loss factor acos_arg = np.exp(-loss_exponent) acos_arg = np.clip(acos_arg, -1.0, 1.0) tip_loss_fac[k] = 2.0/np.pi*np.arccos(acos_arg) @@ -676,11 +676,14 @@ def build_actuator_lines(self, u, inflow_angle, fs, dfd=None): else: if self.lookup_cl is not None: cl, cd, tip_loss_fac, aoa = self.lookup_lift_and_drag(u_rel, twist, self.blade_unit_vec[blade_id]) - + # self.fprint("For blade_id = {4}, the following values come from lookup_lift_and_drag: cl = {0}, cd = {1}, tip_loss_fac = {2}, aoa = {3} \n".format(cl,cd,tip_loss_fac,aoa,blade_id)) # Calculate the lift and drag forces using the relative velocity magnitude - rho = 1.0 + rho = self.params["problem"]["density"] lift = tip_loss_fac*(0.5*cl*rho*chord*self.w*u_rel_mag**2) drag = tip_loss_fac*(0.5*cd*rho*chord*self.w*u_rel_mag**2) + + #self.fprint("For blade_id = {0} with rho = {1}, lift = {2}".format(blade_id,rho,lift)) #DEBUGGING + #self.fprint("drag = {}".format(drag)) #DEBUGGING # Tile the blade coordinates for every mesh point, [numGridPts*ndim x problem.num_actuator_nodes] blade_pos_full = np.tile(self.blade_pos[blade_id], (np.shape(self.coords)[0], 1)) @@ -726,19 +729,23 @@ def build_actuator_lines(self, u, inflow_angle, fs, dfd=None): elif dfd == 'chord': for j in range(self.ndim): dfd_chord[j::self.ndim, k] += vector_nodal_lift[:, j] + vector_nodal_drag[:, j] - + #self.fprint("Blade_id = {} \n".format(blade_id)) # DEBUGGING # Compute the total force vector [x, y, z] at a single actuator node actuator_lift = lift[k]*lift_unit_vec + #self.fprint("actuator_lift = {} \n".format(actuator_lift)) #DEBUGGING actuator_drag = drag[k]*drag_unit_vec + #self.fprint("actuator_drag = {} \n".format(actuator_drag)) #DEBUGGING # Note: since this will be used to define the force (torque) from fluid -> blade # we reverse the direction that otherwise gives the turbine force from blade -> fluid actuator_force = -(actuator_lift + actuator_drag) + #self.fprint("actuator_force = {} \n".format(actuator_force)) #DEBUGGING # actuator_force = -(actuator_lift - actuator_drag) # Find the component in the direction tangential to the blade tangential_actuator_force = np.dot(actuator_force, self.blade_unit_vec[blade_id][:, 2]) + #self.fprint("tangential_actuator_force = {} \n".format(tangential_actuator_force)) #DEBUGGING rotor_plane_force = np.dot(actuator_force, self.blade_unit_vec[blade_id]) # fx.write('%.5f, ' % (rotor_plane_force[0])) # fy.write('%.5f, ' % (rotor_plane_force[1])) @@ -767,6 +774,7 @@ def build_actuator_lines(self, u, inflow_angle, fs, dfd=None): data = np.array(along_blade_quantities[save_val]) nn = self.num_blades*self.num_actuator_nodes + # print(data) # DEBUGGING if np.size(data) == nn: data = data.reshape(1, nn) @@ -776,6 +784,7 @@ def build_actuator_lines(self, u, inflow_angle, fs, dfd=None): global_save_data = np.zeros((self.params.num_procs, nn)) self.params.comm.Allgather(data, global_save_data) + data = np.nanmean(data, axis=0) if self.first_call_to_alm: @@ -864,17 +873,24 @@ def build_actuator_lines(self, u, inflow_angle, fs, dfd=None): def finalize_rotor_torque(self): + # self.fprint("num_procs = {} \n".format(self.params.num_procs)) #DEBUGGING data_in_torque = np.zeros(self.params.num_procs) + #self.fprint("data_in_torque = {} \n".format(data_in_torque)) #DEBUGGING self.params.comm.Gather(self.rotor_torque, data_in_torque, root=0) + #self.fprint("Gathered data_in_torque = {} \n".format(data_in_torque)) #DEBUGGING data_in_torque_count = np.zeros(self.params.num_procs, dtype=int) + #self.fprint("data_in_torque_count = {} \n".format(data_in_torque_count)) #DEBUGGING self.params.comm.Gather(self.rotor_torque_count, data_in_torque_count, root=0) + #self.fprint("Gathered data_in_torque_count = {} \n".format(data_in_torque_count)) #DEBUGGING if self.params.rank == 0: rotor_torque_sum = np.sum(data_in_torque) + #self.fprint("rotor_torque_sum = {} \n".format(rotor_torque_sum)) #DEBUG rotor_torque_count_sum = np.sum(data_in_torque_count) - - # This removes the possibility of a power being doubled or tripled + #self.fprint("rotor_torque_count_sum = {} \n".format(rotor_torque_count_sum)) #DEBUG + + # This removes the possibility of a power being doubled or tripled # if multiple ranks include this turbine and therefore calculate a torque self.rotor_torque = rotor_torque_sum/rotor_torque_count_sum @@ -913,6 +929,7 @@ def turbine_force(self, u, inflow_angle, fs, **kwargs): # Call the function to build the complete mpi_u_fluid array self.get_u_fluid_at_alm_nodes(u) + #self.fprint("mpi_u_fluid = {}".format(self.mpi_u_fluid)) # DEBUGGING # Call the ALM function for this turbine if update: diff --git a/windse/turbine_types/ActuatorLineDolfin.py b/windse/turbine_types/ActuatorLineDolfin.py index 22ed57d1..5f99c338 100644 --- a/windse/turbine_types/ActuatorLineDolfin.py +++ b/windse/turbine_types/ActuatorLineDolfin.py @@ -79,6 +79,8 @@ def load_parameters(self): self.motion_type = self.params["turbines"]["motion_type"] self.use_gauss_vel_probe = self.params["turbines"]["use_gauss_vel_probe"] self.use_ap_linear_interp = self.params["turbines"]["use_ap_linear_interp"] + self.num_blades = self.params["turbines"]["num_blades"] + self.density = self.params["problem"]["density"] def compute_parameters(self): pass @@ -172,8 +174,8 @@ def calc_platform_motion(self,T): def init_alm_calc_terms(self): # compute turbine radius - self.rho = 1.0 - self.num_blades = 3 + self.rho = self.density + self.num_blades = self.num_blades self.radius = 0.5*self.RD self.angular_velocity = 2.0*pi*self.rpm/60.0 @@ -594,7 +596,7 @@ def calculate_tip_loss(self, rdim, aoa): tip_loss_fac = 1.0 else: - loss_exponent = 3.0/2.0*(self.radius-rdim)/(rdim*sin(aoa)) + loss_exponent = self.num_blades/2.0*(self.radius-rdim)/(rdim*sin(aoa)) # DEBUGGING: 3.0 is changed to self.num_blades acos_arg = exp(-loss_exponent) # acos_arg = np.clip(acos_arg, -1.0, 1.0) tip_loss_fac = 2.0/pi*acos(acos_arg) @@ -647,6 +649,7 @@ def build_actuator_node(self, u_k, x_0, x_0_pre_motion, n_0, blade_id, actuator_ w = self.w[actuator_id] twist = self.mtwist[actuator_id] chord = self.mchord[actuator_id] + self.rho = self.density # Density is redefined for calculation of lift and drag # Build a spherical Gaussian with width eps about point x_0 self.gauss_kernel[blade_id][actuator_id] = self.build_gauss_kernel_at_point(x_0)