Linear elasticity using periodic boundary conditions

[dthillaithevan]

Hey,

I have been attempting to use PeriodicSquareMesh to solve a linear elastic problem to compute displacments. My goal is to impose periodic boundary conditions (PBCs) on every exterior surface, apply global strains and then compute the resulting displacements (which should be periodic across every surface).

With direction='both' the solution is appears incorrect with a constant displacement throughout the domain. However, with direction='x' or direction='y' the displacements seem more reasonable (please see attached screenshots of displacements). However, I would like to impose the PBCs on every exterior surface so I believe I need to use the direction='both' option, does anyone have any ideas as to what could be going wrong here? Is there any reason to expect the mathematical formulation to differ from FEniCS?

Thank you!

from firedrake import *

mesh = PeriodicSquareMesh(10,10,1, direction='both', quadrilateral = False) #Constant displacements when direction = 'both'
V = VectorFunctionSpace(mesh, "CG", 1)
U = FunctionSpace(mesh, 'CG', 1)
u = TrialFunction(V)
v = TestFunction(V)# Identity tensor
Id = Identity(mesh.geometric_dimension())
lambda_ = Constant(0.25)
mu = Constant(1)

def epsilon(a):
    return sym(grad(a))

def sigma(a,lmbda_func,mu_func,Eps):
    return lmbda_func*tr(epsilon(a) + Eps)*Id + 2*mu_func*(epsilon(a)+Eps) #Impose additional global strain, Eps

Eps = as_vector([[0,0.5],[0.5,0]]) # Apply unit shear strain
dv = TrialFunction(V)
v_ = TestFunction(V)
v = Function(V)

F = inner(sigma(dv,lambda_,mu, Eps), epsilon(v_))*dx
a = lhs(F)
L = rhs(F)
problem = LinearVariationalProblem(a, L, v, [])
solver = LinearVariationalSolver(problem)
solver.solve()

print (v.vector()[:]) #Constant when direction='both'

direction='y' solution

direction='x' solution

direction='both' solution

[dham]

I think your problem is ill-posed. The elasticity operator with no Dirichlet conditions has a non-trivial nullspace which you'd need to remove. See: https://firedrakeproject.org/solving-interface.html#solving-singular-systems

[dthillaithevan]

Thanks for your quick reply David, I was not aware of this problem.

As a initial starting point I assumed the nullspaces to be the 2 translation vectors and the 1 rotation vector, but this does not seem to be correct. Do you have any suggestions or advice for determining the correct nullspace?

x, y = SpatialCoordinate(mesh)
b0.interpolate(Constant([1, 0]))
b1.interpolate(Constant([0, 1]))
b2.interpolate(as_vector([-y, x]))
nullspace = VectorSpaceBasis([b0, b1, b2])
# nullspace = VectorSpaceBasis(constant=True)
nullspace.orthonormalize()
Eps = as_vector([[0,0.5],[0.5,0]]) # Apply unit x strain
dv = TrialFunction(V)
v_ = TestFunction(V)
v = Function(V)
F = inner(sigma(dv,lambda_,mu, Eps), epsilon(v_))*dx
a = lhs(F)
L = rhs(F)
problem = LinearVariationalProblem(a, L, v, [])
solver = LinearVariationalSolver(problem, nullspace= nullspace)
solver.solve()

[ksagiyam]

Also, I think your residual strain Eps is incompatible with the doubly periodic boundary. Please note that even if strain is periodic, displacement in general is not periodic: we have, e.g., u(y=1) = u(y=0) + cy where cy is some vector constant; consider the following simple shear case:

  ______
 /     /
/_____/

In the above simple example, you might be able to map your problem as ubar = u - cy * y and solve for ubar on the periodic mesh, though.

[dthillaithevan]

Sure, here you go.

import ufl
from firedrake import *

n = 20
mesh = PeriodicSquareMesh(n,n,1, direction='both', quadrilateral = True)
CG = VectorFunctionSpace(mesh, "CG", 1)
elem = ufl.VectorElement("Q", mesh.ufl_cell(), 1, variant="equispaced")
DG = FunctionSpace(mesh, elem)

Eps = as_vector([[0,0.5],[0.5,0]])

outfile = File("DG_macro_EPS.pvd")
outfile.write(project((dot(Eps, SpatialCoordinate(mesh))), DG))

outfile = File("CG_macro_EPS.pvd")
outfile.write(project((dot(Eps, SpatialCoordinate(mesh))), CG))

[ksagiyam]

Sorry again, please use "DQ" instead of "Q".

[dthillaithevan]

That works, thanks a lot for your help Koki, appreciate it!

[dthillaithevan]

Sorry again, please use "DQ" instead of "Q".

Hi Again Koki,

Is there an extension of this code to 3D?

elem = ufl.VectorElement("Q", mesh.ufl_cell(), 1, variant="equispaced")
DG = FunctionSpace(mesh, elem) 

At the moment I get the following error when I run this code on a 3D periodic mesh:

  File "/Users/dee/firedrake/src/ufl/ufl/log.py", line 158, in error
    raise self._exception_type(self._format_raw(*message))

UFLException: Cellname "tetrahedron" invalid for "DQ" finite element.

Thanks!

[ksagiyam]

PeriodicCubeMesh and friends currently do not support hex cells and so use tetrahedral cells, in which case you want "DP" (or "DG") instead of "DQ".