13. Piezoelectric actuation#
Astrid Pechstein
A circular piezoelectric patch actuator is applied to a square aluminum plate. The plate is clamped at one of its lateral faces, the system is free of external mechanical loads. Both patch and plate can be considered as thin structures, the plate being of dimensions \(25\) mm \(\times 25\) mm \(\times 1\) mm, the patch with diameter \(15\) mm and thickness \(0.5\) mm.
See also: Meindlhumer, M., & Pechstein, A. (2018). 3D mixed finite elements for curved, flat piezoelectric structures. International Journal of Smart and Nano Materials.
#%%
from ngsolve import *
import numpy as np
import netgen.occ as ngocc
from ngsolve.webgui import Draw
from netgen.meshing import MeshingParameters
from netgen.meshing import IdentificationType
SetNumThreads(4)
ea = { "euler_angles" : [-50,3,-30] }
#patch geometry
rad=15/2 #radius of patch
hp=0.5 #height of patch
#substrate geometry (rectangle)
len_a = 25 #width
he=1 #height of substrate
BCDisp="fix"
BCPotential="el1|el2"
rect = ngocc.WorkPlane().Rectangle(len_a, len_a).Face()
circ = ngocc.WorkPlane().Circle(len_a/2, len_a/2,rad).Face()
circ.maxh=5
circ.edges.name = "circ"
circ.edges.maxh = 2
plate1 = (rect-circ).Extrude(-he)
plate2 = (circ).Extrude(-he)
patch = (circ).Extrude(hp)
plate1.faces.Min(ngocc.X).name = "fix"
patch.faces.Min(ngocc.Z).name = "el1"
patch.faces.Max(ngocc.Z).name = "el2"
patch.name = "PZT"
plate1.name = "alu"
plate2.name = "alu"
plate1.faces.Min(ngocc.Z).Identify(plate1.faces.Max(ngocc.Z), name="plate_identification", type=IdentificationType.CLOSESURFACES)
plate2.faces.Min(ngocc.Z).Identify(plate2.faces.Max(ngocc.Z), name="plate_identification", type=IdentificationType.CLOSESURFACES)
patch.faces.Min(ngocc.Z).Identify(patch.faces.Max(ngocc.Z), name="patch_identification", type=IdentificationType.CLOSESURFACES)
setup = ngocc.Glue([patch, plate1,plate2])
setup.solids[1].name = "alu"
geo = ngocc.OCCGeometry(setup)
ngmesh = geo.GenerateMesh(MeshingParameters(maxh=5, grading=0.1))
for fd in ngmesh.FaceDescriptors():
# print (fd, fd.bcname)
if fd.bcname == "circ":
fd.domin_singular = True
fd.domout_singular = True
mesh = Mesh(ngmesh)
# mesh.RefineHP(levels=1, factor=0.2)
Draw(mesh, **ea)
BaseWebGuiScene
13.1. Problem setup#
The patch is made from the piezoelectric ceramic PZT-5H and polarized in thickness direction. Both surfaces of the patch are electroded; as an electric voltage is applied to the electrodes, the patch deforms. For an electric field aligned in reverse direction to the polarization direction, an in-plane stretch of the patch is expected.
For a voltage of \(V = 100\) V, a linear mathematical model of the setup is sufficient to capture its behavior.
13.1.1. Unknown fields and constitutive equations#
In linear piezoelasticity, we are concerned with finding two independent unknown fields: the displacement \(\vec u\) and the electric potential \(\varphi\). Strain and electric field are related via
In the current application, we are interested in computing displacements \(\vec u\) of both plate and patch actor, while the electric potential is only relevant in the piezoelectric domain.
The equation of motion presents a relation for the stress tensor,
while Gauss’ law must hold for the dielectric displacement vector with no free charges in the non-conducting piezoelectric domain,
On the boundaries, either the body is fixed (\(\vec u = \vec 0\)) or free of surface stresses (\(\boldsymbol{\sigma} = \boldsymbol{0}\)). The patch is electroded at two of its faces, with voltage prescribed (\(\varphi = V\)) whereas otherwise it is assumed perfectly insulated (\(\vec D \cdot \vec n = 0\)), which is justified as the electric permittivity of surrounding air is negligibly small as compared to the permittivity of piezo-ceramics.
Finally, a linear material law is available for PZT-5H,
Above, \(\mathbf{C}^E\) is a tensor of fourth order, containing stiffnesses at constant electric field. The third order tensor \(\mathbf e\) contains the piezoelectric constants of the piezoelectric material, and \(\boldsymbol{\epsilon}^S\) specifies the material’s electric permittivity. These tensors are in general not isotropic, but show dependence on direction of polarization, which is \(\vec e_z\) in our case. Commonly, Voigt notation is used in implementations, reducing symmetric tensors \(\boldsymbol{\sigma}\) and \(\boldsymbol{\varepsilon}\) to six-dimensional vectors. We comply with this standard, using
For this representation, the tensors/matrices are listed below. Scalings are such that they match computations in mm and V.
# stiffness at constant electric field [N/mm^2]
CEmat = np.array([[127.205, 80.2, 84.6, 0, 0, 0],[80.2, 127.205, 80.2, 0, 0, 0], [84.6, 80.2, 117., 0, 0, 0], \
[0, 0, 0, 23., 0, 0], [0, 0, 0, 0, 23., 0], [0, 0, 0, 0, 0, 23.4]])*1e3
# electric permittivity tensor [C/Vm]
epsilonSmat = np.array( [[15.05, 0, 0],[ 0, 15.05, 0],[ 0, 0, 13.01]])*1e-9
# piezoelectric permittivity [kC/m^2]
emat = np.array([[0,0,0,0, 17.03,0],[0,0,0,17.03,0,0],[-6.62,-6.62,23.24,0,0,0]])*1e-3
voltage = 100
For aluminum, only the stiffness tensor is relevant. For an isotropic material, elastic modulus \(E\) and Poisson ratio \(\nu\) are sufficient.
Emod = 65e3
nu = 0.3
Smat_alu = 1/Emod*np.array([[1, -nu, -nu, 0, 0, 0],
[-nu, 1, -nu, 0, 0, 0],
[-nu, -nu, 1, 0, 0, 0],
[0, 0, 0, 2+2*nu, 0, 0],
[0, 0, 0, 0, 2+2*nu, 0],
[0, 0, 0, 0, 0, 2+2*nu]])
Cmat_alu = np.linalg.inv(Smat_alu)
cf_CE = mesh.MaterialCF({"PZT": CoefficientFunction(tuple(CEmat.reshape(1, -1)[0]), dims=(6,6)), "alu": CoefficientFunction(tuple(Cmat_alu.reshape(1, -1)[0]), dims=(6,6))})
cf_e = CoefficientFunction(tuple(emat.reshape(1, -1)[0]), dims=(3,6))
cf_epsilonS = CoefficientFunction(tuple(epsilonSmat.reshape(1, -1)[0]), dims=(3,3))
cf_rho = mesh.MaterialCF({"alu": 2.5e-9, "PZT": 7.5e-9})
The coupled problem is set up. A compound finite element space, collecting displacement and electric potential, is defined. A GridFunction q collects \(\vec u\) and \(\varphi\).
k_prim = 3
mesh.Curve(k_prim)
fespace_u_prim = VectorH1(mesh, dirichlet=BCDisp, order=k_prim)
# FESpace for el. potential with dirichlet-BC at electrodes
fespace_phi_prim = H1(mesh, order=k_prim, dirichlet=BCPotential, definedon="PZT" )
fespace_prim = fespace_u_prim * fespace_phi_prim
q_prim = GridFunction(fespace_prim)
u_prim, phi_prim = q_prim.components
We use the electric enthalpy \(h\) corresponding to the linear material law,
One easily checks that the variational formulation
corresponds to the partial differential equaitons above.
u_, phi_ = fespace_prim.TrialFunction()
E_ = -grad(phi_)
strain_ = Sym(Grad(u_))
def vec(strain): return CoefficientFunction((strain[0], strain[4], strain[8], (strain[7]+strain[5]), (strain[2]+strain[6]), (strain[1]+strain[3])))
a_prim = BilinearForm(fespace_prim, eliminate_internal=True, symmetric=True)
a_prim += Variation ((0.5*InnerProduct (vec(strain_), cf_CE*vec(strain_) ))*dx).Compile()
a_prim += Variation ((- InnerProduct (cf_e* vec(strain_), E_) - 0.5*InnerProduct(E_, cf_epsilonS*E_))*dx(definedon=mesh.Materials("PZT"))).Compile()
phi_prim.Set(voltage/hp*z, definedon=mesh.Materials("PZT"))
solvers.Newton(a_prim, q_prim, printing=False)
print(u_prim(mesh(len_a,0,-he)))
(-0.00015902960631902638, 3.841709070208746e-05, -0.006152749931929724)
strain_prim = Sym(Grad(u_prim))
E_prim = -Grad(phi_prim)
stress_prim = cf_CE * vec(strain_prim) - cf_e.trans*E_prim
scene = Draw(BoundaryFromVolumeCF(stress_prim[0]), mesh, "s_xx", autoscale=False, min=-2, max=1.5, deformation=1e2*u_prim, **ea)
13.2. Mixed elements#
A variant using mixed TDNNS finite elements may be of advantage in presence of thin structures. Above, a high polynomial order was used to reduce shear locking effects. In the TDNNS method, a different approach is chosen, introducing stresses as additional independent unknowns while reducing the continuity requirements on the displacement.
We derive a formulation for piezoelectric materials sporting stresses \(\boldsymbol{\sigma}\) as independent unknown. To this end, a Legendre transform and derive the total enthalpy \(h\),
Above, \(\mathbf{S}^E\) is the flexibility tensor at constant electric field. The third order tensor \(\mathbf d\) again contains piezoelectric constants, and \(\boldsymbol{\epsilon}^T\) specifies the material’s electric permittivity at constant stress. These constants, rather than \(\mathbf{e}\) and \(\boldsymbol{\epsilon}^S\), are accessible via measurements. Material relations are then given as
In the linear model, they can be computed from each other,
SEmat = np.linalg.inv(CEmat)
dmat = np.matmul(emat, SEmat)
epsilonTmat = epsilonSmat + np.matmul(dmat, emat.T)
cf_SE = mesh.MaterialCF({"PZT": CoefficientFunction(tuple(SEmat.reshape(1, -1)[0]), dims=(6,6)), "alu": CoefficientFunction(tuple(Smat_alu.reshape(1, -1)[0]), dims=(6,6))})
cf_d = CoefficientFunction(tuple(dmat.reshape(1, -1)[0]), dims=(3,6))
cf_epsilonT = CoefficientFunction(tuple(epsilonTmat.reshape(1, -1)[0]), dims=(3,3))
n = specialcf.normal(3)
We set up the finite element spaces for the TDNNS method: tangential-continuous \(H(\operatorname{curl})\) elements for the displacement, normal-normal continuous tensor-valued \(H(\operatorname{div}\operatorname{div})\) elements for the stress, and \(H^1\) elements for the electric field. To hybridize the normal-normal continuity of the stress tensor, an additional Lagrangian multiplier on all interfaces is introduced. It resembles the normal displacement.
k = 2
k_phi = k+1
mesh.Curve(k+1)
fespace_stress = Discontinuous(HDivDiv(mesh, order=k))
fespace_disp = HCurl(mesh, dirichlet=BCDisp, order=k)
fespace_pot = H1(mesh, order=k_phi, dirichlet=BCPotential, definedon="PZT" )
fespace_dispn = HDiv(mesh, order=k, dirichlet=BCDisp, orderinner=0 )
fespace = fespace_disp*fespace_dispn*fespace_stress*fespace_pot
q = GridFunction(fespace)
u, un, sigma, phi = q.components
The piezoelectric variational formulation is now built on the total enthalpy, and reads (with \(\vec u\) and \(\boldsymbol{\sigma}\) both as independent unknowns)
When applying the TDNNS method, a distributional interpretation of the work pair \(\int_\Omega \boldsymbol{\varepsilon} : \boldsymbol{\sigma}\, dV\) is necessary, as displacements are discontinuous. For tangentially continuous displacements \(\vec u\) and stresses with \(\vec n \cdot \boldsymbol{\sigma} \cdot \vec n =: \sigma_{nn}\), the work pair is replaced by
An alternative, equivalent representation needs the tangetial components of stress and displacement vectors, \(\vec u_t = \vec u - \vec u \cdot \vec n \, vec n\), \(\vec \sigma_{nt} = (\boldsymbol{\sigma} \cdot \vec n) - \sigma_{nn} \vec n\),
When hybridizing using \(\tilde u_n\) as independent unknown, the nn-continuity of \(\boldsymbol{\sigma}\) can be dropped, using
This expression is used below.
u_, un_, sigma_, phi_ = fespace.TrialFunction()
E_ = -grad(phi_)
def vec(stress): return CoefficientFunction((stress[0], stress[4], stress[8], 0.5*(stress[7]+stress[5]), 0.5*(stress[2]+stress[6]), 0.5*(stress[1]+stress[3])))
def tang(u): return u - InnerProduct(u,n)*n
def normal(u): return InnerProduct(u,n)*n
a = BilinearForm(fespace, eliminate_internal=True, symmetric=True)
a += Variation ((-0.5*InnerProduct (vec(sigma_), cf_SE*vec(sigma_) ))*dx).Compile()
a += Variation ((- InnerProduct (cf_d* vec(sigma_), E_) - 0.5*InnerProduct(E_, cf_epsilonT*E_))*dx(definedon=mesh.Materials("PZT"))).Compile()
a += Variation ( -InnerProduct(div(sigma_),u_)*dx).Compile()
a += Variation ( InnerProduct(sigma_*n, normal(un_)+tang(u_))*dx(element_boundary=True)).Compile()
phi.Set(voltage/hp*z, definedon=mesh.Materials("PZT"))
solvers.Newton(a, q, printing=None)
print(u(mesh(len_a,0,-he)))
(-0.00015470149595052165, 3.7173061322361094e-05, -0.0059728684199435965)
scene = Draw(BoundaryFromVolumeCF(sigma[0]), mesh, "s_xx", autoscale=False, min=-2, max=1.5, deformation=1e2*u, **ea)
13.3. Computing eigenfrequencies#
In frequency domain, with \(\omega = 2 \pi f\) the (angular) frequency of the system, the variational formulation is extended by inertia terms,
For the case of grounded electrodes, \(\varphi = 0\) on the top and bottom surface of the patch, this translates into a generalized eigenvalue problem,
with \(\mathbf M\) the (generalized) mass matrix of the system, and \(\mathbf K\) the generalized stiffness matrix.
Note that for the piezoelectric case, the mass matrix is not invertible, as inertia affects only the displacements (or accelerations), but not electric potential or stresses. The stiffness matrix \(\mathbf K\) is invertible, but not positive definite. Nevertheless, from physics we expect only positive eigenvalues \(\omega^2\) to the generalized system.
We treat the problem using inverse iteration. Here, care must be taken when choosing the random inital vectors to the iteration, as they must be compatible to the constraint equations where \(\mathbf M\) is zero. We generate according vectors by application of \(\mathbf K^{-1} \mathbf M\).
13.3.1. piezoelectric coupling coefficients#
Different eigenfrequencies are observed when both electrodes are grounded (short circuit) and when the circuit is not closed (open circuit), such that no charges are observed at the electrodes \(\int_{\Gamma} \vec D \cdot \vec n\, ds = 0\), and some constant voltage \(V\) is observed. For \(\omega_{sc,i}\) and \(\omega_{oc,i}\) denoting the different eigenfrequencies, the piezoelectric coupling coefficient \(k_i\) can be computed by
Note that \(k_i\) is real as the open circuit frequencies are always as least as high as the short circuit frequencies.
13.3.2. short circuit case#
In case of short circuit, all electrodes are grounded, and \(\varphi = 0\) on all Dirichlet boundaries. The setup from above can directly be used, only the mass matrix has to be computed. For consistency, the whole setup is repeated below.
k = 2
k_phi = k+1
mesh.Curve(k+1)
fespace_stress = Discontinuous(HDivDiv(mesh, order=k))
fespace_disp = HCurl(mesh, dirichlet=BCDisp, order=k)
fespace_pot = H1(mesh, order=k_phi, dirichlet=BCPotential, definedon="PZT" )
fespace_dispn = HDiv(mesh, order=k, dirichlet=BCDisp, orderinner=0 )
fespace = fespace_disp*fespace_dispn*fespace_stress*fespace_pot
u_, un_, sigma_, phi_ = fespace.TrialFunction()
E_ = -grad(phi_)
def vec(stress): return CoefficientFunction((stress[0], stress[4], stress[8], 0.5*(stress[7]+stress[5]), 0.5*(stress[2]+stress[6]), 0.5*(stress[1]+stress[3])))
def tang(u): return u - InnerProduct(u,n)*n
def normal(u): return InnerProduct(u,n)*n
m = BilinearForm(fespace, symmetric=True)
m += Variation(0.5*cf_rho*InnerProduct(u_, u_)*dx)
a = BilinearForm(fespace, eliminate_internal=True, symmetric=True)
a += Variation ((-0.5*InnerProduct (vec(sigma_), cf_SE*vec(sigma_) ) - InnerProduct (cf_d* vec(sigma_), E_) - 0.5*InnerProduct(E_, cf_epsilonT*E_))*dx).Compile()
a += Variation ( -InnerProduct(div(sigma_),u_)*dx).Compile()
a += Variation ( InnerProduct(sigma_*n, normal(un_)+tang(u_))*dx(element_boundary=True)).Compile()
q.vec[:] = 0
a.AssembleLinearization(q.vec)
m.AssembleLinearization(q.vec)
The eigenvalue solver based on the inverse iteration is listed below,
from scipy.linalg import eigh
from numpy.random import rand
def SolveCondensed(a, res, solver, w):
w[:] = 0
if a.condense:
res.data += a.harmonic_extension_trans * res
w.data = solver * res
w.data += a.harmonic_extension * w
w.data += a.inner_solve * res
else:
w.data = solver * res
## Eigenvalue solver, based on inverse iteration
## more or less, copied from https://ngsolve.org/docu/latest/i-tutorials/unit-2.2-eigenvalues/pinvit.html
def Eigenvalues_InverseIteration(a, m, num, shift):
fes = a.space
u = GridFunction(fes, multidim=num)
inva = a.mat.Inverse(fes.FreeDofs(a.condense))
r = u.vec.CreateVector()
Av = u.vec.CreateVector()
Mv = u.vec.CreateVector()
vecs = []
for i in range(2*num):
vecs.append (u.vec.CreateVector())
for v in u.vecs:
r.FV().NumPy()[:] = rand(fes.ndof)
v.data = Projector(fes.FreeDofs(), True) * r
v.data = m.mat * r
# r.data = inva * v
SolveCondensed(a, r, inva, v)
v.data = r
asmall = Matrix(2*num, 2*num)
msmall = Matrix(2*num, 2*num)
lams = num * [shift]
for i in range(20):
for j in range(num):
vecs[j].data = u.vecs[j]
r.data = a.mat * vecs[j] - lams[j] * m.mat * vecs[j]
# vecs[num+j].data = inva * r
SolveCondensed(a, r, inva, vecs[num+j])
for j in range(2*num):
Av.data = a.mat * vecs[j]
Mv.data = m.mat * vecs[j]
for k in range(2*num):
asmall[j,k] = InnerProduct(Av, vecs[k])
msmall[j,k] = InnerProduct(Mv, vecs[k])
ev,evec = eigh(a=asmall, b=msmall)
lams[:] = ev[:num]
print (i, ":", [lam for lam in lams])
for j in range(num):
u.vecs[j][:] = 0.0
for k in range(2*num):
u.vecs[j].data += float(evec[k,j]) * vecs[k]
return (lams, u)
Compute short circuit frequencies using this solver,
(lami_sc, q_sc) = Eigenvalues_InverseIteration(a, m, 5, 1000**2)
print("\n\nEigenfrequencies [rad/s]")
print([sqrt(l) for l in lami_sc])
print("\n\nEigenfrequencies [1/s]")
print([sqrt(l)/2/np.pi for l in lami_sc])
0 : [12306629963.042618, 116454629502.5949, 159262547041.05182, 539866628133.9658, 1364365944393.2778]
1 : [66596654.62835478, 615315235.7677677, 3450406518.0415573, 7485520594.505121, 15319927163.957657]
2 : [66591498.66065746, 611338830.2790228, 3226381257.874194, 5831935748.3917055, 6370042336.275903]
3 : [66590452.23347678, 611045050.5853229, 3223940793.2687483, 5774486393.630982, 6220554213.1378765]
4 : [66589405.77128258, 610752490.9607207, 3223786931.1281166, 5772628718.415713, 6217798856.895737]
5 : [66588377.25670973, 610713198.8133699, 3223646698.9600143, 5772208144.652463, 6217412029.120544]
6 : [66587352.200832225, 610676022.3542396, 3223524068.094307, 5771821056.821475, 6217084651.193287]
7 : [66586341.704706766, 610652772.1725043, 3223402399.3768563, 5771462218.912892, 6216785634.109495]
8 : [66585336.65466453, 610629997.9164193, 3223282428.91716, 5771112620.2227, 6216492893.566117]
9 : [66584340.36187207, 610607104.0226856, 3223162248.9323487, 5770762873.394534, 6216202840.074675]
10 : [66583353.69564058, 610584615.459571, 3223043285.2512074, 5770417418.90717, 6215915700.823372]
11 : [66582369.86877431, 610561724.4155605, 3222924006.9179373, 5770070175.582284, 6215630021.458786]
12 : [66581397.76797111, 610539168.5486153, 3222805533.2355375, 5769728493.265197, 6215346501.421258]
13 : [66580421.697608665, 610516132.6627824, 3222686227.9135265, 5769380451.445997, 6215064659.042621]
14 : [66579462.018988356, 610493418.0615281, 3222568116.101163, 5769036050.120029, 6214784936.866699]
15 : [66578491.578203894, 610470158.8243984, 3222449172.5181365, 5768686277.244619, 6214506158.503249]
16 : [66577538.71163089, 610447175.6009053, 3222330985.7214093, 5768342084.466382, 6214228848.561713]
17 : [66576569.05031492, 610423600.5627495, 3222211542.960826, 5767988093.793344, 6213952877.344145]
18 : [66575622.596352465, 610400313.503044, 3222093279.9953074, 5767637664.610851, 6213678492.184956]
19 : [66574656.40474142, 610376378.5632991, 3221973865.505077, 5767279391.607135, 6213404728.624552]
Eigenfrequencies [rad/s]
[8159.3294090103645, 24705.79645676899, 56762.433576310636, 75942.60590476952, 78825.15289312512]
Eigenfrequencies [1/s]
[1298.5976077590726, 3932.049629117018, 9034.021885595208, 12086.641121023828, 12545.412722915278]
for i in range(5):
sigma = q_sc.components[2].MDComponent(i)
disp = q_sc.components[0].MDComponent(i)
scene = Draw(BoundaryFromVolumeCF(disp[2]), mesh, "uz", deformation=1e-3*BoundaryFromVolumeCF(disp), **ea)
13.3.2.1. open circuit case#
In this case, an unknown voltage \(V\) is introduced on one of the electrodes, the other electrode is still grounded. We add a single independent unknown value \(V\) using a number finite element space.
k = 2
k_phi = k+1
mesh.Curve(k+1)
fespace_stress = Discontinuous(HDivDiv(mesh, order=k))
fespace_disp = HCurl(mesh, dirichlet=BCDisp, order=k)
fespace_pot = H1(mesh, order=k_phi, dirichlet=BCPotential, definedon="PZT" )
fespace_dispn = HDiv(mesh, order=k, dirichlet=BCDisp, orderinner=0 )
fespace_num = FESpace("number", mesh, definedon=mesh.Materials("PZT"))
fespace = fespace_disp*fespace_dispn*fespace_stress*fespace_pot*fespace_num
q = GridFunction(fespace)
The electric potential is composed from two parts, where one part \(\tilde \varphi\) respects zero voltage on both electrodes, and \(V \varphi_{el}\) models the unknown voltage, with \(\varphi_{el}\) a finite element function that is set to \(\varphi_{el} = 1\) on the electrode in question.
phi_el = GridFunction(fespace_pot)
phi_el.Set(1/hp*z, definedon=mesh.Materials("PZT"))
u, un, sigma, tildephi, V = q.components
phi = tildephi + V*phi_el
E = - Grad(tildephi) - V*Grad(phi_el)
Otherwise, the procedure is unchanged,
u_, un_, sigma_, phi_, V_ = fespace.TrialFunction()
E_ = -grad(phi_)-V_*grad(phi_el)
def vec(stress): return CoefficientFunction((stress[0], stress[4], stress[8], 0.5*(stress[7]+stress[5]), 0.5*(stress[2]+stress[6]), 0.5*(stress[1]+stress[3])))
def tang(u): return u - InnerProduct(u,n)*n
def normal(u): return InnerProduct(u,n)*n
a = BilinearForm(fespace, eliminate_internal=True, symmetric=True)
a += Variation ((-0.5*InnerProduct (vec(sigma_), cf_SE*vec(sigma_) ))*dx).Compile()
a += Variation ((- InnerProduct (cf_d* vec(sigma_), E_) - 0.5*InnerProduct(E_, cf_epsilonT*E_))*dx(definedon=mesh.Materials("PZT"))).Compile()
a += Variation ( -InnerProduct(div(sigma_),u_)*dx).Compile()
a += Variation ( InnerProduct(sigma_*n, normal(un_)+tang(u_))*dx(element_boundary=True)).Compile()
m = BilinearForm(fespace, symmetric=True)
m += Variation(0.5*cf_rho*InnerProduct(u_, u_)*dx)
q.vec[:] = 0
a.AssembleLinearization(q.vec)
m.AssembleLinearization(q.vec)
(lami_oc, q_oc) = Eigenvalues_InverseIteration(a, m, 5, 1000**2)
print("\n\nEigenfrequencies [rad/s]")
print([sqrt(l) for l in lami_oc])
print("\n\nEigenfrequencies [1/s]")
print([sqrt(l)/2/np.pi for l in lami_oc])
0 : [12428362266.618235, 187070325932.51984, 260604599148.79648, 463128333096.19226, 2444689445638.15]
1 : [66950131.82576454, 615820660.7363564, 4061742247.714285, 9764621088.90528, 11233216090.056099]
2 : [66937058.53729574, 611761750.1644193, 3365440889.7156444, 6003211135.816289, 6415804215.896781]
3 : [66935129.17255815, 611742033.616379, 3360019277.594566, 5945113128.688956, 6238872853.254135]
4 : [66933199.708473906, 611723768.9223421, 3359791831.5987916, 5943197687.701218, 6231314816.890937]
5 : [66931453.755083404, 611706184.5042431, 3359591751.7405124, 5942741223.446146, 6230579480.245817]
6 : [66929718.862739526, 611688604.6486713, 3359407109.993349, 5942325296.861118, 6230162793.328171]
7 : [66928169.72966496, 611671288.7292622, 3359226525.537151, 5941928628.983474, 6229763901.742494]
8 : [66926632.14211267, 611653975.2010018, 3359058069.976492, 5941540099.190914, 6229375316.541985]
9 : [66925060.6980679, 611636713.491, 3358888129.068796, 5941147320.875177, 6228994739.100663]
10 : [66923507.71981526, 611619453.8120254, 3358724831.145948, 5940758733.226898, 6228622644.757493]
11 : [66921857.146413945, 611602175.8202469, 3358556397.967463, 5940365850.169574, 6228255422.475521]
12 : [66920226.31864226, 611584901.859938, 3358389643.143736, 5939978647.60568, 6227894389.486891]
13 : [66918451.427142575, 611567565.1648464, 3358216018.7741694, 5939580024.441161, 6227537203.931234]
14 : [66916708.00857649, 611550234.8087065, 3358045101.5144796, 5939184222.615569, 6227185918.284966]
15 : [66914786.38280966, 611532805.6020054, 3357864119.1628, 5938777037.885894, 6226835699.579835]
16 : [66912899.84531207, 611515385.6541488, 3357684238.222408, 5938375460.649728, 6226489346.435639]
17 : [66910790.49661254, 611497834.1468385, 3357492374.358253, 5937954172.600821, 6226143484.352255]
18 : [66908738.325439304, 611480295.0867739, 3357304175.748191, 5937536299.900354, 6225802296.0108595]
19 : [66906418.77030979, 611462592.017859, 3357099513.648421, 5937099374.326156, 6225457665.442678]
Eigenfrequencies [rad/s]
[8179.634390992655, 24727.76965312195, 57940.482511353155, 77052.57539061336, 78901.56947388739]
Eigenfrequencies [1/s]
[1301.829246010946, 3935.546771932121, 9221.51419681137, 12263.298251377046, 12557.574799477774]
We compute the coupling coefficients. For the first, third and fourth mode, a high coupling coefficient indicates that the respective mode can be actuated by the patch actor.
print(" open circuit short circuit coupling coefficient")
for i in range(5):
print(f"omega {i+1}: {sqrt(lami_oc[i])/2/np.pi} \t{sqrt(lami_sc[i])/2/np.pi}\t{np.sqrt(np.abs(lami_oc[i]/lami_sc[i]-1))}")
open circuit short circuit coupling coefficient
omega 1: 1301.829246010946 1298.5976077590726 0.07059258617481691
omega 2: 3935.546771932121 3932.049629117018 0.04218506570649472
omega 3: 9221.51419681137 9034.021885595208 0.2047895955807884
omega 4: 12263.298251377046 12086.641121023828 0.17159668777564988
omega 5: 12557.574799477774 12545.412722915278 0.044043479158252144