Im developing a Elasticity model for underground mining. In fisrt intance i created a 2d de model of plane strain elasticity and i can plot the max and min principal stress like:
from dolfin import *
mesh = Mesh("Brillador.xml")
#plot(mesh)
cd=MeshFunction('size_t',mesh,"Brillador_physical_region.xml")
#plot(cd)
Er=30e3
Ev=15e3
class Young(UserExpression): # UserExpression instead of Expression
def __init__(self, markers, **kwargs):
super().__init__(**kwargs) # This part is new!
self.markers = markers
def eval_cell(self, values, x, cell):
if self.markers[cell.index] == 1:
values[0] = Er
elif self.markers[cell.index] ==2:
values[0] = Ev
elif self.markers[cell.index] ==3:
values[0] = Ev
elif self.markers[cell.index] ==4:
values[0] = Ev
elif self.markers[cell.index] ==5:
values[0] = Ev
elif self.markers[cell.index] ==6:
values[0] = Ev
elif self.markers[cell.index] ==7:
values[0] = Er
elif self.markers[cell.index] ==8:
values[0] = Ev
elif self.markers[cell.index] ==9:
values[0] = Ev
else:
values[0] = Er
E=Young(cd,degree=0)
nu = Constant(0.3)
mu = E/2/(1+nu)
lmbda =E*nu/(1+nu)/(1-2*nu)
x_max = 969.3292
def eps(v):
return sym(grad(v))
def sigma(v):
return lmbda*tr(eps(v))*Identity(2) + 2.0*mu*eps(v)
rho_g = 0.027
f = Constant((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 bottom(x, on_boundary):
return near(x[1], 0.)
def left(x, on_boundary):
return near(x[0], 0.)
def right(x, on_boundary):
return near(x[0], x_max)
bc1 = DirichletBC(V, Constant((0.,0.)), bottom)
bc2 = DirichletBC(V.sub(0), Constant(0.), left)
bc3 = DirichletBC(V.sub(0), Constant(0.), right)
bc=[bc1,bc2,bc3]
u = Function(V, name="Displacement")
solve(a == l, u, bc)
Vsig = TensorFunctionSpace(mesh, "DG",degree=1)
sig = Function(Vsig, name="Cauchy_Stress")
sig.assign(project(-sigma(u), Vsig));
sx = sig.sub(0)
sy = sig.sub(3)
txy = sig.sub(1)
def sigma1(sigma_x, sigma_y, tau_xy): # magnitude of first principal stress
return((sigma_x+sigma_y)/2 + sqrt(((sigma_x-sigma_y)/2)**2 + tau_xy**2))
def sigma3(sigma_x, sigma_y, tau_xy): # magnitude of first principal stress
return((sigma_x+sigma_y)/2 - sqrt(((sigma_x-sigma_y)/2)**2 + tau_xy**2));
S1 = sigma1(sx,sy,txy);
S3 = sigma3(sx,sy,txy)
Vp = TensorFunctionSpace(mesh, "DG",degree=1)
sp = Function(Vsig, name="Principal_Stress")
sp.assign(project(as_tensor(((S1,0.),(0.,S3))), Vp));
Now, im wanna do this in a 3d model, but a have not found the way to get the principal stress max and min in this model. I expect you can help me, the code is:
from dolfin import *
mesh = Mesh("Caserones.xml")
plot(mesh)
#Scaled Variables
E = Constant(20e3) #Módulo de Young
nu = Constant(0.3) #Radio de Poisson
mu = E/(2*(1+nu)) #Módulo de cizalle
lambda_ = E*nu/((1+nu)*(1-2*nu)) #Lamé
rho = Constant(0.027) #densidad
E=Constant(40e3)
nu = Constant(0.3)
mu = E/2/(1+nu)
lmbda =E*nu/(1+nu)/(1-2*nu)
#x_max = 969.3292
def eps(v):
return sym(grad(v))
def sigma(v):
return lmbda*tr(eps(v))*Identity(3) + 2.0*mu*eps(v)
rho_g = 0.027
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 bottom(x, on_boundary):
return (on_boundary and near(x[2],-1142.01835021))
bc1 =DirichletBC(V, Constant((0.,0.,0.)),bottom)
def left(x, on_boundary):
return (on_boundary and near(x[0], 0))
bc2 =DirichletBC(V.sub(1), Constant(0.),left)
def right(x, on_boundary):
return (on_boundary and near(x[0], 1216.8837))
bc3 =DirichletBC(V.sub(1), Constant(0.),right)
def north(x, on_boundary):
return (on_boundary and near(x[1], 0))
bc4 =DirichletBC(V.sub(0), Constant(0.),north)
def south(x, on_boundary):
return (on_boundary and near(x[1], 929.7089))
bc5 =DirichletBC(V.sub(0), Constant(0.),south)
bc=[bc1,bc2,bc3,bc4,bc5]
#Solve
import matplotlib.pyplot as plt
u = Function(V)
solve(a==l, u, bc);
Vsig = TensorFunctionSpace(mesh, "DG",degree=1)
sig = Function(Vsig, name="Cauchy_Stress")
sig.assign(project(-sigma(u), Vsig));