Hi,

Thanks for the explanation. I tried implementing the floating potential BC on the multi-physical piezo model.

The governing equations of the piezo model are shown in the image below:

Currently am currently trying to conduct the stationary analysis of cube structure with piezoelectric material. I have the following boundary conditions:

- fixed displacement - bottom 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 as defined in the figure above.

I have currently implemented the piezoelectric effect with the first and second BCS. Am struggling to implement the floating potential impact on the top boundary marked. I have attached the test code below. I have tried to implement the floating potential effect as you suggested for the previous case. But I notice that the results are not matching with the test case. Looking forward to 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
V0e = FiniteElement("Real", mesh.ufl_cell(),0)
W = FunctionSpace(mesh, MixedElement([Ue, Ve, V0e]))
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)
# 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)
#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)
#Defining the "Trial" functions
#w = Function(W)
soln = Function(W)
(u, v, v0) = split(soln)
#Defining the test functions
(wu, wv, wv0) = TestFunctions(W)
# Mark boundary subdomians
# Create classes for defining parts of the boundaries and the interior of the domain
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
bottom = Bottom()
top = Top()
# Mark subdomains
boundaries = MeshFunction("size_t", mesh, mesh.topology().dim()-1)
boundaries.set_all(0)
bottom.mark(boundaries, 2) # Displacement and Voltage BC
top.mark(boundaries, 3) # Floating potential
ds = Measure('ds',subdomain_data=boundaries)
# Dirichlet boundary condition for fixed support and ground voltage
bcu = DirichletBC(W.sub(0), Constant((0,0,0)), bottom)
bcv1 = DirichletBC(W.sub(1), Constant(0), bottom)
bcs = [bcu, bcv1]
# Stationary analysis
n = FacetNormal(mesh)
Q0 = 1e-6 #C/m^2
meas_Gamma1 = Constant(1) # area of the edge or surface of current flow
D_avg = Constant(-Q0/meas_Gamma1)
Kuuv = inner(sigma(u,v), B(wu))*dx
Kvvu = inner(D_field(u,v),grad(wv))*dx - (dot(D_field(u,v),n)*wv + dot(D_field(wu,wv),n)*(v-v0))*ds(3) + (D_avg - dot (D_field(u,v),n))*wv0*ds(3)
F = Kuuv + Kvvu
# Set up the PDE
solve(F == 0, soln, bcs)
U, V, V0 = soln.split(deepcopy=True)
print(max(abs(U.vector().get_local())))
print(f"Voltage: {V.vector().get_local()}")
```

Results obtained from Fenics:

Displacement: 1.960682369212692e-08 m

Voltage: [ 0. 49.66534269 0. 0. 0. 46.62686436

38.87742427 46.62686436]

Correct results from commercial software:

Displacement: 2.0059095750665296E-8 m

Voltage: [ 0. 48.28872311077018 0. 0. 0. 48.28872311074099

48.28872311077021 48.28872311077022]