In [1]:
from IPython.display import HTML
from math import pi, cos, sqrt
HTML(open("00_custom.css", "r").read())
Out[1]:
The data of our problem, first the numerical constants and later the definitions of the properties depending on the position.
In [2]:
rho = 2500.0
e = 30E9
h = 32.0
m = 80E3
thick = 0.25
Re0 = 1.8
DeltaR = 0.600
def Re(z): return Re0 - DeltaR*z/h
def Ri(z): return Re(z) - thick
def A(z): return pi*(Re(z)**2-Ri(z)**2)
def J(z): return pi*(Re(z)**4-Ri(z)**4)/4
We need to integrate a function, using the trapezoidal approximation
In [3]:
def trap(f,z0,z1,n):
"Approximate the integral of f on interval z0, z1 using trapezoidal rule with n intervals."
dz = (z1-z0)/n
result = (f(z0)+f(z1))/2
for i in range(1,n):
z=z0+i*dz
result = result+f(z)
return result*dz
Let's test the trapezoidal rule function... The volume of the tower is $$V = \int_0^h A(z) \,\text{d}z $$ We have $A(z) = \pi R_e^2 - \pi(R_e^2 -2R_et+t^2)=2\pi R_et-\pi t^2$, substituting $V/\pi=2t\int_0^hR_e\,\text{d}z-ht^2$, but $R_e$ is a linear function of $z$ so we have $$V=2\pi(2R_e(h/2)-t/2) t h$$
In [4]:
vol1 = trap(A,0,h,8)
vol2 = pi*h*thick*(2*Re(h/2)-thick)
print("The volume using trapezes", vol1, "m^3")
print("The vouume, exact formula", vol2, "m^3")
print("While we are at it, the total concrete mass", vol1*rho, "kg")
In [5]:
def compute_gp(psi0, psi2):
mass = trap(lambda z:A(z)*rho*psi0(z)**2, 0, h, 8)
stif = trap(lambda z:J(z)*e*psi2(z)**2, 0, h, 8)
tot_m = m + mass
w2 = stif/tot_m ; w = sqrt(w2)
print("The modal mass, lumped mass excluded, is", mass, "kg")
print("The modal mass of the complete system is", tot_m, "kg")
print("The modal stiffness is ", stif, "N/m")
print("The first eigenvalue is ", w2, "(rad/s)^2")
print("The first circular frequency is ", w, "rad/s")
print("The natural frequency of vibration is ", w/2/pi, "Hz")
print("The natural period of vibration is ", 2*pi/w, "s")
Let's start with the familiar
\begin{align*} \psi(z) &= 1 - \cos\frac{\pi z}{2h}\\ \psi' &= \frac\pi{2h}\sin\frac{\pi z}{2h},\\ \psi'' &= \frac{\pi^2}{4h^2}\cos\frac{\pi z}{2h}. \end{align*}That respects, as requested, the geometrical boundary conditions at $z=0$.
In [6]:
def psi0(z): return 1-cos(pi*z/2/h)
def psi2(z): return pi**2/4/h**2*cos(pi*z/2/h)
print("******** psi = 1 - cos(pi*z/2h) ********")
compute_gp(psi0, psi2)
In [7]:
def psi0(z): return z*z*(z-3*h)/2/h**3
def psi2(z): return 3*(z-h)/h**3
print("******** psi = (z^3 - 3 z^2 h) / (2 h^3)")
compute_gp(psi0, psi2)
In [8]:
def psi0(z): return (z**4 + 6*(h**2*z**2) - 4*(h*z**3)) / (3*h**4)
def psi2(z): return (4*z**2 + 4*h**2 - 8*(h*z)) / h**4
print("******** psi = (z^4 + 6 h^2 z^2 -4 h z^3) / (3 h^4) ********")
compute_gp(psi0, psi2)