Hello,
I am trying to modify the “2D linear elasticity” problem linked:
https://comet-fenics.readthedocs.io/en/latest/demo/elasticity/2D_elasticity.py.html
to make it “1D linear elasticity” problem.
When I compare the 1D displacement extracted from the 2D problem with the displacement obtained directly from 1D problem, there is difference in shape and magnitude.
Guidance in this regard is required.
Below is python script.
Blockquote
from dolfin import *
import matplotlib.pyplot as plt
import numpy as np
2D linear elasticity
L = 25.0
H = 1.0
Nx = 250
Ny = 10
E = Constant(1e5)
nu = Constant(0.3)
mesh = RectangleMesh(Point(0.0, 0.0), Point(L, H), Nx, Ny, “crossed”)
def eps(v):
return sym(grad(v))
model = “plane_stress”
mu = E/2/(1 + nu)
lmbda = Enu/(1 + nu)/(1 - 2nu)
if model == “plane_stress”:
lmbda = 2mulmbda/(lmbda + 2*mu)
def sigma(v):
return lmbdatr(eps(v))Identity(2) + 2.0mueps(v)
rho_g = 1e-3
f = Constant((0.0, -rho_g))
V = VectorFunctionSpace(mesh, ‘Lagrange’, degree=2)
du = TrialFunction(V)
u_ = TestFunction(V)
a = inner(sigma(du), eps(u_))*dx
l = inner(f, u_)*dx
def left(x, on_boundary):
return near(x[0], 0.0)
bc = DirichletBC(V, Constant((0.0, 0.0)), left)
u = Function(V, name =“Displacement”)
solve(a == l, u, bc)
plot(u, mode=“displacement”)
plt.grid(color = ‘k’, linestyle = ‘-’, linewidth = 0.1)
plt.title(“2D Displacemnt”)
plt.xlabel(‘x’)
plt.ylabel(‘Displacement’)
plt.show()
extraction of 1D value from 2D
y_position = H/2
x_min = 0
x_max = L
N_steps = 1000
delta_x = (x_max - x_min)/float(N_steps)
u_vector = np.zeros(N_steps,float)
for i in range(N_steps):
x_position = float(i * delta_x)
u_vector[i] = u(x_position, y_position)[1]
x_values = [i * delta_x for i in range(N_steps)]
plt.xlabel(‘x’)
plt.ylabel(‘Displacement’)
plt.plot(x_values, u_vector)
plt.title(“1D Displacemnt”)
plt.grid(color = ‘k’, linestyle = ‘-’, linewidth = 0.1)
plt.show()
print(“Maximal deflection:”, -u(L, H/2.0)[1])
print(“Beam theory deflection:”, float(3rho_gL4/2/E/H3))
1D linear elasticity
L = 25.0
H = 1.0
Nx = 250
Ny = 10
E = Constant(1e5)
nu = Constant(0.3)
mesh = IntervalMesh(Nx, 0.0, L)
def eps(v):
return sym(grad(v))
model = “plane_stress”
mu = E / (2*(1 + nu))
lmbda = (Enu) / ((1 + nu)(1 - 2*nu))
if model == “plane_stress”:
lmbda = 2mulmbda/(lmbda + 2*mu)
def sigma(v):
return lmbdatr(eps(v))Identity(1) + 2.0mueps(v)
rho_g = 1e-3
f = Constant((-rho_g,))
V = VectorFunctionSpace(mesh, ‘Lagrange’, degree=2)
du = TrialFunction(V)
u_ = TestFunction(V)
a = inner(sigma(du), eps(u_))*dx
l = inner(f, u_)*dx
def left(x, on_boundary):
return near(x[0], 0.0)
bc = DirichletBC(V, Constant((0.0,)), left)
u = Function(V, name =“Displacement”)
solve(a == l, u, bc)
z = np.linspace(0, 25, 100)
u_line = np.array([u(point) for point in z])
plt.plot(z, u_line)
plt.title(“1D Displacemnt”)
plt.grid(color=‘r’, linestyle=‘-’, linewidth=0.1)
plt.legend([‘Displacement’], markerfirst = False, loc=‘best’)
plt.ylabel(‘Displacement’)
plt.show()
print(“Maximal deflection:”, -u(L/2))