Dear everyone:
I am doing an electromechanical coupling analysis of a piezoelectric material beam.
- Boundary conditions for the Mechanical part: prescribed displacement 0 at the left face.
- Boundary conditions for the Electrical part: zero electric potential at the left surface, and floating potential on the right surface (equal voltage at all nodes on this face represent the function of electrodes).
I want to apply the multipoint constraint for the right surface so that nodes in this surface have equal electric potentials.
The way I tried is to pick out a corner point on the right surface as the reference point and make every point inside this surface equal to this reference point and Lagrange multiplier method is adopted.
I have accomplished using Lagrange multiplier method to restrict the surface potential to a known value. But I don’t know how to set the value equal but not a specific constant.
I attached the code below. And it is not working for defining point this way. Could someone give some suggestions about how to realize this target in dolfin2019.2 developed version?
I really appreciate any help you can provide.
from __future__ import print_function
from fenics import *
from dolfin import *
from ufl import nabla_div
import numpy
# Scaled variables
L = 1
W = 0.2 # L is length; W is width
rho = 5550
g = 9.8
# define elasticity tensor C
c11 = 226e9
c12 = 125e9
c13 = 124e9
c33 = 216e9
c44 = 44.2e9 # * 2
c66 = (c11 - c12) / 2 # * 2
# 6x6
C = as_matrix([[c11, c12, c13, 0.0, 0.0, 0.0],
[c12, c11, c13, 0.0, 0.0, 0.0],
[c13, c13, c33, 0.0, 0.0, 0.0],
[0.0, 0.0, 0.0, c44, 0.0, 0.0],
[0.0, 0.0, 0.0, 0.0, c44, 0.0],
[0.0, 0.0, 0.0, 0.0, 0.0, c66]])
# define piezoelectricity tensor e
e15 = 5.8
e31 = -2.2
e33 = 9.3
# 6x3
e1 = as_matrix([[0.0, 0.0, e31],
[0.0, 0.0, e31],
[0.0, 0.0, e33],
[0.0, e15, 0.0],
[e15, 0.0, 0.0],
[0.0, 0.0, 0.0]])
# define permittivity tensor p
p11 = 5.64e-9
p33 = 6.35e-9
# 3x3
p = as_matrix([[p11, 0.0, 0.0],
[0.0, p11, 0.0],
[0.0, 0.0, p33]])
# 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]])
return as_vector([ten[0, 0], ten[1, 1], ten[2, 2], 2 * ten[1, 2], 2 * ten[0, 2], 2 * ten[0, 1]])
def as_voigt_vector2(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], ten[1], ten[2]])
# 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]]])
def from_voigt_vector2(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[1], vec[2]])
# first define the function for the strain tensor in fenics notation
def epsilon(u):
return 0.5 * (nabla_grad(u) + nabla_grad(u).T)
# return sym(nabla_grad(u))
# then define the function for the static electric field in fenics notation
def e_field(phi):
return -nabla_grad(phi)
# now define the stress tensor in fencis notation
def sigma(u, phi):
# calculate sigma with C-matrix in Voigt notation
sig = C * (as_voigt_vector(epsilon(u))) - e1 * as_voigt_vector2(e_field(phi))
# return sigma as symmetric tensor
return from_voigt_vector(sig)
# then define the function for the electric displacement in fenics notation
def e_dis(u, phi):
ed = e1.T * (as_voigt_vector(epsilon(u))) + p * as_voigt_vector2(e_field(phi))
return from_voigt_vector2(ed)
# Create mesh and define function space
mesh = BoxMesh(Point(0, 0, 0), Point(1, 0.2, 0.2), 50, 10, 10)
class RightBoundary(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and near(x[0], 1.0)
class Referencepoint(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and near(x[0], 1.0) and near(x[1], 0.2) and near(x[2], 0.2)
# Marking facets that are on RightBoundary and creating the corresponding submesh
marker = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
marker.set_all(0)
RightBoundary().mark(marker, 1)
submesh = MeshView.create(marker, 1)
pointmarker = MeshFunction("size_t", mesh, mesh.topology().dim() - 2)
Referencepoint().mark(pointmarker, 1)
point = MeshView.create(pointmarker, 1)
V = VectorFunctionSpace(mesh, 'CG', 1)
Q = FunctionSpace(mesh, 'CG', 1)
LM = FunctionSpace(submesh, "CG", 1)
RP = FunctionSpace(point, 'CG', 1)
Z = MixedFunctionSpace(V, Q, LM, RP) #displacement, electric potential, lagrange multiplier, reference point
sol = Function(Z)
(u, phi, L_M, R_P) = TrialFunctions(Z)
(v, q, l_m, r_p) = TestFunctions(Z)
# Define boundary condition
tol = 1E-14
def clamped_boundary(x, on_boundary):
return on_boundary and x[0] < tol
bc1 = DirichletBC(V, Constant((0, 0, 0)), clamped_boundary)
bc2 = DirichletBC(Q, Constant(0), clamped_boundary)
bc = [bc1, bc2]
# Define strain and stress
# Define variational problem
dV = Measure("dx", domain=Z.sub_space(0).mesh())
dV2 = Measure("dx", domain=Z.sub_space(1).mesh())
dL = Measure("dx", domain=Z.sub_space(2).mesh())
d = u.geometric_dimension() # space dimension
f = Constant((0, 0, -rho * g))
T = Constant((0, 0, 0))
cl = Constant(500)
a = inner(sigma(u, phi), epsilon(v)) * dV - dot(e_dis(u, phi), e_field(q)) * dV2
a += L_M * (q - r_p) * dL + l_m * (phi - R_P) * dL
# a += L_M*nabla_grad(q)[1]*dL + l_m*nabla_grad(phi)[1]*dL + L_M*nabla_grad(q)[2]*dL + l_m*nabla_grad(phi)[2]*dL
L = dot(f, v) * dV + dot(T, v) * ds
# solve(a == L, sol, bc)
solve(a == L, sol, bc, solver_parameters={"linear_solver": "direct"})
# solve(a == L, sol, bc, solver_parameters={"linear_solver":"minres"})
u = sol.sub(0)
phi = sol.sub(1)
L_M = sol.sub(2)