both normal and tangential stress boundary conditions for elasticity?

Question

I would like to have both the normal (xx) component and shear (xy) component of a 2D (stress) tensor defined on a boundary (x=const, for instance) for an equation which is of the type

div(sigma) = 0

The normal value of sigma_{xx}=P at x=0,1 can be accounted for via a Neumann boundary condition, in the term obtained by integration by parts in weak form:

inner(sigma, sym(grad(tau))) *dx =   dot(P, tau) *ds

tau is a vector test function for the trial function u, which physically is the displacement.
However, I am having difficulty fixing the shear stress at the boundaries x=0,1 of the unit square, as a simple example of my goal. How can one best accomplish this in FEniCS?

My question is similar to this one, but as far as I can tell, they are distinct.

A minimal running code in FEniCS (in python) which gives incorrect results for sigma is below. Thank you!

from dolfin import *

# Mesh Params ################
E, nu = 1.0, 0.4 #material params
mesh = UnitSquareMesh(100, 100)

print 'Creating Lagrangian function spaces...'
Vv = VectorFunctionSpace(mesh, "Lagrange", 1)
Vf = FunctionSpace(mesh, "Lagrange", 1)
Vt = TensorFunctionSpace(mesh, "Lagrange", 1)

# P is stress vector of stress values normal and tangent to boundary 
P = Expression(('(x[0]-xc)','(x[0]-xc)'),xc=0.5)

class LeftBoundary(SubDomain):
    def inside(self, x, on_boundary):
        tol = 1E-7
        return on_boundary and abs(x[0]) < tol

class RightBoundary(SubDomain):
    def inside(self, x, on_boundary):
        tol = 1E-7
        return on_boundary and abs(x[0] - 1) < tol

left_boundary = LeftBoundary()        
right_boundary = RightBoundary()        

boundary_parts = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
boundary_parts.set_all(0)
right_boundary.mark(boundary_parts, 1)
left_boundary.mark(boundary_parts, 2)
ds = Measure("ds")[boundary_parts]

##############################
# Define variational problem #
##############################
print 'Defining test and trial functions...'
v = TestFunction(Vf)
tau = TestFunction(Vv)
phi = TrialFunction(Vf)
u = TrialFunction(Vv)

# functions for variational problem
d = u.geometric_dimension()
epsilon = sym(nabla_grad(u))
# PLANE STRESS
sigma = E/(1+nu)* (epsilon) + E*nu/((1-nu**2)) * tr(epsilon)*Identity(d)

# Define normal and tangent directions on mesh
norm = FacetNormal(mesh)
tang = as_vector([norm[1], -norm[0]])

tau_n = dot(tau,norm)
tau_t = dot(tau,tang)
P_n = dot(P,norm)
P_t = dot(P,tang)

#################
# SOLVE
#################
au = inner(sigma, sym(nabla_grad(tau)))*dx 
Lu = dot(P_n,tau_n)*ds(1) + dot(P_n,tau_n)*ds(2) + dot(P_t,tau_t)*ds(1) + dot(P_t,tau_t)*ds(2)
u = Function(Vv)
solve(au == Lu, u)

## View the results
sigma = E/(1+nu)* (sym(nabla_grad(u))) + E*nu/((1-nu**2)) * tr(sym(nabla_grad(u)))*Identity(d)
s_V = project( sigma, Vt)
sxx = project(s_V[0,0], Vf); 
sxy = project(s_V[0,1], Vf); 
syx = project(s_V[1,0], Vf);  
syy = project(s_V[1,1], Vf);
plot(sxx, interactive=True, rescale=False,scalarbar=True)
plot(syy, interactive=True, rescale=False,scalarbar=True)
plot(sxy, interactive=True, rescale=False,scalarbar=True)

JeonghunLee · Answer 1 · Oct 4, 2016

Hi.

This question was proposed a long time ago though, let me answer it.

I think the main reason is that your problem does not have boundary conditions on the top and bottom boundaries. You need boundary conditions of displacement to have a well-posed problem. Computer codes usually don't complain it but just give a wrong answer.

And your traction boundary condition can be simply imposed by dot(P, tau)*ds(1), instead of putting normal and tangential components separately.

both normal and tangential stress boundary conditions for elasticity?

1 Answer