### Forcing Elastic modulus to zero when sigma_y < sigma_xx within element

195

views

0

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)

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

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.