diff --git a/regression/scripts/VTOL/test_Stopped_Rotor.py b/regression/scripts/VTOL/test_Stopped_Rotor.py index 2362bb6350..3b888d0423 100644 --- a/regression/scripts/VTOL/test_Stopped_Rotor.py +++ b/regression/scripts/VTOL/test_Stopped_Rotor.py @@ -74,7 +74,7 @@ def main(): # RPM of rotor check during hover RPM = results.segments.climb_1.conditions.propulsion.lift_rotor_rpm[0][0] - RPM_true = 2403.0053609143974 + RPM_true = 2403.004214209376 print(RPM) diff_RPM = np.abs(RPM - RPM_true) print('RPM difference') @@ -83,7 +83,7 @@ def main(): # Battery Energy Check During Transition battery_energy_hover_to_transition = results.segments.transition_1.conditions.propulsion.battery_energy[:,0] - battery_energy_hover_to_transition_true = np.array([3.36372831e+08, 3.34188966e+08, 3.32108809e+08]) + battery_energy_hover_to_transition_true = np.array([3.36372833e+08, 3.34188971e+08, 3.32108816e+08]) print(battery_energy_hover_to_transition) diff_battery_energy_hover_to_transition = np.abs(battery_energy_hover_to_transition - battery_energy_hover_to_transition_true) @@ -266,8 +266,6 @@ def mission_setup(analyses,vehicle): segment.process.iterate.conditions.stability = SUAVE.Methods.skip segment.process.finalize.post_process.stability = SUAVE.Methods.skip segment = vehicle.networks.lift_cruise.add_transition_unknowns_and_residuals_to_segment(segment, - initial_prop_power_coefficient = 0.2, - initial_lift_rotor_power_coefficient = 0.01, initial_throttle_lift = 0.9,) # add to misison diff --git a/regression/scripts/Vehicles/Tiltrotor.py b/regression/scripts/Vehicles/Tiltrotor.py index c2eea1f879..c6e40350ac 100644 --- a/regression/scripts/Vehicles/Tiltrotor.py +++ b/regression/scripts/Vehicles/Tiltrotor.py @@ -412,6 +412,7 @@ def vehicle_setup(): prop.design_thrust = (1.1*Hover_Load)/net.number_of_propeller_engines prop.origin = [[2.,5.7,0.784]] prop.rotation = -1 + prop.sol_tolerance = 1e-4 prop.symmetry = True prop.variable_pitch = True prop.airfoil_geometry = ['../Vehicles/Airfoils/NACA_4412.txt'] diff --git a/regression/scripts/electric_performance/electric_V_h_diagram.py b/regression/scripts/electric_performance/electric_V_h_diagram.py index 3eed3e109b..d5ca2531bf 100644 --- a/regression/scripts/electric_performance/electric_V_h_diagram.py +++ b/regression/scripts/electric_performance/electric_V_h_diagram.py @@ -78,6 +78,7 @@ def main(): [677.94878779, 589.51079809, 503.24817575, 418.7527992 , 335.57780983]] + assert (np.all(np.nan_to_num(np.abs(climb_rate-climb_rate_r)/climb_rate_r) < 1e-6)), "Electric V_h Diagram Regression Failed" return diff --git a/regression/scripts/electric_performance/propeller_single_point.py b/regression/scripts/electric_performance/propeller_single_point.py index 884dc84a7c..9eaaf1cab4 100644 --- a/regression/scripts/electric_performance/propeller_single_point.py +++ b/regression/scripts/electric_performance/propeller_single_point.py @@ -111,11 +111,11 @@ def test_2(): Cp = results.power_coefficient etap = results.efficiency - thrust_r = 645.8624436654369 - torque_r = 127.10005358627888 - power_r = 29281.750271951667 - Cp_r = 0.037380704154712806 - etap_r = 0.22056825075927722 + thrust_r = 645.8643152096909 + torque_r = 127.10023264665021 + power_r = 29281.791524499793 + Cp_r = 0.03738075681718287 + etap_r = 0.22056857916951647 assert (np.abs(thrust - thrust_r) / thrust_r < 1e-6), "Propeller Single Point Regression Failed at Thrust Test" diff --git a/regression/scripts/propeller/propeller_test.py b/regression/scripts/propeller/propeller_test.py index 1580c23233..a54001ef3d 100644 --- a/regression/scripts/propeller/propeller_test.py +++ b/regression/scripts/propeller/propeller_test.py @@ -205,28 +205,28 @@ def main(): plot_results(outputr, rot,'black','-','P') # Truth values for propeller with airfoil geometry defined - F_a_truth = 3352.3646678384716 - Q_a_truth = 978.76084883 - P_a_truth = 202761.66815935 - Cplast_a_truth = 0.10450829 + F_a_truth = 3352.366469630676 + Q_a_truth = 978.76113592 + P_a_truth = 202761.72763161 + Cplast_a_truth = 0.10450832 # Truth values for propeller without airfoil geometry defined - F_truth = 2629.0134410677906 - Q_truth = 787.38493275 - P_truth = 163115.92626346 - Cplast_truth = 0.08407391 + F_truth = 2629.013537561697 + Q_truth = 787.38469662 + P_truth = 163115.87734548 + Cplast_truth = 0.08407389 # Truth values for rotor with airfoil geometry defined - Fr_a_truth = 1499.6757509190838 - Qr_a_truth = 139.10603263 - Pr_a_truth = 28817.42895722 + Fr_a_truth = 1499.6766372165007 + Qr_a_truth = 139.1060306 + Pr_a_truth = 28817.42853679 Cplastr_a_truth = 0.04532838 # Truth values for rotor without airfoil geometry defined - Fr_truth = 1250.1858330726345 - Qr_truth = 121.95427719 - Pr_truth = 25264.24377481 - Cplastr_truth = 0.0397394 + Fr_truth = 1250.1858821890885 + Qr_truth = 121.95416738 + Pr_truth = 25264.22102656 + Cplastr_truth = 0.03973936 # Store errors error = Data() diff --git a/regression/scripts/segments/transition_segment_test.py b/regression/scripts/segments/transition_segment_test.py index dc7ef1695a..cbd13ccee4 100644 --- a/regression/scripts/segments/transition_segment_test.py +++ b/regression/scripts/segments/transition_segment_test.py @@ -60,9 +60,9 @@ def main(): # Truth values - departure_throttle_truth = np.array([0.65161055, 0.65183869, 0.65231038, 0.65255397]) - transition_1_throttle_truth = np.array([0.65435655, 0.64938063, 0.52075798, 0.58035463]) - cruise_throttle_truth = np.array([0.46376496, 0.46410324, 0.46478155, 0.46512157]) + departure_throttle_truth = np.array([0.65161054, 0.65183868, 0.65231039, 0.65255397]) + transition_1_throttle_truth = np.array([0.65434873, 0.64834806, 0.52075861, 0.58035428]) + cruise_throttle_truth = np.array([0.46376169, 0.46409995, 0.4647782 , 0.4651182 ]) transition_y_axis_rotations_truth = np.array([1.36961133, 1.34327318, 1.10250854, 0.06580108]) # Store errors diff --git a/regression/scripts/slipstream/propeller_interactions.py b/regression/scripts/slipstream/propeller_interactions.py index fdb922ed8f..eff0ddf8b0 100644 --- a/regression/scripts/slipstream/propeller_interactions.py +++ b/regression/scripts/slipstream/propeller_interactions.py @@ -64,7 +64,7 @@ def main(): T, Q, P, Cp, outputs , etap = run_downstream_propeller(prop, propeller_wake, conditions, plot_performance=plot_flag) # compare regression results: - T_iso_true, Q_iso_true, P_iso_true, Cp_iso_true, etap_iso_true = 3.22978737601075, 0.07191506, 48.9510291, 0.04572029, 0.58991371 + T_iso_true, Q_iso_true, P_iso_true, Cp_iso_true, etap_iso_true = 3.229745783054581, 0.07191487, 48.95089701, 0.04572016, 0.5899077 assert(abs(np.linalg.norm(T_iso)-T_iso_true)<1e-6) assert(abs(Q_iso-Q_iso_true)<1e-6) @@ -72,7 +72,7 @@ def main(): assert(abs(Cp_iso-Cp_iso_true)<1e-6) assert(abs(etap_iso-etap_iso_true)<1e-6) - T_true, Q_true, P_true, Cp_true, etap_true = 3.4488441181831626,0.07228313,49.20156567,0.04595429,0.62671635 + T_true, Q_true, P_true, Cp_true, etap_true = 3.4500540433063036,0.07229935,49.2126083,0.0459646,0.62679554 assert(abs(np.linalg.norm(T)-T_true)<1e-6) assert(abs(Q-Q_true)<1e-6) diff --git a/regression/scripts/slipstream/slipstream_test.py b/regression/scripts/slipstream/slipstream_test.py index 80ab0bb5c5..491db02e0a 100644 --- a/regression/scripts/slipstream/slipstream_test.py +++ b/regression/scripts/slipstream/slipstream_test.py @@ -208,9 +208,10 @@ def Lift_Rotor_Slipstream(wake_fidelity): def regress_2(results): - CL_truth = 0.41607982 - CDi_truth = 0.00895358 - CM_truth = 0.069561 + CL_truth = 0.41607938 + CDi_truth = 0.00895905 + CM_truth = 0.06956314 + CL = results.CL diff --git a/trunk/SUAVE/Components/Energy/Converters/Rotor.py b/trunk/SUAVE/Components/Energy/Converters/Rotor.py index d42b6771fa..b945d7e2c9 100644 --- a/trunk/SUAVE/Components/Energy/Converters/Rotor.py +++ b/trunk/SUAVE/Components/Energy/Converters/Rotor.py @@ -83,6 +83,7 @@ def __defaults__(self): self.vtk_airfoil_points = 40 self.induced_power_factor = 1.48 # accounts for interference effects self.profile_drag_coefficient = .03 + self.sol_tolerance = 1e-8 self.use_2d_analysis = False # True if rotor is at an angle relative to freestream or nonuniform freestream self.nonuniform_freestream = False diff --git a/trunk/SUAVE/Methods/Propulsion/Rotor_Wake/Fidelity_Zero/fidelity_zero_wake_convergence.py b/trunk/SUAVE/Methods/Propulsion/Rotor_Wake/Fidelity_Zero/fidelity_zero_wake_convergence.py index 9bde2d5610..eb06c23ba7 100644 --- a/trunk/SUAVE/Methods/Propulsion/Rotor_Wake/Fidelity_Zero/fidelity_zero_wake_convergence.py +++ b/trunk/SUAVE/Methods/Propulsion/Rotor_Wake/Fidelity_Zero/fidelity_zero_wake_convergence.py @@ -6,6 +6,7 @@ from SUAVE.Methods.Aerodynamics.Common.Fidelity_Zero.Lift.BET_calculations import compute_airfoil_aerodynamics,compute_inflow_and_tip_loss import numpy as np +import scipy as sp import copy ## @defgroup Methods-Propulsion-Rotor_Wake-Fidelity_Zero @@ -37,7 +38,60 @@ def fidelity_zero_wake_convergence(wake,rotor,wake_inputs): Properties Used: None + """ + + # Unpack some wake inputs + ctrl_pts = wake_inputs.ctrl_pts + Nr = wake_inputs.Nr + Na = wake_inputs.Na + + if wake_inputs.use_2d_analysis: + PSI = np.ones((ctrl_pts,Nr,Na)) + else: + PSI = np.ones((ctrl_pts,Nr)) + + PSI_final,infodict,ier,msg = sp.optimize.fsolve(iteration,PSI,args=(wake_inputs,rotor),xtol=rotor.sol_tolerance,full_output = 1,band=(1,0)) + + if ier!=1: + print("Rotor BEVW did not converge to a solution (Stall)") + + # Calculate the velocities given PSI + va, vt = va_vt(PSI_final, wake_inputs, rotor) + + + return va, vt + + +## @defgroup Methods-Propulsion-Rotor_Wake-Fidelity_Zero +def iteration(PSI, wake_inputs, rotor): + """ + Computes the BEVW iteration. + + Assumptions: + N/A + + Source: + N/A + + Inputs: + B number of rotor blades [-] + beta blade twist distribution [-] + r radius distribution [m] + R tip radius [m] + Wt tangential velocity [m/s] + Wa axial velocity [m/s] + U total velocity [m/s] + Ut tangential velocity [m/s] + Ua axial velocity [m/s] + cos_psi cosine of the inflow angle PSI [-] + sin_psi sine of the inflow angle PSI [-] + piece output of a step in tip loss calculation [-] + + Outputs: + dR_dpsi derivative of residual wrt inflow angle [-] + """ + # Unpack inputs to rotor wake fidelity zero U = wake_inputs.velocity_total Ua = wake_inputs.velocity_axial @@ -59,68 +113,89 @@ def fidelity_zero_wake_convergence(wake,rotor,wake_inputs): a_geo = rotor.airfoil_geometry a_loc = rotor.airfoil_polar_stations cl_sur = rotor.airfoil_cl_surrogates - cd_sur = rotor.airfoil_cd_surrogates - omega = rotor.inputs.omega + cd_sur = rotor.airfoil_cd_surrogates - # Setup a Newton iteration - diff = 1. - tol = 1e-6 # Convergence tolerance - ii = 0 + # Reshape PSI because the solver gives it flat + if wake_inputs.use_2d_analysis: + PSI = np.reshape(PSI,(ctrl_pts,Nr,Na)) + else: + PSI = np.reshape(PSI,(ctrl_pts,Nr)) - if use_2d_analysis: - PSI = np.ones((ctrl_pts,Nr,Na)) + # compute velocities + sin_psi = np.sin(PSI) + cos_psi = np.cos(PSI) + Wa = 0.5*Ua + 0.5*U*sin_psi + Wt = 0.5*Ut + 0.5*U*cos_psi + vt = Ut - Wt + + # compute blade airfoil forces and properties + Cl, Cdval, alpha, Ma, W = compute_airfoil_aerodynamics(beta,c,r,R,B,Wa,Wt,a,nu,a_loc,a_geo,cl_sur,cd_sur,ctrl_pts,Nr,Na,tc,use_2d_analysis) + + # compute inflow velocity and tip loss factor + lamdaw, F, piece = compute_inflow_and_tip_loss(r,R,Wa,Wt,B) + + # compute Newton residual on circulation + Gamma = vt*(4.*np.pi*r/B)*F*(1.+(4.*lamdaw*R/(np.pi*B*r))*(4.*lamdaw*R/(np.pi*B*r)))**0.5 + Rsquiggly = Gamma - 0.5*W*c*Cl + + return Rsquiggly.flatten() + +## @defgroup Methods-Propulsion-Rotor_Wake-Fidelity_Zero +def va_vt(PSI, wake_inputs, rotor): + """ + Computes the inflow velocities + + Assumptions: + N/A + + Source: + N/A + + Inputs: + B number of rotor blades [-] + beta blade twist distribution [-] + r radius distribution [m] + R tip radius [m] + Wt tangential velocity [m/s] + Wa axial velocity [m/s] + U total velocity [m/s] + Ut tangential velocity [m/s] + Ua axial velocity [m/s] + cos_psi cosine of the inflow angle PSI [-] + sin_psi sine of the inflow angle PSI [-] + piece output of a step in tip loss calculation [-] + + Outputs: + dR_dpsi derivative of residual wrt inflow angle [-] + + """ + + # Unpack inputs to rotor wake fidelity zero + U = wake_inputs.velocity_total + Ua = wake_inputs.velocity_axial + Ut = wake_inputs.velocity_tangential + ctrl_pts = wake_inputs.ctrl_pts + Nr = wake_inputs.Nr + Na = wake_inputs.Na + + # Reshape PSI because the solver gives it flat + if wake_inputs.use_2d_analysis: + PSI = np.reshape(PSI,(ctrl_pts,Nr,Na)) else: - PSI = np.ones((ctrl_pts,Nr)) - - PSIold = copy.deepcopy(PSI)*0 - - # BEVW Iteration - while (diff>tol): - # compute velocities - sin_psi = np.sin(PSI) - cos_psi = np.cos(PSI) - Wa = 0.5*Ua + 0.5*U*sin_psi - Wt = 0.5*Ut + 0.5*U*cos_psi - va = Wa - Ua - vt = Ut - Wt - - # compute blade airfoil forces and properties - Cl, Cdval, alpha, Ma, W = compute_airfoil_aerodynamics(beta,c,r,R,B,Wa,Wt,a,nu,a_loc,a_geo,cl_sur,cd_sur,ctrl_pts,Nr,Na,tc,use_2d_analysis) - - # compute inflow velocity and tip loss factor - lamdaw, F, piece = compute_inflow_and_tip_loss(r,R,Wa,Wt,B) - - # compute Newton residual on circulation - Gamma = vt*(4.*np.pi*r/B)*F*(1.+(4.*lamdaw*R/(np.pi*B*r))*(4.*lamdaw*R/(np.pi*B*r)))**0.5 - Rsquiggly = Gamma - 0.5*W*c*Cl - - # use analytical derivative to get dR_dpsi - dR_dpsi = compute_dR_dpsi(B,beta,r,R,Wt,Wa,U,Ut,Ua,cos_psi,sin_psi,piece) - - # update inflow angle - dpsi = -Rsquiggly/dR_dpsi - PSI = PSI + dpsi - diff = np.max(abs(PSIold-PSI)) - PSIold = PSI - - # If omega = 0, do not run BEVW convergence loop - if all(omega[:,0]) == 0. : - break - - # If its really not going to converge - if np.any(PSI>np.pi/2) and np.any(dpsi>0.0): - print("Rotor BEVW did not converge to a solution (Stall)") - break - - ii+=1 - if ii>10000: - print("Rotor BEVW did not converge to a solution (Iteration Limit)") - break + PSI = np.reshape(PSI,(ctrl_pts,Nr)) + # compute velocities + sin_psi = np.sin(PSI) + cos_psi = np.cos(PSI) + Wa = 0.5*Ua + 0.5*U*sin_psi + Wt = 0.5*Ut + 0.5*U*cos_psi + va = Wa - Ua + vt = Ut - Wt + return va, vt ## @defgroup Methods-Propulsion-Rotor_Wake-Fidelity_Zero -def compute_dR_dpsi(B,beta,r,R,Wt,Wa,U,Ut,Ua,cos_psi,sin_psi,piece): +def compute_dR_dpsi(PSI,wake_inputs,rotor): """ Computes the analytical derivative for the BEVW iteration. @@ -148,8 +223,37 @@ def compute_dR_dpsi(B,beta,r,R,Wt,Wa,U,Ut,Ua,cos_psi,sin_psi,piece): dR_dpsi derivative of residual wrt inflow angle [-] """ + # Unpack inputs to rotor wake fidelity zero + U = wake_inputs.velocity_total + Ua = wake_inputs.velocity_axial + Ut = wake_inputs.velocity_tangential + beta = wake_inputs.twist_distribution + r = wake_inputs.radius_distribution + ctrl_pts = wake_inputs.ctrl_pts + Nr = wake_inputs.Nr + Na = wake_inputs.Na + + # Unpack rotor data + R = rotor.tip_radius + B = rotor.number_of_blades + + # Reshape PSI because the solver gives it flat + if wake_inputs.use_2d_analysis: + PSI = np.reshape(PSI,(ctrl_pts,Nr,Na)) + else: + PSI = np.reshape(PSI,(ctrl_pts,Nr)) + + # An analytical derivative for dR_dpsi used in the Newton iteration for the BEVW # This was solved symbolically in Matlab and exported + # compute velocities + sin_psi = np.sin(PSI) + cos_psi = np.cos(PSI) + Wa = 0.5*Ua + 0.5*U*sin_psi + Wt = 0.5*Ut + 0.5*U*cos_psi + + lamdaw, F, piece = compute_inflow_and_tip_loss(r,R,Wa,Wt,B) + pi = np.pi pi2 = np.pi**2 BB = B*B @@ -174,4 +278,10 @@ def compute_dR_dpsi(B,beta,r,R,Wt,Wa,U,Ut,Ua,cos_psi,sin_psi,piece): Utcospsi + Uasinpsi ))/(BBB*pi2*utpUcospsi*utpUcospsi2*((16.*f_wa_2)/(BB*pi2*f_wt_2) + 1.)**(0.5))) dR_dpsi[np.isnan(dR_dpsi)] = 0.1 - return dR_dpsi \ No newline at end of file + + # This needs to be made into a jacobian + dR_dpsi = dR_dpsi.flatten() + L = np.size(PSI) + jac = np.eye(L)*dR_dpsi + + return jac \ No newline at end of file