 # Piezoelectric Analysis

Hi,

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

http://bilenemek.abali.org/

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 ##########################

“”"
\$PhysicalNames
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’)
cells = MeshFunction(‘size_t’, mesh, dim=3)
facets = MeshFunction(‘size_t’, mesh, dim=2)

# 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
#u,phi
W = FunctionSpace(mesh, MixedElement([Ue, Ve]))

# --------------------

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

# 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, vec, vec],
[vec, vec, vec],
[vec, vec, vec]])

# first define the function for the strain tensor in fenics notation

def B(u):

#define the electric field vector
def E_field_vector(Ei):
return as_vector([Ei, Ei, Ei])

#define the electric field in relation to potential
def E_field(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, Di, Di])

# 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*dx

K, M = PETScMatrix(), PETScMatrix()
b = PETScVector()
assemble_system(K_form,L,bcs,A_tensor=K)
assemble_system(M_form,L,bcs,A_tensor=M)
[bc.zero(M) for bc in bcs]

# --------------------Eigen-solver--------------------

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))
eigensolver.solve(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.