Posting
Dear Fenics users, I am trying to assembly a linear system. It should combine linear elasticity and electrostatics to simulate piezoelectric materials.
I've first implemented the equations for the linear elasticity problem, which worked. For the next step i wanted to calculate the effect of a constant electric field.
As you can see i defined it as a 6x1 constant vector of zeros. When i substract that Constant from the defined sigma (as you cann see in the definition of sigma) i recieve the following error.
ArityMismatchTraceback (most recent call last)
in () 86 # Solve the problem 87 result = Function(VV) ---> 88 solve(a == L, result, boundary) 89 90 # Save solution to VTK format
...
ArityMismatch: Adding expressions with non-matching form arguments (Argument(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 6664), VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3)), 1, None),) vs ().
Does somebody has any idea what i'm doing wrong? Without subtracting E it works.
from dolfin import *
dielectric = 28000. * 8.8541878176e-12 * Identity(3)
# derived lame constants
Ey = 100.0 # youngs modulus
nu = 0.3 # poisson ratio
mu = Ey/(2.0*(1.0 + nu))
lmbda = Ey*nu/((1.0 + nu)*(1.0 - 2.0*nu))
# define elasticity tensor
e1 = (1-nu)/(1.0-2.0*nu)
e2 = nu/(1.0-2.0*nu)
# 6x6
c = as_matrix([[e1, e2, e2, 0., 0., 0.],
[e2, e1, e2, 0., 0., 0.],
[e2, e2, e1, 0., 0., 0.],
[0., 0., 0., 1., 0., 0],
[0., 0., 0., 0., 1., 0.],
[0., 0., 0., 0., 0., 1.]])
# rewrite the tensor into a vector using Voigt notation
def as_voigt_vector(ten):
# FEniCS does not know anything about Voigt notation,
# so one need to access the components directly as eps[0, 0] etc.
return as_vector([ten[0,0],ten[1,1],ten[2,2],ten[1,2],ten[0,2],ten[0,1]])
# rewrite a vector into a tensor using Voigt notation
def from_voigt_vector(vec):
# FEniCS does not know anything about Voigt notation,
# so one need to rewrite the vector by hand into a tensor.
return as_tensor([[vec[0], vec[5], vec[4]],
[vec[5], vec[1], vec[3]],
[vec[4], vec[3], vec[2]]])
# first define the function for the strain tensor in fenics notation
def epsilon(u):
# note that it would be shorter to use sym(grad(u)) but the way used now
# represents the formula
return 1./2.*(grad(u) + transpose(grad(u)))
# now define the stress tensor in fencis notation
def sigma(u,E):
# calcuate sigma with C-matrix in Voigt notation
# sigma = 6x1 vector; E = 6x1 vector
sig = (Ey/(1.+nu))*c* (as_voigt_vector(epsilon(u)) - E(u))
# sigma for coupled linear elasicity and electrostatics; d = 3x6 tensor with piezo material constants
#sig = (Ey/(1.+nu))*c* (as_voigt_vector(epsilon(u))-(transpose(d_piezo)*E(u)))
# return sigma as symmetric tensor
return from_voigt_vector(sig)
import numpy
#boundary conditions:
voltage = 10
q = Constant(0)
f = Constant((0., 0., 0.))
t = Constant((0., 0., -1.))
# build the trial and test function spaces over a domain
# Setup domain
grid = 10
Omega = UnitCubeMesh(grid, grid, 5*grid)
x = Omega.coordinates()
x[:, 0] *= 0.1
x[:, 1] *= 0.1
x[:, 2] *= 0.01
Omega.bounding_box_tree().build(Omega)
VV = VectorFunctionSpace(Omega, "Lagrange", 1)
# Get the needed boundaries:
# 1. To apply a predefined stress on it We keep it simple: the top face sees the same stress vector.
# 2. For zero Displacement we choose the bottom face of the cube
class TopBoundary(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and near(x[2], 0.01)
class BottomBoundary(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and near(x[2], 0.)
top_boundary = TopBoundary()
bottom_boundary = BottomBoundary()
# Create the boundary condition for displcement
boundary = []
boundary.append(DirichletBC(VV, Constant((0., 0., 0.)), bottom_boundary))
# As the stress is defined over the outer surface we need to create a
# integral over that surface
boundaries = FacetFunction("size_t", Omega)
boundaries.set_all(0)
top_boundary.mark(boundaries, 1)
ds = Measure("ds", subdomain_data=boundaries)
E = lambda u: as_vector(Constant((0.,0.,0.,0.,0.,0.)))
# Define the variational problem: get all variables and implement function given above
u = TrialFunction(VV)
v = TestFunction(VV)
a = inner(sigma(u,E), grad(v))*dx # typically called a
L = inner(f, v)*dx + inner(t, v)*ds(1) # typically called L
# Solve the problem
result = Function(VV)
solve(a == L, result, boundary)
# Save solution to VTK format
File("displacement.pvd") << result
print "Done"
