In [1]:
%matplotlib inline
import pylab as pl
from math import pi, sin, cos, exp, sqrt
import style ; style.clean()


Out[1]:

Resonant excitation

We want to study the behaviour of an undercritically damped SDOF system when it is subjected to a harmonic force $p(t) = p_o \sin\omega_nt$, i.e., when the excitation frequency equals the free vibration frequency of the system.

Of course, $\beta=1$, $D(\beta,\zeta)|_{\beta=1}=\displaystyle\frac{1}{2\zeta}$ and $\theta=\pi/2$, hence $$\xi(t)=\Delta_{st}\,\frac{1}{2\zeta}\cos\omega_nt.$$

Starting from rest conditions, we have

$$\frac{x(t)}{\Delta_{st}} = \exp(-\zeta\omega_n t)\left( -\frac{\omega_n}{2\omega_D}\sin(\omega_n t) -\frac{1}{2\zeta}\cos(\omega_n t)\right) + \frac{1}{2\zeta}\cos(\omega_n t)$$

and, multiplying both sides by $2\zeta$

\begin{align*} x(t)\frac{2\zeta}{\Delta_{st}} = \bar{x}(t)& = \exp(-\zeta\omega_n t)\left( -\zeta\frac{\omega_n}{\omega_D}\sin(\omega_n t) -\cos(\omega_n t)\right) + \cos(\omega_n t)\\ & = \exp(-\zeta\omega_n t)\left( -\frac{\zeta}{\sqrt{1-z^2}}\sin(\omega_n t) -\cos(\omega_n t)\right) + \cos(\omega_n t). \end{align*}

We have now a normalized function of time that grows, oscillating, from 0 to 1, where the free parameters are just $\omega_n$ and $\zeta$.

To go further, we set arbitrarily $\omega_n=2\pi$ (our plots will be nicer...) and have just a dependency on $t$ and $\zeta$.

Eventually, we define a function of $\zeta$ that returns a function of $t$ only, here it is...


In [2]:
def x_2z_over_dst(z):
    w = 2*pi
    # beta = 1, wn =w
    wd = w*sqrt(1-z*z)
    # Clough Penzien p. 43
    A = z/sqrt(1-z*z)
    def f(t):
        return (cos(wd*t)+A*sin(wd*t))*exp(-z*w*t)-cos(w*t)
    return pl.vectorize(f)

Above we compute some constants that depend on $\zeta$, i.e., the damped frequency and the coefficient in front of the sine term, then we define a function of time in terms of these constants and of $\zeta$ itself.

Because we are going to use this function with a vector argument, the last touch is to vectorize the function just before returning it to the caller.

Plotting our results

We start by using a function defined in the pylab aka pl module to generate a vector whose entries are 1001 equispaced real numbers, starting from zero and up to 20, inclusive of both ends, and assigning the name t to this vector.


In [3]:
t = pl.linspace(0,20,1001)
print(t)


[  0.     0.02   0.04 ...,  19.96  19.98  20.  ]

We want to see what happens for different values of $\zeta$, so we create a list of values and assign the name zetas to this list.


In [4]:
zetas = (.02, .05, .10, .20)
print(zetas)


(0.02, 0.05, 0.1, 0.2)

Now, the real plotting:

  • z takes in turn each of the values in zetas,
    • then we generate a function of time for the current z
    • we generate a plot with a line that goes through the point (a(0),b(0)), (a(1),b(1)), (a(2),b(2)), ... where, in our case, a is the vector t and b is the vector returned from the vectorized function bar_x
    • we make a slight adjustement to the extreme values of the y-axis of the plot
    • we give a title to the plot
    • we FORCE (pl.show()) the plot to be produced.

In [5]:
for z in zetas:
    # call the function of zeta that returns
    # a function of time, assign the name bar_x to this function
    bar_x = x_2z_over_dst(z)
    # do the plotting...
    pl.plot(t,bar_x(t))
    pl.ylim((-1.0, 1.0))
    pl.title(r'$\zeta=%4.2f$'%(z,))
    pl.show()


Wait a minute!

So, after all this work, we have that the greater the damping, the smaller the number of cycles that's needed to reach the maximum value of the response...

Yes, it's exactly like that, and there is a reason. Think of it.

.

.

.

.

.

.

.

.

.

.

We have normalized the response functions to have always a maximum absolute value of one, but in effect the max values are different, and a heavily damped system needs less cycles to reach steady-state because the maximum value is much, much smaller.

Let's plot the unnormalized (well, there's still the $\Delta_{st}$ normalization) responses.

Note the differences with above:

  • we focus on a shorter interval of time and, in each step
    • we don't add a title
    • we don't force the creation of a distinct plot in each cycle,
    • we add a label to each curve

at the end of the cycle,

  • we ask for the generation of a legend that uses the labels we specified to generate a, well, a legend for the curves
  • we ask to plot all the properly labeled curves using pl.plot().

In [6]:
t = pl.linspace(0,5,501)
for z in zetas:
    # call the function of zeta that returns
    # a function of time, assign the name bar_x to this function
    bar_x = x_2z_over_dst(z)
    # do the plotting...
    pl.plot(t,bar_x(t)/2/z, label=r'$\zeta=%4.2f$'%(z,))
pl.legend(ncol=5,loc='lower center', fancybox=1, shadow=1, framealpha=.95)
pl.grid()



In [ ]: