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]:



## Write input file

Set parameters for simulation:

• filename : name of input file
• eltype : type of shell element (see conventions here)
• matname: name of material
• E, $\nu$: elastic properties (Young modulus and Poisson ratio, as TYPE is by default isotropic)
• elname: *Elset assigned name
• pa, pb,ax, $\alpha$: coordinate points, axis and angle used to define local coordinate system (see conventions here)
• intpoints: shell integration points
• F: force at tip of plate [N/mm]


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'

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)

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(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))

outfile.write("**\n")
outfile.write("**\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]:

'../Lab02_abaqus/test_40x8y/test_40x8y.inp'



## Analytical solutions:

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" % σ)




σ max = 34.667 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" % δ)




δ =  329.22 mm



# $\sigma_{11} @\ 350mm$



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)

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]:

array([ 35.3893,  35.3446,  35.3096,  35.2899,  35.2899,  35.3096,
35.3446,  35.3893])



# Tip displacement



In [169]:

Image(filename='../Lab02_abaqus/'+basename+'/Um.png')




Out[169]:




In [ ]:




In [ ]: