In [1]:
import numpy as np
import os
import glob
import shutil
import matplotlib.pyplot as plt
import json
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
Set parameters for simulation:
In [3]:
def write_inp_file(L, z0, nn, eltype, elname, isRiks, matname, A, E, ν, increment, t_data, F):
n_coords = np.zeros((nn,2))
for ni in range(1,nn):
n_coords[ni,:] = ni/(nn-1)*np.array([L, z0])
if isRiks:
basename = 't_ne'+str(nn-1).zfill(2)+'_aR_s'+str(int(t_data[3]*100)).zfill(3)
else:
basename = 't_ne'+str(nn-1).zfill(2)+'_aN_s'+str(int(t_data[3]*100)).zfill(3)
filename = basename+'.inp'
outfile = open(filename, "wt")
outfile.write("** Lab 05 input file test\n")
# NODES section
outfile.write("**\n")
outfile.write("** Nodes\n")
outfile.write("**\n")
outfile.write("*NODE\n")
for i in range(nn):
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(eltype,elname))
for i in range(1,nn):
outfile.write("{0:4d},{1:4d},{2:4d}\n".format(i,i,i+1))
# SOLID section
outfile.write("**\n")
outfile.write("** Solid section\n")
outfile.write("**\n")
outfile.write("*SOLID SECTION, MATERIAL={0}, ELSET={1}\n".format(matname,elname))
outfile.write(str(A)+",\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,\t3\n")
outfile.write("{0},\t1\n".format(nn))
outfile.write("{0},\t3\n".format(nn))
# calculation steps
outfile.write("**\n")
outfile.write("** Step\n")
outfile.write("**\n")
outfile.write("*STEP, NLGEOM, INC={0}\n".format(increment))
if isRiks:
outfile.write("*STATIC, RIKS\n")
else:
outfile.write("*STATIC\n")
outfile.write("{0:8},{1:8},{2:8},{3:8}\n".format(t_data[0], t_data[1], t_data[2], t_data[3]))
# LOADS
outfile.write("**\n")
outfile.write("** Loads\n")
outfile.write("**\n")
outfile.write("*Cload\n")
outfile.write("{0}, {1:2d}, {2}\n".format(nn, 2, -F ))
outfile.write("*OUTPUT,FIELD\n")
outfile.write("*ELEMENT OUTPUT\n")
outfile.write("S,COORD,\n")
outfile.write("*EL PRINT\n")
outfile.write(" S,COORD,\n")
outfile.write("*OUTPUT,FIELD\n")
outfile.write("*NODE OUTPUT\n")
outfile.write("U,COORD\n")
outfile.write("*NODE PRINT\n")
outfile.write("U,COORD\n")
outfile.write("*OUTPUT,HISTORY,FREQUENCY={0}\n".format(increment))
outfile.write("*END STEP\n")
outfile.close()
if not os.path.exists('../Lab05_abaqus/'+basename+'/'):
os.makedirs('../Lab05_abaqus/'+basename+'/')
shutil.move(filename,'../Lab05_abaqus/'+basename+'/'+filename)
In [4]:
L = 2500. # lenght [mm]
z0 = 25. # initial height [mm]
isRiks = True
eltype = 'T3D2H'
elname = 'EALL'
matname = 'material_1'
E = 72000. # modulus [MPa]
ν = 0.33 # Poisson's coefficient
EA = 5e7 # cross sectional properties [N]
A = EA/E # area of cross section [mm^2]
# step parameters
increment = 1000 # max number of calculation steps
# dt0, t_tot, dt_min, dt_max
t_data1 = [0.05, 1.0, 1e-4, 0.05]
t_data2 = [0.25, 1.0, 0.25, 0.25]
ts = [t_data1, t_data2]
# load
F = 15 # load [N]
num_nodes = [2, 3, 6]
for k in num_nodes:
# shall we perform Riks analysis
for isRiks in [True, False]:
# time steps
for td in ts:
write_inp_file(L, z0, k, eltype, elname, isRiks, matname, A, E, ν, increment, td, F)
In [5]:
current_dir = os.getcwd()
os.chdir("../Lab05_abaqus/")
In [6]:
# %%bash
#./lab05.py
In [7]:
os.chdir(current_dir)
$F-w$ analytical relation:
$F = -\frac{EA}{L^3} \cdot ( z^2 \cdot w + \frac{3}{2}z \cdot w^2 + \frac{1}{2}w^3 )$
In [5]:
w_z_max = 2.5
w_z = np.linspace(0,w_z_max,500)
y = w_z - 1.5*w_z**2 + 0.5 * w_z**3
In [6]:
def readDataL05(filename, isRiks):
file=open(filename,'r')
row = file.readlines()
step = []
U2 = []
S11 = []
state = 0
for line in row:
strlist = line.split()
if 'INCREMENT' in strlist and 'SUMMARY' in strlist:
state = 1
elif 'U2' in strlist:
state = 2
elif 'S11' in strlist:
state = 3
elif 'TOTAL' in strlist and state == 1 and not isRiks:
#print(strlist)
step.append(float(strlist[-1]))
state = 0
elif 'FACTOR' in strlist and state == 1 and isRiks:
step.append(float(strlist[4]))
state = 0
elif 'MINIMUM' in strlist and state == 2:
U2.append(float(strlist[2]))
state = 0
elif 'MAXIMUM' == strlist and state == 3:
S11.append(float(strlist[2]))
state = 0
return np.array(step), np.array(U2), np.array(S11)
In [7]:
dirs = glob.glob("../Lab05_abaqus/t_*")
sim_data = {}
for di in dirs:
test = di.split('/')[-1]
if '_aR_' in test:
a, b, c = readDataL05(di+"/"+test+".dat", True)
else:
a, b, c = readDataL05(di+"/"+test+".dat", False)
xy = np.zeros((len(a),2))
xy[:,0] = -b/z0
xy[:,1] = F*a*L**3/(EA*z0**3)
sim_data[test] = xy
In [8]:
#sim_data
In [13]:
w_z_max = 2.5
plt.figure(figsize=(16,10), dpi=300)
plt.plot(w_z,y, '-', lw=3, label='analytical')
for i, test in enumerate(sim_data):
n_el = int(test[4:6])
if '_aR_' in test:
atype = 'Riks'
if n_el == 1:
sym = 'h'
c = ((i+1)/len(sim_data),0.0,0.0)
elif n_el == 2:
sym = '*'
c = ((i+1)/len(sim_data),(i+1)/len(sim_data),0.0)
else:
sym = '^'
c = (0.0,(i+1)/len(sim_data),0.0)
if float(test[-3:])/100 == 0.25:
msize = 12
else:
msize = 6
lb = r"$n_{el}=$"+str(n_el)+" $\Delta l_{max}$="+str(float(test[-3:])/100)+" "+atype
plt.plot(sim_data[test][:,0],sim_data[test][:,1],sym, ms=msize, color=c, label=lb )
plt.xlim([0,w_z_max])
plt.ylim([-0.2,.4])
plt.xticks(np.arange(0.0,w_z_max+.1,0.25))
plt.yticks(np.arange(-0.2,0.4+.1,.1))
plt.title('Force displacement relation - Riks',fontsize=16)
plt.xlabel(r'$-\frac{w}{z}$', fontsize=15)
plt.ylabel(r'$\frac{FL^3}{EAz^3}$', fontsize=15)
plt.legend(loc='upper center', shadow=True, fontsize=14)
plt.grid()
plt.savefig('Lab05_Riks.jpg')
In [17]:
w_z_max =0.5
plt.figure(figsize=(16,10), dpi=300)
plt.plot(w_z,y, '-', lw=3, label='analytical')
for i, test in enumerate(sim_data):
n_el = int(test[4:6])
if '_aN_' in test:
atype = 'NR'
if n_el == 1:
sym = 'h'
c = ((i+1)/len(sim_data),0.0,0.0)
elif n_el == 2:
sym = '*'
c = ((i+1)/len(sim_data),(i+1)/len(sim_data),0.0)
else:
sym = '^'
c = (0.0,(i+1)/len(sim_data),(i+1)/len(sim_data))
if float(test[-3:])/100 == 0.25:
msize = 12
else:
msize = 6
lb = r"$n_{el}=$"+str(n_el)+" $\Delta l_{max}$="+str(float(test[-3:])/100)+" "+atype
plt.plot(sim_data[test][:,0],sim_data[test][:,1],sym, ms=msize, color=c, label=lb )
plt.xlim([0,w_z_max])
plt.ylim([-0.2,.4])
plt.xticks(np.arange(0.0,w_z_max+.1,0.25))
plt.yticks(np.arange(-0.2,0.4+.1,.1))
plt.title('Force displacement relation - Newton Raphson',fontsize=16)
plt.xlabel(r'$-\frac{w}{z}$',fontsize=15)
plt.ylabel(r'$\frac{FL^3}{EAz^3}$',fontsize=15)
plt.legend(loc='upper center', shadow=True,fontsize=14)
plt.grid()
plt.savefig('Lab05_NR.jpg')
In [ ]: