\iffalse

```
from scipy import *
from scipy.integrate import quad

from IPython.core.display import HTML
def css_styling():
    styles = open("custom.css", "r").read()
    return HTML(styles)
css_styling()
```




<style>
    @font-face {
        font-family: "Computer Modern";
        src: url('http://mirrors.ctan.org/fonts/cm-unicode/fonts/otf/cmunss.otf');
    }
    div.cell{
        width:800px;
        margin-left:10% !important;
        margin-right:auto;
    }
    h1 {
        font-family: Candara, Cambria, serif;
    }
    h2 { text-align:center; }
    h4{
        margin-top:12px;
        margin-bottom: 3px;
       }
    div.text_cell_render{
        font-family: Candara, Calibri, Trebuchet, Geneva, sans-serif;
        line-height: 155%;
        font-size: 130%;
        width:800px;
        margin-left:auto;
        margin-right:auto;
    }
    .CodeMirror{
            background-color: #FDFCFF;
            font-family:"Aurulent Sans Mono","Source Code Pro", source-code-pro,Consolas, monospace;
    }
    .prompt{
        display: None;
    }
    .text_cell_render h5 {
        font-weight: 300;
        font-size: 16pt;
        color: #3035A1;
        font-style: italic;
        margin-bottom: .5em;
        margin-top: 0.5em;
        display: block;
    }

    .warning{
        color: rgb( 240, 20, 20 )
        }  

    BODY {background-color: rgb( 252, 255, 220 )}

</style>





    

\fi
# The water tower

The tower is in reinforced concrete,


```
g  = 2500.0
E  = 30E9
g0 = 9.80665 # standard acceleration of gravity
```

the heigth of the tower `H` and the lumped mass `M` at its top (representing the
mass of the tank and of its content) are respectively


```
H = 45.0
M = 1.2E6
```

the cross-section is annular, the outer radius varies linearly from `R0` to `RH`
and the thickness varies linearly from `T0` to `TH`


```
R0 = 3.20 ; RH = 2.40
T0 = 0.25 ; TH = 0.20
```

We can define the radius, the thickness, the cross section area and flexural
inertia of the section as functions of the vertical coordinate `z`,


```
def R(z):
    return R0+(RH-R0)*z/H

def T(z):
    return T0+(TH-T0)*z/H

def A(z):
    t = T(z)
    return pi*(2*t*R(z)-t*t)

def J(z):
    r_ext = R(z)
    r_int = r_ext - T(z)
    return pi/4*(r_ext**4-r_int**4)
```

## Total mass of the beam

We are ready to compute the total mass of the beam and to print aside of the lumped mass, so that we can judge the relevance of the distributed mass on our results.

To compute the beam mass, `mass`, we have to integrate $\gamma A(z)\,\text{d}z$
over the length of the beam. To this purpose we use the library function
`quad` that has 3 mandatory arguments:

  1. a function, library or user defined, to be numerically integrated,
  2. the lower limit of the interval of integration,
  3. the upper limit of the interval of integration.

The integrand is here _defined on the spot_ using the `lambda` syntax
(in `matlab` it is `@` syntax).

Note that the function `quad` returns both the definite integral value
and an estimate of the error, and that we discard the error estimate
by assigning it to the dummy variable `_`.


```
mass, _ = quad(lambda z:A(z)*g, 0, H)
print M, mass
print mass/M
```

    1200000.0 429710.970149
    0.358092475124


The mass of the beam represent a significant fraction of the lumped mass,
we can expect a small, but significant correction from the naive estimate 
$\omega^2 \approx EJ/3ML^3$

# Choice of the shape function

I will simply use the one-minus-quarter-cosine shape function, no discussion, sorry...

\begin{align}
    \psi(z)   & = 1 - \cos\left(\frac\pi2\frac zH\right),\\
    \psi'(z)  & = \frac\pi{2H} \sin\left(\frac\pi2\frac zH\right),
    &\delta(z) &= \int_0^z \psi'^2(s) ds = 
    \frac18\frac\pi H\left(\frac\pi H z- \sin\left(\frac\pi Hz\right)\right),\\
    \psi''(z) & = \frac{\pi^2}{4H^2} \cos\left(\frac\pi2\frac zH\right).
\end{align}

### Definitions

We define here the shape function and its derivatives 
in terms of `Z`, a constant proportional to the wave length.

While we're at it, 
$\int_0^z \psi'^2(s)\,\text{d}s= \frac{\pi}{8H} \left(\frac\pi{H}z 
    - \sin\left(\frac\pi{H}z\right)\right)$
is the increment of vertical displacement at quote `z` for a unit $\delta Z_0$.


```
Z = 2*H/pi

def psi0(z): return 1-cos(z/Z)
def psi1(z): return sin(z/Z)/Z
def psi2(z): return cos(z/Z)/Z/Z
def de_u(z):  return pi*(-sin(pi*z/H) + pi*z/H)/(8*H)
```

# Generalised properties

We compute the generalized mass and stiffness, and also the contributions to the geometric stiffness due to the lumped mass
and the distributed mass


```
mstar = quad(lambda z: psi0(z)**2*A(z)*g, 0, H)[0]
kstar = quad(lambda z: psi2(z)**2*E*J(z), 0, H)[0]
kgeoM = quad(lambda z: psi1(z)**2,        0, H)[0]*M*g0
kgeom = quad(lambda z: de_u(z)*A(z)*g,    0, H)[0]*g0
```

# Results

Using a helper function, we display the results in terms of eigenvalues, frequencies, periods under different assumptions
regarding the geometrical stiffness


```
def pr_res(s,w2):
    wn = sqrt(w2)
    fn = wn/2/pi
    Tn = 1.0/fn
    print "%36s%12.6f%8.5f%8.5f%8.5f" % (s, w2, wn, fn, Tn)
    return None

print
print "                 k*:", kstar,"N/m"
print "                 m*:", mstar,"kg"
print "                  M:    ", M,"kg"
print "           k_G*M*g0:", kgeoM,"N/m"
print "    g0 \int de_u dm: ", kgeom,"N/m"

print
print " "*39+"      w^2 "+"      w "+"      f "+"      T"
pr_res("k*/M ", kstar/M)
pr_res("k*/(M+m*)", kstar/(M+mstar))
pr_res("(k*-kg*M*g0)/(M+m*)", (kstar-kgeoM)/(M+mstar))
pr_res("(k*-kg M g0-\\int(kg dm)g0)/(M+m*)", (kstar-kgeoM-kgeom)/(M+mstar))
```

    
                     k*: 17347995.0148 N/m
                     m*: 81844.8215573 kg
                      M:     1200000.0 kg
               k_G*M*g0: 322625.853333 N/m
        g0 \int de_u dm:  29517.750381 N/m
    
                                                 w^2       w       f       T 
                                   k*/M    14.456663 3.80219 0.60514 1.65252
                               k*/(M+m*)   13.533616 3.67881 0.58550 1.70794
                     (k*-kg*M*g0)/(M+m*)   13.281927 3.64444 0.58003 1.72405
       (k*-kg M g0-\int(kg dm)g0)/(M+m*)   13.258899 3.64128 0.57953 1.72554


## Comments

We have a large dependency of $\omega^2$ on the account of distributed mass, a small but not negligible dependency
on the geometrical stiffness associated with the weight of the lumped mass and a very small dependency on the geometrical stiffness
associated with the distributed weight of the beam

854e4ce70b21eef7ed13af2950b0229c226f0b70
ipynb
Jupyter Notebook
dati_2014/ha01/04_Water_Tower.ipynb
shishitao/boffi_dynamics
365f16d047fb2dbfc21a2874790f8bef563e0947
MIT
dati_2014/ha01/04_Water_Tower.ipynb
shishitao/boffi_dynamics
365f16d047fb2dbfc21a2874790f8bef563e0947
MIT
dati_2014/ha01/04_Water_Tower.ipynb
shishitao/boffi_dynamics
365f16d047fb2dbfc21a2874790f8bef563e0947
MIT
2019-06-23T12:32:39.000Z
2021-08-15T18:33:55.000Z
Qwen/Qwen-72B
1. YES
2. YES