Hi,
Am currently trying to perform the frequency analysis of cube structure with piezoelectric material. I have the following boundary conditions:
- prescribed displacement - left boundary
- Zero Voltage/Grounded - bottom boundary
and I also want to implement the third BC, which is - Floating potential or the surface should have equal potential distribution on the top boundary such as defined in the figure below:
I have currently implemented the piezoelectric effect with the first and second BCS. Am struggling to implement the floating potential effect on the top boundary marked. I have attached the testcode below. Can anyone who has knowledge of floating potential effect help me or suggest me in implementing it or considering it in the analysis. Looking forward for your suggestions.
from dolfin import *
import matplotlib.pyplot as plt
import numpy as np
#units: m, kg, s, V, K
################################## Mesh part ##################################
mesh = UnitCubeMesh(1,1,1)
Ue = VectorElement("CG", mesh.ufl_cell(), 1) # displacement vector element
Ve = FiniteElement("CG", mesh.ufl_cell(), 1) # voltage finite element
W = FunctionSpace(mesh, MixedElement([Ue, Ve]))
print(f"Number of dofs global: {W.dim()}")
# 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
# IEEE notation: x, y, z, yz, xz, xy
# Elasticity tensor c^E
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
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):
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 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 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)
#Defining the "Trial" functions
soln = Function(W)
(u, v) = split(soln)
#Defining the test functions
(wu, wv) = TestFunctions(W)
# Mark boundary subdomians
class Left(SubDomain):
def inside(self, x, on_boundary):
return near(x[1], 0.0)
class Bottom(SubDomain):
def inside(self, x, on_boundary):
return near(x[2], 1)
class Top(SubDomain):
def inside(self, x, on_boundary):
return near(x[2], 0.0)
# Initialize sub-domain instances
left = Left()
bottom = Bottom()
top = Top()
# Mark subdomains
boundaries = MeshFunction("size_t", mesh, mesh.topology().dim()-1)
boundaries.set_all(0)
left.mark(boundaries, 1)
bottom.mark(boundaries, 2)
top.mark(boundaries, 3)
ds = Measure('ds',subdomain_data=boundaries)
# Dirichlet boundary condition for fixed support and ground voltage
bcu = DirichletBC(W.sub(0), Constant((0,0,1e-3)), left)
bcv = DirichletBC(W.sub(1), Constant(0.), bottom)
bcs = [bcu, bcv]
# The frequency Response
freq=50
omega=2*np.pi*freq
Muu = rho_Piezo*(inner(u,wu))*dx
Kuuv = (inner(sigma(u, v), B(wu)))*dx
Kvvu = (inner(D_field(u, v), grad(wv)))*dx
F = - omega**2*Muu + Kuuv + Kvvu
# Set up the PDE
solve(F == 0, soln, bcs)
U, V = soln.split(deepcopy=True)
print(f"Voltage: {V.vector().get_local()}")