We want to know how heavy the ballast of the ground station need to be in order to keep it from tipping over in high winds.
We are assuming the worst-case is that antenna will be covered in ice during a wind storm. The criterion for tipping is when the moment due to wind is equal and opposite to the moment due to gravity.
I'm choosing the ballast weight as the "trim" variable, since it's a lot easier to just swap out some larger blocks than to find out, say, that we need to redesign the ground station for a different mast length.
We also want to know the necessary size and material of the mast, so that it doesn't bend.
The design will be done with an Euler-Bernoulli bending model in this script, but it should be checked in SolidWorks.
In [16]:
import numpy as np
import sys
import matplotlib.pyplot as plt
import sympy as sym
import pandas as pd
import magnitude as mag
from magnitude import mg
mag.new_mag('lbm', mag.Magnitude(0.45359237, kg=1))
mag.new_mag('lbf', mg(4.4482216152605, 'N'))
mag.new_mag('mph', mg(0.44704, 'm/s'))
mag.new_mag('slug', mg(1,'lbf')/mg(1,'ft/s2'))
imperialUnits_typical= {'length':'ft', 'area':'ft2', 'mass':'lbm', 'force':'lbf', 'time':'s'}
bgUnits= {'length':'ft', 'area':'ft2', 'mass':'slug', 'force':'lbm', 'time':'s'}
SIUnits= {'length':'m', 'area':'m2', 'mass':'kg', 'force':'N', 'time':'s'}
myUnits= imperialUnits_typical
from IPython.display import display, Markdown, Latex
def printBig(*message): # like print(), but with header2 formatting
display(Markdown( '## ' + ' '.join([str(x) for x in message]) ))
# drag force:
def dragf(rho, A, Cd, v): return(mg(1/2)*rho*A*Cd*v**2)
# second area moment of an annulus:
def I_annulus(OD, ID): return( mg(np.pi/4)*( (OD/mg(2))**4 - (ID/mg(2))**4 ) )
# area of an annulus:
def annularArea(ID, OD): return(mg(np.pi)*(OD**2/mg(4)-ID**2/mg(4)))
print(sys.version)
%matplotlib inline
sym.init_printing()
In [17]:
FS= mg(1) # factor of safety (NOT SURE IF BUILDING CODES HAVE ONE BUILT IN)
In [18]:
g= mg(9.81,'m/s2') # gravitational acceleration
rho_air= mg(1.225,'kg/m3') # density of air at 0 C
mu_air= mg(1.7e-5,'Pa s') # dynamic viscosity of air at 0 C
nu_air= mu_air/rho_air # kinematic viscosity of air at 0 C
rho_steel= mg(7.8e3,'kg/m3') # density of steel
rho_al= mg(2.7e3,'kg/m3') # density of aluminum
rho_fg= mg(1520,'kg/m3') # density of fiberglass
rho_concrete= mg(2300,'kg/m3') # density of concrete
sigma_y_al= mg(276,'MPa') # tensile yield strength of 6061-T6 aluminum
In [25]:
Cd= mg(1.2) # Drag coefficient of a cylinder/antenna at relevant speeds, (THIS IS A GUESS)
# Fg= mg(300,'lbf') # weight of the ground station (THIS IS A GUESS)
v_wind_actual= mg(100,'mph') # max expexted wind speed (WHAT'S THE JUSTIFICATION HERE?)
v_wind= v_wind_actual*FS # max expexted wind speed (WHAT'S THE JUSTIFICATION HERE?)
# base
Lx_base= mg(5, 'ft')
Ly_base= mg(5, 'ft')
m_base= mg(2,'lbm')
A_base= mg(0,'ft2')
Lz_CoP_base= mg(0,'inch')
Lz_CoM_base= mg(0,'inch')
# antennas
A_70cm= mg(1.0,'ft2') # frontal area of the 70 cm antenna
m_70cm= mg(5,'lbm') # mass of the 70 cm antenna
A_2m= mg(1.2,'ft2') # frontal area of the 2 m antenna
m_2m= mg(6, 'lbm') # mass fo the 2 m antenna
# mast
L_mast= mg(10,'ft') # length of the mast
OD_mast= mg(2+3/8,'inch') # outer diameter of the mast
ID_mast= mg(2,'inch') # inner diameter of the mast
# cross boom
OD_crossBoom= mg(2.0,'inch') # outer diameter of the cross boom
ID_crossBoom= mg(1.5,'inch') # inner diameter of the cross boom
L_crossBoom= mg(60+60+10.75,'inch') # length of the cross boom
rho_crossBoom= rho_fg # density of the cross boom
In [20]:
# frontal areas
A_mast= L_mast*OD_mast
A_crossBoom= L_crossBoom*OD_crossBoom
A_total= A_70cm+A_2m+A_mast+A_crossBoom+A_base
# masses
m_mast= rho_al*L_mast*annularArea(ID=ID_mast, OD=OD_mast)
m_crossBoom= rho_fg*L_crossBoom*annularArea(ID=ID_crossBoom, OD=OD_crossBoom)
m_structure= m_mast+m_crossBoom+m_70cm+m_2m+m_base
# centers of mass
Lz_CoM_70cm= L_mast
Lz_CoM_2m= L_mast
Lz_CoM_mast= L_mast/mg(2)
Lz_CoM_crossBoom= L_mast
In [23]:
print(
'mass contributions:',
'mast:',m_mast,
'base:',m_base,
'cross boom:',m_crossBoom,
'structure (everything but the ballast):',m_structure,
sep='\n' )
print('\n')
print(
'area contributions',
'70cm:', A_70cm.ounit(myUnits['area']),
'2m:', A_2m.ounit(myUnits['area']),
'mast:', A_mast.ounit(myUnits['area']),
'cross boom:', A_crossBoom.ounit(myUnits['area']),
'base:', A_base.ounit(myUnits['area']),
'total:', A_total.ounit(myUnits['area']),
sep='\n')
I'm counting the mass of the feet separately, since they aren't rigidly attached to the rest of the ground station. I assume that if the station tips, it will rotate around the pivots on two of the feet. Meanwhile the third (upstream) foot would dangle from the pivot, keeping its CoM directly below the pivot. So, unlike all the other components, I don't assume that the CoM of the feet is centered on the mast. I take it to be the mass of a single foot on the upstream (-X) leg.
Technically, I'm breaking the right hand rule. I'm using the -Y axis to be the direction of positive rotation. So, on the diagram, a positive rotation would be counter-clockwise. This is just because I'm used to thinking of positive rotation as "out-of-the-page", and the -Y axis is out-of-the-page on the diagram.
In [26]:
# center of pressure for the whole station
# assuming the CoP for each element is at the CoM (no lift on elem, constant density of elem)
Lz_CoP= (
A_mast*Lz_CoM_mast
+A_70cm*Lz_CoM_70cm
+A_2m*Lz_CoM_2m
+A_crossBoom*Lz_CoM_crossBoom
+A_base*Lz_CoP_base
)/A_total
# height of the CoM of the whole station
# (only useful for finding critical tipping angle)
Lz_CoM_total= (
Lz_CoM_mast*m_mast
+Lz_CoM_2m*m_2m
+Lz_CoM_70cm*m_70cm
+Lz_CoM_crossBoom*m_crossBoom
+Lz_CoM_base*m_base
)/m_structure
# drag on the station
#D_wind= mg(1/2)*rho_air*A_total*Cd*v_wind**2
D_wind= dragf(rho=rho_air, A=A_total, Cd=Cd, v=v_wind)
In [29]:
# Moments exerted on the station
M_wind= -D_wind*(Lz_CoP) # negative sign from CW direction
M_structure= m_structure*g*Lx_base/mg(2)
# 0 == M_wind + M_structure + M_foot + M_balast
M_ballast= -M_wind-M_structure
# M_ballast == m_ballast*g*(Lx_base/2)
m_ballast= M_ballast/g/(Lx_base/mg(2))
# m_ballast= M_ballast/g/(Lx_leg) # ballast on the mast
In [ ]:
# # I= pi/4*(r_o^2 - r_i^2)
# I_mast= I_annulus(OD= OD_mast, ID= ID_mast)
# M_bending= D_wind*(Lz_CoP - Lz_leg)
# # sigma_max= M*y/I
# # max stress in the mast, assuming it's anchored at the height of the foot
# sigma_max_mast= M_bending*(OD_mast/mg(2))/I_mast
# print(
# 'tension from bending, if the mast is welded at the height of the feet:',
# sigma_max_mast.ounit('MPa')
# )
# print('tensile yield stress of the mast:', sigma_y_al.ounit('MPa'))
In [31]:
print(Lz_CoP.ounit('ft'))
In [33]:
print('m_ballast:', m_ballast.ounit('lbm'))
print("that's ", (m_ballast/mg(40,'lbm')), 'cinder blocks')
F_p is a point load. F_d is a distributed load.
Drag on the legs is beared equally by the feet and mast (by symmetry) If I ever add in the calculation for the drag on the armpit beams, I'll need to add half their drag to this number for the same reason. And, I'll probably just assume that they have no reaction forces, since they would make things statically indeterminate, which sucks.
Use 2nd law for forces and moments to get the reaction forces... Should probably make a diagram explaining where these loads are applied.
In [ ]:
F_p_legsDrag= mg(1/2)*dragf(rho=rho_air, A=A_legs, Cd=Cd, v=v_wind)
F_p_boom= dragf(
rho=rho_air, Cd=Cd, v=v_wind,
A= A_70cm+A_2m+A_crossBoom
)
F_d_mast= dragf(rho=rho_air, Cd=Cd, A=A_mast, v=v_wind) # distributed over the mast
# about the base:
# sum(M) == 0 == (F_p_legsReact-F_p_legsDrag)*Lz_leg - F_p_boom*L_mast - F_d_mast*L_mast/2
# sum(F) == 0 == -F_p_legsReact + F_p_legsDrag + F_p_boom + F_p_baseReact + F_d_mast
F_p_legsReact= mg(1)/Lz_leg*( F_p_boom*L_mast + F_d_mast*L_mast/mg(2) ) + F_p_legsDrag
F_p_baseReact= F_p_legsReact - F_p_legsDrag - F_p_boom + F_d_mast
In [ ]:
# x = sym.symbols('x')
# print(F_p_boom.ounit('N'))
# # load per length, as a function of length along the mast:
# Fdist_expr_mast= F_d_mast/L_mast \
# + F_p_baseReact*mg(sym.DiracDelta(x),'/m') \
# + (-F_p_legsReact+F_p_legsDrag)*mg(sym.DiracDelta(x-Lz_leg.toval(ounit='m')),'/m') \
# + F_p_boom*mg(sym.DiracDelta(x-L_mast.toval(ounit='m')),'/m')
# print('load distribution expression:\n', Fdist_expr_mast.toval(ounit='N/m'), 'N/m')
# # shear load, as a function of length along the mast:
# F_expr_mast= mg(sym.integrate(Fdist_expr_mast.toval(ounit='N/m'), x),'N')
# print('shear load expression:\n', F_expr_mast.toval(ounit='N'), 'N')
# # bending moment, as a function of length along the mast:
# M_expr_mast= mg(sym.integrate(F_expr_mast.toval(ounit='N'),x), 'N m')
# print('bending load expression:\n', M_expr_mast.toval(ounit='N m'), 'N m')
In [ ]:
# x = sym.symbols('x')
# # load per length, as a function of length along the mast:
# Fdist_expr_mast= F_d_mast.toval('N')/L_mast.toval('m') \
# + F_p_baseReact.toval('N')*sym.DiracDelta(x) \
# + (-F_p_legsReact.toval('N')+F_p_legsDrag.toval('N'))*sym.DiracDelta(x-Lz_leg.toval(ounit='m')) \
# + F_p_boom.toval('N')*sym.DiracDelta(x-L_mast.toval(ounit='m'))
# print('load distribution expression:\n', Fdist_expr_mast.toval(ounit='N/m'), 'N/m')
# # shear load, as a function of length along the mast:
# F_expr_mast= mg(sym.integrate(Fdist_expr_mast.toval(ounit='N/m'), x),'N')
# print('shear load expression:\n', F_expr_mast.toval(ounit='N'), 'N')
# # bending moment, as a function of length along the mast:
# M_expr_mast= mg(sym.integrate(F_expr_mast.toval(ounit='N'),x), 'N m')
# print('bending load expression:\n', M_expr_mast.toval(ounit='N m'), 'N m')
In [ ]:
xs= np.linspace(-1e-6, L_mast.toval(ounit='m')+1e-6, 300)
Fsf= sym.lambdify(x,F_expr_mast.toval(ounit='N'), ['numpy','sympy'])
Fs= mg(np.array([Fsf(x) for x in xs]),'N')
Vs= Fs/annularArea(OD=OD_mast, ID=ID_mast)
plt.plot(xs, Vs.toval(ounit='MPa'))
print(F_d_mast, F_p_baseReact, F_p_legsReact, F_p_legsDrag, F_p_boom, sep='\n')
Msf= sym.lambdify(x, M_expr_mast.toval(ounit='N m'), ['numpy', 'sympy'])
Ms= mg(np.array([Msf(x) for x in xs]),'N m')
plt.figure()
plt.plot(xs, Ms.toval(ounit='N m'))
sigmas= Ms*(OD_mast/mg(2))/I_mast
plt.figure()
plt.plot(xs, sigmas.toval(ounit='MPa'))
In [ ]:
#OD_mast_new=
#ID_mast_new=
In [ ]:
print(D_wind.ounit('lbf'))
print(Lz_CoP.ounit('ft'))
print(m_ballast.ounit('lbm'))
Okay, so, I totally forgot about the "armpit" beams that fix the legs at a particular angle. But, those don't really matter too much, since they're only going to make things worse. The takeaway is that we'd need to add an excessive amount of ballast to the feet (more than doubling their mass), which means we need a different design of ground station (pole on a plate, or something).
It looks like the mast will still need some sort of reinforcement, since it's (currently) way under-built for a 100 mph load.
$ My_{tip} = 0\\ 0 = Fx_{windCR}*Hz_{CoP} - Fg*Lx_{leg} \\ Fx_{windCR} = \frac{Fg*Lx_{leg}}{Hz_{CoP}} $