Matching layers, backing, lenses and finite-element design
The piezoelectric plate of Module 16 is a 34 MRayl resonator trying to push sound into 1.5 MRayl tissue. Everything stacked on it, the matching layers in front, the backing behind and the lens on top, exists to fix that mismatch without losing the pulse. This module builds each piece with one transfer-matrix model you can run in the browser and in MATLAB, then shows where one-dimensional models stop and finite-element tools take over.
You will be able to
- Design one, two and three matching layers and say why \(\sqrt{Z_p Z_L }\), Desilets and binomial impedances differ
- Chain ABCD matrices for plate, backing, layers and load, and read bandwidth, pulse length and insertion loss off the result
- Quote the backing trade-off, \(Q \approx (\pi /2)Z_p /(Z_B + Z_L)\), with numbers
- Size a silicone elevation lens and predict its focus, loss, apodization and spectral downshift
- Explain lateral modes and the w/t rule, and when a PZFlex-class finite-element model is required
Assumes
- Module 16 (Piezoelectrics & KLM): thickness-mode plate, \(k_t\), \(C_0\), the KLM circuit (16.4–16.5)
- Module 15 (Wave physics & harmonics): plane-wave reflection and transmission at an interface (15.2)
- Module 04 (Sampling, Fourier, LTI): the beam as a Fourier transform of the aperture (04.8), for the elevation beam
- Module 27 (Spectral estimation & QUS): attenuation and spectral downshift (27.6)
17.1The mismatch: why a bare plate rings
Acoustic impedance \(Z = \rho c\) (in MRayl = \(10^{6}\) kg m−2 s−1) sets how a plane wave divides at an interface. For a wave going from medium 1 into medium 2 at normal incidence, with pressure amplitudes:
- \(Z_1,\ Z_2\)
- acoustic impedances of the incident and transmitting media, MRayl
- \(R\)
- pressure reflection coefficient
- \(T_p\)
- pressure transmission coefficient
- \(T_I\)
- intensity (power) transmission coefficient
PZT (≈ 30–36 MRayl; this module uses 34.2 MRayl, \(\rho = 7500\) kg/m³, \(v = 4560\) m/s) into tissue (1.5–1.7 MRayl; here 1.5) gives \(R = -0.916\) and \(T_I = 0.161\): −7.9 dB of power one way, and 84 % of each wave reflected back into the plate. The reflected energy bounces between the front and back faces, and every round trip (one period at the half-wave resonance) leaks only a little out. That is ring-down. For a half-wave plate between a backing \(Z_B\) and a load \(Z_L\), the amplitude ratio per round trip is \(r = R_{\mathrm{front}} R_{\mathrm{back}}\), so
- \(r\)
- amplitude ratio per round trip, \(R_{\mathrm{front}} R_{\mathrm{back}}\)
- \(Z_p\)
- plate impedance, MRayl
- \(Z_B,\ Z_L\)
- backing and load impedance, MRayl
- \(Q\)
- quality factor of the loaded plate, dimensionless
- \(\ln\)
- natural logarithm
Bounce diagram: what leaves the plate each round trip
Exact formulaPlane-wave reflection coefficients at normal incidence; ideal \(\lambda /4\) layer at \(f_0\). Bars are drawn in log Z for intuition.Think of the plate as a cavity with two leaky mirrors. Matching layers make the front mirror transparent, backing makes the back mirror absorbing. Both shorten the pulse; only the first sends the energy where you want it.
17.2The quarter-wave layer is an impedance transformer
A layer of impedance \(Z_1\) and thickness d in front of a load \(Z_L\) presents the input impedance of a transmission line \((\theta = kd = 2\pi fd/c_1)\):
- \(Z_1\)
- layer impedance, MRayl
- \(Z_L\)
- load impedance, MRayl
- \(\theta\)
- phase thickness \(kd = 2\pi f d/c_1\), rad
- \(d,\ \lambda_1\)
- layer thickness and wavelength in the layer, m
- \(Z_{\mathrm{in}}\)
- impedance seen looking into the layer, MRayl
- \(j\)
- imaginary unit
- \(\tan\)
- tangent (argument in rad)
The wave picture says the same thing. The reflection back into the plate is the sum of the reflection at the plate/layer face and the reflection from the layer/load face, delayed by the round trip \(2\theta\):
- \(\Gamma\)
- pressure reflection seen from the plate
- \(r_1,\ r_2\)
- reflection coefficients of the plate–layer and layer–load interfaces
- \(Z_1,\ Z_L\)
- as in (17.3)
- \(Z_p\)
- plate impedance, MRayl
- \(\theta\)
- one-way phase thickness of the layer, as in (17.3), rad
- \(j\)
- imaginary unit
The band over which transmission stays above one half \((\lvert \Gamma \rvert \le 1/\sqrt{2})\) has the closed form of the quarter-wave transformer (Pozar, Microwave Engineering):
- \(\Delta f\)
- band over which \(\lvert\Gamma\rvert \le \Gamma_m\), Hz
- \(f_0\)
- design (quarter-wave) frequency, Hz
- \(\Gamma_m\)
- largest tolerated reflection magnitude
- \(Z_p,\ Z_L\)
- plate and load impedance, MRayl
- \(\arccos\)
- inverse cosine, result in rad
One layer: transmission band and the two echoes that cancel
Exact formulaClosed-form \(\Gamma\) of one lossless layer between a semi-infinite plate and load, plane wave at normal incidence; no resonator, backing or electrical port."Quarter wave" means a quarter of the wavelength in the layer, \(\lambda_1 = c_1 /f_0\), not in the plate or in tissue. A 2500 m/s epoxy layer at 5 MHz is 125 µm thick; at 15 MHz it is 41.7 µm, and a ±5 µm lapping tolerance becomes ±12 % of the thickness.
17.3Broadband matching and the ABCD stack model
Section 17.2 matched a semi-infinite plate at one frequency. A real transducer is a resonant plate with a backing and an electrical port, and the goal is a short pulse, not perfect transmission at \(f_0\). Three design rules are quoted (impedances in any consistent unit, exponents on \(Z_p\) and \(Z_L\)):
| Rule | One layer | Two layers (\(Z_1\) next to the plate) | What it optimizes |
|---|---|---|---|
| Quarter-wave (geometric mean) | \(Z_p^{1/2} Z_L^{1/2} = 7.16\) | binomial: \(Z_p^{3/4} Z_L^{1/4} = 15.65\), \(Z_p^{1/4} Z_L^{3/4} = 3.28\) | Perfect transmission at \(f_0\); binomial = maximally flat transmission for a semi-infinite source (small-reflection theory) |
| Desilets, Fraser & Kino (1978) | \(Z_p^{1/3} Z_L^{2/3} = 4.25\) | \(Z_p^{4/7} Z_L^{3/7} = 8.95\), \(Z_p^{1/7} Z_L^{6/7} = 2.34\) | Derived from the KLM model of the resonant plate for a broadband, short impulse response |
| Binomial, three layers | \(Z_p^{7/8} Z_L^{1/8}\), \(Z_p^{1/2} Z_L^{1/2}\), \(Z_p^{1/8} Z_L^{7/8} = 23.1\), 7.16, 2.22 | Steps in ln Z weighted 1 : 3 : 3 : 1 | |
Desilets values sit closer to the load than the quarter-wave values because the plate itself is a resonator: it already provides a band-pass response, and a lower-impedance layer damps it more evenly across the band. Souquet and colleagues (1979) gave a related single-layer criterion, \((2Z_p Z_L^2)^{1/3}\), from an equal-quality-factor condition between the electrical and mechanical branches of the equivalent circuit.
The stack as a chain of 2 × 2 matrices. Every uniform layer relates force and particle velocity on its two faces by
- \(F,\ u\)
- force (N) and particle velocity (m/s) at each face
- \(\gamma\)
- propagation constant, 1/m
- \(\alpha\)
- attenuation, Np/m
- \(d\)
- layer thickness, m
- \(Z\)
- mechanical impedance of the layer, kg/s
- \(A\)
- element area, m²
- \(F_{\mathrm{in}},\ u_{\mathrm{in}}\)
- force and velocity at the input face
- \(F_{\mathrm{out}},\ u_{\mathrm{out}}\)
- force and velocity at the output face
- \(\rho,\ c\)
- density (kg/m³) and sound speed (m/s) of the layer
- \(\omega,\ j\)
- angular frequency (rad/s) and imaginary unit
- \(\cosh,\ \sinh\)
- hyperbolic functions of the complex argument \(\gamma d\)
The KLM model of Module 16 puts the plate's electrical port at its centre. With \(C_0 = \varepsilon^S A/t\), \(h^2 = k_t^2 c^D /\varepsilon^S\), \(Z_0 = \rho vA\) and \(\theta_p = \omega t/v\):
- \(h\)
- piezoelectric constant \(h_{33}\), V/m
- \(\theta_p\)
- plate phase thickness \(\omega t/v\), rad
- \(Z_0\)
- plate mechanical impedance \(\rho v A\), kg/s
- \(\varphi\)
- KLM transformer ratio (this chain's convention)
- \(\mathbf{L}_p,\ \mathbf{L}_i\)
- line matrices of the front half-plate and of layer \(i\)
- \(Z_{\mathrm{back}}\)
- input impedance of the rear half-plate on the backing, kg/s
- \(V,\ I\)
- electrical-port voltage (V) and current (A)
- \(\omega,\ j\)
- angular frequency (rad/s) and imaginary unit
- \(C_0\)
- clamped capacitance, F
- \(t\)
- plate thickness, m
- \(\mathbf{L}_1,\dots,\mathbf{L}_N\)
- line matrices of the \(N\) matching layers
- \(F_L,\ u_L\)
- force (N) and velocity (m/s) at the load
- \(X_1\)
- series reactance, Ω
- \(\sin\)
- sine (argument in rad)
With source resistance \(R_s\) and load \(F_L = Z_L A u_L\), the chain gives \(u_L /V_s\), and from it
- \(\eta\)
- one-way transducer power gain
- \(R_s\)
- source resistance, Ω
- \(Z_L A\)
- load mechanical impedance, kg/s
- \(u_L,\ F_L\)
- velocity (m/s) and force (N) at the load
- \(V_s\)
- source voltage, V
- \(H_{\mathrm{pe}}\)
- pulse-echo voltage ratio
- \(f\)
- frequency, Hz
- \(A\)
- element area, m² (as in (17.6))
- \(\lvert\cdot\rvert\)
- magnitude
Reference transducer for all numbers below: PZT-5H-like plate (\(k_t = 0.50\), \(\varepsilon^S /\varepsilon_0 = 1470\)), 0.456 mm thick (\(\lambda /2\) at 5 MHz), area 20 mm² (a 5 mm disk, \(C_0 = 571\) pF), 50 Ω source and receiver, backing 5 MRayl, load 1.5 MRayl, lossless layers \(\lambda /4\) thick in a 2500 m/s material. The pulse length is the −20 dB envelope duration of the pulse-echo impulse response after a receive low-pass (flat to 9 MHz, cos² taper to 12 MHz) that removes the third-harmonic plate resonance.
Matching-layer designer (KLM plate + ABCD layers)
Approx. model1D KLM plate + ABCD layers: thickness mode, plane waves at normal incidence, semi-infinite backing and load, 50 Ω, layer loss \(\alpha \propto f\), no bond lines or lateral modes.| Design (backing 5 MRayl) | Z (MRayl) | \(f_c\) (MHz) | −6 dB FBW | −20 dB pulse (µs) | Two-way loss (dB) |
|---|---|---|---|---|---|
| No layer | – | 4.65 | 27 % | 1.25 | −12.8 |
| One layer, \(\sqrt{Z_p Z_L }\) | 7.16 | 4.77 | 56 % | 1.03 | −5.6 |
| One layer, Desilets | 4.25 | 4.71 | 44 % | 0.78 | −5.4 |
| Two layers, Desilets | 8.95 / 2.34 | 4.67 | 66 % | 0.81 | −5.8 |
| Two layers, binomial | 15.65 / 3.28 | 4.76 | 79 % | 1.00 | −6.2 |
| Three layers, binomial | 23.1 / 7.16 / 2.22 | 4.39 | 84 % | 1.06 | −7.3 |
| Two layers, grid search for the shortest pulse | 8.0 / 2.4 | 4.60 | 63 % | 0.73 | −5.8 |
Two lessons sit in this table. First, the −6 dB fractional bandwidth and the pulse length disagree: the single layer with the widest −6 dB band (63 % at \(Z_1 = 10\) MRayl in the grid of m17_matching_layers.m) has a passband dip that just stays above −6 dB, while the shortest single-layer pulse occurs at 4.25–4.5 MRayl, the Desilets value. Second, the grid optimum for two layers (8.0 / 2.4) lies within about 11 % of Desilets (8.95 / 2.34). The closed forms are good starting points; the final numbers come from a model of your plate, your electrical port and your materials.
Bandwidth here is what Module 10 called the axial-resolution budget and what Module 27 needs for spectral QUS: a −6 dB FBW of 66 % at 4.7 MHz spans about 3.1–6.2 MHz. Harmonic imaging wants the receive band to reach \(2f_0\), which is why harmonic-capable probes push past 80 % FBW with two or three layers.
17.4Material realities: what you can actually build
The design rules return a number; the materials catalogue returns a range. Typical values (state ranges, they depend on formulation and frequency):
| Material | Z (MRayl) | Role | Notes |
|---|---|---|---|
| PZT ceramics | ≈ 30–36 | active plate | Module 16; 1-3 composites are ≈ 10–25 MRayl depending on ceramic volume fraction, which eases matching |
| Unfilled epoxy | ≈ 2.8–3.2 | outer layer, bonding | \(c \approx 2500\text{–}2700\) m/s, \(\rho \approx 1.1\text{–}1.2\) g/cm³ |
| Particle-filled epoxy (alumina, tungsten and others) | ≈ 3 to above 10 | inner layers, backing | Z set by filler fraction; attenuation and scattering rise with particle size and frequency. Alumina/epoxy nanocomposites: 2.8–5.1 MRayl, about 15 dB/mm at 40 MHz reported |
| Parylene C | ≈ 2.6 | outer layer at high frequency | vapour-deposited, thickness controlled to micrometres |
| RTV silicone (filled) | ≈ 1.0–1.6 | lens | \(c \approx 930\text{–}1050\) m/s; ≈ 1–4.5 dB/mm at 5 MHz depending on fillers (one measured RTV: 1.44 dB/mm at 6 MHz, 1.11 MRayl; Section 17.15), rising faster than linearly with frequency |
| Tissue / water | 1.5–1.7 / 1.48 | load | \(c \approx 1450\) (fat) to 1600 m/s |
Four effects decide whether a design survives production. Thickness sets frequency: a layer 10 % too thick moves its transmission peak to 0.91 \(f_0\) but costs only 0.5 dB at \(f_0\). Attenuation in the layer is enhanced by the standing wave inside it: at 0.5 dB/mm/MHz a \(\lambda /4\) layer absorbs 0.73 dB one way, 2.3 times its single-pass loss, and because \(\alpha \propto f\) while \(d \propto 1/f\) the loss in dB is the same at every design frequency. Bondlines: a thin glue film between two stiff materials acts as a compliance; 5 µm of 3 MRayl epoxy costs 0.5 dB at 5 MHz and 4.7 dB at 20 MHz. Availability: the 15–23 MRayl inner layers of binomial designs are hard to make as castable epoxies, which is one practical reason the lower Desilets values are popular.
17.5Backing: damping against sensitivity
The backing absorbs the energy the plate radiates rearward. A backing equal to \(Z_p\) removes the back mirror entirely (one-pass emission, very broadband, but at least half of the energy goes into the backing each way: half when the layers present \(Z_p\) to the plate, about 96 % for a bare face into tissue); air backing sends everything forward but rings. The KLM model quantifies the trade:
- \(Z_{L,\mathrm{eff}}\)
- load seen by the plate through the matching layers, MRayl
- \(Q_m\)
- mechanical quality factor of the loaded plate
- \(Z_p,\ Z_B\)
- plate and backing impedance, MRayl (as in (17.2))
For two Desilets layers, raising the backing from air to 34 MRayl widens the −6 dB band from 63 % to 83 % and raises the two-way loss from 2.0 dB to 17.5 dB. Without matching layers, the same sweep runs from a 16 % band at −3.6 dB (on air backing the 50 Ω electrical port, not the backing, does the damping) to 127 % at −31.2 dB: backing alone buys bandwidth only by spending sensitivity. Production arrays typically sit at moderate backing impedances (a few MRayl) and get most of their bandwidth from the matching layers.
The backing must also absorb what it takes. Energy that reaches the back of the backing and reflects returns to the plate as a late echo at \(2d_b /c_b\). The round-trip attenuation \(2\alpha_b d_b\) plus the transmission losses must push that echo below the displayed dynamic range. Heavily filled epoxies (for example tungsten in epoxy) provide both a tunable Z and high scattering loss; some designs also shape or angle the rear face so the reflection does not return coherently.
Backing trade-off: bandwidth, loss, pulse length and the back echo
Approx. model1D KLM plate + ABCD layers: thickness mode, plane waves at normal incidence, semi-infinite backing and load, 50 Ω, layer loss \(\alpha \propto f\), no bond lines or lateral modes. Back echo from \(2\alpha_b d_b\) with \(c_b = 2000\) m/s.17.6The acoustic lens: elevation focus, loss and downshift
A 1D array focuses electronically in azimuth but needs a fixed focus in elevation. The usual solution is a cylindrical lens of RTV silicone, slower than tissue (\(c_l \approx 950\text{–}1050\) m/s against \(c_t \approx 1540\) m/s). A slower material delays the wave, so the lens must be thickest at the centre: convex. Paraxially, a lens of thickness \(T(y) = T_c - y^2 /(2R)\) adds a delay \(T(y)(1/c_l - 1/c_t)\); focusing at F requires a delay that falls as \(y^2 /(2Fc_t)\) away from the centre. Equating the \(y^2\) terms:
- \(F\)
- elevation focal length from the lens apex, mm
- \(R\)
- radius of curvature of the lens, mm
- \(c_l,\ c_t\)
- sound speed in the lens and in tissue, m/s
Example: \(F = 40\) mm needs \(R = 21.6\) mm. Over a 6 mm elevation aperture the sag is 0.209 mm; with a 0.3 mm edge the centre is 0.509 mm thick. The exact Snell ray trace focuses the marginal ray 0.92 mm short of the paraxial focus (spherical aberration, 2.3 % of F), small against the elevation depth of field at F# 6.7.
Lens attenuation (≈ 2–4.5 dB/mm at 5 MHz for filled silicones) has three effects:
- Loss: two-way through the centre at 5 MHz with 3 dB/mm, 2 × 3 × 0.509 = 3.1 dB.
- Inverse apodization: the centre is thickest, so it is attenuated most (two-way 3.1 dB at the centre, 1.8 dB at the edge). The aperture is weighted towards its edges, which narrows the main lobe slightly and raises sidelobes (−26.5 dB → −25.7 dB two-way in this example).
- Downshift: attenuation that rises with frequency multiplies the pulse-echo spectrum by a low-pass factor. For a Gaussian spectrum of standard deviation \(\sigma\) and a linear law (\(\beta\) Np/MHz two-way), the centre moves by \(-\beta \sigma^2\). With \(\alpha \propto f^{1.5}\) and the aperture-mean thickness (0.439 mm), a 70 % band at 5 MHz shifts its centroid to 4.81 MHz.
The lens impedance should match tissue (filled silicones reach 1.4–1.6 MRayl; unfilled ones are nearer 1.0–1.1 MRayl, \(R \approx -0.2\), a −14 dB reverberation source between lens face and matching stack).
| Typical probe | \(\lambda\) (mm) | H (mm) | \(F_{\mathrm{el}}\) (mm) | R (mm) | Slice \(0.886\lambda F/H\) (mm) | \(H^2 /(4\lambda)\) (mm) |
|---|---|---|---|---|---|---|
| 3.5 MHz abdominal | 0.44 | 12 | 60 | 32.4 | 1.95 | 81.8 |
| 7 MHz linear | 0.22 | 4.5 | 20 | 10.8 | 0.87 | 23.0 |
Lens calculator: focus, rays, apodization and downshift
Approx. modelExact Snell ray trace; elevation beam as a CW Fraunhofer sum at the focus; lens loss \(\alpha = \alpha_5 (f/5)^{n}\) applied to a Gaussian spectrum.A measured spectrum through the lens is already downshifted. When you compare a probe's measured centre frequency with a KLM prediction of the bare stack, include the lens (as one more lossy ABCD layer) or you will "find" a centre-frequency error (0.19 MHz in the example above) that is pure lens loss. The same applies to QUS calibration: the reference phantom spectrum carries the lens filter, which cancels only if the probe and settings are unchanged.
17.7Kerf, sub-dicing and lateral modes
An array element is not a plate. It is a bar of width w (azimuth), thickness t and length equal to the elevation height, cut free from its neighbours by kerfs. When w is comparable with t, the thickness motion couples through Poisson's ratio into width motion, and the bar has resonances that a 1D model cannot know about.
- w/t ≪ 1: the bar is a slender column. Its thickness resonance drops from the plate value v/(2t) to the bar value; in plane strain \(f\,t/c_L = \tfrac{1}{2}\sqrt{1 - 2\nu }/(1 - \nu) = 0.441\) for \(\nu = 0.32\). For real piezoceramics the bar mode also has a higher coupling (\(k_{33}^{\prime}\) above \(k_t\)), one reason narrow elements are used.
- \(w/t \approx 0.8\text{–}2\): the first width mode, near f \(t/c_L \approx 0.441/(w/t)\), crosses the thickness mode. The modes repel and share energy: the spectrum splits into two or three peaks inside the band.
- w/t ≫ 1: the thickness mode returns to 0.5 with a comb of weak width harmonics beside it.
The common design guideline keeps w/t below about 0.6 (the exact value depends on the material and the analysis); an element that must be wider is sub-diced: cut into two or more sub-elements that are wired together, so each piece is slender. Composites carry the same logic to pillars, whose aspect ratio and lateral period must keep the lateral (stop-band) resonances of the periodic structure well above the thickness mode. Kerf filler and the shared backing and matching layers couple neighbours mechanically; designs commonly target adjacent-element crosstalk around −30 dB or lower.
The demo below is a small 2D elastic finite-difference model (velocity-stress staggered grid, Virieux 1986) of a free bar, plane strain, isotropic with \(\nu = 0.32\), driven by a uniform thickness-direction eigenstress the way an electric field between two electrodes drives a poled plate. It runs in about 50 ms in your browser; the MATLAB version is m17_fdtd_plate_modes.m.
Finite-difference bar: spectrum, mode chart and mode shape
SimulationSolver: 2D elastic velocity-stress staggered-grid FDTD (Virieux 1986), written for this page in JS; plane strain, isotropic, \(\nu = 0.32\), free bar, 16 cells per thickness, Courant 0.5, light damping, eigenstress drive. The mode chart dots were precomputed with the same JS code (spot-checked: identical peaks at \(w/t = 0.5\), 1.0, 1.5) and agree with m17_fdtd_plate_modes.m. Not PZFlex or any commercial FE code.Lateral modes show up in channel data as a ringing tail or a second spectral peak on some elements, and in the beam as unexpected sidelobes at the lateral-mode frequency. When a measured pulse-echo spectrum of an array element shows a notch or a double peak that the KLM model of the stack does not predict, the first suspect is a lateral mode, the second a bondline.
17.8Finite-element design and the workflow
A 1D model (KLM, Mason, ABCD) is exact for an infinite plate with plane layers. It fails as soon as the structure has edges: lateral modes and kerfs (Section 17.7), crosstalk through filler, backing and lens, element directivity and baffle effects, curved lenses, non-uniform poling, and composites whose effective properties depend on pillar geometry. These need a 2D or 3D model of the actual structure.
| Approach | Examples | Strength | Cost and limits |
|---|---|---|---|
| 1D equivalent circuit / transfer matrix | KLM, Mason, this module's ABCD chain | Milliseconds; ideal for layer impedances, thicknesses, backing, electrical tuning; easy to fit to measured impedance | No edges, no lateral modes, no crosstalk, no beam |
| Explicit time-domain FE / FD | PZFlex (Weidlinger Associates, from the late 1980s; its solver was later commercialised as OnScale, which Ansys agreed to acquire in 2022), in-house FDTD codes | Large transient wave problems, broadband in one run, scales to many elements and long propagation paths | Time step capped by stability (CFL); mesh must resolve the shortest wavelength everywhere |
| Implicit / frequency-domain FE | COMSOL Multiphysics, Ansys Mechanical (harmonic and modal analyses) | Exact eigenmodes and admittance at chosen frequencies; strong multiphysics coupling | One solve per frequency; large 3D models need large memory |
No PZFlex, OnScale, COMSOL or Ansys model was run for this module, and no figure reproduces their output. The 1D results come from the KLM + ABCD chain (JS and Octave), the lateral-mode results from the small 2D FDTD of Section 17.7 (JS and Octave). That FDTD has the core ingredients of an explicit time-domain code (staggered grid, CFL-limited step, free boundaries) but no piezoelectric coupling, anisotropy, electrodes, kerf filler, backing, lens or absorbing boundaries, which are exactly what a commercial FE model adds.
Meshing and stability. Second-order schemes (central-difference FD, and lumped-mass linear FE on a uniform mesh, which reduce to the same stencil in 1D) under-estimate the phase velocity: −1.24 % at 10 points per wavelength and Courant number 0.5, −0.31 % at 20. That is 4.5° of phase per wavelength travelled at 10 points, which accumulates to more than a full cycle over 100 \(\lambda\). The mesh is set by the shortest wavelength anywhere in the model, which is usually a shear wave in the epoxies of the matching layers, backing and kerf filler, not the longitudinal wave in the ceramic (very soft silicone is usually modelled as a lossy, nearly fluid medium). An explicit scheme needs \(c_{\mathrm{max}} \Delta t/\Delta x \le 1/\sqrt{d}\) on this stencil (\(d = 2\) or 3 dimensions), so halving \(\Delta x\) doubles the number of steps: cost grows as (points per wavelength)\(^{d+1}\).
Boundaries. Real models are truncated with absorbing boundaries (perfectly matched layers, Berenger 1994, or graded damping zones) or by extending a lossy backing until its echoes are below the noise. Unchecked boundary reflections look exactly like ringing.
Outputs that are compared with measurements: electrical impedance magnitude and phase versus frequency (impedance analyzer), pulse-echo waveform and spectrum from a flat reflector (centre frequency, −6 dB bandwidth, loop sensitivity), surface displacement (laser vibrometry), crosstalk to neighbouring elements in dB versus frequency, and element directivity. Module 21 covers the measurements.
Workflow. KLM/ABCD to choose impedances, thicknesses and backing; FE to check the real geometry (kerf, sub-dicing, lens, crosstalk) and to refine; prototype; measure impedance and pulse-echo; then fit the model's material constants (\(k_t\), \(\varepsilon^S\), losses, layer speeds) to the measurement and iterate. A model you have fitted to your own measured impedance is worth more than any datasheet.
17.9Visual summary
Each card is a question an expert answers in five seconds, the answer with its number, and a figure computed from the formula or from this module's models.
17.10How experts think about it
An expert sees the stack as a cascade of transmission lines with one resonator in it. Impedances decide where energy goes (forward, backward, reflected), thicknesses decide at which frequency, and losses decide how much survives. Everything one-dimensional is a 2 × 2 matrix product and takes milliseconds; anything with an edge is a finite-element problem. The numbers below are the ones quoted in design reviews.
\(T_I = 4Z_1 Z_2 /(Z_1 + Z_2)^{2}\) depends only on the ratio \(r = Z_p /Z_L\): \(T_I = 4r/(1 + r)^{2}\). At \(r = 22.8\) (34.2/1.5) that is 0.161. A 1-3 composite at 15 MRayl \((r = 10)\) already gets 0.33 across. The curve is symmetric in log r: going from tissue into PZT is just as bad.
Each round trip (one period at resonance) multiplies the amplitude by \(r = R_{\mathrm{front}} R_{\mathrm{back}} \approx 1 - 2(Z_B + Z_L)/Z_p\), so the energy lost per cycle is ≈ \(4(Z_B + Z_L)/Z_p\) and \(Q = 2\pi/(\text{that})\). The exact form \(\pi /\ln (1/\lvert r\rvert)\) (blue) and the approximation (pink, dashed) agree within 3 % up to \(Z_B = 10\) MRayl; beyond that the approximation overestimates Q (by 14 % at 20 MRayl).
The transmission-line input impedance \(Z_1 (Z_L + jZ_1 \tan \theta)/(Z_1 + jZ_L \tan \theta)\) becomes \(Z_1^2 /Z_L\) when \(\tan \theta \to \infty\). The figure plots \(\operatorname{Re} Z_{\mathrm{in}}\) and \(\operatorname{Im} Z_{\mathrm{in}}\) versus \(f/f_0\) for the √ layer: \(\operatorname{Re} Z_{\mathrm{in}}\) rises from 1.5 to 34.2 MRayl at \(f_0\) and \(\operatorname{Im} Z_{\mathrm{in}}\) crosses zero there.
At \(\theta = \pi /2\) the reflection is ∝ \(r_1 - r_2\), zero when \(Z_1^2 = Z_p Z_L\). Off \(f_0\) the two echoes rotate apart; the band where \(1 - \lvert \Gamma \rvert ^{2}\) ≥ ½ follows the quarter-wave-transformer formula, 57.7 % here.
The plate is already a band-pass resonator. A lower-impedance layer under-matches at \(f_0\) (77 % transmission) and damps the response more evenly across the band. In this module's KLM chain the shortest −20 dB pulse occurs for \(Z_1 = 4.25\text{–}4.5\) MRayl, right on the Desilets value, while \(\sqrt{Z_p Z_L }\) gives 1.03 µs. Figure: −20 dB pulse length versus \(Z_1\) (backing 5 MRayl), with the two values marked.
Every rule is a staircase in log Z from the plate to the load. Binomial steps are symmetric (weights 1 : 2 : 1); Desilets steps are skewed towards the load because the plate's own resonance supplies part of the transformation. A grid search on the full model puts the shortest two-layer pulse at 8.0 / 2.4 MRayl, within about 11 % of Desilets.
Binomial designs with 0, 1, 2 and 3 layers (backing 5 MRayl) give −6 dB FBW of 27, 56, 79 and 84 %: the increments shrink while every layer adds attenuation, bondlines and a thickness tolerance. Three-layer designs are used when harmonic imaging or very wide bands justify the extra process steps.
The transmission peak sits where the layer is \(\lambda /4\), so it moves to \(f_0 /s\) for a thickness scaled by s (0.909 and 1.111 \(f_0\) for \(s = 1.1\) and 0.9), while the value at \(f_0\) drops only to 0.887. A thickness error looks like a centre-frequency error, not a sensitivity error. At 15 MHz a ±5 µm lapping tolerance on a 42 µm layer is ±12 %.
Inside a \(\lambda /4\) layer between \(Z_p\) and \(Z_L\) the wave bounces several times before leaving, so the path is longer than d. The single pass is \(\alpha d = 0.5 \times 5 \times 0.125 = 0.31\) dB at 5 MHz; the standing wave makes it 0.73 dB. Since \(d = c/(4f)\) and \(\alpha = \alpha_1 f\), \(\alpha d = \alpha_1 c/4\) is frequency independent. Figure: two-way loss of the Desilets-2 transducer versus layer loss (0 → 2 dB/mm/MHz costs 5.2 dB).
A thin low-impedance film between two stiff layers behaves like a series spring of compliance \(d/(\rho c^2)\). Its effect grows with \(kd = 2\pi fd/c\), so a bond that is invisible at 3 MHz dominates at 20–40 MHz. The figure plots transmission at \(f_0\) through a 3 MRayl bond under a √ layer for 2, 5 and 10 µm.
Energy sent into the backing is lost. With two Desilets layers the bandwidth is already 63 % on air backing, so heavy backing adds 20 points of FBW for 15.5 dB of two-way loss. Without layers the same sweep gives 16 % → 127 % FBW for −3.6 → −31.2 dB of loss. Figure: FBW (blue, %) and two-way loss (pink, dB) versus \(Z_B\) for Desilets-2.
A wave entering the backing returns from its rear face after \(2d_b /c_b\). Its level relative to the main pulse is roughly the rear-face reflection times \(10^{-2\alpha_b d_b/20}\) (\(\alpha_b\) in dB/mm). With 3 dB/mm at 5 MHz, 10 mm gives 60 dB, 3 mm only 18 dB. Since \(\alpha_b\) rises with frequency, the low end of the band sets the thickness. Figure: round-trip attenuation versus thickness for 2, 5 and 10 dB/mm.
Delay through thickness T is \(T(1/c_l - 1/c_t)\); a convex profile \(T_c - y^2 /(2R)\) gives the parabolic delay of a focus at F when \(1/(Fc_t) = (1/R)(1/c_l - 1/c_t)\). Because \(c_t /c_l - 1 = 0.54\) is small, a 6 % drop in tissue speed lengthens F by 20 %. Figure: F/R versus \(c_l\) for tissue at 1540 and 1450 m/s.
Multiplying a Gaussian spectrum \(\exp (-(f - f_0)^{2}/2\sigma^2)\) by \(e^{-\beta f}\) completes the square: the centre moves to \(f_0 - \beta \sigma^2\). For 0.88 mm of round-trip silicone at 3 dB/mm (5 MHz), \(\alpha \propto f^{1.5}\), and a 70 % band, the centroid moves 0.19 MHz. Same physics as tissue attenuation in Module 27.
Weight \(w(y) = 10^{-\alpha T(y)/20}\) (\(\alpha\) in dB/mm, one way) is smallest where the lens is thickest. Edge-heavy weighting narrows the main lobe a little (1.82 → 1.80 mm at −6 dB, two-way) and raises sidelobes (−26.5 → −25.7 dB). The effect grows with frequency and aperture: at 10 MHz with a 12 mm aperture the centre-to-edge difference exceeds 10 dB two-way. Figure: thickness profile and one-way weight across the elevation aperture.
The first symmetric width mode of a bar resonates near \(f \approx c_p /(2w)\), with \(c_p = 0.882 c_L\) in plane strain for \(\nu = 0.32\). It meets the thickness mode at \(w/t \approx 0.9\); around there the FDTD mode chart (dots, size ∝ strength) shows split, repelling branches. Below \(w/t \approx 0.6\) a single mode remains; above 2 the thickness mode is back near 0.5.
A plate is laterally clamped by its own extent, so it vibrates with the plate-wave stiffness \(c_L^2\). A slender bar can expand sideways, which softens it to the extensional stiffness \(E/(1 - \nu^2)\) in plane strain: \(c_p /c_L = \sqrt{1 - 2\nu }/(1 - \nu) = 0.882\). A 5.0 MHz plate design diced into \(w/t = 0.5\) elements resonates at 4.3 MHz. KLM (pink line) predicts 0.5 for every width; the FDTD's strongest mode (dots) does not.
For the second-order staggered scheme, \(\sin (\omega \Delta t/2) = S \sin (k\Delta x/2)\) with Courant number S, so the numerical phase speed falls short of c by an amount that shrinks as 1/ppw². At 10 points per wavelength and \(S = 0.5\) the phase drifts 4.5° per wavelength, 450° over 100 \(\lambda\). Epoxy shear waves (≈ 1000–1200 m/s) are about four times shorter than the ceramic's longitudinal wave (≈ 4500 m/s) and set \(\Delta x\); soft silicone, whose shear speed is only tens of m/s, is normally modelled as a lossy, nearly fluid medium so that its shear wave does not dictate the mesh.
17.11Cheat sheet
| Concept | Formula | Ultrasound use |
|---|---|---|
| Interface | \(R = (Z_2 - Z_1)/(Z_2 + Z_1)\), \(T_I = 4Z_1 Z_2 /(Z_1 + Z_2)^{2}\) | PZT → tissue: −0.916, 0.161 |
| Plate Q | \(Q = \pi /\ln (1/\lvert R_f R_b \rvert) \approx (\pi /2)Z_p /(Z_B + Z_L)\) | Ring-down, backing choice |
| \(\lambda /4\) transformer | \(Z_{\mathrm{in}} = Z_1^2 /Z_L\), \(Z_1 = \sqrt{Z_p Z_L }\) | Narrowband match, 7.16 MRayl |
| Quarter-wave band | \(\Delta f/f_0 = 2 - (4/\pi)\arccos[2\Gamma_m\sqrt{Z_p Z_L}/(\lvert Z_L - Z_p\rvert\sqrt{1 - \Gamma_m^2})]\) | 57.7 % half-power |
| Desilets | \((Z_p Z_L^2)^{1/3}\); \(Z_p^{4/7} Z_L^{3/7}\), \(Z_p^{1/7} Z_L^{6/7}\) | 4.25; 8.95 / 2.34 MRayl |
| Binomial | ln Z steps ∝ binomial coefficients | 2 layers 15.65 / 3.28; 3 layers 23.1 / 7.16 / 2.22 |
| Layer ABCD | \(\begin{bmatrix} \cosh\gamma d & Z\sinh\gamma d \\ \sinh\gamma d/Z & \cosh\gamma d \end{bmatrix}\), \(\gamma = \alpha + j\omega /c\) | Chain plate, layers, lens, load |
| Transducer gain | \(\eta = 4R_s Z_L A\lvert u_L /V_s \rvert ^{2}\), two-way loss \(20 \log _{10}\eta\) | Insertion loss, sensitivity |
| Layer thickness | \(d = c_1 /(4f_0)\) | 125 µm at 5 MHz, 41.7 µm at 15 MHz (\(c_1 = 2500\) m/s) |
| Lens focus | \(F = R/(c_t /c_l - 1)\) | 1.85 R for silicone in tissue |
| Lens sag | \(s = R - \sqrt{R^2 - (H/2)^{2}} \approx H^2 /(8R)\) | 0.209 mm for \(R = 21.6\), \(H = 6\) mm |
| Downshift | \(\Delta f = -\beta \sigma^2\) (Gaussian spectrum, linear law) | Lens and tissue filtering |
| Bar vs plate | \(c_p /c_L = \sqrt{1 - 2\nu }/(1 - \nu)\) | 0.882 for \(\nu = 0.32\) |
| Width mode | \(f \approx m c_p /(2w)\), m odd | Keep w/t ≲ 0.6 or sub-dice |
| FD dispersion | \(\sin (\omega \Delta t/2) = S \sin (k\Delta x/2)\) | 15–20 ppw at \(S \approx 0.5\) |
| Explicit stability | \(c_{\mathrm{max}} \Delta t/\Delta x \le 1/\sqrt{d}\) | \(\mathrm{cost} \propto \mathrm{ppw}^{d+1}\) |
17.12Interview and qualifying-exam questions
I1Why does a transducer need a matching layer, and what impedance would you pick for one layer on PZT into tissue? · screen
PZT at ≈ 34 MRayl against tissue at 1.5 MRayl transmits only 16 % of the power per pass and reflects 92 % of the pressure, so the plate rings. A \(\lambda /4\) layer transforms the load to \(Z_1^2 /Z_L\); \(\sqrt{Z_p Z_L } \approx 7.2\) MRayl makes that equal to \(Z_p\) at \(f_0\). For a short pulse I would start lower, at the Desilets value \((Z_p Z_L^2)^{1/3} \approx 4.3\) MRayl, and refine with a KLM model. Follow-up: "Why lower?" Because the plate is already resonant; under-matching at \(f_0\) flattens the band.
I2Two matching layers: give the impedances and say where the numbers come from. · onsite
Desilets, Fraser and Kino (1978): \(Z_1 = Z_p^{4/7} Z_L^{3/7} \approx 9\) MRayl next to the plate and \(Z_2 = Z_p^{1/7} Z_L^{6/7} \approx 2.3\) MRayl outside, from a KLM analysis of the resonant plate. The microwave binomial transformer gives 15.7 / 3.3 MRayl, which maximizes flat transmission from a semi-infinite source and in a transducer gives a wider −6 dB band but a ringier pulse. Follow-up: "What would you actually build?" Whatever filled epoxy or ceramic gets near 9 MRayl with acceptable loss, then re-optimize the second layer in the model with the real first-layer properties.
I3What does the backing do, and why not always use a backing matched to the ceramic? · screen
It absorbs the rearward wave and damps the plate: \(Q \approx (\pi /2)Z_p /(Z_B + Z_L)\). A matched backing gives the shortest pulse but sends at least half of the energy backwards on transmit and loses as much on receive (half each way only if the layers present \(Z_p\) to the plate; a bare PZT face into tissue sends about 96 % backwards). Without matching layers, going from air to a matched backing adds about 28 dB of two-way loss in the reference design (−3.6 to −31.2 dB) for 16 % → 127 % bandwidth. With two good matching layers the band is already about 63 % on air backing, so designers use moderate backing (a few MRayl) and let the layers do the broadbanding. Follow-up: the backing must also be thick and lossy enough that its rear echo is below the dynamic range.
I4Silicone lens: convex or concave, and what radius for a 40 mm elevation focus? · screen
Convex, because silicone (≈ 1000 m/s) is slower than tissue: the thick centre delays the central rays. \(F = R/(c_t /c_l - 1)\), so \(R = 40 \times 0.54 = 21.6\) mm. Follow-up: "What happens in a patient with a thick fat layer?" With 1450 m/s the same lens focuses at 48 mm; the elevation focus is set by the lens and the tissue, not by the scanner.
I5What does lens attenuation do to the image besides costing sensitivity? · onsite
Two things. It low-pass filters the pulse (centre frequency down by about \(\beta \sigma^2\), 0.2–0.3 MHz is typical at 5 MHz), which lowers axial resolution and biases spectral QUS unless the reference has the same filter. And it weights the elevation aperture towards the edges because the convex lens is thickest at the centre: an inverse apodization that narrows the elevation main lobe slightly and raises sidelobes, more so at high frequency and wide aperture. Follow-up: "How would you reduce it?" Lower-loss lens materials, thinner lens with a shorter radius at the cost of more aberration, or elevation focusing in a 1.5D/1.75D array instead.
I6Your measured centre frequency is 10 % below the KLM design. List the suspects in order. · research
(1) The lens: it downshifts the spectrum and is often left out of the model. (2) Element geometry: a slender element resonates at the bar frequency, about 12 % below the plate value for \(\nu \approx 0.3\) (FDTD: 0.440 against 0.5). (3) Layer thickness or speed: a layer 10 % thick moves its band 10 % down. (4) Material constants: \(v^D\), \(k_t\) and \(\varepsilon^S\) differ from the datasheet after dicing and poling. (5) Bondlines at high frequency. The test that separates them: measure the electrical impedance of a bare element (\(f_r\), \(f_a\) give \(k_{\mathrm{eff}}\) and v), then of the stack without lens, then with lens, and fit each step.
I7What are lateral modes and how do you avoid them in an array? · onsite
Width resonances of the element bar, near \(f \approx c_p /(2w)\), that couple to the thickness mode through Poisson's ratio when w is comparable with t. They split the spectrum and add ringing and spurious beam lobes. The guideline keeps w/t below about 0.6; wider elements are sub-diced into two or more slender pieces wired in parallel. In composites the pillar aspect ratio and the lateral period play the same role. Follow-up: "How would you check?" A 2D/3D FE model of the element with kerf and filler, then laser vibrometry or the impedance spectrum of a diced sample.
I8When is a KLM model enough, and when do you need finite elements? · onsite
KLM/ABCD is exact for a laterally infinite plate with plane layers: it is the right tool to choose layer impedances and thicknesses, backing and electrical tuning, and it fits measured impedance well for large single elements. FE is needed for anything with edges: array elements with w/t near 1, kerf and sub-dicing, crosstalk through filler, backing and lens, curved lenses, element directivity, composites with non-trivial pillar geometry. Follow-up: "Explicit or implicit?" Explicit time-domain for broadband transient and large models; implicit frequency-domain for exact modes, admittance at a few frequencies and strongly coupled multiphysics.
I9You build an explicit FE model of a 7.5 MHz probe with a silicone lens. How do you choose the mesh and time step? · research
Find the shortest wavelength anywhere at the highest frequency of interest, including shear waves in the epoxies (treat the soft silicone lens as a lossy, nearly fluid medium) and the upper band edge (say 12 MHz). Use 15–20 elements per that wavelength for second-order schemes (−0.3 to −0.55 % phase speed error at \(S = 0.5\)), more if the waves travel many wavelengths. The time step follows from the fastest wave (ceramic longitudinal) and the smallest element through the CFL limit, \(c_{\mathrm{max}} \Delta t/\Delta x \le 1/\sqrt{d}\) for this stencil (codes add safety factors). Then check convergence by refining once and comparing impedance and pulse-echo. Follow-up: cost grows as \(\mathrm{ppw}^4\) in 3D, so exploit symmetry and model one element plus neighbours.
I10Design a matching stack for a 15 MHz probe on a 1-3 composite at 15 MRayl. What changes compared with PZT? · research
The mismatch is smaller (\(T_I = 0.33\) bare), so lower layer impedances: Desilets one layer \((15 \times 2.25)^{1/3} = 3.2\) MRayl, two layers \(15^{4/7}\,1.5^{3/7} = 5.6\) and \(15^{1/7}\,1.5^{6/7} = 2.1\) MRayl, all reachable with filled epoxies and parylene. At 15 MHz a 2500 m/s layer is 42 µm thick, so lapping tolerance and bondlines (keep them to 1–2 µm) dominate, and layer loss in dB stays the same as at 5 MHz for \(\alpha \propto f\). The composite's lower Z and higher \(k_t\) give bandwidth with lighter backing. Follow-up: check the composite's lateral (inter-pillar) resonances against 2 × 15 MHz if harmonic imaging is planned.
17.13MATLAB lab
Base MATLAB only (no toolboxes); each script runs as-is in MATLAB R2019b+ and GNU Octave 8. Paste into the editor and press Run.
m17_matching_layers.mKLM plate + ABCD layers: table of \(f_c\), FBW, pulse length and loss for seven designs, single-layer sweep · change ZB, Rs, alpha_dB, the design list · runtime ≈ 8 s
%% m17_matching_layers.m -- KLM + ABCD model of a backed plate with 0-3 matching layers
% Purpose : compare matching-layer design rules (quarter-wave sqrt, Desilets
% 1978, binomial) with one self-contained KLM/ABCD model and print
% centre frequency, -6 dB fractional bandwidth, -20 dB pulse length
% and two-way insertion loss (perfect reflector at the face).
% Model : PZT-5H-like plate (rho 7500 kg/m^3, v 4560 m/s, kt 0.50,
% eps_S/eps0 1470), lambda/2 at 5 MHz, area 20 mm^2, 50 ohm
% source and receiver, backing 5 MRayl, load 1.5 MRayl, layers
% lambda/4 at 5 MHz in a c = 2500 m/s material, lossless.
% Change : ZB (backing), Rs, the design list, layer loss alpha_dB.
% Expect : none 27 % / 1.25 us / -12.8 dB; sqrt 56 % / 1.03 us / -5.6 dB;
% Desilets-1 44 % / 0.78 us / -5.4 dB; Desilets-2 66 % / 0.81 us;
% binomial-2 79 % / 1.00 us; shortest single-layer pulse at
% 4.25-4.5 MRayl (the Desilets value). Runtime ~ 8 s (Octave).
% Note : fc, FBW and loss are read in the fundamental band (peak searched
% below 9 MHz), so the 3rd-harmonic resonance near 15 MHz is never
% mistaken for the main band (it matters for light backings).
clear; close all; rng(1);
%% Parameters
rho = 7500; v = 4560; kt = 0.50; epsS = 1470*8.854e-12; % piezo
f0 = 5e6; t = v/(2*f0); A = 20e-6; % plate thickness (m), area (m^2)
Zp = rho*v/1e6; % 34.2 MRayl
ZB = 5; ZL = 1.5; Rs = 50; % MRayl, MRayl, ohm
cLay = 2500; alpha_dB = 0; % layer speed (m/s), loss (dB/mm/MHz)
N = 2048; fs = 100e6; df = fs/N; f = (1:N/2-1)*df; % frequency grid (Hz)
% one ABCD for a (lossy) line: Z (Rayl*m^2), length d (m), speed c, loss a (Np/m)
lineM = @(Z, d, c, fr, a) [cosh((a + 1i*2*pi*fr/c)*d), Z*sinh((a + 1i*2*pi*fr/c)*d); ...
sinh((a + 1i*2*pi*fr/c)*d)/Z, cosh((a + 1i*2*pi*fr/c)*d)];
zin = @(Mx, Zt) (Mx(1,1)*Zt + Mx(1,2)) / (Mx(2,1)*Zt + Mx(2,2));
% -20 dB envelope duration with linear interpolation of both crossings (us)
dur20 = @(e, i1, i2, th) ((i2 + (e(i2) - th)/(e(i2) - e(i2+1))) - (i1 - 1 + (th - e(i1-1))/(e(i1) - e(i1-1))))/fs*1e6;
%% Designs (impedances in MRayl)
R = @(e) Zp^e * ZL^(1-e);
names = {'none', 'sqrt (QW)', 'Desilets-1', 'Desilets-2', 'binomial-2', 'binomial-3', 'grid best-2'};
Zs = {[], R(1/2), R(1/3), [R(4/7) R(1/7)], [R(3/4) R(1/4)], [R(7/8) R(1/2) R(1/8)], [8.0 2.4]};
C0 = epsS*A/t; Z0 = rho*v*A; h2 = kt^2*rho*v^2/epsS; % clamped capacitance, line impedance, h^2
Hpe = zeros(numel(Zs), N/2-1); eta = Hpe;
for q = 1:numel(Zs)
for k = 1:numel(f)
w = 2*pi*f(k); th = w*t/v;
X1 = h2*sin(th)/(w^2*Z0); phi = (w*Z0/(2*sqrt(h2)))/sin(th/2);
half = lineM(Z0, t/2, v, f(k), 0);
Zb = zin(half, ZB*1e6*A); % back branch seen from the centre node
Mf = half;
for L = 1:numel(Zs{q})
Mf = Mf*lineM(Zs{q}(L)*1e6*A, cLay/(4*f0), cLay, f(k), alpha_dB*f(k)/1e6*1000/8.686);
end
T = [1, 1i*(X1 - 1/(w*C0)); 0, 1] * [1/phi, 0; 0, phi] * [1, 0; 1/Zb, 1] * Mf;
ZLa = ZL*1e6*A;
u = 1/((T(1,1)*ZLa + T(1,2)) + Rs*(T(2,1)*ZLa + T(2,2))); % load velocity per source volt
F = ZLa*u;
eta(q,k) = 4*Rs*abs(u)^2*ZLa; % one-way transducer power gain
Hpe(q,k) = 2*Rs/ZLa*F^2; % pulse-echo Vout/Vs (reciprocity, det = 1)
end
end
%% Metrics
tvec = (0:N-1)/fs*1e6; tau0 = 2e-6; % time axis (us), display delay (s)
% receive low-pass for the time-domain view: 1 below 9 MHz, cos^2 taper to 0 at 12 MHz
% (removes the 3rd-harmonic thickness resonance near 14 MHz; the fundamental band is untouched)
Glp = (f < 9e6) + (f >= 9e6 & f < 12e6).*cos(pi/2*(f - 9e6)/3e6).^2;
fprintf('%-12s %-20s %7s %7s %9s %8s\n', 'design', 'Z (MRayl)', 'fc MHz', 'FBW %', 'L-20 us', 'IL2 dB');
rf = zeros(numel(Zs), N);
for q = 1:numel(Zs)
m = abs(Hpe(q,:)); [pk, ip] = max(m.*(f < 9e6)); % fundamental band only
lo = ip; while lo > 1 && m(lo) > pk/2, lo = lo - 1; end
hi = ip; while hi < numel(m) && m(hi) > pk/2, hi = hi + 1; end
fl = interp1(m(lo:lo+1), f(lo:lo+1), pk/2); fh = interp1(m(hi-1:hi), f(hi-1:hi), pk/2);
fc = (fl + fh)/2; fbw = (fh - fl)/fc;
S = zeros(1, N); S(2:N/2) = 2*Hpe(q,:).*Glp.*exp(-1i*2*pi*f*tau0); % analytic-signal spectrum
a = ifft(S); env = abs(a); rf(q,:) = real(a);
idx = find(env > 0.1*max(env)); L20 = dur20(env, idx(1), idx(end), 0.1*max(env));
zstr = sprintf('%.2f ', Zs{q}); if isempty(Zs{q}), zstr = '-'; end
fprintf('%-12s %-20s %7.2f %7.0f %9.2f %8.1f\n', names{q}, zstr, fc/1e6, 100*fbw, L20, 20*log10(2*pk));
end
%% Single layer: which Z gives the shortest pulse?
Zgrid = 2:0.25:12; L20g = zeros(size(Zgrid)); fbwg = L20g;
for g = 1:numel(Zgrid)
H = zeros(1, numel(f));
for k = 1:numel(f)
w = 2*pi*f(k); th = w*t/v; X1 = h2*sin(th)/(w^2*Z0); phi = (w*Z0/(2*sqrt(h2)))/sin(th/2);
half = lineM(Z0, t/2, v, f(k), 0); Zb = zin(half, ZB*1e6*A);
Mf = half*lineM(Zgrid(g)*1e6*A, cLay/(4*f0), cLay, f(k), 0);
T = [1, 1i*(X1 - 1/(w*C0)); 0, 1] * [1/phi, 0; 0, phi] * [1, 0; 1/Zb, 1] * Mf;
ZLa = ZL*1e6*A; u = 1/((T(1,1)*ZLa + T(1,2)) + Rs*(T(2,1)*ZLa + T(2,2)));
H(k) = 2*Rs/ZLa*(ZLa*u)^2;
end
S = zeros(1, N); S(2:N/2) = 2*H.*Glp.*exp(-1i*2*pi*f*tau0); env = abs(ifft(S));
idx = find(env > 0.1*max(env)); L20g(g) = dur20(env, idx(1), idx(end), 0.1*max(env));
m = abs(H); [pk, ip] = max(m.*(f < 9e6)); lo = ip; while lo > 1 && m(lo) > pk/2, lo = lo - 1; end
hi = ip; while hi < numel(m) && m(hi) > pk/2, hi = hi + 1; end
fbwg(g) = (f(hi) - f(lo))/((f(hi) + f(lo))/2);
end
[Lmin, ib] = min(L20g);
fprintf('\nSingle layer: shortest -20 dB pulse %.2f us at Z1 = %.2f MRayl (Desilets %.2f, sqrt %.2f)\n', Lmin, Zgrid(ib), R(1/3), R(1/2));
fprintf(' pulse length at Z1 = 4.25: %.2f us, at 7.25: %.2f us (the minimum is flat from 4 to 5 MRayl)\n', L20g(Zgrid == 4.25), L20g(Zgrid == 7.25));
[~, iw] = max(fbwg);
fprintf('Single layer: widest -6 dB band %.0f %% at Z1 = %.2f MRayl (ripple counts as band)\n', 100*fbwg(iw), Zgrid(iw));
%% Figures
figure('Position', [80 80 980 640]);
subplot(2,2,1); plot(f/1e6, 20*log10(eta([1 2 3 4],:))', 'LineWidth', 1.3); grid on; xlim([0 10]); ylim([-40 0]);
xlabel('frequency (MHz)'); ylabel('20 log_{10} \eta(f) (dB, two-way)'); title('Pulse-echo loss, backing 5 MRayl');
legend(names([1 2 3 4]), 'Location', 'south');
subplot(2,2,2); plot(f/1e6, 20*log10(eta([4 5 6 7],:))', 'LineWidth', 1.3); grid on; xlim([0 10]); ylim([-40 0]);
xlabel('frequency (MHz)'); ylabel('dB'); title('Two and three layers'); legend(names([4 5 6 7]), 'Location', 'south');
subplot(2,2,3); sel = [1 3 4 5];
for k = 1:4, plot(tvec - tau0*1e6, rf(sel(k),:)/max(abs(rf(sel(k),:))) - 2.2*(k-1), 'LineWidth', 1.1); hold on; end
xlim([-0.3 3]); grid on; xlabel('time (\mus)'); ylabel('normalised, offset'); title('Pulse-echo impulse responses');
legend(names(sel), 'Location', 'southeast');
subplot(2,2,4); [ax, h1, h2] = plotyy(Zgrid, L20g, Zgrid, 100*fbwg); grid on;
xlabel('single-layer Z_1 (MRayl)'); ylabel(ax(1), '-20 dB pulse length (\mus)'); ylabel(ax(2), '-6 dB FBW (%)');
title('One layer: shortest pulse vs widest band');
What to look for: the printed table matches Section 17.3 (none 27 % / −12.8 dB, √ 56 % / −5.6 dB, Desilets-1 44 % / 0.78 µs, Desilets-2 66 %, binomial-2 79 % / 1.00 µs); the single-layer sweep finds the shortest pulse at 4.25–4.5 MRayl and the widest −6 dB band at 10 MRayl, where the passband dip just clears −6 dB.
m17_quarter_wave_tolerance.mOne \(\lambda /4\) layer: band against the closed form, thickness error, resonant layer absorption, glue bondlines · change Zp, ZL, bond thickness · runtime ≈ 1 s
%% m17_quarter_wave_tolerance.m -- one quarter-wave layer: bandwidth, tolerance, loss, bondline
% Purpose : semi-infinite piezo (Zp) radiating into tissue (ZL) through one
% layer. Checks Z1 = sqrt(Zp ZL), the -3 dB transmission band against
% the closed-form quarter-wave-transformer formula, and what thickness
% error, layer attenuation and a thin glue bondline do.
% Change : Zp (a 1-3 composite is 10-25 MRayl), ZL, bond thickness, alpha.
% Expect : T(no layer) = 0.161 (-7.9 dB); -3 dB band 0.711-1.289 f0 = 57.7 %
% (formula 57.7 %); +-10 % thickness error: T(f0) = 0.887, band
% centre moves by -+10 %; alpha = 0.5 dB/mm/MHz absorbs 0.73 dB one
% way at ANY f0 (alpha ~ f, d ~ 1/f); 5 um bond (3 MRayl): -0.5 dB at
% 5 MHz, -4.7 dB at 20 MHz. Runtime ~ 1 s.
clear; close all; rng(1);
%% Media (MRayl) and helpers
Zp = 34.2; ZL = 1.5; Z1 = sqrt(Zp*ZL);
f = linspace(0.001, 2, 4001); % frequency / f0
% ABCD of a line with electrical angle th (rad), impedance Z, loss factor a*d (Np)
lineM = @(Z, th, ad) [cosh(ad + 1i*th), Z*sinh(ad + 1i*th); sinh(ad + 1i*th)/Z, cosh(ad + 1i*th)];
Tpow = @(Zin) 1 - abs((Zin - Zp)./(Zin + Zp)).^2; % power transmission from Zp into the stack
%% 1. No layer and the ideal quarter-wave layer
T0 = 4*Zp*ZL/(Zp + ZL)^2;
th = pi/2*f; % electrical length, pi/2 at f0
Zin = Z1*(ZL + 1i*Z1*tan(th))./(Z1 + 1i*ZL*tan(th));
T1 = Tpow(Zin);
band = f(T1 >= 0.5);
Gm = sqrt(0.5); % |Gamma| at the -3 dB transmission edge
bwFormula = 2 - 4/pi*acos(2*Gm*sqrt(Zp*ZL)/(abs(ZL - Zp)*sqrt(1 - Gm^2)));
fprintf('No layer: T = %.3f (%.1f dB), R_p = %.3f\n', T0, 10*log10(T0), (ZL - Zp)/(ZL + Zp));
fprintf('Z1 = sqrt(Zp ZL) = %.2f MRayl; Zin(f0) = Z1^2/ZL = %.2f MRayl\n', Z1, Z1^2/ZL);
fprintf('-3 dB transmission band: %.3f to %.3f f0 = %.1f %% (closed form %.1f %%)\n', band(1), band(end), 100*(band(end) - band(1)), 100*bwFormula);
r1 = (Z1 - Zp)/(Z1 + Zp); r2 = (ZL - Z1)/(ZL + Z1);
fprintf('Two interface reflections at Z1 = sqrt: r1 = %.3f, r2 = %.3f (equal, cancel at f0)\n', r1, r2);
%% 2. Thickness error: scale the electrical length by s
sv = [0.9 0.95 1.05 1.1]; Ts = zeros(numel(sv), numel(f));
for k = 1:numel(sv)
tk = pi/2*f*sv(k); Ts(k,:) = Tpow(Z1*(ZL + 1i*Z1*tan(tk))./(Z1 + 1i*ZL*tan(tk)));
[~, ip] = max(Ts(k,:));
[~, i0] = min(abs(f - 1));
fprintf('thickness x %.2f: T(f0) = %.3f (%.2f dB), peak at %.3f f0\n', sv(k), Ts(k, i0), 10*log10(Ts(k, i0)), f(ip));
end
fprintf('quarter-wave thickness in a 2500 m/s layer: %.0f um at 5 MHz, %.1f um at 15 MHz\n', 2500/(4*5e6)*1e6, 2500/(4*15e6)*1e6);
%% 3. Layer attenuation (alpha in dB/mm/MHz, linear in f), loss at f0
c1 = 2500; al = [0 0.5 1 2];
for f0 = [5e6 15e6]
d = c1/(4*f0); % quarter-wave thickness (m)
for a = al
ad = a*(f0/1e6)*d*1e3/8.686; % Np through the layer
M = lineM(Z1, pi/2, ad); v = M*[1; 1/ZL]; % [p_in; u_in] for p_L = 1, u_L = 1/ZL
Zi = v(1)/v(2); absorb = (1/ZL)/real(v(1)*conj(v(2))); % P_load / P_into_layer
fprintf('f0 = %4.1f MHz, alpha = %.1f dB/mm/MHz: absorption %.2f dB, total one-way %.2f dB\n', f0/1e6, a, -10*log10(absorb), -10*log10(absorb*Tpow(Zi)));
end
end
%% 4. Glue bondline between piezo and layer (Zb = 3 MRayl, cb = 2500 m/s)
f0s = linspace(1e6, 40e6, 400); bonds = [2 5 10]*1e-6; Tb = zeros(numel(bonds), numel(f0s));
for b = 1:numel(bonds)
for k = 1:numel(f0s)
Mb = lineM(3, 2*pi*f0s(k)*bonds(b)/2500, 0); Ml = lineM(Z1, pi/2, 0); M = Mb*Ml;
Tb(b,k) = Tpow((M(1,1)*ZL + M(1,2))/(M(2,1)*ZL + M(2,2)));
end
end
for fq = [5e6 20e6 30e6]
[~, k] = min(abs(f0s - fq));
fprintf('bondline at %2.0f MHz: 2 um %.2f dB, 5 um %.2f dB, 10 um %.2f dB\n', fq/1e6, 10*log10(Tb(:,k)));
end
%% Figures
figure('Position', [80 80 1000 330]);
subplot(1,3,1); plot(f, T1, 'LineWidth', 1.6); hold on; plot(f, T0 + 0*f, '--'); plot([band(1) band(end)], [0.5 0.5], 'k-o');
grid on; xlabel('f / f_0'); ylabel('power transmission'); title('One \lambda/4 layer, Z_1 = 7.16 MRayl'); legend('\lambda/4 layer', 'no layer', '-3 dB band', 'Location', 'south');
subplot(1,3,2); plot(f, Ts', 'LineWidth', 1.2); hold on; plot(f, T1, 'k', 'LineWidth', 1.6); grid on; xlim([0.5 1.5]); ylim([0.4 1.02]);
xlabel('f / f_0'); ylabel('power transmission'); title('Thickness error x0.9 ... x1.1'); legend('0.90', '0.95', '1.05', '1.10', 'nominal', 'Location', 'south');
subplot(1,3,3); plot(f0s/1e6, 10*log10(Tb'), 'LineWidth', 1.4); grid on; xlabel('design frequency f_0 (MHz)'); ylabel('transmission at f_0 (dB)');
title('Glue bondline (3 MRayl) under the layer'); legend('2 \mum', '5 \mum', '10 \mum', 'Location', 'southwest');
What to look for: band 0.712–1.288 \(f_0\) = 57.7 % (closed form 57.7 %); \(r_1 = r_2 = -0.654\); ±10 % thickness gives \(T(f_0) = 0.887\) with the peak at 1.111 / 0.909 \(f_0\); 0.5 dB/mm/MHz absorbs 0.73 dB at both 5 and 15 MHz; a 5 µm bond costs 0.52 dB at 5 MHz and 4.73 dB at 20 MHz.
m17_lens_design.mSilicone lens: paraxial vs exact ray-trace focus, inverse apodization, elevation beam, spectral downshift · change cl, ct, H, alpha5, nexp · runtime ≈ 2 s
%% m17_lens_design.m -- silicone elevation lens: focal length, aberration, apodization, downshift
% Purpose : design a plano-convex RTV-silicone lens for a 40 mm elevation focus,
% check the paraxial focal length F = R/(c_t/c_l - 1) against an exact
% Snell ray trace, and quantify what lens attenuation does: a
% centre-heavy loss (inverse apodization) and a downshift of the
% pulse-echo spectrum.
% Change : cl (900-1100 m/s), ct, H (elevation aperture), alpha5, nexp, te.
% Expect : R = 21.6 mm, sag 0.209 mm, centre 0.509 mm; marginal ray crosses
% the axis 0.92 mm short of the paraxial focus; two-way loss at 5 MHz
% 3.06 dB (centre) vs 1.80 dB (edge); -6 dB elevation width at focus
% 1.82 mm uniform, 1.80 mm lens-weighted, first sidelobe -26.5 ->
% -25.7 dB; spectrum centroid 5.00 -> 4.81 MHz; F = 48.0 mm in fat-like
% 1450 m/s tissue.
% Runtime < 2 s.
clear; close all; rng(1);
%% Parameters
cl = 1000; ct = 1540; % lens and tissue sound speed (m/s)
Ftarget = 40e-3; % elevation focus (m)
H = 6e-3; % elevation aperture (m)
te = 0.3e-3; % lens thickness at the elevation edge (m)
f0 = 5e6; % centre frequency (Hz)
alpha5 = 3.0; nexp = 1.5; % lens attenuation: alpha5 dB/mm at 5 MHz, alpha ~ f^nexp
lam = ct/f0;
%% Geometry
R = Ftarget*(ct/cl - 1); % radius of curvature for the target focus
sag = R - sqrt(R^2 - (H/2)^2); Tc = te + sag; % sagitta and centre thickness
T = @(y) Tc - (R - sqrt(R^2 - y.^2)); % lens thickness at elevation y
fprintf('R = %.2f mm, sag = %.3f mm, centre thickness = %.3f mm, F# = %.2f\n', R*1e3, sag*1e3, Tc*1e3, Ftarget/H);
%% Exact ray trace (rays leave the face normally, refract at the spherical surface)
ys = linspace(0.05e-3, H/2, 40); zc = zeros(size(ys));
for k = 1:numel(ys)
y = ys(k); zs = T(y); % surface point (z measured from the face)
nrm = [y, zs - (Tc - R)]/R; % outward unit normal (y, z), centre at z = Tc - R
thi = asin(abs(nrm(1))); % incidence angle (ray direction is +z)
tht = asin(ct/cl*sin(thi)); % Snell: sin(tht)/ct = sin(thi)/cl
dev = tht - thi; % ray turns toward the axis by dev
zc(k) = zs + y/tan(dev); % axis crossing (from the face)
end
Fpar = Tc + R/(ct/cl - 1); % paraxial focus measured from the face
fprintf('Paraxial focus %.2f mm from the face (F = %.2f mm from the apex)\n', Fpar*1e3, (Fpar - Tc)*1e3);
fprintf('Marginal ray (y = %.1f mm) crosses at %.2f mm: spherical aberration %.2f mm\n', H/2*1e3, zc(end)*1e3, (zc(end) - Fpar)*1e3);
for ctt = [1450 1540 1580]
fprintf(' tissue %4d m/s -> F = %.1f mm\n', ctt, R/(ctt/cl - 1)*1e3);
end
%% Attenuation weighting and elevation beam at the focus
alpha = @(fr) alpha5*(fr/5e6).^nexp; % dB/mm
y = linspace(-H/2, H/2, 601); dy = y(2) - y(1);
w = 10.^(-alpha(f0)*T(y)*1e3/20); % one-way amplitude weight
fprintf('Two-way loss at %.0f MHz: centre %.2f dB, edge %.2f dB\n', f0/1e6, 2*alpha(f0)*Tc*1e3, 2*alpha(f0)*te*1e3);
x = linspace(-4e-3, 4e-3, 801);
Bu = abs(exp(-1i*2*pi*x(:)*y/(lam*Ftarget))*ones(numel(y),1))*dy;
Bw = abs(exp(-1i*2*pi*x(:)*y/(lam*Ftarget))*w(:))*dy;
P2u = (Bu/max(Bu)).^2; P2w = (Bw/max(Bw)).^2; % two-way (transmit = receive)
wid = @(p) 2*interp1(p(x >= 0), x(x >= 0), 0.5); % -6 dB full width (two-way amplitude)
fprintf('Two-way -6 dB elevation width: uniform %.2f mm (0.886 lam F/H = %.2f), lens-weighted %.2f mm\n', ...
wid(P2u)*1e3, 0.886*lam*Ftarget/H*1e3, wid(P2w)*1e3);
pk = @(p) max(p(x > lam*Ftarget/H*0.9 & x < 2.2*lam*Ftarget/H));
fprintf('First sidelobe (two-way): uniform %.1f dB, lens-weighted %.1f dB\n', 20*log10(pk(P2u)), 20*log10(pk(P2w)));
%% Spectral downshift (Gaussian pulse-echo spectrum, 70 % -6 dB bandwidth)
fr = linspace(0.5e6, 10e6, 2000); sig = 0.70*f0/(2*sqrt(2*log(2)));
S0 = exp(-(fr - f0).^2/(2*sig^2));
Tbar = mean(T(y)); % aperture-mean thickness
S1 = S0.*10.^(-2*alpha(fr)*Tbar*1e3/20); % two-way through the mean thickness
cen = @(S) sum(fr.*S)/sum(S);
[~, i1] = max(S1);
fprintf('Mean thickness %.3f mm; centroid %.2f -> %.2f MHz; peak %.2f -> %.2f MHz\n', Tbar*1e3, cen(S0)/1e6, cen(S1)/1e6, f0/1e6, fr(i1)/1e6);
b = 2*alpha5/5*Tbar*1e3/8.686; % Np/MHz if alpha were linear (nexp = 1)
fprintf('Linear-law estimate beta*sigma^2 = %.2f MHz (beta = %.3f Np/MHz)\n', b*(sig/1e6)^2, b);
%% Figures
figure('Position', [60 60 1000 620]);
subplot(2,2,1); hold on;
for k = 1:4:numel(ys)
plot([0 T(ys(k))]*1e3, [ys(k) ys(k)]*1e3, 'b'); plot([T(ys(k)) zc(k)]*1e3, [ys(k) 0]*1e3, 'b');
plot([0 T(ys(k))]*1e3, -[ys(k) ys(k)]*1e3, 'b'); plot([T(ys(k)) zc(k)]*1e3, -[ys(k) 0]*1e3, 'b');
end
plot(Fpar*1e3, 0, 'ro', 'MarkerFaceColor', 'r'); grid on; xlabel('depth z (mm)'); ylabel('elevation y (mm)');
title(sprintf('Exact rays, R = %.1f mm, c_l = %d m/s', R*1e3, cl)); xlim([0 45]);
subplot(2,2,2); plot(y*1e3, T(y)*1e3, 'LineWidth', 1.5); hold on; plot(y*1e3, w, 'LineWidth', 1.5); grid on;
xlabel('elevation y (mm)'); legend('thickness (mm)', 'one-way weight at f_0', 'Location', 'south'); title('Lens profile and attenuation weight');
subplot(2,2,3); plot(x*1e3, 20*log10(P2u + 1e-6), x*1e3, 20*log10(P2w + 1e-6), 'LineWidth', 1.3); grid on; ylim([-60 0]);
xlabel('elevation x at focus (mm)'); ylabel('two-way (dB)'); legend('uniform', 'lens-weighted'); title('Elevation beam at 40 mm');
subplot(2,2,4); plot(fr/1e6, S0, fr/1e6, S1/max(S1), 'LineWidth', 1.3); grid on; xlabel('frequency (MHz)'); ylabel('normalised');
legend('before lens', 'after lens (two-way, normalised)'); title('Lens loss is a low-pass filter');
What to look for: \(R = 21.6\) mm, sag 0.209 mm; marginal ray 0.92 mm short of the paraxial focus; \(F = 48.0\) mm at 1450 m/s; two-way loss 3.06 dB (centre) and 1.80 dB (edge); −6 dB width 1.82 mm (uniform, matching \(0.886\lambda F/H\)) and 1.80 mm (lens-weighted); centroid 5.00 → 4.81 MHz.
m17_fdtd_plate_modes.m2D elastic staggered-grid FDTD of a free bar for 31 widths: face-velocity spectra and the mode chart · change nu, ppt, the w/t list · runtime ≈ 15 s
%% m17_fdtd_plate_modes.m -- 2D elastic FDTD of a free piezo-like bar: lateral modes vs w/t
% Purpose : show why 1D (KLM) models fail for array elements. A bar of width w
% and thickness t (plane strain: the element is long in elevation) is
% driven by a uniform thickness-direction eigenstress, like the field
% between two electrodes. The spectrum of the mean face velocity
% shows the thickness mode, the width modes, and their coupling.
% Method : velocity-stress staggered grid (Virieux 1986), isotropic solid,
% nu = 0.32 (ceramic-like), free surfaces by the vacuum formulation,
% 16 cells per thickness, Courant number 0.5. Normalised units:
% t = 1, cL = 1, rho = 1, so f t / cL = 0.5 is the plate thickness
% mode (5 MHz for a 0.456 mm plate with cL = 4560 m/s).
% Change : wt list, nu, ppt (cells per thickness), Tend.
% Expect : w/t = 0.2 -> single resonance at f t/cL ~ 0.44 (thin bar, theory
% 0.441); w/t ~ 0.8-2 -> split, coupled modes; w/t = 6 -> ~0.5 with
% width harmonics. First width mode follows 0.441/(w/t).
% Runtime ~ 15 s in Octave.
clear; close all; rng(1);
nu = 0.32; ppt = 16; Tend = 60; damp = 0.002;
cL = 1; cS = cL*sqrt((1 - 2*nu)/(2*(1 - nu))); mu = cS^2; lam = cL^2 - 2*mu;
dx = 1/ppt; dt = 0.5*dx/cL; nt = round(Tend/dt); g = 1 - damp;
sig = 0.35/0.5; t0 = 4*sig; gp = @(u) exp(-0.5*((u - t0)/sig).^2); % Gaussian eigenstress drive
fr = 0.02:0.005:1.0; % analysis frequencies (f t / cL)
win = 0.5 - 0.5*cos(2*pi*(0:nt-1)/nt); % Hann window on the record
E = exp(-1i*2*pi*(0:nt-1)'*dt*fr); % DFT matrix (nt x nf)
wts = [0.2 0.3:0.1:3.0 4 6];
chart = []; specs = zeros(numel(wts), numel(fr));
tic;
for q = 1:numel(wts)
wt = wts(q); nxb = max(2, round(wt*ppt)); nzb = ppt; pad = 2;
nx = nxb + 2*pad; nz = nzb + 2*pad;
S = false(nz, nx); S(pad+1:pad+nzb, pad+1:pad+nxb) = true; % solid cells (rows = z, cols = x)
Sx = [S(:, 2:end), false(nz, 1)]; Sz = [S(2:end, :); false(1, nx)]; Sxz = [Sz(:, 2:end), false(nz, 1)];
rx = S + Sx; bx = zeros(nz, nx); bx(rx > 0) = 2./rx(rx > 0); % buoyancy at vx nodes
rz = S + Sz; bz = zeros(nz, nx); bz(rz > 0) = 2./rz(rz > 0); % buoyancy at vz nodes
Mxz = mu*(S & Sx & Sz & Sxz); % shear modulus at sxz nodes
L2M = (lam + 2*mu)*S; LL = lam*S;
vx = zeros(nz, nx); vz = vx; sxx = vx; szz = vx; sxz = vx; rec = zeros(1, nt);
I = 2:nx-1; J = 2:nz-1; top = pad + nzb; % row of the top face vz nodes
for n = 1:nt
tt = (n - 1)*dt; sv = gp(tt) - gp(tt - dt);
vx(J,I) = g*vx(J,I) + dt/dx*bx(J,I).*(sxx(J,I+1) - sxx(J,I) + sxz(J,I) - sxz(J-1,I));
vz(J,I) = g*vz(J,I) + dt/dx*bz(J,I).*(sxz(J,I) - sxz(J,I-1) + szz(J+1,I) - szz(J,I));
ex = (vx(J,I) - vx(J,I-1))/dx; ez = (vz(J,I) - vz(J-1,I))/dx;
sxx(J,I) = sxx(J,I) + dt*(L2M(J,I).*ex + LL(J,I).*ez);
szz(J,I) = szz(J,I) + dt*(LL(J,I).*ex + L2M(J,I).*ez) + sv*S(J,I);
sxz(J,I) = sxz(J,I) + dt*Mxz(J,I).*((vx(J+1,I) - vx(J,I)) + (vz(J,I+1) - vz(J,I)))/dx;
rec(n) = mean(vz(top, pad+1:pad+nxb));
end
sp = abs((rec.*win)*E); sp = sp/max(sp); specs(q,:) = sp;
pk = find(sp(2:end-1) > sp(1:end-2) & sp(2:end-1) >= sp(3:end) & sp(2:end-1) > 0.1) + 1;
chart = [chart; [wt*ones(numel(pk), 1), fr(pk)', sp(pk)']]; %#ok<AGROW>
band = pk(fr(pk) > 0.35 & fr(pk) < 0.62); [~, ib] = max(sp(band));
if any(abs(wt - [0.2 0.6 1 1.5 2 3 6]) < 1e-9)
fprintf('w/t = %.1f: peaks at f t/cL = %s (strongest in 0.35-0.62: %.3f)\n', wt, sprintf('%.3f ', fr(pk)), fr(band(ib)));
end
end
fprintf('thin-bar limit (plane strain) f t/cL = 0.5*sqrt((1-2nu))/(1-nu) = %.3f; plate = 0.500\n', 0.5*sqrt(1 - 2*nu)/(1 - nu));
fprintf('first width mode estimate 0.441/(w/t): w/t = 1.5 -> %.3f, 2 -> %.3f, 3 -> %.3f\n', 0.441/1.5, 0.441/2, 0.441/3);
fprintf('total FDTD time %.1f s for %d bars\n', toc, numel(wts));
%% Figures
figure('Position', [60 60 1000 400]);
subplot(1,2,1); sel = find(ismember(round(wts*10), [2 6 10 15 30]));
for k = 1:numel(sel), plot(fr, specs(sel(k),:) - 1.1*(k-1), 'LineWidth', 1.3); hold on; end
grid on; xlabel('f t / c_L (0.5 = 5 MHz for t = 0.456 mm)'); ylabel('normalised |v_z| (offset)');
legend(arrayfun(@(v) sprintf('w/t = %.1f', v), wts(sel), 'UniformOutput', false), 'Location', 'east');
title('Mean face velocity spectrum');
subplot(1,2,2); scatter(chart(:,1), chart(:,2), 8 + 40*chart(:,3), chart(:,3), 'filled'); hold on;
ww = linspace(0.15, 6, 300); plot(ww, 0.441./ww, 'k--', ww, 3*0.441./ww, 'k:', [0 6], [0.5 0.5], 'r-', [0 6], [0.441 0.441], 'b-');
ylim([0 1]); xlim([0 6.2]); grid on; xlabel('width / thickness w/t'); ylabel('resonance f t / c_L');
title('Mode chart: KLM (red line) only knows the plate'); colorbar;
What to look for: \(w/t = 0.2\) resonates at f \(t/c_L = 0.440\) (thin-bar theory 0.441); \(w/t = 1.0\) splits into 0.365, 0.435, 0.560; the first width mode sits at or just below the estimate 0.441/(w/t) (FDTD 0.275, 0.215, 0.145 against 0.294, 0.221, 0.147 at \(w/t = 1.5\), 2, 3; the largest gap is nearest the crossing with the thickness branch); at \(w/t = 6\) the strongest peak is back near 0.5. The browser demo uses the same scheme and reproduces these peaks.
m17_stack_dimensions.mTypical array stack at 3.5 and 7 MHz: plate, Desilets layers, pitch, kerf, w/t, lens R, sag, slice thickness; draws both stacks to scale · change f0 list, layer speed, H, F · runtime ≈ 1 s
%% m17_stack_dimensions.m -- typical array acoustic stack at 3.5 and 7 MHz
% Purpose : compute the dimensions drawn in Figures 17.C-17.E: active-layer
% thickness (lambda/2 plate), Desilets two-layer impedances and
% lambda/4 thicknesses, pitch, kerf, element w/t (sub-dicing check),
% cylindrical lens radius / sag / centre thickness and loss, and the
% elevation slice thickness and near-field length.
% Model : rules of thumb only (no wave solver). PZT-5H-like plate
% (7500 kg/m^3, 4560 m/s), layers c = 2500 m/s, tissue 1540 m/s,
% RTV lens 1000 m/s with 3 dB/mm at 5 MHz and alpha ~ f^1.5,
% linear array pitch = lambda. Plate speed is used for t; a slender
% element resonates lower (m17_fdtd_plate_modes.m).
% Change : f0 list, H and F per frequency, kerf, cLay, the array pitch rule.
% Expect : 3.5 MHz: t 0.651 mm, layers 8.95/2.34 MRayl at 179 um, pitch
% 0.440 mm, w/t 0.63, R 32.4 mm, sag 0.560 mm, slice 1.95 mm.
% 7 MHz: t 0.326 mm, layers 89 um, pitch 0.220 mm, w/t 0.61,
% R 10.8 mm, sag 0.237 mm, slice 0.87 mm. Runtime < 1 s.
clear; close all; rng(1);
%% Parameters
f0s = [3.5 7]; % MHz
Hs = [12 4.5]; % elevation height (mm)
Fs = [60 20]; % elevation focus (mm)
kerf = [0.030 0.020]; % dicing kerf (mm)
vP = 4560; rhoP = 7500; Zp = rhoP*vP/1e6; % plate speed (m/s), density, MRayl
ZL = 1.5; cLay = 2500; % tissue MRayl, layer speed (m/s)
ct = 1540; cl = 1000; te = 0.30; % tissue / lens speeds (m/s), lens edge (mm)
a5 = 3; nexp = 1.5; % lens loss dB/mm at 5 MHz, exponent
Z1 = Zp^(4/7)*ZL^(3/7); Z2 = Zp^(1/7)*ZL^(6/7); % Desilets two-layer
figure('Position', [80 80 900 420]);
for q = 1:numel(f0s)
f0 = f0s(q); lam = ct/(f0*1e3); % mm
t = vP/(2*f0*1e3); % plate lambda/2 (mm)
d = cLay/(4*f0*1e3); % lambda/4 layer (mm)
pitch = lam; w = pitch - kerf(q);
R = Fs(q)*(ct/cl - 1); sag = R - sqrt(R^2 - (Hs(q)/2)^2); Tc = te + sag;
alpha = a5*(f0/5)^nexp; % dB/mm at f0
slice = 0.886*lam*Fs(q)/Hs(q); nf = Hs(q)^2/(4*lam);
fprintf('--- %.1f MHz (lambda %.3f mm in tissue)\n', f0, lam);
fprintf('plate t = %.3f mm (Zp %.1f MRayl); ML1 %.2f MRayl, ML2 %.2f MRayl, each %.0f um\n', t, Zp, Z1, Z2, 1e3*d);
fprintf('pitch %.3f mm, kerf %.0f um, width %.3f mm, w/t = %.2f', pitch, 1e3*kerf(q), w, w/t);
if w/t > 0.6, fprintf(' -> sub-dice: halves at w/t = %.2f\n', (w - kerf(q))/2/t); else, fprintf(' (ok)\n'); end
fprintf('lens R = %.1f mm, sag %.3f mm, centre %.3f mm, two-way centre loss %.1f dB\n', R, sag, Tc, 2*alpha*Tc);
fprintf('elevation F %.0f mm (F# %.2f): slice %.2f mm, near field H^2/(4 lambda) = %.1f mm\n', Fs(q), Fs(q)/Hs(q), slice, nf);
% stack to scale (azimuth section, 3 elements)
subplot(1, 2, q); hold on;
z = 0; lay = {t, d, d, Tc}; col = {[0.85 0.55 0.3], [0.3 0.65 0.4], [0.6 0.85 0.65], [0.75 0.7 0.95]};
for e = 0:2
x0 = e*pitch;
patch(x0 + [0 w w 0], [0 0 t t], col{1});
patch(x0 + [0 w w 0], t + [0 0 d d], col{2});
patch(x0 + [0 w w 0], t + d + [0 0 d d], col{3});
end
patch([0 3*pitch-kerf(q) 3*pitch-kerf(q) 0], t + 2*d + [0 0 Tc Tc], col{4});
patch([0 3*pitch-kerf(q) 3*pitch-kerf(q) 0], [-0.2 -0.2 0 0], [0.8 0.8 0.8]);
axis equal; xlim([-0.05 3*pitch+0.05]); ylim([-0.25 t+2*d+Tc+0.1]);
xlabel('azimuth (mm)'); ylabel('thickness (mm)');
title(sprintf('%.1f MHz: plate, 2 layers, lens (to scale)', f0));
end
What to look for: the numbers of Figure 17.C and 17.E. At 3.5 MHz: PZT 0.651 mm, layers 8.95 / 2.34 MRayl at 179 µm, linear pitch 0.440 mm with \(w/t = 0.63\) (sub-dice), \(R = 32.4\) mm, sag 0.560 mm, slice 1.95 mm, near field 81.8 mm. At 7 MHz: PZT 0.326 mm, layers 89 µm, pitch 0.220 mm, \(w/t = 0.61\), \(R = 10.8\) mm, sag 0.237 mm, slice 0.87 mm, near field 23.0 mm.
17.14Self-check
Q1A 1-3 composite at 20 MRayl radiates into water (1.48 MRayl). What fraction of the power crosses a bare interface, and what single \(\lambda /4\) layer gives a perfect match at \(f_0\)?
\(T_I = 4 \times 20 \times 1.48/(21.48)^{2} = 0.257\) (−5.9 dB). \(Z_1 = \sqrt{20 \times 1.48} = 5.44\) MRayl.
Q2Estimate the mechanical Q of a PZT plate (34.2 MRayl) with a 3 MRayl backing radiating into tissue (1.5 MRayl), no matching layer. How many cycles to ring down 20 dB?
\(Q \approx (\pi /2)(34.2/4.5) = 11.9\). Exactly, \(r = 0.916 \times 0.839 = 0.768\), \(Q = \pi /\ln (1/0.768) = 11.9\); per cycle −2.29 dB, so about 8.7 cycles to −20 dB.
Q3Compute the Desilets two-layer impedances for PZT (34.2 MRayl) into tissue (1.5 MRayl) and the thickness of each layer at 7 MHz if the speeds are 2800 and 2300 m/s.
\(Z_1 = 34.2^{4/7}\,1.5^{3/7} = 8.95\) MRayl, \(Z_2 = 34.2^{1/7}\,1.5^{6/7} = 2.34\) MRayl. \(d_1 = 2800/(4 \times 7 \times 10^{6}) = 100\) µm, \(d_2 = 2300/(28 \times 10^{6}) = 82\) µm.
Q4A lens of RTV with \(c_l = 980\) m/s should focus at 25 mm in 1540 m/s tissue. What radius, and what sag over a 4 mm elevation aperture?
\(c_t /c_l - 1 = 0.571\), \(R = 25 \times 0.571 = 14.3\) mm. Sag = \(14.3 - \sqrt{14.3^{2} - 2^{2}} = 0.140\) mm (\(H^2 /8R = 0.140\) mm).
Q5The matching layer of a 5 MHz design came out 8 % too thick. Roughly where does its transmission peak move, and how much sensitivity do you lose at 5 MHz?
The peak moves to \(f_0 /1.08 = 4.63\) MHz. At 5 MHz the loss is small, a few tenths of a dB (0.14 dB at 5 %, 0.52 dB at 10 % for the ideal single layer). The symptom is a low centre frequency, not low sensitivity.
Q6Why does the −6 dB fractional bandwidth reward a matching design that rings?
The −6 dB band is measured between the outermost crossings; a passband with a dip that stays above −6 dB counts as one wide band. In time, a spectral ripple is a pair of echoes, so the pulse has a tail. Judge designs on both the band and the −20 dB (or −40 dB) pulse length.
Q7An element of a 5 MHz array is 0.25 mm wide and 0.40 mm thick. Is it in the safe range for lateral modes, and what thickness resonance should a plane-strain bar model predict relative to the plate value?
\(w/t = 0.63\), at the edge of the usual guideline (≈ 0.6); FDTD at \(w/t = 0.6\) still shows a single mode. Its resonance is near f \(t/c_L \approx 0.42\) (FDTD at 0.6) instead of 0.5, about 16 % below the plate value: design the plate thickness with the bar velocity, not the plate velocity.
17.15Figures from the literature
Measured data to hold the models of this module against. Images are not reproduced; open the linked figure. The two plots are re-computed here from numbers or equations printed in the papers, with the solver named.
Single λ/4 layer against a graded layer, measured pulse-echo
Z. Li, D.-Q. Yang, S.-L. Liu et al., "Broadband gradient impedance matching using an acoustic metamaterial for ultrasonic transducers," Scientific Reports 7, 42863 (2017). doi:10.1038/srep42863 (CC BY 4.0). Fig. 5.
Fig. 5(a–d) compares the echo and spectrum of a 4 MHz ceramic transducer with one conventional \(\lambda /4\) layer and with a graded (cone-structured) layer, on a polystyrene block. The paper reports a measured −6 dB bandwidth of 58.4 % for the \(\lambda /4\) layer (58.2 % simulated) and about 107 % for the graded layer. Look at the single-layer band next to this module's one-layer \(\sqrt{Z_p Z_L }\) prediction (56 % in Demo 17.3, a different ceramic, backing and electrical port): the same order, and a factor of two short of what a graded profile reaches. Going much wider than one resonant layer allows needs more layers or a graded impedance.
Tapered (graded) matching layers on a 1-3 composite
S. N. Ramadas, R. L. O'Leary, A. J. Mulholland, G. Hayward, A. Mackintosh, A. Troge, R. A. Pethrick, "Tapered transmission line technique based graded matching layers for thickness mode piezoelectric transducers," Proc. 2009 IEEE International Ultrasonics Symposium, 1695–1698. doi:10.1109/ULTSYM.2009.5441463 (author preprint: Strathprints 29323). Fig. 2 and Fig. 6.
Fig. 2 shows exponential and triangular impedance tapers that match a 16 MRayl composite to water and their modelled input reflection; Fig. 6 is the measured pulse-echo impulse response and spectrum of the resulting 1 MHz device into water with a 6 MRayl backing, about 80 % fractional bandwidth. The plot below recomputes the exponential taper of their Eq. (2) with this module's ABCD code. Look at the trade: the taper is a high-pass match that needs a total length of about half a wavelength at the lowest frequency, while a single \(\lambda /4\) layer is a narrow band-pass that repeats at odd multiples of \(f_0\).
Lens materials: speed, attenuation and impedance, measured
M. Audoin, A. Salari, B. G. Tomov, K. F. Pedersen, J. A. Jensen, E. V. Thomsen, "Diverging polymer acoustic lens design for high-resolution row-column array ultrasound transducers," IEEE Trans. Ultrason. Ferroelectr. Freq. Control 71(1), 202–213 (2024). doi:10.1109/TUFFC.2023.3327567 (accepted manuscript: DTU Orbit). Table II, Fig. 4, Fig. 6.
Table II lists sound speed, attenuation at 6 MHz and impedance for the lens candidates (RTV 1020 m/s, 1.44 dB/mm, 1.11 MRayl; PDMS 930 m/s, 3.10 dB/mm, 1.30 MRayl; and faster polymers for diverging lenses). Fig. 4 shows the measured elevation-plane transmit fields with each lens on a 6 MHz, 0.2 mm pitch row-column array, and Fig. 6 the measured transmit impulse responses in time and frequency for each lens (the no-lens reference is quoted at 6.25 MHz with 103 % −3 dB bandwidth), where the spectral filtering by lens attenuation computed in Demo 17.5 can be read off a real probe. Note the spread: this RTV sits below the 2–4.5 dB/mm assumed in this module's lens examples, so the lens loss and downshift there are on the high side for an unfilled RTV.
KLM-designed transducers at 100–300 MHz, measured
C. Fei, C. T. Chiu, X. Chen, Z. Chen, J. Ma, B. Zhu, K. K. Shung, Q. Zhou, "Ultrahigh frequency (100 MHz–300 MHz) ultrasonic transducers for optical resolution medical imagining," Scientific Reports 6, 28360 (2016). doi:10.1038/srep28360 (CC BY 4.0). Fig. 2.
Per its caption, Fig. 2 shows the measured electrical impedance (top row) and pulse-echo responses with normalised spectra (bottom row) of LiNbO3 transducers designed with a KLM-based tool, with a vapour-deposited parylene layer of about 1.5 µm as matching and protective layer. All show more than 40 % −6 dB bandwidth; the 300 MHz device measured 275 MHz against 315 MHz simulated. Look at the gap: at these frequencies a micrometre of thickness or bond error is a large fraction of a wavelength, the regime where Section 17.4 says bondlines and thickness tolerance dominate.
When λ/4 layers stop being practical: mass-spring matching
J. Chen, C. Fei, D. Lin, P. Gao, J. Zhang, Y. Quan, D. Chen, D. Li, Y. Yang, "A review of ultrahigh frequency ultrasonic transducers," Frontiers in Materials 8, 733358 (2022). doi:10.3389/fmats.2021.733358. Figs. 9, 11 and 12.
Fig. 9 is the pulse-echo response and spectrum of a 200 MHz LiNbO3 transducer; Figs. 11 and 12 sketch metal–polymer matching structures based on mass-spring and transmission-line networks, used when no material with the required impedance can be made thin and uniform enough. Read them as the thin-layer limit of the ABCD chain in Section 17.3: for kd ≪ 1 a dense layer contributes mainly its \(Z \sinh \gamma d \approx j\omega \rho d\) term (a mass), a soft layer mainly its \(\sinh (\gamma d)/Z \approx j\omega d/(\rho c^2)\) term (a compliance).
17.16Go deeper
- C. S. Desilets, J. D. Fraser, G. S. Kino, "The design of efficient broad-band piezoelectric transducers," IEEE Trans. Sonics Ultrason. SU-25(3), 115–125, 1978. The 1/3 and 4/7, 1/7 exponents.
- J. Souquet, P. Defranould, J. Desbois, "Design of low-loss wide-band ultrasonic transducers for noninvasive medical application," IEEE Trans. Sonics Ultrason. SU-26(2), 75–81, 1979.
- R. Krimholtz, D. A. Leedom, G. L. Matthaei, "New equivalent circuits for elementary piezoelectric transducers," Electronics Letters 6(13), 398–399, 1970.
- S. Sherrit, S. P. Leary, B. P. Dolgin, Y. Bar-Cohen, "Comparison of the Mason and KLM equivalent circuits for piezoelectric resonators in the thickness mode," Proc. IEEE Ultrasonics Symposium, 1999.
- G. S. Kino, Acoustic Waves: Devices, Imaging, and Analog Signal Processing, Prentice-Hall, 1987. Transducer equivalent circuits and matching.
- T. L. Szabo, Diagnostic Ultrasound Imaging: Inside Out, Academic Press (2nd ed. 2014). Transducer design, lenses, arrays.
- D. M. Pozar, Microwave Engineering, Wiley. Quarter-wave, binomial and Chebyshev transformers.
- W. A. Smith, B. A. Auld, "Modeling 1-3 composite piezoelectrics: thickness-mode oscillations," IEEE Trans. Ultrason. Ferroelectr. Freq. Control 38(1), 40–47, 1991.
- J. Virieux, "P-SV wave propagation in heterogeneous media: velocity-stress finite-difference method," Geophysics 51(4), 889–901, 1986. The staggered scheme used in Demo 17.6.
- J.-P. Berenger, "A perfectly matched layer for the absorption of electromagnetic waves," J. Comput. Phys. 114, 185–200, 1994.
- R. Lerch, "Simulation of piezoelectric devices by two- and three-dimensional finite elements," IEEE Trans. Ultrason. Ferroelectr. Freq. Control 37(3), 233–247, 1990.
- Explicit time-domain finite elements for transducers: the PZFlex papers by G. L. Wojcik and colleagues (Weidlinger Associates) in the IEEE Ultrasonics Symposium proceedings of the 1990s.