In this example, we demonstrate how to use python to preprocess data by creating a frame structure. This example requires some Python familiarity. However, as we are developing a frame structure which contains mainly frame elements, the procedure is relatively straightforward.
Objective
We are going to create a three--storey frame as shown below.
Prerequisites
We are going to find section properties from AISC table, which is available online.
Database v15.0
Type EDI_Std_Nomenclature AISC_Manual_Label T_F W A d ddet \
0 W W44X335 W44X335 F 335.0 98.5 44 44
1 W W44X290 W44X290 F 290.0 85.4 43.6 43.625
2 W W44X262 W44X262 F 262.0 77.2 43.3 43.25
3 W W44X230 W44X230 F 230.0 67.8 42.9 42.875
4 W W40X655 W40X655 T 655.0 193.0 43.6 43.625
Ht h ... rts ho PA PA2 PB PC PD T WGi WGo
0 – – ... 4.24 42.2 132 – 148 104 120 38 5.5 –
1 – – ... 4.2 42 131 – 147 103 119 38 5.5 –
2 – – ... 4.17 41.9 131 – 147 102 118 38 5.5 –
3 – – ... 4.13 41.7 130 – 146 102 118 38 5.5 –
4 – – ... 4.71 40.1 132 – 149 104 121 34 7.5 –
[5 rows x 84 columns]
Section properties such as area and moment of inertia can be extracted from this table using section designations. It is possible to define such a function now.
def from_table(designation: str):
index = section_table.index[section_table['AISC_Manual_Label'] == designation].tolist()
assert len(index) == 1
a = section_table.at[index[0], 'A'] # area
sx = section_table.at[index[0], 'Sx'] # elastic modulus
ix = section_table.at[index[0], 'Ix'] # moment of inertia
zx = section_table.at[index[0], 'Zx'] # plastic modulus
return a, sx, ix, zx
Then we can use this function to extract section properties from the table. For example,
print(from_table('W44X230'))
(67.8, 971.0, 20800.0, 1100.0)
Geometry of the Frame Structure
For simplicity, we assume the frame structure has the same column/beam section for all the columns/beams on the same floor. Under such a condition, two lists of section designations can be provided so that elements can be created. Similarly, geometry information such as floor height, bay span, as well as floor mass, can be provided in the same way.
For example, we can define several lists as follows.
The similar procedure can be used to generate fibre based elements with slight modifications. Here we use F21 element. It relies on sections, and designations can be directly used to create sections. In this example, we use US2D category.
with open('fibre_frame.sp', 'w') as f:
f.write(f'material Bilinear1D 1 {e} {fy} {hardening}\n\n')
for v in section_pool.values():
f.write(f'section US2D {v.name} {v.tag} 1\n')
f.write('\n')
for v in element_pool.values():
f.write(f'element F21 {v.tag} {v.node_i.tag} {v.node_j.tag} {v.section.tag}\n')
with open('fibre_frame.sp', 'r') as f:
print(f.read())
material Bilinear1D 1 29 0.05 0.01
section US2D W21X68 1 1
section US2D W21X68 2 1
section US2D W21X68 3 1
section US2D W14X193 4 1
section US2D W14X159 5 1
section US2D W14X159 6 1
element F21 1 2 6 1
element F21 2 6 10 1
element F21 3 3 7 2
element F21 4 7 11 2
element F21 5 4 8 3
element F21 6 8 12 3
element F21 7 1 2 4
element F21 8 5 6 4
element F21 9 9 10 4
element F21 10 2 3 5
element F21 11 6 7 5
element F21 12 10 11 5
element F21 13 3 4 6
element F21 14 7 8 6
element F21 15 11 12 6
Generate Mass
For simplicity, we only apply horizontal mass to each node, the procedure is similar to that of frame elements.
To keep code modular, we define a mass grid to store mass at each node.
column_num = len(x_coor)
if 1 == column:
portions = 1.
else:
portions = 2. * (column_num - 1)
# unit conversion
portions *= 12 * 1000
# calculate tributary mass by assuming uniform density
# thus exterior nodes take half of the mass of interior nodes
mass_grid = np.zeros((len(y_coor), column_num))
for i in range(1, len(y_coor)):
for j in range(column_num):
mass_grid[i, j] = mass[i - 1] / portions
for j in range(1, column_num - 1):
mass_grid[i, j] *= 2.
print(mass_grid)
del column_num, portions
with open('mass.sp', 'w') as f:
for node_tag, mass in zip(node_grid.flatten(), mass_grid.flatten()):
if mass == 0.:
continue
f.write(f'element Mass {element_tag} {node_tag} {mass :.4e} 1 2\n')
element_tag += 1
with open('mass.sp', 'r') as f:
print(f.read())
element Mass 16 2 2.0833e-04 1 2
element Mass 17 6 4.1667e-04 1 2
element Mass 18 10 2.0833e-04 1 2
element Mass 19 3 2.0833e-04 1 2
element Mass 20 7 4.1667e-04 1 2
element Mass 21 11 2.0833e-04 1 2
element Mass 22 4 2.0833e-04 1 2
element Mass 23 8 4.1667e-04 1 2
element Mass 24 12 2.0833e-04 1 2
Analysis Settings
The basic geometry of the model is defined in node.sp, beam.sp, column.sp and mass.sp. We import those files in the main script.
file_list = ['node.sp', 'beam.sp', 'column.sp', 'mass.sp']
main_file = open('three-storey-frame.sp', 'w')
for file in file_list:
main_file.write(f'file {file}\n')
The boundary conditions can be assigned by simply fixing all bottom nodes.
main_file.write(f'fix 1 P 1 {" ".join(str(tag) for tag in node_grid[0, :])}\n')
pass
For results, we can, for example, record the nodal displacement history.
recorder_tag = 1
main_file.write(f'\nhdf5recorder {recorder_tag} Node U ' + ' '.join([str(tag) for tag in node_pool.keys()]) + '\n')
pass
For illustration, we can also record element yield flag at both ends so that the corresponding plastic hinge distribution can be generated. The yield flag is not generally available in other elements.
recorder_tag += 1
main_file.write(
f'\nhdf5recorder {recorder_tag} Element YF ' + ' '.join([str(tag) for tag in element_pool.keys()]) + '\n')
pass
We keep writing dynamic analysis step and the corresponding settings into the file.
Finally, let's pack everything into an archive so that it can be downloaded.
import os
from zipfile import ZipFile
file_list.append('fibre_frame.sp')
file_list.append('three-storey-frame.sp')
file_list.append('ELNS')
with ZipFile('three-storey-frame.zip', 'w') as f:
for file in file_list:
f.write(file)
Run Analysis
The analysis can be run by calling the executable.
The online documentation is not calling the executable, but it is possible to run the analysis locally with the application available.
from shutil import which
def run_analysis():
if os.name == 'nt':
exe_name = 'suanPan.exe'
else:
exe_name = 'suanpan'
if which(exe_name):
os.system(f'{exe_name} -np -f three-storey-frame.sp')
run_analysis()
Roof Displacement History
The result file will be generated and stored in .h5 file. We can read the file and plot the results.
import h5py
import matplotlib.pyplot as plt
import glob
def plot_result():
h5_file = glob.glob('*U.h5')
if len(h5_file) == 0:
return
plt.figure(figsize=(10, 7))
with h5py.File(h5_file[0], 'r') as result_file:
for group in result_file.values():
for key, data in group.items():
plt.plot(data[:, 0], data[:, 1] / max(y_coor) * 100., label=key)
plt.xlabel('Time (s)')
plt.ylabel('Drift (%)')
plt.legend()
plt.show()
plot_result()
Clean up the files to end this example.
for file in file_list:
if file == 'ELNS':
continue
if os.path.exists(file):
os.remove(file)
Let's create some TikZ commands to be used to plot plastic hinge distribution.
tikz: list = []
# generate coordinates to represent nodal positions
for node in node_pool.values():
tikz.append(f'\\coordinate(N{node.tag})at({node.x / 1E2},{node.y / 1E2});')
# draw frame elements
for element in element_pool.values():
node_i = element.node_i
node_j = element.node_j
if node_i.x == node_j.x:
# vertical columns
tikz.append(
f'\\draw[line width=1mm](N{element.node_i.tag})--(N{element.node_j.tag})node[midway,fill=white,font=\\tiny,rotate=90]{{{element.section.name}}};')
else:
# horizontal beams
tikz.append(
f'\\draw[line width=1mm](N{element.node_i.tag})--(N{element.node_j.tag})node[midway,fill=white,font=\\tiny]{{{element.section.name}}};')
hinge_label = '[fill=red,circle,draw,inner sep=0,minimum size=3mm]'
h5_file = glob.glob('*YF.h5')
if len(h5_file) > 0:
with h5py.File(h5_file[0], 'r') as h_file:
h5_prefix = h5_file[0].split('.')[0]
for g in h_file.values():
for k, v in g.items():
element_tag = int(k.replace(h5_prefix, ''))
node_i = element_pool[element_tag].node_i.tag
node_j = element_pool[element_tag].node_j.tag
if v[:, 1].max() > .5:
tikz.append(f'\\node{hinge_label}at($(N{node_i})!.1!(N{node_j})$){{}};')
if v[:, 2].max() > .5:
tikz.append(f'\\node{hinge_label}at($(N{node_i})!.9!(N{node_j})$){{}};')
with open(f'DIST.tex', 'w') as f:
f.write('\\begin{tikzpicture}[scale=2]\n')
f.write('\n'.join(tikz))
f.write('\n\\end{tikzpicture}\n')
with open(f'DIST.tex', 'r') as f:
print(f.read())
The ground motion shall be applied. The ELNS file contains one of the accelerograms of the 1940 El Centro Earthquake. It is normalised so that the maximum amplitude is unity, we assign a PGA of 0.5g.
