2
$\begingroup$

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

$$ \nabla \cdot \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:

$$ \int_\Omega \sigma : \nabla \tau = \int_{\partial\Omega} P \cdot \tau \, d s . $$

$\tau$ is a vector test function for the trial function $\vec{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)
$\endgroup$
  • $\begingroup$ Is your question whether imposing the tangential component of the stress is mathematically well-posed (and/or how to formulate a mathematical method to impose it), or is it about how to implement it in FEniCS? $\endgroup$ – Wolfgang Bangerth Oct 9 '15 at 11:40
  • $\begingroup$ My concern is implementing it in FEniCS. I've edited the question to reflect this. $\endgroup$ – NPMitchell Oct 9 '15 at 15:12
  • $\begingroup$ In that case I believe that your question is better suited to the FEniCS help forums or mailing lists. $\endgroup$ – Wolfgang Bangerth Oct 11 '15 at 1:48
  • $\begingroup$ FEniCS has its own Q&A forum: fenicsproject.org/qa (where this question has just been asked, I see: fenicsproject.org/qa/8335/…). $\endgroup$ – Christian Clason Oct 11 '15 at 12:42