RasmahRasmah

Pyroelectricity

Overview

A pyroelectric crystal generates an electric field when its temperature changes. The linearized constitutive law (Nye 1985; Newnham 2005) couples mechanics, electricity and temperature:

\[\sigma = C\varepsilon - e^\top E - \beta\,\Delta T, \qquad D = e\,\varepsilon + \epsilon\,E + p\,\Delta T\]

with σ stress, D electric displacement, E = −∇φ, ΔT the temperature change, p the pyroelectric coefficient (along a fixed axis), and β the thermal-stress coefficient. Rasmah assembles the coupled (u, φ, T) block system and solves it in one shot.

using Rasmah

A temperature gradient across a slab

Take a pyroelectric slab, clamp the displacement, impose a linear temperature field T = z, and short the z = 0 and z = 1 electrodes (φ = 0). The temperature gradient alone must induce a depolarising potential.

m = tetrahedralize_box(1.0, 1.0, 1.0, 4, 4, 4)
n = size(m.nodes, 2)

mat = Material(youngs_modulus = 1e3, poisson_ratio = 0.3, relative_permittivity = 2.0,
               pyroelectric_coefficient = 1.0, thermal_conductivity = 1.0)

fixed_u = reduce(vcat, [[3i - 2, 3i - 1, 3i] for i in 1:n])   # clamp all displacement dofs
fixed_T = collect(1:n)
gT = [m.nodes[3, i] for i in 1:n]                        # T = z
fixed_φ = [i for i in 1:n if (m.nodes[3, i] ≈ 0.0 || m.nodes[3, i] ≈ 1.0)]
gφ = zeros(length(fixed_φ))

res = solve(Pyroelectricity(), m, mat;
            fixed_u = fixed_u, gu = zeros(length(fixed_u)),
            fixed_φ = fixed_φ, gφ = gφ,
            fixed_T = fixed_T, gT = gT)
(u = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0  …  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], φ = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0  …  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], T = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0  …  1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0])

Verification — the pyroelectric potential

The 1D solution is φ(z) = −p ΔT z(1−z)/(2ε), so the mid-slab potential is φ(0.5) = −p ΔT L/(8ε) = −1/16:

center = argmin(sum(abs2, m.nodes[:, i] .- [0.5, 0.5, 0.5]) for i in 1:n)
println("φ(centre) = ", round(res.φ[center]; digits = 4), "  (analytic −1/16 = ", round(-1.0 / 16.0; digits = 4), ")")
println("max |u|   = ", round(maximum(abs, res.u); digits = 2), "  (clamped ⇒ 0)")
φ(centre) = -0.0625  (analytic −1/16 = -0.0625)
max |u|   = 0.0  (clamped ⇒ 0)

Next steps