hello,I was trying to do a simulation of a biweekly compression test (2D), and the part that added confining pressure to the left and right sides of the model stumped me. I don’t know how to properly add these pressures, I tried to write a part, but it did not work
from dolfin import *
L = 0.3
H = 0.6
Nx = 50
Ny = 25
mesh = RectangleMesh(Point(0., 0.), Point(L, H), Nx, Ny, "crossed")
def eps(v):
return sym(grad(v))
E = Constant(1e5)
nu = Constant(0.3)
model = "plane_stress"
mu = E/2/(1+nu)
lmbda = E*nu/(1+nu)/(1-2*nu)
if model == "plane_stress":
lmbda = 2*mu*lmbda/(lmbda+2*mu)
def sigma(v):
return lmbda*tr(eps(v))*Identity(2) + 2.0*mu*eps(v)
rho_g = 10.
f = Constant((0, -rho_g))
class Top(SubDomain):
def inside(self,x,on_boundary):
return near(x[1],H) and on_boundary
class Left(SubDomain):
def inside(self,x,on_boundary):
return near(x[0],0) and on_boundary
class Right(SubDomain):
def inside(self,x,on_boundary):
return near(x[0],0.3) and on_boundary
facets = MeshFunction('size_t',mesh,1)
Top().mark(facets,1)
Left().mark((facets,2))
Right().mark(facets,3)
ds = Measure('ds',subdomain_data=facets)
V = VectorFunctionSpace(mesh, 'Lagrange', degree=2)
#Bottom fixing
def bottom(x,on_boundary):
return near(x[1],0) and on_boundary
bc = DirichletBC(V,Constant((0,0,0)),bottom)
T1 = Constant(0.,-5000.) #the loading pressure
T2 = Constant(6000.,0.) #the confine pressure of the left boundary
T3 = Constant(-6000.,0.) #the confine pressure of the right boundary
du = TrialFunction(V)
u_ = TestFunction(V)
a = inner(sigma(du), eps(u_))*dx
l = inner(f, u_)*dx + dot(T1,u_)*ds(1)+dot(T2,u_)*ds(2)+dot()+dot(T3,u_)*ds(3)
u = Function(V, name="Displacement")
solve(a == l, u, bc)