Hello,
I am trying to solve for stresses in a domain consisting of two (connected) subdomains with different material properties (See figure below). The top subdomain material is linear elastic while the bottom subdomain material experiences plastic deformation and stresses.
Currently I am able to solve for stresses in the full domain assuming both subdomain materials to be linear elastic. The working code for the same is attached below:
Stresses in full domain assuming both subdomain materials linear elastic:
from __future__ import print_function
from dolfin import *
import matplotlib.pyplot as plt
from matplotlib import cm
from fenics import *
from mshr import *
import numpy as np
from ufl import nabla_div,sinh,exp
import random
############################# Create main mesh ##################################
mesh = RectangleMesh(Point(0.0,-100e-6),Point(100e-6,100e-6),100,100)
dFE = FiniteElement("DG", mesh.ufl_cell(), 0)
FE = FiniteElement("Lagrange", mesh.ufl_cell() ,1)
vFE = VectorElement(FE)
tFE = TensorElement(dFE)
U = FunctionSpace(mesh, dFE)
V = FunctionSpace(mesh, vFE)
W = FunctionSpace(mesh, tFE)
#################### Subdomains & Material Properties #########################
# Top subdomain
E0 = Constant(11.0e3)
nu0 = Constant(0.29)
mu_0 = E0/(2*(1+nu0))
lambda_0 = E0*nu0/((1+nu0)*(1-2*nu0))
sig0_0 = 400
# Bottom subdomain
E1 = Constant(9.0e3)
nu1 = Constant(0.42)
mu_1 = E1/(2*(1+nu1))
lambda_1 = E1*nu1/((1+nu1)*(1-2*nu1))
sig0_1 = 0.4
tol = 1E-14
class Omega_0(SubDomain):
def inside(self, x, on_boundary):
return x[1] >= 0.0 - tol
class Omega_1(SubDomain):
def inside(self, x, on_boundary):
return x[1] < 0.0 + tol
materials = MeshFunction('size_t',mesh,2)
subdomain_0 = Omega_0()
subdomain_1 = Omega_1()
materials.set_all(0)
subdomain_1.mark(materials, 1)
class M(UserExpression):
def __init__(self, materials, mu_0, mu_1, **kwargs):
super().__init__(**kwargs)
self.materials = materials
self.mu_0 = mu_0
self.mu_1 = mu_1
def eval_cell(self, values, x, cell):
if self.materials[cell.index] == 0:
values[0] = self.mu_0
else:
values[0] = self.mu_1
mu = M(materials, mu_0, mu_1, degree=0)
class L(UserExpression):
def __init__(self, materials, lambda_0, lambda_1, **kwargs):
super().__init__(**kwargs)
self.materials = materials
self.lambda_0 = lambda_0
self.lambda_1 = lambda_1
def eval_cell(self, values, x, cell):
if self.materials[cell.index] == 0:
values[0] = self.lambda_0
else:
values[0] = self.lambda_1
lambda_ = L(materials, lambda_0, lambda_1, degree=0)
class Y(UserExpression):
def __init__(self, materials, sig0_0, sig0_1, **kwargs):
super().__init__(**kwargs)
self.materials = materials
self.sig0_0 = sig0_0
self.sig0_1 = sig0_1
def eval_cell(self, values, x, cell):
if self.materials[cell.index] == 0:
values[0] = self.sig0_0
else:
values[0] = self.sig0_1
sig0 = Y(materials, sig0_0, sig0_1, degree=0)
class El(UserExpression):
def __init__(self, materials, E0, E1, **kwargs):
super().__init__(**kwargs)
self.materials = materials
self.E0 = E0
self.E1 = E1
def eval_cell(self, values, x, cell):
if self.materials[cell.index] == 0:
values[0] = self.E0
else:
values[0] = self.E1
E = El(materials, E0, E1, degree=0)
########################### Boundary Conditions ############################
class TopBoundary(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and near(x[1], 100e-6, tol)
class BottomBoundary(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and near(x[1],-100e-6, tol)
class Interface(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and near (x[1], 0.0, tol) and not near(x[0], 0.) and not near(x[0], 100e-6)
class LeftBoundary(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and near(x[0], 0.0, tol)
class RightBoundary(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and near(x[0], 100e-6, tol)
top_boundary = TopBoundary()
bottom_boundary = BottomBoundary()
interface = Interface()
left_boundary = LeftBoundary()
right_boundary = RightBoundary()
boundary_parts = MeshFunction("size_t", mesh, mesh.topology().dim()-1)
boundary_parts.set_all(2)
top_boundary.mark(boundary_parts, 1)
bottom_boundary.mark(boundary_parts, 3)
left_boundary.mark(boundary_parts, 4)
right_boundary.mark(boundary_parts, 4)
ds = Measure("ds")[boundary_parts]
#Left boundary
bcL = DirichletBC(V.sub(0), Constant(0.), left_boundary)
#Right boundary
bcR = DirichletBC(V.sub(0), Constant(0.), right_boundary)
#Bottom boundary
bcBottom = DirichletBC(V, Constant((0.,0.)), bottom_boundary)
bcs = [bcL, bcR, bcBottom]
######################## Stress & Strain Operators #############################
def eps(u):
return sym(grad(u))
def sigma(u):
return lambda_*tr(eps(u))*Identity(2) + 2.0*mu*eps(u)
############### Solving stress field in full domain #################
u = TrialFunction(V)
v = TestFunction(V)
f = Constant((0, 0))
pressure = Constant((0, -10.0))
a = inner(sigma(u), eps(v))*dx
L = dot(f, v)*dx + dot(pressure, v)*ds(1)
u = Function(V)
solve(a == L, u, bcs)
sig = sigma(u)
fileResults = XDMFFile("Stress_linearelastic.xdmf")
fileResults.parameters["flush_output"] = True
fileResults.parameters["functions_share_mesh"] = True
disp = Function(V, name='Displacement')
disp.assign(u)
fileResults.write(disp,0.)
stress = Function(W, name='Stress')
stress.assign(project(sig, W))
fileResults.write(stress,0.)
I am also using the following reference to solve for plastic stresses: https://comet-fenics.readthedocs.io/en/latest/demo/2D_plasticity/vonMises_plasticity.py.html
However, I am facing difficulty in figuring out how to solve for the full domain stress field with linear elastic stresses in the top subdomain and plastic stresses in the bottom subdomain. Would using a mixed element formulation work? Any leads on this would be much appreciated. Thanks!