Wind Loading Analysis for the

Stand-alone EB Ground Station

Intro

Explanations

Resources

TODO

  • confirm the dimensions from the discussion with Glenn
    • antenna dimensions (frontal area with ice!)
    • rotator
    • L-braces
    • base plate
    • weight of a single cinder block
  • tidy-up the code and annotations to match the independent ground station
  • make a non-shitty version of the loading diagram

Problem Description

Failure Mode: Tipping

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.

Failure Mode: Bending

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.

Solution

Design Process
  1. choose arbitrary dimensions for the ground station (based on parts)
  2. find the worst-case loading
  3. determine if the ground station will tip over
    • If it will, add weight to the feet.
  4. determine if any part of the ground station will yield
    • If it will, go to 1, using a larger mast.
Information I/O:
  • Empirical C_d data: (typical element diameter, ripped plots from Munson book, wind speed) -> approximate C_d for the whole antenna (probably about 1.2)
  • Drag equation: (wind speed, air density, C_d, frontal area) -> drag on element
  • (element dimensions, material) -> element CoM
  • 2nd law for moments: (CoMs of all elements, masses of all elements, drag on elements) -> necessary ballast mass to prevent tipping

Setup

Imports


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()


3.4.3 (default, Nov 17 2016, 01:08:31) 
[GCC 4.8.4]

Diagram

input parameters


In [17]:
FS= mg(1)                           # factor of safety (NOT SURE IF BUILDING CODES HAVE ONE BUILT IN)
physical quantities

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
ground station parameters

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

secondary parameters


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')


mass contributions:
6.8414 kg
2.0000 lbm
4.4763 kg
17.2144 kg


area contributions
70cm:
1.0000 ft2
2m:
1.2000 ft2
mast:
1.9792 ft2
cross boom:
1.8160 ft2
base:
0.0000 ft2
total:
5.9951 ft2

Calculations

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.

Tipping

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.

find drag, CoP, CoM


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)

balance moments and find required ballast


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

bending on a welded mast

(probably not relevant)


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'))

report results


In [31]:
print(Lz_CoP.ounit('ft'))


8.3494 ft

In [33]:
print('m_ballast:', m_ballast.ounit('lbm'))
print("that's ", (m_ballast/mg(40,'lbm')), 'cinder blocks')


m_ballast: 576.0760 lbm
that's  14.4019 cinder blocks

Bending

determine the loads applied to the mast

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

symbolically find loading on the 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}} $