<jupyter_start><jupyter_text>Bathe Model Problem March 2020, By <NAME>This is an OpenSeesPy simulation of a simple linear 3 DOF system introduced by [Bathe and Noh](https://www.sciencedirect.com/science/article/abs/pii/S0045794912000107). 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]](https://link.springer.com/article/10.1007/s00366-020-00960-w). ![ModelProblem](imgs/ModelProblem.PNG)<jupyter_code>import numpy as np
import openseespy.opensees as ops
import matplotlib.pyplot as plt<jupyter_output><empty_output><jupyter_text>No units were defined in the original paper. So, I just assume the base units as follows:<jupyter_code>## Units
m = 1               # Meters
KN = 1              # KiloNewtons
s = 1               # Seconds

kg = KN*(s**2)/m    # mass unit (derived)<jupyter_output><empty_output><jupyter_text>Model specifications are defined as follows:<jupyter_code># 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
# Dynamic Analysis Parameters
dt = 0.2618
time = 10<jupyter_output><empty_output><jupyter_text>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]](https://link.springer.com/article/10.1007/s00366-020-00960-w).<jupyter_code>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)<jupyter_output><empty_output><jupyter_text>Model ConstructionI 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.<jupyter_code>def build_model():
    
    ops.wipe()
    ops.model('basic','-ndm',1,'-ndf',1)

    # Adding nodes
    [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)

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

    # Adding Nodal masses
    [ops.mass(n+1, mass) for n, mass in enumerate([m_1,m_2,m_3])];
    print('Model built successfully!')<jupyter_output><empty_output><jupyter_text>Dynamic AnalysisTo 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](http://web.mit.edu/kjb/www/Principal_Publications/On_a_Composite_Implicit_Time_Integration_Procedure_for_Nonlinear_Dynamics.pdf), named as 'TRBDF2' (with asymptotic annihilation property) are used.<jupyter_code>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)}<jupyter_output><empty_output><jupyter_text>Let's perform analysis:<jupyter_code>newmark = do_dynamic_analysis(time, dt,['Newmark',0.5,0.25])
bathe = do_dynamic_analysis(time, dt,['TRBDF2'])<jupyter_output>Model built successfully!
Done with  Newmark !
Model built successfully!
Done with  TRBDF2 !<jupyter_text>Visualization<jupyter_code>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);
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);<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1665175909361956, 'fraction_numerical': 0.0511329755515802, 'mean_word_length': 4.0330082520630155, 'pattern_counts': {'":': 0, '<': 33, '<?xml version=': 0, '>': 33, 'https://': 3, 'lorem ipsum': 0, 'www.': 2, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '8', 'dir': 'jupyter-structured-clean-dedup', 'id': '558254', 'license': 'MIT', 'n_tokens_mistral': 2705, 'n_tokens_neox': 2471, 'n_words': 528, 'path': 'Bathe Model Problem.ipynb', 'repo_name': 'mgrubisic/OpenSeesPy-Examples', 'size': '6723'}