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'] = 20
from IPython.display import Image
%matplotlib inline



# Lab 7 - Buckling behavior under biaxial compression

Module for computations



In [2]:

import FEM_utilities as FEM



## Write input file

Set parameters for simulation:

• a, b, h: geometry of plate
• nx, ny, number of elements in x and y directions
• E, $\nu$: elastic properties (Young modulus and Poisson ratio, as TYPE is by default isotropic)
• $\gamma$: scale factor of $\sigma_y$ wrt $\sigma_x$
• num_modes: buckling modes to be computed (by default 1)

Function write_inp_file defines:

• filename : name of input file
• eltype : type of shell element (see conventions here)
• matname: name of material
• elname: *Elset assigned name
• bc, load: names for boundary condition, external and internal nodes loaded
• pa, pb,ax, $\alpha$: coordinate points, axis and angle used to define local coordinate system (see conventions here)
• intpoints: shell integration points


In [3]:

def write_inp_file(a, b, h, nx, ny, E, ν, γ, num_modes=1):
"""
function to write inp file:
a        : lenght of plate
b        : height
h        : thickness
nx       : number of elements in x direction
ny       : number of elements in y direction
E        : Young's modulus
nu       : Poisson ratio
gamma    : sy = gamma*sx
num_modes: number of eigenvalues to be estimated
"""

elname = "EALL"
matname = "mat1"

bcs = ["x0_nodes", "xL_nodes", "y0_nodes", "yW_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 = 'p_nx'+str(nx).zfill(3)+'_ny'+str(ny).zfill(3)+'_g'+str(int(γ*100)).zfill(3)

filename = basename+'.inp'

outfile = open(filename, "wt")

outfile.write("** Lab 07 input file - buckling of plate under biaxial load\n")

nodes, elements, y0, yL, x0, xW = FEM.NodesLab07(nx, ny, a, b)

bc_nodes = [y0, yL, x0, xW]

# 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")

for i, bc in enumerate(bcs):
outfile.write("*NSET, nset = {0}\n".format(bc))

n_lines = np.ceil(len(bc_nodes[i])/16)

for i_line in range(int(n_lines)):
bc_str = "{0:4d}".format(int(bc_nodes[i][16*i_line]))

for j in range(16*i_line+1,min(len(bc_nodes[i]),16*(i_line+1)) ):
bc_str+=",{0:4d}".format(int(bc_nodes[i][j]))
bc_str+="\n"
outfile.write(bc_str)

# il seguente funziona se <= 16
#bc_str = "{0:4d}".format(int(bc_nodes[i][0]))
#for j in range(1,len(bc_nodes[i])):
#     bc_str+=",{0:4d}".format(int(bc_nodes[i][j]))
#bc_str+="\n"
#outfile.write(bc_str)

# SIGMA x
if len(yL > 2):
outfile.write("*Nset, nset = {0}\n".format(sx_int))

yLint = yL[1:-1]

n_lines = np.ceil(len(yLint)/16)

for i_line in range(int(n_lines)):
il_str = "{0:4d}".format(int(yLint[16*i_line]))
for j in range(16*i_line+1,min(len(yLint),16*(i_line+1)) ):
il_str+=",{0:4d}".format(int(yLint[j]))
il_str+="\n"
outfile.write(il_str)

outfile.write("*Nset, nset = {0}\n".format(sx_ext))
el_str = "{0:4d},{1:4d}\n".format(int(yL[0]),int(yL[-1]))
outfile.write(el_str)

# SIGMA y
if len(xW > 2):
outfile.write("*Nset, nset = {0}\n".format(sy_int))

xWint = xW[1:-1]

n_lines = np.ceil(len(xWint)/16)

for i_line in range(int(n_lines)):
il_str = "{0:4d}".format(int(xWint[16*i_line]))
for j in range(16*i_line+1,min(len(xWint),16*(i_line+1)) ):
il_str+=",{0:4d}".format(int(xWint[j]))
il_str+="\n"
outfile.write(il_str)

outfile.write("*Nset, nset = {0}\n".format(sy_ext))
el_str = "{0:4d},{1:4d}\n".format(int(xW[0]),int(xW[-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(h,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")
for bc in bcs:
outfile.write("{0}, 3\n".format(bc))

outfile.write("{0}, 1\n".format(bcs[0]))
outfile.write("{0}, 2\n".format(bcs[2]))

outfile.write("**\n")
outfile.write("**\n")
#outfile.write("{0}, {1:2d}, {2}\n".format(loadnodes, 2, -0.5 ))
if len(bcs[1]) > 2:
outfile.write("{0}, {1:2d}, {2}\n".format(sx_int, 1, -1.0*(b*h)/ny ))
outfile.write("{0}, {1:2d}, {2}\n".format(sx_ext, 1, -0.5*(b*h)/ny ))

if len(bcs[3]) > 2:
outfile.write("{0}, {1:2d}, {2}\n".format(sy_int, 2, -1.0*(a*h)*γ/nx ))
outfile.write("{0}, {1:2d}, {2}\n".format(sy_ext, 2, -0.5*(a*h)*γ/nx ))

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('../Lab07_abaqus/'+basename+'/'):
os.makedirs('../Lab07_abaqus/'+basename+'/')

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



Parameters:



In [4]:

a = 700. # lenght [mm]
b = 350.  # width  [mm]
h = 1.3   # thickness [mm]



Compute node coordinates, elements, constrained and loaded nodes



In [5]:

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

nx = [7,14,14,50,100]
ny = [7,7,14,50,50]

#test many nodes
#nx = [50,100]
#ny = [50,50]

# list of gammas
γ_v = [0.01, 0.05, 0.1, 0.2, 0.25, 0.3, 0.4, 0.5, 0.75, 1.0, 1.5, 2.0, 2.5, 4.0, 6.0, 7.5, 9.0]

#test
#γ_v = [1.0]

if len(nx) == len(ny):
for i in range(len(nx)):
for γ_act in γ_v:
write_inp_file(a, b, h, nx[i], ny[i], E, ν, γ_act)
else:
print("nx and ny vector lenghts don't match")




In [6]:

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("/"):
#relativer 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 [7]:

# commented out as already computed:
#generateCases("../Lab07_abaqus/","p_")



## Analytical solution:

$\sigma_x = \frac{\pi^2}{b^2} \cdot \frac{1}{h} \cdot Dk$

$k = \frac{\left[ \left( \frac{m}{r} \right)^2 + n^2 \right]^2}{\left( \frac{m}{r} \right)^2 + \gamma n^2} \quad m,n = 1 \ldots \infty$

$\sigma_{y} = \gamma \cdot \sigma_{x}$

$r = \frac{a}{b} \quad D = \frac{E h^3}{12\left(1-\nu^2\right)}$



In [8]:

r = a/b

D = E*h**3/(12*(1-ν**2))

γ_a = np.linspace(0.0,10.0, 10000)

σ_x = np.zeros_like(γ_a)
σ_y = np.zeros_like(γ_a)

k_a = np.zeros_like(γ_a)
m_min = np.zeros_like(γ_a)
n_min = np.zeros_like(γ_a)

for i, γ in enumerate(γ_a):
m_min[i] = 0
n_min[i] = 0
k_a[i] = np.inf
for m in range(1,15):
for n in range(1,15):
k_tmp = ((m/r)**2 + n**2)**2/((m/r)**2+γ*n**2)
if k_tmp < k_a[i]:
k_a[i] = k_tmp
m_min[i] = m
n_min[i] = n

σ_x[i] = np.pi**2/b**2*D*k_a[i]/h
σ_y[i] = γ * σ_x[i]



# Simulations Data



In [9]:

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

dirs = glob.glob("../Lab07_abaqus/p_*")

p_eig = {}

for di in dirs:
test = di.split('/')[-1]
if a != -1.0:
if test[2:13] not in p_eig:
p_eig[test[2:13]] = [[],[]]
p_eig[test[2:13]][0].append(a)
p_eig[test[2:13]][1].append(a*float(test[-3:])/100)




In [11]:

plt.figure(figsize=(16,10), dpi=300)
plt.plot(σ_x, σ_y, lw=3., label='analytical')

for act_test in p_eig:
lb = r'$n_x=$'+str(int(act_test[2:5]))+' $n_y=$'+str(int(act_test[8:11]))
plt.plot(p_eig[act_test][0],p_eig[act_test][1],'h',ms=8,label=lb)

#plt.xlim([0.,4])
#plt.ylim([0.0,1.6])
plt.xticks(np.arange(0.,4.1,0.5))
plt.yticks(np.arange(0.0,1.61,.2))
plt.title('Buckling behavior under biaxial load', fontsize=20, fontweight="bold")
plt.xlabel(r'$\sigma_x [MPa]$', fontsize=20)
plt.ylabel(r'$\sigma_y [MPa]$', fontsize=20)
plt.grid()
plt.savefig('Lab07.jpg');







In [12]:

plt.figure(figsize=(16,10), dpi=300)
plt.plot(γ_a, m_min, lw=3.,label='m')
plt.plot(γ_a, n_min, lw=3.,label='n')




Out[12]:

[<matplotlib.lines.Line2D at 0x7fe7f02a4f60>]




In [ ]:




In [ ]: