# Piezoelectric Analysis

**URL:** <https://fenicsproject.discourse.group/t/piezoelectric-analysis/6760>\
**Category:** Uncategorized\
**Created:** [October 7, 2021, 2:17am UTC](https://fenicsproject.discourse.group/t/piezoelectric-analysis/6760 "2021-10-07T02:17:07Z")\
**Posts on this page:** 8\
**Page:** 1

<div class="post-metadata">

**Author:** ![Gunasheela](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/gunasheela/32/2573_2.png) [@Gunasheela](https://fenicsproject.discourse.group/u/Gunasheela)\
**Post date:** [October 7, 2021, 2:17am UTC](https://fenicsproject.discourse.group/t/piezoelectric-analysis/6760/1 "2021-10-07T02:17:07Z")

</div>

Hi,

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

 ![image](https://global.discourse-cdn.com/free1/uploads/fenicsproject1/original/2X/d/d05c3390d43640717bd8dabd469279799ce70f21.png)

In my model, i would also need to include the “voltage coupling” on the top surface as one is available in ANSYS:  
 ![Capture](https://global.discourse-cdn.com/free1/uploads/fenicsproject1/original/2X/0/042bacdcf00106329b4a35076da1e996162e2f79.png)

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?

---

<div class="post-metadata">

**Author:** ![bilenemek](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/bilenemek/32/32_2.png) [@bilenemek](https://fenicsproject.discourse.group/u/bilenemek)\
**Post date:** [October 7, 2021, 1:32pm UTC](https://fenicsproject.discourse.group/t/piezoelectric-analysis/6760/2 "2021-10-07T13:32:35Z")

</div>

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/](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

---

<div class="post-metadata">

**Author:** ![Gunasheela](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/gunasheela/32/2573_2.png) [@Gunasheela](https://fenicsproject.discourse.group/u/Gunasheela)\
**Post date:** [October 11, 2021, 9:48pm UTC](https://fenicsproject.discourse.group/t/piezoelectric-analysis/6760/3 "2021-10-11T21:48:32Z")

</div>

## 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’)  
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  
#u,phi  
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()  
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)

---

<div class="post-metadata">

**Author:** ![bilenemek](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/bilenemek/32/32_2.png) [@bilenemek](https://fenicsproject.discourse.group/u/bilenemek)\
**Post date:** [October 12, 2021, 11:41am UTC](https://fenicsproject.discourse.group/t/piezoelectric-analysis/6760/4 "2021-10-12T11:41:21Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![Gunasheela](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/gunasheela/32/2573_2.png) [@Gunasheela](https://fenicsproject.discourse.group/u/Gunasheela)\
**Post date:** [November 16, 2021, 1:22am UTC](https://fenicsproject.discourse.group/t/piezoelectric-analysis/6760/5 "2021-11-16T01:22:23Z")

</div>

Dear Sir,

I was able to simulate the piezo model successfully. Now I also want to include the floating potential on one of the surfaces, which is similar to voltage coupling in ansys. How can that be implemented ?.The equation is:  
∇(ϵ∇φ)=0  
, and boundary conditions φ=0 on Γ0 φ=1 on Γ1 and φ=const on Γ2, we want to find the value of this constant as a result of solution of the equation.

---

<div class="post-metadata">

**Author:** ![bilenemek](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/bilenemek/32/32_2.png) [@bilenemek](https://fenicsproject.discourse.group/u/bilenemek)\
**Post date:** [November 16, 2021, 9:33am UTC](https://fenicsproject.discourse.group/t/piezoelectric-analysis/6760/6 "2021-11-16T09:33:42Z")

</div>

Hi there,

I guess you want to have div( epsilon grad(phi) ) = 0 on the boundary, instead. If so, it is the Neumann boundary condition and you already get it (since it is zero), whenever you skip implementing a Dirichlet on one boundary. Of course for other boundaries like gamma0 and 1,2 you can apply phi=0,1,const by using DirichletBC() function in FEniCS, but I am sure you already did this.

If I could not understand what you were asking for, please pardon me, and rewrite your question.

Best, Emek

---

<div class="post-metadata">

**Author:** ![Gunasheela](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/gunasheela/32/2573_2.png) [@Gunasheela](https://fenicsproject.discourse.group/u/Gunasheela)\
**Post date:** [November 16, 2021, 7:07pm UTC](https://fenicsproject.discourse.group/t/piezoelectric-analysis/6760/7 "2021-11-16T19:07:24Z")

</div>

Dear Sir,

I need to implement the effect, which is briefed in the link below:

[https://fenicsproject.org/qa/7536/specifying-constant-potential-boundary-condition/](https://fenicsproject.org/qa/7536/specifying-constant-potential-boundary-condition/)

Am not sure, whether the Neumann boundary condition implies the above effect.

---

<div class="post-metadata">

**Author:** ![bilenemek](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/bilenemek/32/32_2.png) [@bilenemek](https://fenicsproject.discourse.group/u/bilenemek)\
**Post date:** [November 17, 2021, 9:33am UTC](https://fenicsproject.discourse.group/t/piezoelectric-analysis/6760/8 "2021-11-17T09:33:22Z")

</div>

Hi Gunasheela

I read the link and did not understand. The equation there is unknown to me, since it means grad( epsilon grad(phi) ) = 0 meaning

\nabla ( \epsilon \nabla \phi ) = 0 \ ,

equivalently

( \epsilon \phi\_{,i} )\_{,j} = 0 \ .

and this is a matrix not a scalar. Here comma means a partial space derivative, you use .dx in FEniCS for that. In this case, I do not know how to derive a weak form.

Therefore, I have assumed that you want to write div( epsilon grad(phi) ) = 0, means

\nabla \cdot (\epsilon \nabla \phi ) = 0 \ , \\ (\epsilon \phi\_{,i} )\_{,i} = 0

I wonder if you simply want to solve a surface, on top of that a gel kind of electrolytic paste and then a metal plate to control the voltage? Then simply do this and solve everywhere, what you have been solving in the above code (you use v that is \phi), but with different permittivity values. I guess this is meant by a floating electrode.

Best, Emek