RasmahRasmah

A Piezoelectric Stack Actuator (u + φ)

using Rasmah
using LinearAlgebra

Overview

Piezoelectric materials turn voltage into displacement and displacement into voltage. They are the active element in every ultrasound transducer, the “buzzer” in a smoke alarm, and the spark that lights a barbecue igniter. The coupling is directional: a field applied along the material’s poling axis stretches it by ε = d₃₃ E, and a mechanical strain pushes charge back.

This tutorial models a PZT-5A stack actuator — a 10 × 10 × 10 mm ceramic block, silvered on its flat faces to form the electrodes, driven by 100 V across its thickness. Piezoelectricity() couples the displacement :u and the electric potential :φ through the piezoelectric strain matrix d and the permittivity, both read from the material. The block is the building unit of real piezo stacks (a column of such blocks, poled along their height, is the actuator of a fuel injector or a nano-positioning stage). We draw it as a single box, clamp its backing face, and check the thickness-mode displacement against the textbook Δt = d₃₃ V.

The workflow is the canonical one, with face selection on the exact BRep before meshing:

shape → boundary sets → mesh → boundary conditions → solve

The physics

For a ceramic poled along z (the 6mm symmetry class of PZT), the piezoelectric strain matrix has three independent coefficients; the one we use here is the longitudinal charge constant d₃₃, the strain per volt per metre along the poling axis:

\[ \varepsilon_{zz} = d_{33}\,\frac{V}{t} \qquad\Longrightarrow\qquad \Delta t = d_{33}\, V . \]

The coupled system is the symmetric saddle

\[ \begin{bmatrix} K_{uu} & K_{u\varphi} \\ K_{u\varphi}^{\mathsf T} & -K_{\varphi\varphi} \end{bmatrix} \begin{bmatrix} u \\ \varphi \end{bmatrix} = \begin{bmatrix} f \\ q \end{bmatrix}, \]

solved with LU. A second constant, d₃₁, describes the lateral contraction that accompanies the thickness stretch — free here, so the block gets slightly narrower as it stretches.

Drawing the block

A stack element is a box with silvered top and bottom faces:

a, t = 0.010, 0.010                       # 10 × 10 × 10 mm (m)

block = box(a, a, t)
b = brep(block)
BRep{Vertex{Float64}, Edge{LineCurve{Float64}}, Loop, Face{PlaneSurface{Float64}}, Shell, Solid}(Vertex{Float64}[Vertex{Float64}([-0.005, -0.005, -0.005]), Vertex{Float64}([0.005, -0.005, -0.005]), Vertex{Float64}([0.005, 0.005, -0.005]), Vertex{Float64}([-0.005, 0.005, -0.005]), Vertex{Float64}([-0.005, -0.005, 0.005]), Vertex{Float64}([0.005, -0.005, 0.005]), Vertex{Float64}([0.005, 0.005, 0.005]), Vertex{Float64}([-0.005, 0.005, 0.005])], Edge{LineCurve{Float64}}[Edge{LineCurve{Float64}}(1, 2, LineCurve{Float64}([-0.005, -0.005, -0.005], [0.005, -0.005, -0.005])), Edge{LineCurve{Float64}}(2, 3, LineCurve{Float64}([0.005, -0.005, -0.005], [0.005, 0.005, -0.005])), Edge{LineCurve{Float64}}(3, 4, LineCurve{Float64}([0.005, 0.005, -0.005], [-0.005, 0.005, -0.005])), Edge{LineCurve{Float64}}(4, 1, LineCurve{Float64}([-0.005, 0.005, -0.005], [-0.005, -0.005, -0.005])), Edge{LineCurve{Float64}}(5, 6, LineCurve{Float64}([-0.005, -0.005, 0.005], [0.005, -0.005, 0.005])), Edge{LineCurve{Float64}}(6, 7, LineCurve{Float64}([0.005, -0.005, 0.005], [0.005, 0.005, 0.005])), Edge{LineCurve{Float64}}(7, 8, LineCurve{Float64}([0.005, 0.005, 0.005], [-0.005, 0.005, 0.005])), Edge{LineCurve{Float64}}(8, 5, LineCurve{Float64}([-0.005, 0.005, 0.005], [-0.005, -0.005, 0.005])), Edge{LineCurve{Float64}}(1, 5, LineCurve{Float64}([-0.005, -0.005, -0.005], [-0.005, -0.005, 0.005])), Edge{LineCurve{Float64}}(2, 6, LineCurve{Float64}([0.005, -0.005, -0.005], [0.005, -0.005, 0.005])), Edge{LineCurve{Float64}}(3, 7, LineCurve{Float64}([0.005, 0.005, -0.005], [0.005, 0.005, 0.005])), Edge{LineCurve{Float64}}(4, 8, LineCurve{Float64}([-0.005, 0.005, -0.005], [-0.005, 0.005, 0.005]))], Loop[Loop([-1, -4, -3, -2]), Loop([5, 6, 7, 8]), Loop([1, 10, -5, -9]), Loop([3, 12, -7, -11]), Loop([9, -8, -12, 4]), Loop([2, 11, -6, -10])], Face{PlaneSurface{Float64}}[Face{PlaneSurface{Float64}}(PlaneSurface{Float64}([-0.005, -0.005, -0.005], [0.01, 0.0, 0.0], [0.0, 0.01, 0.0]), [1]), Face{PlaneSurface{Float64}}(PlaneSurface{Float64}([-0.005, -0.005, 0.005], [0.01, 0.0, 0.0], [0.0, 0.01, 0.0]), [2]), Face{PlaneSurface{Float64}}(PlaneSurface{Float64}([-0.005, -0.005, -0.005], [0.01, 0.0, 0.0], [0.0, 0.0, 0.01]), [3]), Face{PlaneSurface{Float64}}(PlaneSurface{Float64}([-0.005, 0.005, -0.005], [0.01, 0.0, 0.0], [0.0, 0.0, 0.01]), [4]), Face{PlaneSurface{Float64}}(PlaneSurface{Float64}([-0.005, -0.005, -0.005], [0.0, 0.01, 0.0], [0.0, 0.0, 0.01]), [5]), Face{PlaneSurface{Float64}}(PlaneSurface{Float64}([0.005, -0.005, -0.005], [0.0, 0.01, 0.0], [0.0, 0.0, 0.01]), [6])], Shell[Shell([1, 2, 3, 4, 5, 6])], Solid[Solid([1])])
vtk(b)

Boundary sets and meshing

The backing face (bottom) is bonded to the load frame, so it is clamped and grounded; the radiating face (top) carries the drive voltage:

mesh = tetrahedralize(b, 8)

bottom = faces_in_box(b, [-0.006, -0.006, -0.0051], [0.006, 0.006, -0.0049])
top    = faces_in_box(b, [-0.006, -0.006,  0.0049], [0.006, 0.006,  0.0051])

println("nodes: ", size(mesh.nodes, 2), ", elements: ", size(mesh.elements, 2))
nodes: 729, elements: 3072

Driving the block

Ground the backing face, drive the radiating face with 100 V, and clamp the backing mechanically:

V = 100.0
res = solve(Piezoelectricity(), mesh, pzt;
            fixed_u = bottom => [0.0, 0.0, 0.0],
            fixed_φ = [bottom => 0.0, top => V])

println("potential range: ", round(minimum(ustrip.(res.φ)); digits = 1), " … ",
        round(maximum(ustrip.(res.φ)); digits = 1), " V")
potential range: 0.0 … 100.0 V

The potential ramps linearly from 0 to 100 V through the thickness (the electrostatic problem is a 1-D Laplace field), so the electric field is a uniform E = V/t. The mechanical response is a uniform thickness stretch:

d33 = piezoelectric_d33(pzt)

println("FEM thickness change: ", round(maximum(ustrip.(res.u[3:3:end])); digits = 10), " m")
println("hand, Δt = d₃₃ V:    ", round(d33 * V; digits = 10), " m")
println("relative error:      ", round((maximum(ustrip.(res.u[3:3:end])) - d33 * V) / (d33 * V) * 100; digits = 2), " %")
FEM thickness change: 3.67e-8 m
hand, Δt = d₃₃ V:    3.74e-8 m
relative error:      -2.0 %

XX classic thickness-mode actuation, 37.4 nm at 100 V. (Driven at its thickness resonance — ≈ 200 kHz for a 10 mm block — the same stack radiates acoustic power, which is exactly how a piezo “buzzer” works.)

vtk(mesh; field = ustrip.(res.φ), fieldname = "electric potential (V)")
umag = [norm(ustrip.(res.u)[(3i-2):(3i)]) for i in 1:size(mesh.nodes, 2)]
vtk(mesh; field = umag, fieldname = "|u| (thickness mode)")

Notes

  • The d (3×6) matrix and permittivity are taken from the material via piezoelectric_strain_matrix and relative_permittivity; piezoelectric_d33 reads the single coefficient used above.
  • The default solver is LU (the coupled system is an indefinite saddle).
  • The converse effect — strain producing charge — is the same matrix, so the identical model covers a force-sensing element: load the face, read the potential.

Design (topology optimization)

PiezoelectricTopOpt() minimizes the coupled electro-mechanical compliance with the density method — the density scales the full piezoelectric block (elasticity + coupling + permittivity):

φ, hist = solve(PiezoelectricTopOpt(), mesh, pzt;
                fixed_u = bottom => [0.0, 0.0, 0.0], loads_u = top => [1e6, 0.0, 0.0],
                fixed_φ = [bottom => 0.0, top => V],
                volfrac = 0.5, max_iter = 40)

println("final electro-mechanical compliance = ", round(hist[end]; digits = 4))
println("final volume fraction               = ", round(sum(φ) / length(φ); digits = 4))
final electro-mechanical compliance = 0.0008
final volume fraction               = 0.5
vtk(mesh; field = φ, fieldname = "piezoelectric density")

Next steps