Gravitational Neumann BC

1

I want to create a Neuann bondaty condition in a model from a underground mining cavity. This conditions represent a compress stress condition which depend from the gravity. An illustrative figure from my model:

targets_3
targets_33102×2222 244 KB

The rigt and left edge have a Dirichlet boundary condition, this is 0 cisplacement in y x axis and free in y.

How i can create this neumann bc in mi model?? The code is:

from dolfin import *
mesh = Mesh("p1_2vetas.xml")
#plot(mesh)
cd=MeshFunction('size_t',mesh,"p1_2vetas_physical_region.xml")
#plot(cd)

Ev = 15 #Young Veta
Em = 30 #Young Macizo

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] = Em
        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] = Em
        elif self.markers[cell.index] ==5:
            values[0] = Ev
        elif self.markers[cell.index] ==6:
            values[0] = Ev
        else:
            values[0] = Em
E=Young(cd,degree=0)
nu = Constant(0.3)
mu = E/2/(1+nu)
lmbda =E*nu/(1+nu)/(1-2*nu)
rho_g = 0.027

def eps(v):
  return sym(grad(v))

def sigma(v):
  return lmbda*tr(eps(v))*Identity(2) + 2.0*mu*eps(v)

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;

x_max = 914.5695 
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)