Piezoelectric Analysis


I want to perform the modal analysis of the device with a piezoelectric effect. The governing equations of my model are shown below:

In my model, i would also need to include the “voltage coupling” on the top surface as one is available in ANSYS:

How is it possible to achieve the coupling of DOFs over a surface in FEniCS?, If am not wrong coupling and Constrain are two different effects?
Also, I would like to know the implementation of the Piezoelectric effect if one is already an expert in it?

If you want to do a modal analysis, then you need to start learning dolfinx because there you can setup a problem in complex numbers.

Instead, if you want to solve the coupled equations in time, dolfin is pretty good at this.

You can have a look in #029 at my web site


the first application: 5.1. Piezoelectric transducer and pyroelectric energy harvester is a transient solution of thermomechanics and electromagnetism. The code in FEniCS is also on the web site (tab: supply code, #029).

If something is unclear, simply ask.

Best, Emek

Thanks for the suggestion, I tried to follow your supply code and i have implemented the above equations, But I get an error stating that :

*** Error: Unable to assemble system.
*** Reason: expected a bilinear form for a.
*** Where: This error was encountered inside SystemAssembler.cpp.
*** Process: 0

*** DOLFIN version: 2019.2.0.dev0
*** Git changeset: unknown

I have attached my code below, It would be appreciable to point out the mistake:
from dolfin import *
import matplotlib.pyplot as plt
import numpy as np
import timeit
#units: m, kg, s, V, K

################################## Mesh part ##########################

2D 2 “top”
2D 3 “bottom”
3D 1 “cube”

load the mesh data from Gmsh

mesh_file = ‘cube.h5’
mesh = Mesh()
hdf = HDF5File(mesh.mpi_comm(), mesh_file, ‘r’)
hdf.read(mesh, ‘/mesh’, False)
cells = MeshFunction(‘size_t’, mesh, dim=3)
hdf.read(cells, ‘/subdomains’)
facets = MeshFunction(‘size_t’, mesh, dim=2)
hdf.read(facets, ‘/facets’)

JZ: use the Measure function

dx = Measure(‘dx’, domain=mesh, subdomain_data=cells)
ds = Measure(‘ds’, domain=mesh, subdomain_data=facets)

Function spaces

N = FacetNormal(mesh)
Ue = VectorElement(“CG”, mesh.ufl_cell(), 2) # displacement vector element
Ve = FiniteElement(“CG”, mesh.ufl_cell(), 1) # voltage finite element
W = FunctionSpace(mesh, MixedElement([Ue, Ve]))

Parameters of Piezomaterial (PIC-255)


c_11 = 1.23e11
c_12 = 7.67e10
c_13 = 7.025e10
c_33 = 9.711e10
c_44 = 2.226e10

eps_11 = 8.234
eps_33 = 7.588

e_15 = 11.91
e_13 = -7.15
e_33 = 13.75

rho_Piezo = Constant(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

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., e_15],
[0., 0., 0., 0., e_15, 0.],
[e_13, e_13, e_33, 0., 0., 0.]])

transpose form of piezo tensor

e_piezo = as_tensor([[0., 0., e_13],
[0., 0., e_13],
[0., 0., e_33],
[0., 0, 0.],
[0., e_15, 0.],
[e_15, 0., 0.]])

Permittivitatstensor 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],ten[0,1],ten[1,2],ten[0,2]])

rewrite a vector into a tensor using Voigt notation

def voigt3stress(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[3], vec[5]],
[vec[3], vec[1], vec[4]],
[vec[5], vec[4], 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)
(U, V) = split(w)

#Defining the test functions
(utest, vtest) = TestFunctions(W)

dirichlet boundary condition for fixed support

bcs = [DirichletBC(W.sub(0), Constant((0.,0.,0.)), facets, 2)]

Forms & matrices

U_form = inner(sigma(U, V), B(utest)) * dx
V_form = inner(D_field(U,V), grad(vtest)) * dx
K_form = U_form + V_form
M_form = rho_Piezo*inner(U, utest)*dx
L = utest[0]*dx

K, M = PETScMatrix(), PETScMatrix()
b = PETScVector()
[bc.zero(M) for bc in bcs]


eigensolver = SLEPcEigenSolver(K, M)
eigensolver.parameters[‘problem_type’] = ‘gen_hermitian’
eigensolver.parameters[‘spectral_transform’] = ‘shift-and-invert’
eigensolver.parameters[‘spectral_shift’] = 0.
N_eig = 6 # number of eigenvalues
print(“Computing {} first eigenvalues…”.format(N_eig))

I never use Voigt notation to do calculations. You may try to isolate the term causing the issue. Using so many python def in the UFL is not a good idea at first sight.