Forcing Elastic modulus to zero when sigma_y < sigma_xx within element


85
views
0
4 months ago by
hsk  
Hi, I am solving linear elasticity problem. I want to force a constraint as sigma_xx > sigma_y then elastic modulus within element set to zero for the new step.
 small help can improve my solution.

Thanks in advance

Here I am attaching my code:

File attached: testElasticity.py (1.21 KB)

 

Community: FEniCS Project

1 Answer


3
4 months ago by
pf4d  
Convert your simga_xx to a numpy array, then find the indices where it is less than sigma_y, and finally set the 'E' vector to zero at these indicies:

from fenics import * 

mesh = UnitSquareMesh(2,2)
V    = VectorFunctionSpace(mesh,"Lagrange",1)
P1   = FunctionSpace(mesh, "CG", 1)

# elasticity is now a vector :
E = Function(P1)
E.interpolate(Constant(10.0)) 

# boundaries :
class Bottom(SubDomain): 
  def inside(self,x,on_boundary):
    return (near (x[1],0.0) and on_boundary)
bottom = Bottom()

class Top(SubDomain):
  def inside(self , x, on_boundary ):
    return (near (x[1],1.0) and on_boundary)
top = Top()

# boundary conditions :
bc_fixed = DirichletBC(V, Constant ((0,0)),  bottom)
bcT      = DirichletBC(V.sub(1), Constant(1.0),  top)
bcs      = [bc_fixed, bcT]

# Lame parameters :
nu    = 0.3
lmbda = E*nu/((1.0 + nu )*(1.0-2.0*nu))
mu    = E/(2*(1+nu))

# stress tensor :
def sigma(v): 
  return 2.0*mu*sym(grad(v)) + lmbda*tr(sym(grad(v)))*Identity(2)

# variational form :
u  = TrialFunction(V)
v  = TestFunction(V)
du = Function(V)
a  = inner(grad(v), sigma(u))*dx

# solution :
problem_u = LinearVariationalProblem(lhs(a), rhs(a), du, bcs)
solver_u  = LinearVariationalSolver(problem_u)
solver_u.solve()

point = (1.0,1.0)
print 'numerical u at the corner point:', du(point)

(u1, u2) = du.split(True)
plot(u1, interactive = True)
plot(u2, interactive = True)

# Now I want to set E ==0 when sigma[0][0] > sigma_y
sigma_y     = 0.5
sigma_xx    = project(sigma(du)[0,0], P1)
sigma_xx_a  = sigma_xx.vector().array()
E_a         = E.vector().array()

E_a[sigma_xx_a > sigma_y] = 0.0

E.vector().set_local(E_a)

plot(E, interactive = True)
Please login to add an answer/comment or follow this question.

Similar posts:
Search »