In [1]:

import numpy as np
import os
import glob
import shutil
import matplotlib.pyplot as plt
import json
import subprocess
from matplotlib import rcParams
rcParams['font.family'] = 'serif'
rcParams['font.size'] = 16
from IPython.display import Image
%matplotlib inline



Module for computations (if needed)



In [2]:

import FEM_utilities as FEM



# Lab 06

## Write input file - Shell



In [3]:

def write_shell_inp_file(L, h, b, num_e, y_e, E, ν, num_modes):
"""
function to write inp file:
L        : lenght of cantilever
h        : height
b        :  thickness
num_e    : number of elements
y_e      : number of elements in y direction
E        : Young's modulus
nu       : Poisson ratio
num_modes: number of eigenvalues to be estimated
"""

elname = "EALL"
matname = "mat1"
bcname = "fixed_nodes"

pa = [1.0,0.0,0.0]
pb = [0.0,1.0,0.0]
# pc = [0.0,0.0,0.0]
ax = 3
α = 0.0

intpoints = 5

basename = 'b1_ne'+str(num_e).zfill(3)+'_ye'+str(y_e).zfill(2)+'_aS'

filename = basename+'.inp'

outfile = open(filename, "wt")

outfile.write("** Lab 06 input file buckling - shell\n")

nodes, elements, fixed, loaded = FEM.Nodes(num_e, y_e, L, h)

# NODES section
outfile.write("**\n")
outfile.write("** Nodes\n")
outfile.write("**\n")
outfile.write("*NODE, nset = nglobal\n")

for i in range(nodes.shape[0]):
nodestring = "{0:4d}".format(int(nodes[i,0]))
for j in range(1,nodes.shape[1]):
nodestring+=",{0:8}".format(nodes[i,j])
nodestring+="\n"
outfile.write(nodestring)

# ELEMENTS section
outfile.write("**\n")
outfile.write("** Elements\n")
outfile.write("**\n")
outfile.write("*ELEMENT, type = S4R\n")

for i in range(elements.shape[0]):
elstring = "{0:4d}".format(int(elements[i,0]))
for j in range(1,elements.shape[1]):
elstring+=",{0:4d}".format(int(elements[i,j]))
elstring+="\n"
outfile.write(elstring)

# MATERIAL section
outfile.write("**\n")
outfile.write("** Materials\n")
outfile.write("**\n")
outfile.write("*MATERIAL, name = {0}\n".format(matname))
outfile.write("*Elastic\n")
outfile.write("{0},{1:6}\n".format(E,ν))

# SETS section
# NODES
outfile.write("**\n")
outfile.write("** Sets\n")
outfile.write("**\n")
outfile.write("*NSET, nset = {0}\n".format(bcname))

fix_str = "{0:4d}".format(int(fixed[0]))

for i in range(1,len(fixed)):
fix_str+=",{0:4d}".format(int(fixed[i]))
fix_str+="\n"
outfile.write(fix_str)

il_str+="\n"
outfile.write(il_str)

outfile.write(el_str)

# ELEMENTS
outfile.write("*Elset, elset = {0}, generate\n".format(elname))
outfile.write("{0:4d},{1:4d},{2:4d}\n".format(1,elements.shape[0],1))

# ORIENTATION
outfile.write("**\n")
outfile.write("** LOCAL ORIENTATION\n")
outfile.write("**\n")
outfile.write("*orientation, name = local_orientation\n")
outfile.write("".join(str(pa+pb))[1:-1])
outfile.write("\n")
outfile.write("{0},{1:4}\n".format(ax,α))

# SHELL PROPERTIES
outfile.write("**\n")
outfile.write("** SHELL PROPERTIES\n")
outfile.write("**\n")
outfile.write("*Shell Section, elset = {0}, material = {1}, orientation = local_orientation\n" \
.format(elname,matname))
outfile.write("{0:4},{1:4d}\n".format(b,intpoints))

# calculation steps
outfile.write("**\n")
outfile.write("** STEP\n")
outfile.write("**\n")

outfile.write("*STEP\n")
outfile.write("*BUCKLE\n")
outfile.write("{0},\n".format(num_modes))

# BOUNDARY CONDITIONS
outfile.write("**\n")
outfile.write("** BOUNDARY CONDITIONS\n")
outfile.write("**\n")
outfile.write("*Boundary\n")
outfile.write("{0}, ENCASTRE\n".format(bcname))

outfile.write("**\n")
outfile.write("**\n")
#outfile.write("{0}, {1:2d}, {2}\n".format(loadnodes, 2, -0.5 ))
outfile.write("{0}, {1:2d}, {2}\n".format(loadintnodes, 2, -1.0/y_e ))
outfile.write("{0}, {1:2d}, {2}\n".format(loadextnodes, 2, -0.5/y_e ))

outfile.write("*NODE FILE, LAST MODE={0}\n".format(num_modes))
outfile.write("U, RF\n")
#outfile.write("*EL FILE, ELSET=PRINT, LAST MODE={0}\n".format(num_modes))
#outfile.write("ENER,\n")
#outfile.write("ELEN,\n")
outfile.write("*END STEP\n")
outfile.close()

if not os.path.exists('../Lab06_abaqus/'+basename+'/'):
os.makedirs('../Lab06_abaqus/'+basename+'/')

shutil.move(filename,'../Lab06_abaqus/'+basename+'/'+filename)



## Write input file - Beam



In [4]:

def write_beam_inp_file(L, h, b, num_e, E, ν, num_modes):
"""
function to write inp file:
L        : lenght of cantilever
h        : height
b        :  thickness
num_e    : number of elements
E        : Young's modulus
nu       : Poisson ratio
num_modes: number of eigenvalues to be estimated
"""

basename = 'b1_ne'+str(num_e).zfill(3)+'_aB'

filename = basename+'.inp'

outfile = open(filename, "wt")

outfile.write("** Lab 06 input file buckling - beam\n")

elname = "EALL"
matname = "mat1"

n_coords = np.zeros((num_e+1,2))

for ni in range(1,num_e+1):
n_coords[ni,0] = ni/(num_e)*L

# NODES section
outfile.write("**\n")
outfile.write("** Nodes\n")
outfile.write("**\n")
outfile.write("*NODE\n")

for i in range(num_e+1):
nodestring = "{0:4d},{1:8},{2:8}\n".format(i+1,n_coords[i,0],n_coords[i,1])
outfile.write(nodestring)

# ELEMENTS section
outfile.write("**\n")
outfile.write("** Elements\n")
outfile.write("**\n")
outfile.write("*ELEMENT, TYPE={0}, ELSET={1}\n".format('B31',elname))
for i in range(1,num_e+1):
outfile.write("{0:4d},{1:4d},{2:4d}\n".format(i,i,i+1))

# BEAM section

outfile.write("**\n")
outfile.write("** Beam section\n")
outfile.write("**\n")
outfile.write("*BEAM SECTION, SECTION=RECT, MATERIAL={0}, ELSET={1}\n".format(matname,elname))
outfile.write("{0},{1}\n".format(b,h))
# default orientation?
outfile.write("0,0,-1\n")

# MATERIAL section
outfile.write("**\n")
outfile.write("** Materials\n")
outfile.write("**\n")
outfile.write("*MATERIAL, name = {0}\n".format(matname))
outfile.write("*ELASTIC\n")
outfile.write("{0},{1:6}\n".format(E,ν))

# BOUNDARY CONDITIONS
outfile.write("**\n")
outfile.write("** Boundary conditions\n")
outfile.write("**\n")
outfile.write("*BOUNDARY\n")
outfile.write("1,\t1,\t6\n")

# calculation steps

outfile.write("**\n")
outfile.write("** Step\n")
outfile.write("**\n")

outfile.write("*STEP\n")
outfile.write("*BUCKLE\n")
outfile.write("{0},\n".format(num_modes))

outfile.write("**\n")
outfile.write("**\n")
outfile.write("{0}, {1:2d}, {2}\n".format(num_e+1, 2, -1. ))

outfile.write("*NODE FILE, LAST MODE={0}\n".format(num_modes))
outfile.write("U, RF\n")
#outfile.write("*EL FILE, ELSET=PRINT, LAST MODE={0}\n".format(num_modes))
#outfile.write("ENER,\n")
#outfile.write("ELEN,\n")
outfile.write("*END STEP\n")

outfile.close()

if not os.path.exists('../Lab06_abaqus/'+basename+'/'):
os.makedirs('../Lab06_abaqus/'+basename+'/')

shutil.move(filename,'../Lab06_abaqus/'+basename+'/'+filename)



Set parameters for simulation:

• L, b, h : cantilever dimensions
• eltype : type of beam or shell element
• eig_num : maximum eigenvalue to be estimated
• E, $\nu$: elastic properties (Young modulus and Poisson ratio, as TYPE is by default isotropic)
• num_elems: number of elements in x direction
• yel_shell: number of elements in y direction for shell case


In [5]:

L = 200. # lenght [mm]
h = 15.  # height [mm]
b = 1.5  # thickness [mm]

eltype = ['S4R','B31']

eig_num = 6

E = 72000. # modulus [MPa]
ν = 0.33   # Poisson's coefficient

# number of elements

num_elems = [4, 5, 10, 20, 25, 50, 100,200]

yel_shell = [1, 2, 5,15]

for ne in num_elems:
for et in eltype:
# time steps
if et[0] == 'S':
for ye in yel_shell:
write_shell_inp_file(L, h, b, ne, ye, E, ν, eig_num)
elif et[0] == 'B':
write_beam_inp_file(L, h, b, ne, E, ν, eig_num)
else:
print("Element type not recognized")




In [6]:

#old NOT used
def generateCases_old(loc, prefix, filename):
"""
generates all cases from specified root directory and starting with prefix
loc : root directory
"""
dirs = glob.glob(loc+"/"+prefix+"*")

command = ['abaqus']

options = [['datacheck','interactive'],['interactive','continue']]

outfile = open(loc+"/"+filename, "wt")

outfile.write("#!/bin/bash\n\n")

for di in dirs:
test = di.split('/')[-1]
outfile.write("cd {0}\n\n".format(di))
outfile.write("/usr/simulia/abaqus/Commands/abaqus j={0} datacheck interactive\n".format(test))
outfile.write("/usr/simulia/abaqus/Commands/abaqus j={0} interactive continue\n\n".format(test))
outfile.write("cd ..\n\n")

outfile.close()
os.chmod(loc+"/"+filename, 0o744)




In [7]:

def generateCases(loc, prefix):
"""
generates all cases from specified root directory and starting with prefix
loc : root directory
"""
curr_dir = os.getcwd()

if not loc.startswith("/"):
#relative path
base_dir = os.path.abspath(loc)
else:
base_dir = loc

#print(curr_dir)
dirs = glob.glob(loc+"/"+prefix+"*")

command = ['/usr/simulia/abaqus/Commands/abaqus']

options = [['datacheck','interactive'],['interactive','continue']]

for di in dirs:
test = di.split('/')[-1]
os.chdir(di)
job="j={0}".format(test)

for opt in options:
subprocess.call(command+[job]+opt)

os.chdir(base_dir)

os.chdir(curr_dir)



# Generate Abaqus cases



In [8]:

generateCases("../Lab06_abaqus/","b1_")



# Analytical solution:

$P_{cr}$ analytical relation:

$P_{cr} = 4.0129 \cdot \frac{1}{L^2} \cdot \sqrt{E I_{22} G J_{t}}$

where:

• $G = \frac{E}{2 \cdot (1+\nu)}$
• $I_{22} = \frac{1}{12} \cdot hb^3$
• $J_{t} = \beta \cdot hb^3$
• $\beta = 0.312 \quad (\frac{h}{b} = 10)$


In [9]:

G = E/(2*(1+ν))
I22 = h*b**3/12
Jt = .312*h*b**3

Pcr = 4.0129/L**2*np.sqrt(E*I22*G*Jt)

print("Pcr = ",Pcr)




Pcr =  36.1527829995



# Simulations data



In [10]:

file=open(filename,'r')

eig = -1.0

for line in row:
strlist = line.split()
if len(strlist) == 2:
if strlist[0] == '1':
eig = np.abs(float(strlist[1]))
break
file.close()
return eig




In [18]:

dirs = glob.glob("../Lab06_abaqus/b1_*")

#print(dirs)

bl_el = []
bl_eig = []

shell_data = {}

for di in dirs:
test = di.split('/')[-1]

if '_aS' in test:
#print("shell")
yel = int(test[11:13])
if yel not in shell_data:
shell_data[yel] = [[],[]]
if a != -1.0:
shell_data[yel][0].append(int(test[5:8]))
shell_data[yel][1].append(a)
elif '_aB' in test:
#print("beam")
if a != -1.0:
bl_el.append(int(test[5:8]))
bl_eig.append(a)
else:
print("Element type not recognized")

beam_el = np.array(bl_el)
beam_eig = np.array(bl_eig)

sort_index = np.argsort(beam_el)
beam_el = beam_el[sort_index]
beam_eig = beam_eig[sort_index]

#print(beam_el)
#print(beam_eig)

for ye in shell_data:
tmp_el = np.array(shell_data[ye][0])
tmp_eig = np.array(shell_data[ye][1])
sort_index = np.argsort(tmp_el)
shell_data[ye] =[tmp_el[sort_index],tmp_eig[sort_index]]

print(shell_data)




{2: [array([  4,   5,  10,  20,  25,  50, 100, 200]), array([ 41.981,  40.72 ,  38.896,  38.409,  38.348,  38.266,  38.244,
38.239])], 1: [array([  4,   5,  10,  20,  25,  50, 200]), array([ 41.058,  39.722,  37.525,  36.476,  36.278,  35.964,  35.871])], 5: [array([  4,   5,  10,  20,  25,  50, 100, 200]), array([ 42.18 ,  40.911,  39.03 ,  38.378,  38.24 ,  37.915,  37.747,
37.689])], 15: [array([  4,   5,  10,  20,  25,  50, 100, 200]), array([ 42.209,  40.939,  39.052,  38.386,  38.239,  37.868,  37.663,
37.59 ])]}




In [17]:

max_beam_el = np.max(beam_el)

max_shell_el = 0

for ye in shell_data:
tmp_max = np.max(shell_data[ye][0])
if tmp_max > max_shell_el:
max_shell_el = tmp_max

max_el = max(max_shell_el, max_beam_el)

x_pcr = np.linspace(0,max_el,10)
y_pcr = Pcr*np.ones_like(x_pcr)

plt.figure(figsize=(16,10), dpi=300)

plt.plot(x_pcr,y_pcr, '-', lw=3, label='$P_{cr}$ analytical')

plt.plot(beam_el,beam_eig,'-', lw=2, color=(0.5,0.0,0.0), label="$P_{cr}$ beam")
plt.plot(beam_el,beam_eig,'h', ms=8, color=(0.5,0.0,0.0))

for i, ye in enumerate(shell_data):
plt.plot(shell_data[ye][0],shell_data[ye][1],'-', lw=2, color=(0.0,(i+1)/len(shell_data),0.0), \
label="$P_{cr}$ shell $n_y$="+str(ye)  )
plt.plot(shell_data[ye][0],shell_data[ye][1],'h', ms=8, color=(0.0,(i+1)/len(shell_data),0.0) )

plt.xlim([0,max_el])
plt.ylim([np.floor(Pcr-1),np.ceil(Pcr+6)])
plt.xticks([0.0,5,10,20,25,50,100,200])
plt.yticks(np.arange(np.floor(Pcr-1),np.ceil(Pcr+6)+.1,.5))
plt.xlabel(r'number of elements',fontsize=14)
plt.ylabel(r'Critical load $[N]$',fontsize=14)
plt.grid()
plt.savefig('Lab06a.jpg')







In [ ]: