In [152]:
import numpy as np
import os
import shutil
import matplotlib.pyplot as plt
import json
from matplotlib import rcParams
rcParams['font.family'] = 'serif'
rcParams['font.size'] = 14
from IPython.display import Image
Module for computations
In [153]:
import FEM_utilities as FEM
Parameters:
In [154]:
a = 1000. # lenght [mm]
b = 200. # width [mm]
t = 1.5 # thickness [mm]
Compute node coordinates, elements, constrained and loaded nodes
In [155]:
nx = 40
ny = 8
nodes, elements, fixed, loaded = FEM.Nodes(nx, ny, a, b)
In [156]:
#nodes
In [157]:
#elements
In [158]:
#fixed
In [159]:
#loaded
Set parameters for simulation:
In [160]:
basename = 'test_'+str(nx)+'x'+str(ny)+'y'
filename = basename+'.inp'
eltype = 'S4R' #'SC8R'
matname = 'material_1'
E = 72000. # modulus [MPa]
ν = 0.33 # Poisson's coefficient
elname = 'plate_1'
bcname = 'fixed_nodes'
loadextnodes = 'loaded_nodes_external'
loadintnodes = 'loaded_nodes_internal'
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
F = 0.02 # load [N/mm]
In [161]:
outfile = open(filename, "wt")
outfile.write("** Lab 02 input file test\n")
# 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 = {0}\n".format(eltype))
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)
if len(loaded > 2):
outfile.write("*Nset, nset = {0}\n".format(loadintnodes))
il_str = "{0:4d}".format(int(loaded[1]))
for i in range(2,len(loaded)-1):
il_str+=",{0:4d}".format(int(loaded[i]))
il_str+="\n"
outfile.write(il_str)
outfile.write("*Nset, nset = {0}\n".format(loadextnodes))
el_str = "{0:4d},{1:4d}\n".format(int(loaded[0]),int(loaded[-1]))
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(t,intpoints))
# calculation steps
outfile.write("**\n")
outfile.write("** STEP\n")
outfile.write("**\n")
outfile.write("*Step, name = step_1\n")
outfile.write("*Static\n")
# BOUNDARY CONDITIONS
outfile.write("**\n")
outfile.write("** BOUNDARY CONDITIONS\n")
outfile.write("**\n")
outfile.write("*Boundary\n")
outfile.write("{0}, ENCASTRE\n".format(bcname))
# LOADS
outfile.write("**\n")
outfile.write("** LOADS\n")
outfile.write("**\n")
outfile.write("*Cload\n")
outfile.write("{0}, {1:2d}, {2}\n".format(loadextnodes, 3, F*b/ny/2 ))
outfile.write("{0}, {1:2d}, {2}\n".format(loadintnodes, 3, F*b/ny ))
# FIELD OUTPUT
outfile.write("**\n")
outfile.write("** FIELD OUTPUT\n")
outfile.write("**\n")
outfile.write("*Output, field, variable=PRESELECT\n")
# HISTORY OUTPUT
outfile.write("**\n")
outfile.write("** HISTORY OUTPUT\n")
outfile.write("**\n")
outfile.write("*Output, history, variable=PRESELECT\n")
# PRINT DISPLACEMENT AND STRESS DATA
outfile.write("**\n")
outfile.write("** PRINT RESULTS TO FILE\n")
outfile.write("**\n")
outfile.write("*NODE FILE\n")
outfile.write("U,\n")
outfile.write("RF,\n")
outfile.write("*EL FILE\n")
outfile.write("S,\n")
outfile.write("*End Step\n\n")
outfile.close()
copy file to analysis directory
In [162]:
if not os.path.exists('../Lab02_abaqus/'+basename):
os.makedirs('../Lab02_abaqus/'+basename)
shutil.copy(filename,'../Lab02_abaqus/'+basename)
Out[162]:
Maximum stress at $x = 350 mm$:
$\sigma_{xx} = \frac{(a-x) \cdot F}{\frac{1}{12}\cdot bt^3} \cdot \frac{t}{2} = \frac{6 \cdot f \cdot (a-x)}{t^2}$
In [163]:
σ = 6*F*(a-350)/t**2
print("σ max = %6.3f MPa" % σ)
Tip displacement:
$$\delta = \frac{F \cdot a^3}{3 \cdot EI} = \frac{4 \cdot f \cdot l^3}{E \cdot t^3}$$
In [164]:
δ = 4*F*a**3/(E*t**3)
print("δ = %7.2f mm" % δ)
In [165]:
Image(filename='../Lab02_abaqus/'+basename+'/S11.png')
Out[165]:
In [166]:
S11_file = open('../Lab02_abaqus/'+basename+'/'+basename+'_S11.rpt')
In [167]:
S11_array = np.zeros(ny)
nodes = [int(y*nx+7./20.*nx) for y in range(ny)]
#print(nodes)
rows = S11_file.readlines()
i_array = 0
for line in rows:
data = line.split()
if len(data)==4 and data[0].isdigit():
if int(data[0]) in nodes:
S11_array[i_array] = float(data[2])
i_array += 1
#print(S11_array)
In [168]:
S11_array
Out[168]:
In [169]:
Image(filename='../Lab02_abaqus/'+basename+'/Um.png')
Out[169]:
In [ ]:
In [ ]: