231 lines
No EOL
8.3 KiB
Python
231 lines
No EOL
8.3 KiB
Python
# William Johnstone
|
|
# Apologies to anyone reading this, it became a mess. I wanted to tidy it but we don't always get what we want.
|
|
|
|
import matplotlib.pyplot as plt
|
|
import math
|
|
|
|
def raymers_iterate(current_max_weight, payload_weight, fuel_weight, a_constant, c_constant, variable_sweep_constant):
|
|
empty_max_weight_ratio = (a_constant * current_max_weight**c_constant * variable_sweep_constant)
|
|
|
|
new_max_weight = payload_weight / (1 - (fuel_weight / current_max_weight) - (empty_max_weight_ratio))
|
|
return new_max_weight
|
|
|
|
payload_mass_metric = 125 # Grams
|
|
payload_mass_imperial = payload_mass_metric * 0.002204623 # Converting grams to pounds
|
|
|
|
a_constant = 0.86 # Sailplane - Unpowered
|
|
c_constant = -0.05 # Sailplane - Unpowered
|
|
|
|
max_weight_guess = 6 # Initial weight guess of 6 lbs
|
|
current_estimate = max_weight_guess
|
|
|
|
imperial_estimates = []
|
|
|
|
for i in range(0,11):
|
|
current_estimate = raymers_iterate(current_estimate, payload_mass_imperial, 0, a_constant, c_constant, 1)
|
|
imperial_estimates.append(current_estimate)
|
|
|
|
metric_estimates = []
|
|
|
|
for imperial_estimate in imperial_estimates:
|
|
metric_estimates.append(imperial_estimate / 0.002204623)
|
|
|
|
plt.plot(metric_estimates, linestyle='--', marker='o')
|
|
|
|
# This section has been commented out because jordan (senft) does not believe the customer wants to see the numbers
|
|
# Add labels to each point
|
|
#for index, weight in enumerate(metric_estimates):
|
|
# plt.text(index, weight + 5, f'{index}: {weight:.2f}', fontsize=9, ha='right')
|
|
|
|
plt.title("Aircraft Raymer's Method Weight Estimation")
|
|
plt.xlabel("Iteration Number")
|
|
plt.ylabel("Weight Estimation (g)")
|
|
plt.grid(True)
|
|
#plt.show()
|
|
|
|
print("Raymer's Aircraft Calculations (IMPERIAL UNITS)\n")
|
|
print("Weight Estimation (mass): %0.4f lbs" % imperial_estimates[-1])
|
|
|
|
# Wing loading
|
|
typical_weight_area_ratio = 6 # Historical sailplane ratio
|
|
print("Typical Weight Area Ratio: %0.4f lbs/ft²" % typical_weight_area_ratio)
|
|
|
|
# To find area of the wings we must work backwards (ratio = weight/area)
|
|
# Area of wings
|
|
wing_area = imperial_estimates[-1] / typical_weight_area_ratio
|
|
print("Main Wing Area: %0.4f ft²" % wing_area)
|
|
|
|
# L/D Ratio
|
|
ld_ratio = 7.5 / 1.5
|
|
print("Lift to Drag Ratio: %0.4f" % ld_ratio)
|
|
|
|
# Aspect Ratio (Wing span) - Need L/D
|
|
aspect_ratio = 4.464 * (ld_ratio**0.69)
|
|
print("Aspect Ratio: %0.4f" % aspect_ratio)
|
|
|
|
# Wing span
|
|
wing_span = math.sqrt(aspect_ratio*wing_area)
|
|
print("Main Wing Span %0.4f ft" % wing_span)
|
|
|
|
# Chord Length
|
|
chord_length = wing_area / wing_span
|
|
print("Main Wing Chord Length %0.4f ft" % chord_length)
|
|
|
|
# Fuselage Length
|
|
fuselage_length = 0.86 * imperial_estimates[-1]**0.48
|
|
print("Fuselage Length %0.4f ft" % fuselage_length)
|
|
|
|
# Tail Positioning (Moment arm)
|
|
tail_position = fuselage_length * 0.65
|
|
print("Tail Positioning (Moment arm location) %0.4f ft" % tail_position)
|
|
|
|
# Vertical Tail Area
|
|
C_VT = 0.02 # Assumed constant (sail plane, raymers)
|
|
vertical_tail_area = C_VT * ((wing_span * wing_area) / tail_position)
|
|
print("Vertical Tail Area: {:.4f} ft²".format(vertical_tail_area))
|
|
|
|
# Horizontal Tail Area
|
|
C_HT = 0.5 # Assumed constant (sail plane, raymers)
|
|
horizontal_tail_area = C_HT * ((chord_length * wing_area) / tail_position)
|
|
print("Horizontal Tail Area: {:.4f} ft²\n".format(horizontal_tail_area))
|
|
|
|
|
|
# Convert everything to metric
|
|
|
|
print("End of imperial numbers (THANK THE LORD!)\n")
|
|
|
|
print("Weight Estimation (mass): %0.4f kg" % (metric_estimates[-1]*(10**-3)))
|
|
|
|
# Total Lift Force
|
|
lift_force_total = (metric_estimates[-1]*(10**-3)) * 9.81
|
|
#lift_force_per_ft = lift_force_total / wing_span
|
|
#lift_force_per_wing_per_ft = lift_force_per_ft / 2
|
|
print("Total Lift Force: {:.4f} N\n".format(lift_force_total))
|
|
#print("Lift Force per unit length: {:.4f} lbf/ft".format(lift_force_per_ft))
|
|
#print("Lift Force per unit length (1 wing): {:.4f} lbf/ft".format(lift_force_per_wing_per_ft))
|
|
|
|
# Bending Moments
|
|
|
|
number_of_wings = 2
|
|
|
|
print("BENDING MOMENTS\n")
|
|
wing_span_metric = wing_span * 0.3048
|
|
print("Wing Span: {:.4f} m".format(wing_span_metric))
|
|
|
|
# Bending Moment at Wing Root
|
|
bending_moment_wing = (lift_force_total / number_of_wings) * (wing_span_metric / 4)
|
|
print(f"Bending Moment at of 1 Wing from root: {bending_moment_wing:.4f} Nm\n")
|
|
|
|
# Max Stress Calcs
|
|
|
|
B = 4 * 0.001 # Width in m
|
|
D = 4 * 0.001 # Height or depth in m
|
|
|
|
# 4mm Spars
|
|
area_moment_inertia = ((B) * (D**3))/12
|
|
print(f"Second moment of inertia for 4mm: {area_moment_inertia} m^4")
|
|
|
|
max_stress = bending_moment_wing * ((D/2) / area_moment_inertia)
|
|
print(f"Max Stress for 4mm spar: {max_stress:.4f} Pa or {max_stress/1000000:.4f} MPa")
|
|
|
|
B = 6 * 0.001 # Width in m
|
|
D = 6 * 0.001 # Height or depth in m
|
|
|
|
# 6mm Spars
|
|
area_moment_inertia = ((B) * (D**3))/12
|
|
print(f"Second moment of inertia for 6mm: {area_moment_inertia} m^4")
|
|
|
|
max_stress = bending_moment_wing * ((D/2) / area_moment_inertia)
|
|
print(f"Max Stress for 6mm spar: {max_stress:.4f} Pa or {max_stress/1000000:.4f} MPa")
|
|
|
|
|
|
B = 12 * 0.001 # Width in m
|
|
D = 2 * 0.001 # Height or depth in m
|
|
|
|
# 6mm Spars
|
|
area_moment_inertia = ((B) * (D**3))/12
|
|
print(f"Second moment of inertia for new test: {area_moment_inertia} m^4")
|
|
|
|
max_stress = bending_moment_wing * ((D/2) / area_moment_inertia)
|
|
print(f"Max Stress for new test spar: {max_stress:.4f} Pa or {max_stress/1000000:.4f} MPa")
|
|
|
|
# Aerofoil Spars
|
|
|
|
B = 18 * 0.001 # Width in m
|
|
D = 3.5 * 0.001 # Height or depth in m
|
|
|
|
# 6mm Spars
|
|
area_moment_inertia = ((B) * (D**3))/12
|
|
print(f"Second moment of inertia for 6x3mm: {area_moment_inertia} m^4")
|
|
|
|
max_stress = bending_moment_wing * ((D/2) / area_moment_inertia)
|
|
print(f"Max Stress for 6x3mm: {max_stress:.4f} Pa or {max_stress/1000000:.4f} MPa")
|
|
|
|
B = 6 * 0.001 # Width in m
|
|
D = 0.5 * 0.001 # Height or depth in m
|
|
|
|
# 6mm Spars
|
|
area_moment_inertia = ((B) * (D**3))/12
|
|
print(f"Second moment of inertia for 6x3.5mm: {area_moment_inertia} m^4")
|
|
|
|
max_stress = bending_moment_wing * ((D/2) / area_moment_inertia)
|
|
print(f"Max Stress for 6x3.5mm: {max_stress:.4f} Pa or {max_stress/1000000:.4f} MPa")
|
|
|
|
|
|
# Bending Moment for tail (Rohan's numbers )
|
|
bending_moment_wing = (lift_force_total / number_of_wings) * (wing_span_metric / 4)
|
|
print(f"Bending Moment at of 1 Wing from root: {bending_moment_wing:.4f} Nm\n")
|
|
|
|
|
|
|
|
######################################################################
|
|
# 27/11/24
|
|
# Updated stress calcs with new weight (from detailed CAD model, including payload mass)
|
|
|
|
print("#############################################################")
|
|
print("Everything after this is based on numbers from final CAD model. Thus more relevant")
|
|
|
|
print("Weight (Final, from CAD Model) (mass): %0.4f kg" % (262*(10**-3)))
|
|
|
|
# Total Lift Force
|
|
lift_force_total = (262*(10**-3)) * 9.81
|
|
#lift_force_per_ft = lift_force_total / wing_span
|
|
#lift_force_per_wing_per_ft = lift_force_per_ft / 2
|
|
print("Total Lift Force: {:.4f} N\n".format(lift_force_total))
|
|
|
|
# Bending Moments
|
|
|
|
number_of_wings = 2
|
|
|
|
print("BENDING MOMENTS (LATEST, THESE ARE THE ONLY ONES RELEVANT)\n")
|
|
wing_span_metric = wing_span * 0.3048
|
|
print("Wing Span: {:.4f} m".format(wing_span_metric))
|
|
|
|
# Bending Moment at Wing Root
|
|
bending_moment_wing = (lift_force_total / number_of_wings) * (wing_span_metric / 4)
|
|
print(f"Bending Moment at of 1 Wing from root: {bending_moment_wing:.4f} Nm\n")
|
|
|
|
# Max bending stress for Balsa (https://woodworkly.com/is-balsa-wood-strong/)
|
|
bending_strength_psi_lower = 2550
|
|
bending_strength_psi_upper = 3170
|
|
bending_strength_megapascal_lower = bending_strength_psi_lower * 0.006894757 # (1 psi = 0.006894757 MPa)
|
|
bending_strength_megapascal_upper = bending_strength_psi_upper * 0.006894757 # (1 psi = 0.006894757 MPa)
|
|
|
|
# Max Stress Calcs
|
|
|
|
B = 18 * 0.001 # Width in m
|
|
D = 2.5 * 0.001 # Height or depth in m
|
|
|
|
# 4mm Spars
|
|
area_moment_inertia = ((B) * (D**3))/12
|
|
print(f"Second moment of inertia for final model: {area_moment_inertia} m^4")
|
|
|
|
max_stress = bending_moment_wing * ((D/2) / area_moment_inertia)
|
|
|
|
max_stress_megapascal = max_stress/1000000
|
|
|
|
print(f"Max Stress for final model spar: {max_stress:.4f} Pa or {max_stress_megapascal:.4f} MPa")
|
|
|
|
print(f"Balsa Wood Bending Stress Upper Limit: {bending_strength_megapascal_upper:.4f}MPa, Lower Limit: {bending_strength_megapascal_lower:.4f} MPa")
|
|
|
|
print(f"Safety Factor: {bending_strength_megapascal_lower/max_stress_megapascal:.1f} - {bending_strength_megapascal_upper/max_stress_megapascal:.1f}") |