

In [1]:

    
%pylab inline









    



Populating the interactive namespace from numpy and matplotlib

Constant acceleration algorithm

We are, again and again, integrating the equation of motion for the same example, that is a damped SDOF with period $T=0.6\\,{}$s, with a triangular loading with a peak value of 40 kN.

The loading function



In [2]:

    
def p(t):
    if t < t0 :
        return 54*(t0-t)**2*t/t0**3
    return 0.00

The SDOF system parameters

and some derived values



In [3]:

    
mass = 3.
k    = 1200.
zeta = 0.10
wn =sqrt(k/mass)
wd   = wn * sqrt(1.00-zeta**2)
damp = 2*zeta*mass*wn
t0 = 0.12

print mass, damp, k
print wn, wd









    



3.0 12.0 1200.0
20.0 19.8997487421

Initialization of the CA algorithm



In [4]:

    
# time step duration
h = t0/30.

# we require the response from 0 to 6 s, it is convnient
# to define a slightly modified duration
duration = 10*t0 + h/2

# The constants used by the algorithm
k_ = k + 2*damp/h + 4*mass/h/h
cv = 2*damp + 4*mass/h
ca = 2*mass

# We'll use these three containers to store our results
x = [] ; v = [] ; t = []

Initial conditions of the system



In [5]:

    
T = 0.00
X = 0.00
V = 0.00
P = p(T)
A = (P - V*damp - X*k)/mass

Iteration of the elementary, step-wise incremental solution



In [6]:

    
while T < duration:
    x.append(X) ; v.append(V) ; t.append(T)
    # print "%6.3f   %+12.10f %+12.10f" % (t, X, V)
    T = T+h
    Ph = p(T)
    dp_ = (Ph-P) + cv*V + ca*A
    dx  = dp_/k_
    dv  = 2*dx/h - 2*V
    X = X+dx ; V = V+dv
    P = Ph ; A = (P - damp*V - k*X)/mass

Plotting the resulting displacements



In [7]:

    
plot(t, [1000*X for X in x]) ; xlim((0,10*t0)) ; xlabel('t/s') ; ylabel('x/mm') ; grid() ; xticks(linspace(0,10*t0,11));