---
jupyter:
  description: 'Chladni figures with planex2d and NGSolve: the modes of a square plate free at all edges, its natural
    frequencies against the tabulated frequency coefficients, and the nodal lines where sand gathers.'
---

# Chladni figures

What do the sand patterns on a bowed plate show? Ernst Chladni strewed sand on metal plates and bowed their edges. At
each natural frequency, the plate vibrates in one of its modes, and the sand gathers on the nodal lines, where the plate
stands still. The frequencies of a square plate free at all its edges are tabulated as the coefficients
$\lambda^2 = \omega a^2 \sqrt{\rho h / D}$, with $a$ the side, $h$ the thickness, $\rho$ the density and
$D = E h^3 / (12 (1 - \nu^2))$ the bending stiffness. For $\nu = 0.3$, the first are 13.463, 19.596, 24.270, 34.788
(twice) and 61.090 (twice), as Apaipong and Sompornjaroensuk list them, in agreement with Leissa's classic tables.

Here, NGSolve finds the modes of a thin square plate, the frequencies are checked against the coefficients, and the
nodal lines are drawn.

![Nine modes of a free square plate, their deflection in color and their nodal lines](../_static/gallery/chladni-figures.png){.result}

```python
import math

import ngsolve as ngs
import numpy as np
import planex as px
import scipy.sparse as sp
from scipy.sparse.linalg import eigsh
```

```python nbsphinx="hidden"
from docs_figures import MAP, Row
```

## The plate

The plate is a square of side 1, 0.002 thick, with $E = 1$, $\rho = 1$ and $\nu = 0.3$, so the coefficients come out
directly. It bends as a Reissner-Mindlin plate: the deflection $w$ and the rotations $\theta$ are unknowns of their
own, coupled by the shear $\nabla w - \theta$. For a plate this thin, the model approaches the classic thin plate.
Nothing holds the edges, so they are free:

```python
a, h, E, nu, rho = 1.0, 0.002, 1.0, 0.3, 1.0
D, G = E * h**3 / (12 * (1 - nu**2)), E / (2 * (1 + nu))
plate = px.QuickSketch().RectangleC(a, a).Face()
mesh = ngs.Mesh(px.OCCGeometry(plate).GenerateMesh(maxh=a / 12))

X = ngs.H1(mesh, order=6) * ngs.VectorH1(mesh, order=5)
(w, th), (v, et) = X.TnT()
bend = ngs.Sym(ngs.Grad(th)), ngs.Sym(ngs.Grad(et))
stiff = D * ((1 - nu) * ngs.InnerProduct(*bend) + nu * ngs.div(th) * ngs.div(et)) \
    + 5 / 6 * G * h * (ngs.grad(w) - th) * (ngs.grad(v) - et)
mass = rho * h * w * v + rho * h**3 / 12 * th * et
stiffness, inertia = ngs.BilinearForm(stiff * ngs.dx).Assemble(), ngs.BilinearForm(mass * ngs.dx).Assemble()
K = sp.csr_matrix(stiffness.mat.CSR(), shape=(X.ndof, X.ndof))  # the forms stay alive: CSR() shares their memory
M = sp.csr_matrix(inertia.mat.CSR(), shape=(X.ndof, X.ndof))
```

## The frequencies

SciPy's `eigsh` finds the lowest modes. The first three have no frequency: the plate moves up and down or tilts as a
whole, without bending. The next ones are the Chladni modes, and their coefficients agree with the table to about four
digits:

```python
values, vectors = eigsh(K, k=12, M=M, sigma=-1e-6, which="LM")
order = np.argsort(values)
values, vectors = values[order], vectors[:, order]
table = [13.463, 19.596, 24.270, 34.788, 34.788, 61.090, 61.090, 63.655, 69.241]
for value, listed in zip(values[3:], table):
    print("λ²", round(math.sqrt(value) * a**2 * math.sqrt(rho * h / D), 3), "| table", listed)
```

## The nodal lines

On a grid over the plate, the deflection $w$ of a mode changes its sign across the nodal lines, which are drawn in
dark, the deflection in color (Figure 1). Where two modes share a frequency, any combination of them vibrates at it as
well, and the solver returns one of them:

```python
mode = ngs.GridFunction(X)
grid = np.linspace(-a / 2, a / 2, 121) * (1 - 1e-9)
gh, gv = np.meshgrid(grid, grid)


def deflection(k):
    mode.vec.FV().NumPy()[:] = vectors[:, k]
    shape = np.asarray(mode.components[0](mesh(gh.ravel(), gv.ravel()))).reshape(gh.shape)
    return shape / np.abs(shape).max()
```

```python nbsphinx="hide-input"
Row(*[(plate, f"λ² = {math.sqrt(values[k]) * a**2 * math.sqrt(rho * h / D):.2f}",
       {"image": {"values": deflection(k), "extent": (-a / 2, a / 2, -a / 2, a / 2), "cmap": MAP.reversed(),
                  "clim": (-1, 1), "levels": [0]}}) for k in range(3, 9)],
    caption="Six modes of a free square plate and their nodal lines")
```

```python nbsphinx="hidden"
import sys
sys.path.insert(0, "..")
import matplotlib.pyplot as plt
from docs_figures import INK
from gallery import save_figure
fig, axes = plt.subplots(3, 3, figsize=(6, 5.4))
for ax, k in zip(axes.ravel(), range(3, 12)):
    shown = ax.imshow(deflection(k), origin="lower", extent=(-a / 2, a / 2, -a / 2, a / 2), cmap=MAP.reversed(),
                      vmin=-1, vmax=1)
    ax.contour(gh, gv, deflection(k), levels=[0], colors=INK, linewidths=1.6)
    px.Plot(plate, ax=ax, face_col=None, line_col=INK)
fig.colorbar(shown, ax=axes, label="deflection", fraction=0.04, pad=0.03)
save_figure(fig, "chladni-figures", bare=True, bars={9: "deflection"})
```