Pure Neumann boundary conditions in mixed formulation of 2d elasticity

Question

I am trying to implement a 2d elasticity problem using mixed formulation, since I want to simulate a material with a Poisson's ratio of 0.45 or larger. I would like to impose a constant stress at the boundary of the square. For this, I try to transform the BoundarySource class from here to the 2d problem.

However, I do not understand the meaning/order of the "values" array. Using the current configuration, I would expect a constant stress of 1.0 at one of the boundaries, but I get a strange periodic displacement pattern. What am I doing wrong and where does the pattern come from? Do I use the correct function spaces for this issue? I would be very grateful, if you could help me.

from dolfin import *

# Create mesh
mesh = UnitSquareMesh(32, 32)

# Define function spaces and mixed (product) space
TCG = TensorFunctionSpace(mesh, "CG", 2)
VCG = VectorFunctionSpace(mesh, "CG", 2)
W = TCG * VCG

nu = 0.10
E = 1.0
mu = Constant(E / (2.0 * (1.0 + nu)))
lmbda = Constant(E * nu / ((1.0 + nu) * (1.0 - 2.0 * nu)))

def compliance(sigma):
    return 1.0 / 2.0 * mu * (sigma - lmbda / (2 * mu + 2 * lmbda) * tr(sigma) * Identity(2))

# Define trial and test functions
(sigma, u) = TrialFunctions(W)
(tau, v) = TestFunctions(W)

# Define source function
f = Constant((0.0, 0.0))

# Define variational form
a = (inner(compliance(sigma), tau) + inner(div(tau), u) + inner(div(sigma), v))*dx
L = - inner(f, v)*dx# + inner(u0, dot(n, tau))*ds(1) + inner(u0, dot(n, tau))*ds(2)

# Define function G such that G \cdot n = g
class BoundarySource(Expression):
    def __init__(self, mesh):
        self.mesh = mesh
    def eval_cell(self, values, x, ufc_cell):
        cell = Cell(self.mesh, ufc_cell.index)
        n = cell.normal(ufc_cell.local_facet)
        values[0] = 1.0
        values[1] = 0.0
        values[2] = 0.0
        values[3] = 0.0
    def value_shape(self):
        return (2, 2)

G = BoundarySource(mesh)

# Define essential boundary
def boundary(x):
    return x[1] < DOLFIN_EPS or x[1] > 1.0 - DOLFIN_EPS

bc = DirichletBC(W.sub(0), G, boundary)

# Compute solution
w = Function(W)
solve(a == L, w, bc)
(sigma, u) = w.split()

# Plot sigma and u
plot(u)
interactive()

Comments

Hi, I would be surprised if CG2-CG2 element satisfies the inf-sup condition in your mixed problem. Perhaps element like Taylor-Hood (sigma in CG2 and u in CG1) is better. As for the boundary condition; values defines in a given point a 2x2 matrix, say G, such that G.n is the vector you want to prescribe on the boundary. In the following G is always an identity

class BoundarySource(Expression):
    def __init__(self, mesh):
        self.mesh = mesh
    def eval_cell(self, values, x, ufc_cell):
        cell = Cell(self.mesh, ufc_cell.index)
        n = cell.normal(ufc_cell.local_facet)
        values[0] = 1.0 
        values[1] = 0.0
        values[2] = 0.0
        values[3] = 1.0
    def value_shape(self):
        return (2, 2)

so the surface stress is 1*n.

---

Hi MiroK,

wow, the example works now and is definitely something I can build upon!

Thanks!