Hi,
I am performing frequency domain analysis of a piezoelectric cube model. The boundary conditions for the mechanical part: prescribed displacement at the left face,
electrical part: grounded or zero voltage at the bottom face, and floating potential(equal voltage at all nodes on this face). I am currently struggling to apply the multipoint constraint or the floating potential BC so that nodes in this region have equal voltage values. I have written the multipoint constraint part in the code below, but I doubt the implementation part to achieve what I wanted. can someone kindly go through the implementation and suggest to me the necessary changes.
import dolfinx
import numpy as np
from mpi4py import MPI
import dolfinx.fem as fem
from petsc4py import PETSc
from ufl import (Circumradius, FacetNormal, SpatialCoordinate, TrialFunction,
TestFunction, div, ds, grad, inner,dot, nabla_grad, nabla_div,Identity,dot, dx)
from dolfinx.cpp.mesh import CellType, locate_entities_boundary
from dolfinx.mesh import MeshTags, locate_entities
from dolfinx.mesh import locate_entities_boundary
from dolfinx.fem import (Constant, DirichletBC, Function,NonlinearProblem, LinearProblem, FunctionSpace, VectorFunctionSpace,
locate_dofs_topological,assemble_scalar)
import ufl
from ufl import (Form, SpatialCoordinate, VectorElement, FiniteElement, TensorElement, MixedElement,as_tensor,as_vector,dot,TestFunction, TrialFunction,TestFunctions, TrialFunctions,dx, ds, grad, FacetNormal,inner, max_value,nabla_grad, nabla_div, Identity)
from scipy.optimize import root
from petsc4py.PETSc import ScalarType
import dolfinx.io
from petsc4py import PETSc
from dolfinx.fem import (assemble_matrix, assemble_vector, form, apply_lifting, locate_dofs_geometrical, set_bc)
import dolfinx_mpc
import dolfinx_mpc.utils
from dolfinx.generation import BoxMesh
from dolfinx.common import Timer, list_timings, TimingType
from dolfinx.io import XDMFFile
L, B, H = 1,1,1
mesh = BoxMesh(MPI.COMM_WORLD, [np.array([0,0,0]), np.array([L, B, H])], [1,1,1], cell_type=CellType.tetrahedron)
U = VectorElement("CG", mesh.ufl_cell(), 1)# # displacement vector element
V = FiniteElement("CG", mesh.ufl_cell(), 1) # voltage finite element
U1, V1 = FunctionSpace(mesh, U), FunctionSpace(mesh, V)
W = FunctionSpace(mesh, MixedElement([U, V]))
num_dofs_global = W.dofmap.index_map.size_global * W.dofmap.index_map_bs
print(f"Number of dofs global: {num_dofs_global}")
# Parameters of Anisotropic Piezomaterial (PIC-255)
c_11 = (1.23e11)
c_12 = (7.67e10)
c_13 = (7.023e10)
c_33 = (9.711e10)
c_44 = (2.226e10)
# Relative permitivity matrix parameters
eps_0 = (8.854e-12) # relative permitivity of free space
eps_11 = (1649 * eps_0)
eps_33 = (1750 * eps_0)
# Coupling matrix parameters
e_13 = (-7.15)
e_33 = (13.75)
e_15 = (11.91)
rho_Piezo = (7800) #in kg/m3
# Order of matrices similar to ANSYS notation: x, y, z, xy, yz,xz
# IEEE notation: x, y, z, yz, xz, xy
# Elasticity tensor c^E
#IEEE Notation
C = as_tensor([[c_11, c_12, c_13, 0., 0., 0.],
[c_12, c_11, c_13, 0., 0., 0.],
[c_13, c_13, c_33, 0., 0., 0.],
[0., 0., 0., c_44, 0., 0],
[0., 0., 0., 0., c_44, 0.],
[0., 0., 0., 0., 0., (c_11-c_12)/2]])
# piezoelectric coupling tensor e
et_piezo = as_tensor([[0.0, 0.0, 0.0, 0.0, e_15, 0.0],
[0.0, 0.0, 0.0, e_15, 0.0, 0.0],
[e_13, e_13, e_33, 0.0, 0.0, 0.0]])
# transpose form of piezo tensor
e_piezo = as_tensor([[0.0, 0.0, e_13],
[0.0, 0.0, e_13],
[0.0, 0.0, e_33],
[0.0, e_15, 0.0],
[e_15, 0.0, 0.0],
[0.0, 0.0, 0.0]])
# Permittivitats tensor epsilon^S
eps_s = as_tensor([[eps_11, 0., 0.],
[0., eps_11, 0.],
[0., 0., eps_33]])
# rewrite the tensor into a vector using Voigt notation
def strain3voigt(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],2*ten[1,2],2*ten[0,2],2*ten[0,1]])
# rewrite a vector into a tensor using Voigt notation
def voigt3stress(vec):
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 B(u):
return 0.5*(nabla_grad(u) + nabla_grad(u).T)
#define the electric field vector
def E_field_vector(Ei):
return as_vector([Ei[0], Ei[1], Ei[2]])
#define the electric field in relation to potential
def E_field(v):
return -grad(v)
def Sigma(u):
return voigt3stress(dot(C, strain3voigt(B(u))))
#definition of the different components of my Stress tensor equation (sigma=C^E:Bu-e_piezo.E)
def sigma(u, v):
return voigt3stress(dot(C, strain3voigt(B(u))) - e_piezo * E_field_vector(E_field(v)))
# the electric displacements vector
def disp_D(Di):
return as_vector([Di[0], Di[1], Di[2]])
# definition of the different components of my electric displacements equation (D=e^s:Bu+eps^S.E)
# elect_disp_D for coupled linear elasicity and electrostatics
def D_field(u, v):
D = et_piezo * strain3voigt(B(u)) + eps_s * E_field(v)
return disp_D(D)
# Mechanical: Prescribed Displacement Dirichlet BC
def clamped_boundary(x):
return np.isclose(x[1], 0)
def f(x):
vals = np.zeros(x.shape)
vals[2] = 1e-3
return vals
fdim = mesh.topology.dim - 1
boundary_facets = locate_entities_boundary(mesh, fdim, clamped_boundary)
u_D = Function(U1)
u_D.interpolate(f)
bc1 = DirichletBC(u_D, locate_dofs_topological((W.sub(0), U1), fdim, boundary_facets), W.sub(0))
# Electrical: Voltage Dirichlet BC
def Ground(x):
return np.isclose(x[2], 1)
ground_facets = locate_entities_boundary(mesh, fdim, Ground)
v_D = Function(V1)
with v_D.vector.localForm() as loc:
loc.set(0)
bc2 = DirichletBC(v_D, locate_dofs_topological((W.sub(1),V1), fdim, ground_facets), W.sub(1))
bcs = [bc1,bc2]
# Floating Potential
def Voltage(x):
return np.isclose(x[2], 0)
v_facets = locate_entities_boundary(mesh, fdim, Voltage)
v_dofs = locate_dofs_topological(V1, mesh.topology.dim-1, v_facets)
def create_rigid_constraint(mpc, dofs, subspace=None):
dof_coordinates = mpc.V.tabulate_dof_coordinates()[dofs, :]
for i in range(len(dof_coordinates) - 1):
mpc.create_general_constraint(
{dof_coordinates[i + 1, :].tobytes(): {dof_coordinates[0, :].tobytes(): 1}},
subspace,
subspace,
)
with Timer("~RIGID: Initialize MPC"):
mpc = dolfinx_mpc.MultiPointConstraint(V1)
create_rigid_constraint(mpc, v_dofs)
mpc.finalize()
# Trial and Test functions
(u, v) = TrialFunctions(W)
(wu, wv) = TestFunctions(W)
T = Constant(mesh, ScalarType((0,0,0)))
Tv = Constant(mesh, ScalarType((0)))
# The frequency domain is : [0,130] with interval of 10.
start = 10; d_freq = 10; freq=40
omega=2*np.pi*freq
# Matrices
Muu = rho_Piezo*(inner(u,wu))*dx
Kuuv = (inner(sigma(u, v), B(wu)))*dx
Kvvu = (inner(D_field(u, v), grad(wv)))*dx
soln = Function(W)
a = - omega**2*Muu + Kuuv + Kvvu
L = inner(T, wu) * ds
# Set up the PDE
problem = dolfinx_mpc.LinearProblem(a, L, mpc, bcs)
soln = problem.solve()
Ur = soln.sub(0).collapse()
Vr = soln.sub(1).collapse()