# Bathe Model Problem¶

#### March 2020, By Amir Hossein Namadchi¶

This is an OpenSeesPy simulation of a simple linear 3 DOF system introduced by Bathe and Noh. The model was used to demonstrate the ability of Bathe method in filtering out the unwanted artificial high-frequency responses. The effectiveness of numerical dissipation in a specific time integration algorithm could be evaluated by analyzing this benchmark model problem [1].

In [4]:
import numpy as np
import openseespy.opensees as ops
import matplotlib.pyplot as plt


No units were defined in the original paper. So, I just assume the base units as follows:

In [5]:
## Units
m = 1               # Meters
KN = 1              # KiloNewtons
s = 1               # Seconds

kg = KN*(s**2)/m    # mass unit (derived)


Model specifications are defined as follows:

In [6]:
# Nodal mass values
m_1, m_2, m_3 = 0*kg, 1*kg, 1*kg

# Modulus of Elasticity
E_1, E_2 = (10.0**7)*(KN/m**2) , 1.0*(KN/m**2)

# Prescribed displacment @ node 1
d_p = lambda t: np.sin(1.2*t)*m

In [7]:
# Dynamic Analysis Parameters
dt = 0.2618
time = 10


Here's the exact solution for comparison. It can be obtained using modal decomposition technique, considering only the lowest frequency plus the static correction [1].

In [8]:
d_2_ex = lambda t : (2.7273e-7)*np.sin(t) + np.sin(1.2*t)
d_3_ex = lambda t : (2.7273)*np.sin(t) - (2.2727)*np.sin(1.2*t)

v_2_ex = lambda t : (2.7273e-7)*np.cos(t) + 1.2*np.cos(1.2*t)
v_3_ex = lambda t : (2.7273)*np.cos(t) - 1.2*(2.2727)*np.cos(1.2*t)

a_2_ex = lambda t : -(2.7273e-7)*np.sin(t) - 1.2*1.2*np.sin(1.2*t)
a_3_ex = lambda t : -(2.7273)*np.sin(t) + 1.2*1.2*(2.2727)*np.sin(1.2*t)


### Model Construction¶

I use zeroLength element to model the spings, then add nodal masses via mass command. Also,I wrapped the whole modeling steps in build_model() function to use it later.

In [21]:
def build_model():

ops.wipe()
ops.model('basic','-ndm',1,'-ndf',1)

[ops.node(n+1, 0.0) for n in range(3)]

# Defining Material (spring stiffness)
ops.uniaxialMaterial('Elastic', 1, E_1)
ops.uniaxialMaterial('Elastic', 2, E_2)

ops.element('zeroLength', 1, *[1,2], '-mat', 1, '-dir', 1)
ops.element('zeroLength', 2, *[2,3], '-mat', 2, '-dir', 1)

[ops.mass(n+1, mass) for n, mass in enumerate([m_1,m_2,m_3])];
print('Model built successfully!')


### Dynamic Analysis¶

To use different time integration algorithms, the function do_dynamic_analysis is defined herein. Just pass the time, dt and integrator_params and run it multiple times with different integrators. Here, Newmark with $\gamma=0.5$ and $\beta=0.25$ (a non-dissipative algorithm) and the Bathe method, named as 'TRBDF2' (with asymptotic annihilation property) are used.

In [22]:
def do_dynamic_analysis(time, dt, integrator_params):

build_model()
time_domain = np.arange(0,time,dt)

# Constraints
ops.timeSeries('Path', 1,
'-dt', dt,
'-values', *np.vectorize(d_p)(time_domain),
'-time', *time_domain)
ops.pattern('Plain', 1, 1)
ops.sp(1, 1, 1)

# Analysis
ops.constraints('Transformation')
ops.numberer('Plain')
ops.system('ProfileSPD')
ops.algorithm('Linear')
ops.integrator(*integrator_params)
ops.analysis('Transient')

time_lst =[0]           # list to hold time stations for plotting
res_node_2 = [[0,0,0]]  # response params of node 2
res_node_3 = [[0,0,0]]  # response params of node 3

for i in range(len(time_domain)):
ops.analyze(1, dt)
time_lst.append(ops.getTime())
res_node_2.append([ops.nodeDisp(2,1),
ops.nodeVel(2,1),
ops.nodeAccel(2,1)])
res_node_3.append([ops.nodeDisp(3,1),
ops.nodeVel(3,1),
ops.nodeAccel(3,1)])

print('Done with ', integrator_params[0],'!')

return {'time_list':np.array(time_lst),
'Node 2': np.array(res_node_2),
'Node 3': np.array(res_node_3)}


Let's perform analysis:

In [23]:
newmark = do_dynamic_analysis(time, dt,['Newmark',0.5,0.25])
bathe = do_dynamic_analysis(time, dt,['TRBDF2'])

Model built successfully!
Done with  Newmark !
Model built successfully!
Done with  TRBDF2 !


### Visualization¶

In [32]:
plt.figure(figsize=(12,4))
plt.plot(np.arange(0, time, 0.2*dt),
np.vectorize(v_2_ex)(np.arange(0,time,0.2*dt)),
color = '#7D3C98', linewidth = 2, label = 'Exact')
plt.plot(newmark['time_list'], newmark['Node 2'][:,1],
color = '#d62d20', linewidth=1.5, ls='--', label = 'Newmark (CAAM)' )
plt.plot(bathe['time_list'], bathe['Node 2'][:,1],
color = '#008000', linewidth=1.75, label = 'Bathe (TRBDF2)', marker = 'o',ls=':')

plt.ylabel('Velocity of Node 2 (m/s)', {'fontname':'Cambria', 'fontstyle':'italic','size':14})
plt.xlabel('Time (sec)', {'fontname':'Cambria', 'fontstyle':'italic','size':14})
plt.xlim([0.0, time])
plt.ylim([-3, 3])
plt.grid()
plt.legend()
plt.yticks(fontname = 'Cambria', fontsize = 14)
plt.xticks(fontname = 'Cambria', fontsize = 14);

In [33]:
plt.figure(figsize=(12,4))
plt.plot(np.arange(0, time, 0.2*dt),
np.vectorize(a_2_ex)(np.arange(0,time,0.2*dt)),
color = '#7D3C98', linewidth = 2, label = 'Exact')
plt.plot(newmark['time_list'], newmark['Node 2'][:,2],
color = '#d62d20', linewidth=1.5, ls='--', label = 'Newmark (CAAM)' )
plt.plot(bathe['time_list'], bathe['Node 2'][:,2],
color = '#008000', linewidth=1.75, label = 'Bathe (TRBDF2)', marker = 'o',ls=':')

plt.ylabel('Acceleration of Node 2 m/s^2', {'fontname':'Cambria', 'fontstyle':'italic','size':14})
plt.xlabel('Time (sec)', {'fontname':'Cambria', 'fontstyle':'italic','size':14})
plt.xlim([0.0, time])
#plt.ylim([-30, 30])
plt.grid()
plt.legend()
plt.yticks(fontname = 'Cambria', fontsize = 14)
plt.xticks(fontname = 'Cambria', fontsize = 14);


### Closure¶

As expected, non-dissipative methods like CAAM are unable to represent the exact response due to spurious high frequency oscillations. On the other hand, TRBDF2 scheme can effectively damp out these oscillations via algorithmic damping.

### References¶

• Bathe, K.J. and Noh, G., 2012. Insight into an implicit time integration scheme for structural dynamics. Computers & Structures, 98, pp.1-6.
• Namadchi, A.H., Jandaghi, E. and Alamatian, J., 2020. A new model-dependent time integration scheme with effective numerical damping for dynamic analysis. Engineering with Computers, pp.1-16.