Buckling in Compression Model (Linear Elastic)

1

While I’m studying the 2D linear elasticity model (with some small changes in the geometry), I’ve tried to recalculate the Young’s Modulus, as can be seen at the end of the code. But, I couldn’t get the same elastic modulus that I’ve determined as a material property at the beginning of the code -which may be expected. In general, there are some points needed to be clear such as:

  1. How should the finite element type and degree be chosen during the projection of the stress (DG or P, 0 or 1)?
  2. To find the average max stress and strain on the part, is assembling the stress and strain vectors the appropriate way?
  3. In the tutorial of the example, is is said that it is assumed there is no buckling: but, it seems there is a buckling in the displacement in the X direction. So, why is this happening?
from fenics import *
import matplotlib.pyplot as plt
import numpy as np

def bottom(x, on_boundary):
    return (on_boundary and near(x[1], 0.0))

# Strain function
def epsilon(u):
    return sym(grad(u))

# Stress function
def sigma(u):
    return lambda_*div(u)*Identity(2) + 2*mu*epsilon(u)

# Density (kg/m3)
rho = Constant(1240.0)

# Young's modulus (Pa) and Poisson's ratio
E = 33e6 # from Essentium TPU74D TDS
nu = 0.39 #from https://www.hindawi.com/journals/amse/2019/5838361/

# Lame's constants
lambda_ = E*nu/(1+nu)/(1-2*nu)
mu = E/2/(1+nu)

l_x, l_y = 5.0, 5.0  # Domain dimensions
n_x, n_y = 20, 20  # Number of elements

# Load
g_z = -200000
b_z = -10
g = Constant((0.0, g_z))
b = Constant((0.0, b_z))

# Model type
model = "plane_strain"
if model == "plane_stress":
    lambda_ = 2*mu*lambda_/(lambda_+2*mu)

mesh = RectangleMesh(Point(0.0, 0.0), Point(l_x, l_y), n_x, n_y)

# Definition of Neumann boundary condition domain
boundaries = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
boundaries.set_all(0)
# Neumann boundary conditions (top)
top = AutoSubDomain(lambda x: near(x[1], l_y))

top.mark(boundaries, 1)
ds = ds(subdomain_data=boundaries)

V = VectorFunctionSpace(mesh, "CG", 1) # displacement function space
u_tr = TrialFunction(V)
u_test = TestFunction(V)

bc = DirichletBC(V, Constant((0.0, 0.0)), bottom)

a = inner(sigma(u_tr), epsilon(u_test))*dx
l = rho*dot(b, u_test)*dx + inner(g, u_test)*ds(1)

u = Function(V)
A_ass, L_ass = assemble_system(a, l, bc)
solve(A_ass, u.vector(), L_ass)
print(np.amax(u.vector()[:]))

# --------------------
# Plot Displacement in Y dir.
# --------------------
plt.figure()
c = plot(u.sub(1))
plt.colorbar(c)
plt.title('Displacement in y-direction')
plt.show()

# --------------------
# Plot Displacement in X dir.
# --------------------
plt.figure()
c = plot(u.sub(0))
plt.colorbar(c)
plt.title('Displacement in x-direction')
plt.show()

stress = sigma(u)
# --------------------
# Project Stress
# --------------------
def local_project(fce, space):
    lp_trial, lp_test = TrialFunction(space), TestFunction(space)
    lp_a = inner(lp_trial, lp_test)*dx
    lp_L = inner(fce, lp_test)*dx
    local_solver = LocalSolver(lp_a, lp_L)
    local_solver.factorize() # caches the lhs
    lp_f = Function(space)
    local_solver.solve_local_rhs(lp_f)
    return lp_f

V0 = TensorFunctionSpace(mesh, "P", 1)
stress_2 = local_project(stress, V0)

strain = epsilon(u)

# --------------------
# Calculate the average strain and stress
# --------------------
strain_xx = assemble(strain[0,0]*dx)

strain_yy = assemble(strain[1,1]*dx)

stress_xx = assemble(stress_2[0,0]*dx)

stress_yy = assemble(stress_2[1,1]*dx)

YM = round(np.abs(stress_yy)/np.abs(strain_yy)/1e6, 2) # in MPa

disp_x