© The Author(s), Published by EDP Sciences, 2026

1 Introduction

In-plane electrostatic ultrasound transducers – a particular kind of CMUT¹ that is based on lateral motion – are microfabricated devices that emit waves by squeezing the gas in a narrow cavity between a flexible electrode (“beam”) and a static electrode (“stator”). These electrodes are carved into the depth of the chip, such that, upon electrostatic actuation, the beam is deflected in the in-plane direction whilst the gas is expelled in the out-of-plane direction, through a microslit, into the atmosphere. We have previously reported a proof of concept of this device in [1]. Figure 1 presents a scheme of said transducer, where its main components are identified: the cavity between beam and stator, the clearances that enable the motion of the beam, and the slits through which the gas is expelled. These components result in a significant damping force that determines the performance of the device near its resonance frequency. The purpose of this work is to develop analytic expressions that describe the propagation of acoustic waves in an “L-CMUT”, comparing them against finite-element simulations and experiments.

Figure 1. (A) Photograph of a fabricated L-CMUT and scheme of (B) its top view and (C) its cross section. The device consists of a beam that is pulled towards a stator by an electrostatic force. As this occurs, the gas in the cavity between beam and stator is expelled towards the slits and then departs from the chip, emitting ultrasound waves.

Given that the dimensions of an L-CMUT are in the same order of magnitude as the viscous and thermal boundary layers, the application of the Full Linearised Navier–Stokes (FLNS) equations is required; nevertheless, an approximate form of these, known as the “Low Reduced Frequency” (LRF) model [2], is applicable in this case. In contrast, the applicability of the Reynolds equation is limited to low frequencies in which this device does not operate, and the classical Helmholtz wave equation does not consider the visco-thermal effects. The condition for applying the LRF approximation is that the reduced frequency be much smaller than unity, that is, k = ω l/c ₀ ≪ 1. Taking a typical microslit aperture of l = 30 μm as the reference length, this requirement is met in the case of air (c = 343 m/s) up to frequencies of about 180 kHz, which is enough for the typical operation range of an L-CMUT (10–100 kHz). A direct consequence of the LRF assumption is that the pressure field in the direction of the reference length is constant. In this case of a Cartesian coordinate system, this enables a convenient separation of the cross-section coordinate (i.e. along the reference length) from the propagation coordinates. By applying the LRF equations to the geometry of an L-CMUT, we were able to develop an approximate analytical model for the prediction of the frequency-dependent damping force that is exerted onto the beam and the corresponding gas volume that is expelled from the device.

2 The LRF model for the L-CMUT

By scaling the state variables of the gas (pressure p, temperature T and density ρ) with respect to the values of the quiescent state (P ₀, T ₀, ρ ₀) and the velocity of the particles ($\overrightarrow{v}$) with respect to the speed of sound (c ₀), the Navier–Stokes equations can be formulated in terms of the following dimensionless variables:

• S² = ɭ²ρ0ω/μ: The shear wave number squared (a Reynolds number).

• k = ωɭ/c₀: The reduced frequency.

• σ² = μcp/κ: The Prandtl number.

• γ = c_p/c_v : The ratio of specific heats (adiabatic constant).

Here, μ denotes the dynamic viscosity and κ the thermal conductivity. Let the critical dimension be aligned to the y-axis (l = l _y ), leaving the the x- and z-axes for the propagation directions. The dimensionless coordinates are then $\widehat{x}=x/l_y$, $\widehat{y}=y/l_y$, and $\widehat{z}=z/l_y$. For small harmonic amplitudes at frequency ω, the LRF approximation results in the following set of linear partial differential equations for the time harmonic variations² p _{a c} , T _{a c} , ρ _{a c} and ${\overrightarrow{v}}_{\mathit{ac}}$ of the normalised state variables and particle velocity, respectively:

$$\begin{array}{cc}& \frac{\partial p_{\mathit{ac}}}{\partial \widehat{y}}=0\hfill \end{array}$$(1)

$$\begin{array}{cc}& j{s}^2{v}_{ac,\widehat{x}}=-\frac{s^2}{k\gamma }\frac{\partial p_{\mathit{ac}}}{\partial \widehat{x}}+\frac{{\partial }^2{v}_{ac,\widehat{x}}}{\partial {\widehat{y}}^2}\hfill \end{array}$$(2)

$$\begin{array}{cc}& j{s}^2{v}_{ac,\widehat{z}}=-\frac{s^2}{k\gamma }\frac{\partial p_{\mathit{ac}}}{\partial \widehat{z}}+\frac{{\partial }^2{v}_{ac,\widehat{z}}}{\partial {\widehat{y}}^2}\hfill \end{array}$$(3)

$$\begin{array}{cc}& \frac{\partial v_{ac,\widehat{x}}}{\partial \widehat{x}}+\frac{\partial v_{ac,\widehat{y}}}{\partial \widehat{y}}+\frac{\partial v_{ac,\widehat{z}}}{\partial \widehat{z}}+jk{\rho }_{\mathit{ac}}=0\hfill \end{array}$$(4)

$$\begin{array}{cc}& j{s}^2{\sigma }^2{T}_{\mathit{ac}}=\frac{{\partial }^2{T}_{\mathit{ac}}}{\partial {\widehat{y}}^2}+j\left(\frac{\gamma -1}{\gamma }\right)s^2{\sigma }^2{p}_{\mathit{ac}}\hfill \end{array}$$(5)

$$\begin{array}{cc}& p_{\mathit{ac}}={\rho }_{\mathit{ac}}+T_{\mathit{ac}}.\hfill \end{array}$$(6)

The two main acoustic components of the L-CMUT structure, the beam-stator cavity and the micro-slits, are depicted in Figure 2. Let us denote the boundaries of each structure in the respective axes by $\{{\widehat{x}}_{-},{\widehat{x}}_{+}\}$, $\{{\widehat{y}}_{-},{\widehat{y}}_{+}\}$ and $\{{\widehat{z}}_{-},{\widehat{z}}_{+}\}$, such that $\mathrm{\Delta }\widehat{y}={\widehat{y}}_{+}-{\widehat{y}}_{-}$ indicates its dimensionless thickness. The reference length for the case of the slit is denoted l_y = h_s, and that of the cavity is l_y = g₀. The major difference between the slit and the cavity is that the y-boundaries of the cavity are not plane but a function of x, turning the thickness, $\mathrm{\Delta }\widehat{y}(\widehat{x})$, into a function of $\widehat{x}$ as well. The ${\widehat{y}}_{-}(\widehat{x})$ curve of the cavity corresponds to the deflected shape of the beam, which can be calculated by a static force balance at the biasing point (see [3]), whereas the ${\widehat{y}}_{+}(\widehat{x})$ curve corresponds to the shape of the stator, which is freely determined by the designer. Schenk et al. [3] show why the spatial zeroth mode of a clamped-clamped Euler–Bernoulli beam (Eq. (9)) is sufficient to model very accurately the entire electrostatic bifurcation pattern of a beam under electrostatic load. This results in the general deflection profile described in equation (7), with z _DC = a ₀/g ₀ as the dimensionless static deflection. The arguments of [3] extend into the stiffness-dominated dynamic range relevant to this paper. Moreover, we argue in [1] that the stator profile ought to have the same shape as the beam deflection profile in order to avoid additional electrostatic instabilities. Combining these findings leads to a convenient expression for the thickness $\mathrm{\Delta }\widehat{y}(\widehat{x})$ that only depends on two design variables: the static deflection (z _DC) and the stator shape factor (σ _s ).

Figure 2. Illustration of the dimensions and coordinates used to describe the geometry of the slits and the beam-stator cavity of an L-CMUT.

$$\begin{array}{cc}\hfill {\widehat{y}}_{-}\left(\widehat{x}\right)& =z_{\mathrm{DC}}{\varphi }_0(g_0\widehat{x}/L)\hfill \end{array}$$(7)

$$\begin{array}{cc}\hfill {\widehat{y}}_{+}\left(\widehat{x}\right)& =1+{\sigma }_s{\varphi }_0(g_0\widehat{x}/L).\hfill \end{array}$$(8)

Here ${\varphi }_0(\xi )$ denotes the dimensionless and normalized Euler–Bernoulli zero mode (cf. [3]) which is defined on the interval $\xi \in [-1/2,+1/2]$,

$$\begin{array}{cc}\hfill {\varphi }_0(\xi )& =\frac{cos({\beta }_0/2)cosh({\beta }_0\xi )-cosh({\beta }_0/2)cos({\beta }_0\xi )}{cos({\beta }_0/2)-cosh({\beta }_0/2)},\hfill \\ \hfill {\beta }_0& \approx 4.73.\hfill \end{array}$$(9)

Due to the no-slip condition, we may assume a zero particle velocity component parallel to the silicon walls of the L-CMUT. The particle velocity component perpendicular to the static silicon walls vanishes because there is no mass transport across these boundaries. Likewise, we assume the silicon walls to have no temperature difference with the environment, and due to the relatively high heat capacity of the semiconductor material, that temperature fluctuations vanish at the silicon boundaries. By inserting the ansatz in equation (10) into equation (4) (as suggested in [2]), the differential equation for the profile shape function $C(\widehat{x},\widehat{y})$ in equation (11) emerges.

$$\begin{array}{cc}& T_{\mathit{ac}}=-\left(\frac{\gamma -1}{\gamma }\right)C(\widehat{x},\widehat{y})p_{\mathit{ac}}\hfill \end{array}$$(10)

$$\begin{array}{cc}& \frac{{\partial }^{2}C}{\partial {\widehat{y}}^2}-j{s}^2{\sigma }^{2}C(\widehat{x},\widehat{y})=j{s}^2{\sigma }^2.\hfill \end{array}$$(11)

Similarly, inserting the ansatz in equations (12) and (13) results in a nearly identical differential equation for the profile shape function of the particle velocity, $A(\widehat{x},\widehat{y})$, save for the absence of the Prandtl number (σ ²). We conclude, then, that these two functions, which are governed by differential equations with the same form and have the same boundary conditions, are identical, save for the absence of the Prandtl number in the solution of $A(\widehat{x},\widehat{y})$. This coincides with the observations of Stinson [4] and Beltman [2].

$$\begin{array}{cc}& v_{ac,\widehat{x}}=-\frac{j}{k\gamma }A(\widehat{x},\widehat{y})\frac{\partial p_{\mathit{ac}}}{\partial \widehat{x}}\hfill \end{array}$$(12)

$$\begin{array}{cc}& v_{ac,\widehat{z}}=-\frac{j}{k\gamma }A(\widehat{x},\widehat{y})\frac{\partial p_{\mathit{ac}}}{\partial \widehat{z}}\hfill \end{array}$$(13)

$$\begin{array}{cc}& \frac{{\partial }^{2}A}{\partial {\widehat{y}}^2}-j{s}^{2}A(\widehat{x},\widehat{y})=j{s}^2.\hfill \end{array}$$(14)

The solution of equation (11) for $C({\widehat{y}}_{-})=C({\widehat{y}}_{+})=0$ is the cosine function in equation (15), and the solution of equation (14) for $A({\widehat{y}}_{-})=A({\widehat{y}}_{+})=0$ is the cosine function in equation (16).

$$\begin{array}{cc}\hfill C(\widehat{x},\widehat{y})& =\frac{cos\left[\sqrt{\frac{s^2{\sigma }^2}{j}}(\widehat{y}-\frac{1}{2}({\widehat{y}}_{+}(\widehat{x})+{\widehat{y}}_{-}(\widehat{x})))\right]}{cos\left[\sqrt{\frac{s^2{\sigma }^2}{j}}\frac{1}{2}({\widehat{y}}_{+}(\widehat{x})-{\widehat{y}}_{-}(\widehat{x}))\right]}-1\hfill \end{array}$$(15)

$$\begin{array}{c}\hfill \\ \hfill A(\widehat{x},\widehat{y})& =\frac{cos\left[\sqrt{\frac{s^2}{j}}(\widehat{y}-\frac{1}{2}({\widehat{y}}_{+}(\widehat{x})+{\widehat{y}}_{-}(\widehat{x})))\right]}{cos\left[\sqrt{\frac{s^2}{j}}\frac{1}{2}({\widehat{y}}_{+}(\widehat{x})-{\widehat{y}}_{-}(\widehat{x}))\right]}-1.\hfill \end{array}$$(16)

The profile shape function for the temperature, $C(\widehat{x},\widehat{y})$, also serves to calculate the density fluctuations. If equation (10) is inserted into the ideal gas equation (Eq. (5)), we obtain the following relation:

$$\begin{array}{c}\hfill {\rho }_{\mathit{ac}}(\widehat{x},\widehat{y},\widehat{z})=[1+\left(\frac{\gamma -1}{\gamma }\right)C(\widehat{x},\widehat{y})]p_{\mathit{ac}}(\widehat{x},\widehat{z}).\end{array}$$(17)

The previous results can then be incorporated into the continuity equation (Eq. (3)) to obtain a differential equation for the pressure waves:

$$\begin{array}{cc}& \frac{j}{k\gamma }\frac{\partial }{\partial \widehat{x}}\left(A(\widehat{x},\widehat{y})\frac{\partial p_{\mathit{ac}}}{\partial \widehat{x}}\right)+\frac{j}{k\gamma }A(\widehat{x},\widehat{y})\frac{{\partial }^2{p}_{\mathit{ac}}}{\partial {\widehat{z}}^2}\hfill \\ \hfill & -jk[1+\left(\frac{\gamma -1}{\gamma }\right)C(\widehat{x},\widehat{y})]p_{\mathit{ac}}=\frac{\partial v_{ac,\widehat{y}}}{\partial \widehat{y}}·\hfill \end{array}$$(18)

We proceed by integrating along $\widehat{y}$. To this end it is helpful to define the following auxiliary functions,

$$\begin{array}{cc}\hfill B(\widehat{x})& ={\int }_{{\widehat{y}}_{-}(\widehat{x})}^{{\widehat{y}}_{+}(\widehat{x})}A(\widehat{x},\widehat{y})\mathrm{d}\widehat{y}\hfill \\ \hfill & =2\sqrt{\frac{j}{s^2}}tan\left[\frac{1}{2}\sqrt{\frac{s^2}{j}}\mathrm{\Delta }\widehat{y}\left(\widehat{x}\right)\right]-\mathrm{\Delta }\widehat{y}(\widehat{x})\hfill \\ \hfill D(\widehat{x})& ={\int }_{{\widehat{y}}_{-}(\widehat{x})}^{{\widehat{y}}_{+}(\widehat{x})}C(\widehat{x},\widehat{y})\mathrm{d}\widehat{y}\hfill \\ \hfill & =2\sqrt{\frac{j}{s^2{\sigma }^2}}tan\left[\frac{1}{2}\sqrt{\frac{s^2{\sigma }^2}{j}}\mathrm{\Delta }\widehat{y}(\widehat{x})\right]-\mathrm{\Delta }\widehat{y}(\widehat{x})\hfill \end{array}$$(19)

$$\begin{array}{cc}\hfill G(\widehat{x})& :=\frac{1}{n(\widehat{x})}={\int }_{{\widehat{y}}_{-}(\widehat{x})}^{{\widehat{y}}_{+}(\widehat{x})}[1+\left(\frac{\gamma -1}{\gamma }\right)C(\widehat{x},\widehat{y})]\mathrm{d}\widehat{y}\hfill \\ \hfill & =\mathrm{\Delta }\widehat{y}(\widehat{x})+\left(\frac{\gamma -1}{\gamma }\right)D(\widehat{x}).\hfill \end{array}$$(20)

This leads to the equation

$$\begin{array}{cc}& \frac{\partial }{\partial \widehat{x}}\left(B(\widehat{x})\frac{\partial p_{\mathit{ac}}}{\partial \widehat{x}}\right)+B(\widehat{x})\frac{{\partial }^2{p}_{\mathit{ac}}}{\partial {\widehat{z}}^2}-\frac{k^2\gamma }{n(\widehat{x})}{p}_{\mathit{ac}}(\widehat{x},\widehat{z})\hfill \\ \hfill & =-jk\gamma (v_{\widehat{y}+}(\widehat{x})-v_{\widehat{y}-}(\widehat{x})).\hfill \end{array}$$(21)

Here, $v_{\widehat{y}+}$ and $v_{\widehat{y}-}$ are the velocities with which the boundary surfaces at ${\widehat{y}}_{+}$ and ${\widehat{y}}_{-}$ oscillate in the normal direction, respectively. This term on the right hand side acts as a source term, which models in particular the air displacement caused by the motion of the actuator. Given that no flow can take place at the $\{{\widehat{x}}_{-},{\widehat{x}}_{+}\}$ boundaries, we further require ${\partial p_{\mathit{ac}}/\partial \widehat{x}|}_{\widehat{x}={\widehat{x}}_{\pm }}=0$ as a relevant boundary condition for this coordinate. A suitable solution strategy begins by expanding the pressure in terms of the Hilbert space base $\{{\psi }_m(\widehat{x})\}$

$$\begin{array}{cc}\hfill {\psi }_m(\widehat{x})& =\{\begin{array}{cc}\sqrt{\frac{{\kappa }_x}{2\pi }}\hfill & \text{if}m=0\hfill \\ \sqrt{\frac{{\kappa }_x}{\pi }}cos\left(m{\kappa }_x\widehat{x}\right)\hfill & \text{if}m>0\hfill \end{array};\hfill \\ \hfill {\kappa }_x& =\frac{2\pi }{\mathrm{\Delta }\widehat{x}}·\hfill \end{array}$$(22)

This set of functions $\{{\psi }_m(\widehat{x})\}$ is adapted to the boundary conditions, it is orthonormal and complete [5] on the interval $[{\widehat{x}}_{-},{\widehat{x}}_{+}]$, i.e. we have,

$$\begin{array}{cc}\hfill p_{\mathit{ac}}(\widehat{x},\widehat{z})& =\sum _{m=0}^{\infty }{\pi }_m(\widehat{z}){\psi }_m(\widehat{x})\hfill \end{array}$$(23)

$$\begin{array}{cc}\hfill {\pi }_m(\widehat{z})& ={\int }_{{\widehat{x}}_{-}}^{{\widehat{x}}_{+}}{p}_{\mathit{ac}}(\widehat{x},\widehat{z}){\psi }_m(\widehat{x})\mathrm{d}\widehat{x}.\hfill \end{array}$$(24)

The expansion of equation (24) leads in the general case of large beam deflections to an infinite system of coupled linear ordinary differential equations with constant coefficients for the functions ${\pi }_n(\widehat{z})$ as outlined in Section 2.3 (cf. Eq. (77)). The solution of the general case requires resorting to Fredholm theory [5] and, numerically speaking, to the spectral theory of nonnormal matrices and operators [6]. This is due to the non-elliptic nature of the Navier–Stokes setting, reflecting the fact that various kinds of thermo-viscous wave patterns are potentially created inside a heavily driven L-CMUT. While proceeding along this line is very well possible and certainly provides deep insights, we will rather focus in this paper on two special cases. In Section 2.1 we will first dwell on the fact that the L-CMUT designs discussed here are typically operated at beam deflection amplitudes that are small compared to the slit aperture. In this case the approximation $\mathrm{\Delta }\widehat{y}(\widehat{x})\approx 1$ holds and, as we will see, all of the above mentioned equations for the ${\pi }_n(\widehat{z})$ decouple. In Section 2.2 we will focus on the fact that the L-CMUT uses an in-plane actuator emitting the ultrasound in the direction perpendicular to a thin wafer stack. This implies that the dimensions of the cavity in the $\widehat{x}$-direction are very large compared to the sound propagation length (z direction). In this case, the pressure gradient in the $\widehat{x}$-direction can be neglected in equation (22), which allows for an analytical solution in closed form.

2.1 Special case I: small deflections

When the difference between the static deflection and the stator recess is small compared to the slit aperture (|z _DC − σ _s |≪1) we can approximate $\mathrm{\Delta }\widehat{y}(\widehat{x})\approx 1$. This means that the coefficient functions $B(\widehat{x})$ and $D(\widehat{x})$ in equation (22) essentially become constants:

$$\begin{array}{cc}\hfill B(\widehat{x})& \approx B=2\sqrt{\frac{j}{s^2}}tan\left[\frac{1}{2}\sqrt{\frac{s^2}{j}}\right]-1\hfill \end{array}$$(25)

$$\begin{array}{cc}\hfill D(\widehat{x})& \approx D=2\sqrt{\frac{j}{s^2{\sigma }^2}}tan\left[\frac{1}{2}\sqrt{\frac{s^2{\sigma }^2}{j}}\right]-1=\frac{1}{n}·\hfill \end{array}$$(26)

As a consequence, equation (22) is reduced to the form laid out in [2]:

$$\begin{array}{cc}\hfill \frac{{\partial }^2{p}_{\mathit{ac}}}{\partial {\widehat{x}}^2}& +\frac{{\partial }^2{p}_{\mathit{ac}}}{\partial {\widehat{z}}^2}-\frac{k^2\gamma }{\mathit{nB}}{p}_{\mathit{ac}}(\widehat{x},\widehat{z})\hfill \\ \hfill & =-j\frac{k\gamma }{B}(v_{\widehat{y}+}(\widehat{x})-v_{\widehat{y}-}(\widehat{x})).\hfill \end{array}$$(27)

We now insert the expansion of p _{a c} in terms of cosine functions (Eq. (24)), bearing in mind that the following eigenfunction property holds:

$$\begin{array}{c}\hfill \frac{{\partial }^2{\psi }_m}{\partial {\widehat{x}}^2}=-m^2{\kappa }_x^2{\psi }_m(\widehat{x}).\end{array}$$(28)

This leads to a decoupled system of equations for each mode, m, given the linear independence of the eigenfunctions.

$$\begin{array}{cc}& \frac{{\mathrm{d}}^2{\pi }_m}{\mathrm{d}{\widehat{z}}^2}-(m^2{\kappa }_x^2+\frac{k^2\gamma }{\mathit{nB}}){\pi }_m(\widehat{z})=-j\frac{k\gamma }{B}{R}_m\hfill \end{array}$$(29)

$$\begin{array}{cc}& R_m={\int }_{{\widehat{x}}_{-}}^{{\widehat{x}}_{+}}(v_{\widehat{y}+}(\widehat{x})-v_{\widehat{y}-}(\widehat{x})){\psi }_m(\widehat{x})\mathrm{d}\widehat{x}.\hfill \end{array}$$(30)

The key point is here that this infinite system of ordinary differential equations for the pressure amplitudes ${\pi }_m(\widehat{z})$ is completely decoupled, i.e. no energy is transferred between the respective spatial modes in this case. We will now apply this result to the micro-slit and to the cavity.

2.1.1 Solution for the slit

In the case of the slit, the pressure waves agitated in the cavity enter the channel through the opening at ${\widehat{z}}_{-}$ and are propagated towards the opposite end at ${\widehat{z}}_{+}$, where a reflection is likely to take place. Given that the walls are static, we set R _m  = 0 ∀m. We can solve equation (30) simply by inserting a forward- and a backward-travelling exponential with a propagation constant ±Γ _m .

$$\begin{array}{cc}\hfill {\pi }_m(\widehat{z})& =a_m^{+}{e}^{-{\mathrm{\Gamma }}_m\widehat{z}}+a_m^{-}{e}^{+{\mathrm{\Gamma }}_m\widehat{z}}\hfill \end{array}$$(31)

$$\begin{array}{cc}\hfill {\mathrm{\Gamma }}_m& =\sqrt{{\left(2\pi m\frac{h_s}{W_s}\right)}^2+\frac{k^2\gamma }{\mathit{nB}}}·\hfill \end{array}$$(32)

Notice that equation (32) is not a one-dimensional treatment of the slits but the calculation of how the ${\psi }_m(\widehat{x})$ propagation modes are modulated along the $\widehat{z}$-axis. According to equation (30) we can focus on the case where the reflected wave returns through the original propagation mode, not exciting any neighbouring modes. In this particular case, the distributed volume velocity at each port, $q(\widehat{x})$, (that is, the average z-velocity along the y-direction multiplied by the thickness) is given by the sum of the individual contributions of each propagation mode. We observe this by integrating equation (13) in $\widehat{y}$ and evaluating the pressure gradient at either the inlet (${\widehat{z}}_{-}$) or the outlet (${\widehat{z}}_{+}$) – see equations (34) and (35). Given that the pressure gradient is proportional to the pressure (by virtue of the exponential form of π _m ), a mode-specific impedance, ${\overline{z}}_m$, can be defined as the proportionality constant between the pressure and the distributed volume velocity of each propagation mode.

$$\begin{array}{cc}\hfill q_{\mathrm{in}}(\widehat{x})& =-\frac{\mathit{jB}}{k\gamma }(c_0{h}_s)\sum _{m=0}^{\infty }{\psi }_m(\widehat{x}){\frac{\mathrm{d}{\pi }_m}{\mathrm{d}\widehat{z}}|}_{\widehat{z}={\widehat{z}}_{-}}\hfill \\ \hfill & =\sum _{m=0}^{\infty }\frac{p_{\mathrm{in},m}}{{\overline{z}}_{\mathrm{in},m}}{\psi }_m(\widehat{x})\hfill \end{array}$$(33)

$$\begin{array}{cc}\hfill q_{\mathrm{out}}(\widehat{x})& =-\frac{\mathit{jB}}{k\gamma }(c_0{h}_s)\sum _{m=0}^{\infty }{\psi }_m(\widehat{x}){\frac{\mathrm{d}{\pi }_m}{\mathrm{d}\widehat{z}}|}_{\widehat{z}={\widehat{z}}_{+}}\hfill \\ \hfill & =\sum _{m=0}^{\infty }\frac{p_{\mathrm{out},m}}{{\overline{z}}_{\mathrm{out},m}}{\psi }_n(\widehat{x})\hfill \end{array}$$(34)

$$\begin{array}{cc}\hfill p_{\mathrm{in}}(\widehat{x})& =P_0\sum _{m=0}^{\infty }{\pi }_m({\widehat{z}}_{-}){\psi }_m(\widehat{x})=\sum _{m=0}^{\infty }{p}_{\mathrm{in},m}{\psi }_n(\widehat{x})\hfill \end{array}$$(35)

$$\begin{array}{cc}\hfill p_{\mathrm{out}}(\widehat{x})& =P_0\sum _{m=0}^{\infty }{\pi }_m({\widehat{z}}_{+}){\psi }_m(\widehat{x})\hfill \\ \hfill & =\sum _{m=0}^{\infty }{p}_{\mathrm{out},m}{\psi }_m(\widehat{x}).\hfill \end{array}$$(36)

Let r _m  = a _m ⁻/a _m ⁺ be the reflection coefficient for the m-th mode. This reflection coefficient can be expressed in terms of the respective impedance ${\overline{z}}_{\mathrm{out},m}$ at the outlet:

$$\begin{array}{c}\hfill r_m=\frac{{\mathrm{\Gamma }}_m{\overline{z}}_{\mathrm{out},m}-{\overline{z}}_s}{{\mathrm{\Gamma }}_m{\overline{z}}_{\mathrm{out},m}+{\overline{z}}_s}{e}^{-2{\mathrm{\Gamma }}_m{L}_s/h_s},{\overline{z}}_s=\frac{{\rho }_0{c}_{0}k}{jB{h}_s}·\end{array}$$(37)

The input impedance ${\overline{z}}_{\mathrm{in},m}$ is related to this reflection coefficient by

$$\begin{array}{c}\hfill {\overline{z}}_{\mathrm{in},m}=\frac{1}{{\mathrm{\Gamma }}_m}\left(\frac{1+r_m}{1-r_m}\right){\overline{z}}_s.\end{array}$$(38)

For the case of a sound-soft boundary condition (p _out = 0), equation (39) can be expressed in an approximate form by inserting tanh(z)/z ≈ 1 − z ²/3 at low frequencies (s ² ≪ 1), which leads to

$$\begin{array}{c}\hfill {\overline{z}}_{\mathrm{in},0}\approx 12\mu \left(\frac{L_s}{h_s^3}\right)[1-4j\omega \frac{\mu }{P_0}{\left(\frac{L_s}{h_s}\right)}^2]·\end{array}$$(39)

The real part corresponds to a well-known expression in acoustics textbooks [7], but the imaginary part obtained by this coarse approximation deviates from the corresponding low-frequency expression.

A further quantity of interest is the frequency-dependent ratio, η _slit, of the total volume flow at the outlet of the slit, Q _out, to the total volume flow at the inlet, Q _in. The flows are calculated by integrating the respective velocities over $\widehat{x}$, and noticing that only the zeroth spatial mode, ${\psi }_0(\widehat{x})$, has a non-zero integral and contributes to the total volume flow.

$$\begin{array}{cc}\hfill Q_{\mathrm{in}}& =h_s{\int }_{{\widehat{x}}_{-}}^{{\widehat{x}}_{+}}{q}_{\mathrm{in}}(\widehat{x})\mathrm{d}\widehat{x}\hfill \\ \hfill & =-\frac{\mathit{jB}}{k\gamma }\left(c_0{h}_s{W}_s\right)\sqrt{\frac{{\kappa }_x}{2\pi }}{\pi }_0^{\prime }\left({\widehat{z}}_{-}\right)\hfill \end{array}$$(40)

$$\begin{array}{cc}\hfill Q_{\mathrm{out}}& =h_s{\int }_{{\widehat{x}}_{-}}^{{\widehat{x}}_{+}}{q}_{\mathrm{out}}(\widehat{x})\mathrm{d}\widehat{x}\hfill \\ \hfill & =-\frac{\mathit{jB}}{k\gamma }(c_0{h}_s{W}_s)\sqrt{\frac{{\kappa }_x}{2\pi }}{\pi }_0^{\prime }\left({\widehat{z}}_{+}\right)\hfill \end{array}$$(41)

$$\begin{array}{cc}\hfill {\eta }_{\text{slit}}:=\frac{Q_{\mathrm{out}}}{Q_{\mathrm{in}}}& =\frac{{\pi }_0^{\prime }\left({\widehat{z}}_{+}\right)}{{\pi }_0^{\prime }\left({\widehat{z}}_{-}\right)}=\frac{r_0{e}^{{\mathrm{\Gamma }}_0{L}_s/h_s}-e^{-{\mathrm{\Gamma }}_0{L}_s/h_s}}{r_0-1}·\hfill \end{array}$$(42)

2.1.2 Solution for the cavity

The excitation in this case does not occur due to an incoming wave, but due to the moving fluid-beam interface at ${\widehat{y}}_{-}(\widehat{x})$. This boundary oscillates with a velocity amplitude _vb and has a spatial velocity profile given by the Euler–Bernoulli zero mode ${\varphi }_0(g_0\widehat{x}/L)$ shown in equation (9)

$$\begin{array}{c}\hfill v_{\widehat{y}+}(\widehat{x})=0,v_{\widehat{y}-}(\widehat{x})=\frac{v_b}{c_0}{\varphi }_0(g_0\widehat{x}/L).\end{array}$$(43)

We now proceed to compute the source terms R _m defined in equation (31), which are essentially the inner product between the Euler–Bernoulli curve, ${\varphi }_0(g_0\widehat{x}/L)$, and the respective spatial mode, ${\psi }_m(\widehat{x})$, multiplied by a constant. After some algebraic manipulation, this integral can be conveniently expressed in the following form:

$$\begin{array}{cc}\hfill {\varphi }_m& ={\int }_{{\widehat{x}}_{-}}^{{\widehat{x}}_{+}}{\varphi }_0\left(\frac{g_0}{L}\widehat{x}\right){\psi }_m(\widehat{x})\mathrm{d}\widehat{x}\hfill \end{array}$$(44)

$$\begin{array}{cc}\hfill {\varphi }_0& =\sqrt{\frac{L}{g_0}}\frac{2}{{\beta }_0}(csc\left(\frac{{\beta }_0}{2}\right)-\text{csch}\left(\frac{{\beta }_0}{2}\right))\hfill \end{array}$$(45)

$$\begin{array}{cc}\hfill {\varphi }_m& =\frac{\sqrt{2}{(-1)}^m{\beta }_0^4}{{\beta }_0^4-{(2m\pi )}^4}{\varphi }_0\text{for}m\ge 1\hfill \end{array}$$(46)

$$\begin{array}{cc}\hfill R_0& =-\frac{v_b}{c_0}{\varphi }_0\hfill \end{array}$$(47)

$$\begin{array}{cc}\hfill R_m& =\frac{{\varphi }_m}{{\varphi }_0}\text{for}m\ge 1.\hfill \end{array}$$(48)

According to equation (47) the absolute value of the ratio R _m /R ₀ rapidly decays with $O{[m]}^{-4}$. This reveals that the beam motion mostly excites the zeroth and first spatial modes, with negligible contributions from the higher modes to the generated sound pressure field.

Having calculated the effect of the oscillating ${\widehat{y}}_{-}$ boundary, we examine now the relevant boundary conditions at the bottom (${\widehat{z}}_{-}$) and the top side (${\widehat{z}}_{+}$) of the cavity. In the case of an ideal clearance, no leakage flow should take place at the bottom of the cavity, so we could expect $\partial p_{\mathit{ac}}/\partial \widehat{z}{|}_{\widehat{z}={\widehat{z}}_{-}}=0$. The typical dimensions of the clearance (h _c  = 1 μ m), the frequencies of interest (10–100 kHz), and the 1/h _c ³ behaviour of the impedance of this narrow slit are good reasons to neglect the leakage on the bottom side. On the top side we introduce a mode-specific output impedance ${\overline{z}}_{\mathrm{out},m}$, as in the case of the slit. In summary, our boundary conditions are:

$$\begin{array}{cc}& {\frac{\mathrm{d}{\pi }_m}{\mathrm{d}\widehat{z}}|}_{\widehat{z}=0}=0\hfill \end{array}$$(49)

$$\begin{array}{cc}& {\pi }_m({\widehat{z}}_{+})=-\frac{{\overline{z}}_{\mathrm{out},m}}{{\overline{z}}_c}{\frac{\mathrm{d}{\pi }_m}{\mathrm{d}\widehat{z}}|}_{\widehat{z}={\widehat{z}}_{+}},{\overline{z}}_c=\frac{{\rho }_0{c}_{0}k}{jB{g}_0}·\hfill \end{array}$$(50)

With these boundary conditions, the solution to equation (52) can be expressed in terms of hyperbolic functions as follows:

$$\begin{array}{cc}& {\pi }_m(\widehat{z})=a_{m}cosh({\mathrm{\Gamma }}_m\widehat{z})+b_{m}sinh({\mathrm{\Gamma }}_m\widehat{z})+c_m\hfill \end{array}$$(51)

$$\begin{array}{cc}& b_m=0,c_m=\frac{\mathit{jk}}{{\mathrm{\Gamma }}_m^{2}B}{R}_m\hfill \end{array}$$(52)

$$\begin{array}{cc}& \frac{c_m}{a_m}=-\frac{{\overline{z}}_{\mathrm{out},m}}{{\overline{z}}_c}{\mathrm{\Gamma }}_{m}sinh({\mathrm{\Gamma }}_{m}h/g_0)-cosh({\mathrm{\Gamma }}_{m}h/g_0)\hfill \end{array}$$(53)

$$\begin{array}{cc}& {\mathrm{\Gamma }}_m=\sqrt{{\left(2\pi m\frac{g_0}{L}\right)}^2+\frac{k^2\gamma }{\mathit{nB}}}·\hfill \end{array}$$(54)

Having solved the pressure field, we can calculate the frequency-dependent force (F _b ) exerted by the pressure field onto the beam. We are only interested in the Galerkin projection onto the Euler–Bernoulli zero mode (Eq. (9)), as this is the one that dominates the beam dynamics of the L-CMUT.

$$\begin{array}{cc}\hfill F_b& =P_0{g}_0^2\sum _{m=0}^{\infty }{\int }_{{\widehat{x}}_{-}}^{{\widehat{x}}_{+}}{\psi }_m(\widehat{x}){\varphi }_0\left(\frac{g_0\widehat{x}}{L}\right)\mathrm{d}\widehat{x}\hfill \\ \hfill & \times {\int }_{{\widehat{z}}_{-}}^{{\widehat{z}}_{+}}{\pi }_m(\widehat{z})\mathrm{d}\widehat{z}.\hfill \end{array}$$(55)

Note that for comparison with our FEM results, ${\varphi }_0\left(\frac{x}{L}\right)$ is normalised in formula equation (56) with respect to the maximum value. Using the integrals

$$\begin{array}{c}\hfill {\mathrm{\Pi }}_m={\int }_0^{h/g_0}{\pi }_m(\widehat{z})d\widehat{z}=\frac{a_m}{{\mathrm{\Gamma }}_m}sinh\left({\mathrm{\Gamma }}_m\frac{h}{g_0}\right)+c_m\frac{h}{g_0}\end{array}$$(56)

we find the relation

$$\begin{array}{c}\hfill F_b=P_0{g}_{0}L\sum _{m=0}^{\infty }{\mathrm{\Pi }}_m{\varphi }_m.\end{array}$$(57)

We can also compute the ratio of the net output to input volume velocity, which will only depend on the contribution of the zeroth spatial mode ${\psi }_0(\widehat{x})$, as was observed previously. A calculation similar to the case of the slit yields

$$\begin{array}{cc}\hfill Q_{\mathrm{out}}& =-\frac{\mathit{jB}}{k\gamma }\left(c_0{g}_{0}L\right)\sqrt{\frac{{\kappa }_x}{2\pi }}{\pi }_0^{\prime }\left({\widehat{z}}_{+}\right)\hfill \\ \hfill & =-\frac{\mathit{jB}}{k\gamma }\left(c_0{g}_0^{3/2}{L}^{1/2}\right){\mathrm{\Gamma }}_0{a}_{0}sinh\left({\mathrm{\Gamma }}_0\frac{L}{g_0}\right)·\hfill \end{array}$$(58)

The introduced volume velocity, Q _in, is directly proportional to the integral of the ϕ ₀ mode along $\widehat{x}$ (a constant that we denote by χ):

$$\begin{array}{cc}\hfill Q_{\mathrm{in}}& =g_{0}h{v}_b{\int }_{{\widehat{x}}_{-}}^{{\widehat{x}}_{+}}{\varphi }_0(g_0\widehat{x}/L)\mathrm{d}\widehat{x}\hfill \\ \hfill & =v_{b}hL\chi ,\chi \approx 0.523.\hfill \end{array}$$(59)

Hence, the following ratio of output to input volume velocities is obtained:

$$\begin{array}{cc}\hfill {\eta }_{\text{cav}}:=& \frac{Q_{\mathrm{out}}}{Q_{\mathrm{in}}}\hfill \\ \hfill =& -\frac{\mathit{jB}}{k\gamma }\left(\frac{c_0{g}_0^{3/2}}{\chi v_{b}h{L}^{1/2}}\right){\mathrm{\Gamma }}_0{a}_{0}sinh\left({\mathrm{\Gamma }}_0\frac{h}{g_0}\right)·\hfill \end{array}$$(60)

2.2 Special case II: thin actuator

In case of the L-CMUT, the dimensions of the cavity in the $\widehat{x}$-direction are usually large compared to the distance in $\widehat{z}$-direction, which the sound wave has to travel in order to leave the MEMS chip. In this case, the pressure gradient in the $\widehat{x}$-direction can be neglected as compared to the pressure gradient in $\widehat{z}$-direction. This substantially simplifies equation (22), which becomes

$$\begin{array}{cc}& (\frac{{\partial }^2}{\partial {\widehat{z}}^2}-\frac{k^2\gamma }{n(\widehat{x})B(\widehat{x})})p_{\mathit{ac}}(\widehat{x},\widehat{z})\hfill \\ \hfill & =-\frac{jk\gamma }{B(\widehat{x})}(v_{\widehat{y}+}(\widehat{x})-v_{\widehat{y}-}(\widehat{x})).\hfill \end{array}$$(61)

This simplification is only interesting for the beam-stator cavity but not for the slits, which are elongated in the $\widehat{z}$-direction. In case of the cavity, we can propose:

$$\begin{array}{cc}\hfill p_{\mathit{ac}}(\widehat{x},\widehat{z})& =a(\widehat{x})cosh[\mathrm{\Gamma }(\widehat{x})\widehat{z}]+c(\widehat{x})\hfill \end{array}$$(62)

$$\begin{array}{cc}\hfill \mathrm{\Gamma }(\widehat{x})& =\sqrt{\frac{k^2\gamma }{n(\widehat{x})B(\widehat{x})}}\hfill \end{array}$$(63)

$$\begin{array}{cc}\hfill c(\widehat{x})& =-\frac{jn(\widehat{x})}{k}\frac{v_b}{c_0}{\varphi }_0\left(\frac{g_0}{L}\widehat{x}\right).\hfill \end{array}$$(64)

The coefficient function $a(\widehat{x})$ can the be expressed in terms of the boundary value functions $p_{\mathrm{out}}(\widehat{x})=p_{\mathit{ac}}(\widehat{x},{\widehat{z}}_{+})$ and $q_{\mathrm{out}}(\widehat{x})=q(\widehat{x},{\widehat{z}}_{+})$. The result is

$$\begin{array}{cc}\hfill \frac{c(\widehat{x})}{a(\widehat{x})}& =-\frac{1}{{\overline{z}}_c}\frac{p_{\mathrm{out}}(\widehat{x})}{q_{\mathrm{out}}(\widehat{x})}\frac{B(\widehat{x})}{B}\mathrm{\Gamma }(\widehat{x})sinh\left(\mathrm{\Gamma }(\widehat{x})\frac{h}{g_0}\right)\hfill \\ \hfill & -cosh\left(\mathrm{\Gamma }(\widehat{x})\frac{h}{g_0}\right)\hfill \end{array}$$(65)

where we have used the definitions in equations (26) and (51) for B and ${\overline{z}}_c$ respectively. Notice that neglecting the pressure gradient in the $\widehat{x}$-direction in equation (22) does not mean that pressure is constant along the $\widehat{x}$-axis. From the solution of the pressure field (Eq. (63)), we can calculate the force F _b exerted on the beam and the expelled volume velocity Q _out.

$$\begin{array}{cc}\hfill F_b& =P_0{g}_0^2{\int }_{{\widehat{x}}_{-}}^{{\widehat{x}}_{+}}[\frac{a(\widehat{x})}{\mathrm{\Gamma }(\widehat{x})}sinh\left(\mathrm{\Gamma }(\widehat{x})\frac{L}{g_0}\right)+c(\widehat{x})\frac{L}{g_0}]\hfill \\ \hfill \\ \hfill & \times {\varphi }_0\left(\frac{g_0}{L}\widehat{x}\right)\mathrm{d}\widehat{x}\hfill \\ \hfill Q_{\mathrm{out}}& =-\frac{j}{k\gamma }(c_0{g}_0^2){\int }_{{\widehat{x}}_{-}}^{{\widehat{x}}_{+}}B(\widehat{x})\mathrm{\Gamma }(\widehat{x})a(\widehat{x})\hfill \\ \hfill & \times sinh\left(\mathrm{\Gamma }(\widehat{x})\frac{h}{g_0}\right)\mathrm{d}\widehat{x}.\hfill \end{array}$$(66)

In Section 3 we evaluate this integral numerically for the case of a sond-soft boundary condition at the cavity outlet, i.e. $p_{\mathrm{out}}(\widehat{x})=0$:

$$\begin{array}{cc}\hfill F_b& =g_0{\int }_{{\widehat{x}}_{-}}^{{\widehat{x}}_{+}}{w}_{\mathit{vt}}(\widehat{x})\mathrm{d}\widehat{x},\hfill \\ \hfill w_{\mathit{vt}}(\widehat{x})& =\frac{jn(\widehat{x})}{k}\left(\frac{P_0{g}_0{v}_b}{c_0}\right)\frac{h}{g_0}{\varphi }_0\left(\frac{g_0}{L}\widehat{x}\right)\hfill \\ \hfill & \times [1-\frac{tanh[\mathrm{\Gamma }(\widehat{x})\frac{h}{g_0}]}{\mathrm{\Gamma }(\widehat{x})\frac{h}{g_0}}]·\hfill \end{array}$$(67)

In the low-frequency limit, equation (69) should reproduce the predictions obtained by Melnikov et al. [8] and Wall et al. [9] by an analysis of the Reynolds equation. Since $n(\widehat{x})B(\widehat{x})$ has the expansion

$$\begin{array}{c}\hfill n(\widehat{x})B(\widehat{x})=\frac{g_0^2{\rho }_0\mathrm{\Delta }\widehat{y}{(\widehat{x})}^2}{12\mu j}\omega +O{[\omega ]}^2\end{array}$$(68)

we find

$$\begin{array}{c}\hfill \frac{tanh\left[\mathrm{\Gamma }(\widehat{x})\frac{h}{g_0}\right]}{\mathrm{\Gamma }(\widehat{x})\frac{h}{g_0}}=1-j\frac{4\mu \gamma }{c_0^2{\rho }_0}{\left(\frac{h}{g_0}\right)}^2\omega +O{[\omega ]}^2\end{array}$$(69)

which, once inserted into equation (70), results in the same expression for the distributed load obtained by Melnikov et al. [8]:

$$\begin{array}{c}\hfill w_{\mathit{vt}}(\widehat{x})=4\mu {\left[\frac{h}{g_0\mathrm{\Delta }y(\widehat{x})}\right]}^3{v}_b{\varphi }_0\left(\frac{g_0\widehat{x}}{L}\right)+O[\omega ].\end{array}$$(70)

The corresponding force projected onto the ϕ ₀ mode is obtained by performing the integral in equation (69). Wall et al. [9] analysed this case and found a very accurate algebraic approximation for the Θ(z) function that emerges from this integral.

$$\begin{array}{cc}\hfill F_b& =2\mu L{\left(\frac{h}{g_0}\right)}^3{v}_b\mathrm{\Theta }(z_{\mathrm{DC}}-{\sigma }_s)+O[\omega ]\hfill \end{array}$$(71)

$$\begin{array}{cc}\hfill \mathrm{\Theta }(z)& ={\int }_{-1/2}^{1/2}\frac{2{\varphi }_0^2(\xi )}{{[1-z{\varphi }_0(\xi )]}^3}\mathrm{d}\xi \hfill \end{array}$$(72)

$$\begin{array}{cc}& \approx -\frac{1}{162}+\frac{1}{82}{(1-z)}^{-1/2}-\frac{2}{109}{(1-z)}^{-3/2}\hfill \\ \hfill & +\frac{62}{77}{(1-z)}^{-5/2}.\hfill \end{array}$$(73)

Along this line and using the Chebyshev–Edgeworth expansion (cf. Eq. (76)) of reference [3], it is possible to compute higher order frequency corrections and to determine the frequency domain for which the Reynolds theory yields accurate results within accepted tolerance margins for F _b .

2.3 Outlook: general case

To enhance the sound pressure generated by L-CMUTs it may be attempted to design for harder driving conditions or to generate more volume flow per stroke by fabricating actuators with larger beam height h. Both strategies lead beyond the realm of the special cases studied in Sections 2.1 and 2.2. For the general case we have to replace equation (30) by the following infinite set of coupled ordinary differential equations with constant coefficients for the amplitudes ${\pi }_m(\widehat{z})$,

$$\begin{array}{cc}& \sum _{n=0}^{\infty }[-B_{m,n}{\pi }_n^{″}(\widehat{z})+(B_{m,n}^{\prime }+k^2\gamma G_{m,n}){\pi }_n(\widehat{z})]\hfill \\ \hfill & =jk\gamma R_m,\hfill \\ \hfill & B_{m,n}={\int }_{{\widehat{x}}_{-}}^{{\widehat{x}}_{+}}{\psi }_m(\widehat{x})B(\widehat{x}){\psi }_n(\widehat{x})\mathrm{d}\widehat{x},\hfill \\ \hfill & B_{m,n}^{\prime }={\int }_{{\widehat{x}}_{-}}^{{\widehat{x}}_{+}}{\psi }_m^{\prime }(\widehat{x})B(\widehat{x}){\psi }_n^{\prime }(\widehat{x})\mathrm{d}\widehat{x},\hfill \\ \hfill & G_{m,n}={\int }_{{\widehat{x}}_{-}}^{{\widehat{x}}_{+}}{\psi }_m(\widehat{x})G(\widehat{x}){\psi }_n(\widehat{x})\mathrm{d}\widehat{x}.\hfill \end{array}$$(74)

Notice that in a classical mechanical setting the “mass matrix” B _m, n and the “stiffness matrix” $B{\prime }_{m,n}+k^2\gamma G_{m,n}$ would be real-valued, positive definite matrices. This is not the case in the Navier–Stokes setting studied here, where B _m, n , $B{\prime }_{m,n}$ and G _m, n are infinite-dimensional, complex-valued objects, causing the solution theory of equation (77) to be much more involved than the mechanical case. Nevertheless, it is possible to cast equation (77) into a form only containing compact operators, i.e. into a form that allows for well-defined approximations by finite-dimensional matrices (finite rank approximation). Moreover, the theory of finite-rank approximations of compact operators [5] asserts that larger matrices give conceptually higher numerical accuracy. Implementing finite rank approximations requires in our case analysing how quickly the off-diagonal elements of, say, B _m, n vanish. Using trigonometric identities it is readily seen that these matrix elements are linear combinations of the generalized Fourier coefficients B _m as defined in,

$$\begin{array}{c}\hfill B(\widehat{x})=\sum _{m=0}^{\infty }{B}_m{\psi }_m(\widehat{x}).\end{array}$$(75)

This simple observation enables the required analysis of the off-diagonal matrix elements B _m, n . Finally, going for higher volume flows requires using the full Rayleigh-Sommerfeld integrals [10] to model the sound radiation, to compute the relevant impedance matrix, and to quantify the mode mixing caused inside the L-CMUT by radiating sound into the ambient.

2.4 Approximate model for the aggregate parameters of the complete L-CMUT structure

By inspection of Figure 1, one can identify that the L-CMUT structure is composed of one narrow cavity coupled to the bottom slit and one wide cavity coupled to the top slit. These two cavities are interconnected through the clearances above and below the beam, through which no appreciable flow is expected to take place because of their narrow opening. Nonetheless, the exact nature of the coupling that takes place as waves depart from the cavity towards the respective slit is a matter that requires a deeper analysis. There is, unfortunately, an abrupt mismatch of the cross-sections of the cavity and the slit. As a first approximation, we assume here that each propagation mode of the cavity, ψ _m, cav, is losslessly transmitted to the respective propagation mode of the slit, ψ _m, slit. This is equivalent to saying that both the pressure and volume velocity are preserved at the point where the cross-section changes, which requires compensating the value of the particle velocity abruptly. If this is the case, then the mode-specific, acoustic impedance perceived at the input of the slit (Eq. (39)) turns out to be the same mode-specific impedance imposed as a load at the top side of the cavity (Eq. (51)). This is the simplest kind of coupling between the two elements that can be calculated. In reality we may expect energy losses to take place in this transition, and perhaps also a more complex coupling between the propagation modes of slit and cavity, which would lead to deviations from this model. It is also worth noting that every slit is shared by two cavities, or, equivalently, the impedance of “half a slit” is imposed to the cavity (whereby the same pressure exercised by the full slit is offered to half of the flow). As a result, the total acoustic parameters of an L-CMUT device composed of N _b beams – the equivalent friction coefficient (c _{e q}  = F _b /v _b ) and the equivalent piston area (S _{e q}  = Q _out/v _b ) of the respective outlet – can be calculated as follows:

$$\begin{array}{cc}\hfill c_{\mathit{eq}}& =N_b(c_{\text{cav,bot}}[2{\overline{z}}_{\text{slit,bot}}]+c_{\text{cav,top}}\left[2{\overline{z}}_{\text{slit,top}}\right])\hfill \end{array}$$(76)

$$\begin{array}{cc}\hfill S_{eq,\text{bot}}& =N_b\chi hL{\eta }_{\text{cav,bot}}{\eta }_{\text{slit,bot}}\hfill \end{array}$$(77)

$$\begin{array}{cc}\hfill S_{eq,\text{top}}& =N_b\chi hL{\eta }_{\text{cav,top}}{\eta }_{\text{slit,top}}.\hfill \end{array}$$(78)

3 Finite-element model

The analytic solutions developed in the previous sections were compared against equivalent, finite-element models (FEM) that solve the Full Linearised Navier–Stokes (FLNS) equations. These models were programmed in ANSYS ©2023R1, activating the FLNS option for the FLUID220 element. Although the calculation of the pressure field did not impose stringent requirements on the mesh quality, the calculation of volume flow at the respective surface (i.e. the surface integral of the perpendicular component of the velocity) required a very fine mesh indeed. First, the proposed LRF solution of the idealised geometries (rectangular slit and curved cavity) will be compared against the equivalent ANSYS simulation, imposing a sound-soft boundary condition at the outlet (p_out = 0). Then, the approximate L-CMUT model will be compared against a two- and three-dimensional finite-element model of the complete structure.

3.1 Simulation of isolated components

Figure 3 shows a comparison of the calculation of the acoustic parameters of the rectangular slit, as obtained by both the analytic and finite-element models. Only the zeroth propagation mode was taken into consideration for the simulation. For this propagation mode, in which pressure is uniform in the x-axis, a two-dimensional model (i.e. consisting of only one element in x) suffices to capture the physics of the slit. A spatial resolution of 40 elements across h_s was required to compute the volume flow with enough accuracy. The results of the LRF equations and of the FEM simulation match perfectly. A further comparison with the approximate curve for the Reynolds regime (s² ≪ 1) is also shown in the picture to illustrate the realm in which the real part of the acoustic parameters can be approximated by this simpler expression. The imaginary part, however, ought not to be calculated with this approximation. It is noticeable that the real part of the impedance remains nearly constant for a wide frequency range, but at high enough frequencies it rises very sharply. In the same frequency range where the real part rises, the imaginary part, which would otherwise increase linearly with frequency, falls suddenly. This is due to a reflection taking place at the open end of the slit, which can be recognised by an inspection of the velocity fields. The volume transfer ratio of the slit remains very close to unity for a wide frequency range. It is only when said reflections become relevant that its phase is shifted and it magnitude rises above unity. This does not mean that the slit works more effectively at high frequencies; instead, the sharp increase in the impedance makes the overall ratio of output volume flow to input pressure lower than at low frequencies.

Figure 3. Calculation of (A) the acoustic impedance and (B) the volume transfer ratio of a rectangular slit (L s  = 400 μm, h s  = 30 μm) according to the LRF model, the FLNS finite-element model, and the simplified expressions for the Reynolds regime (s 2 ≪ 1).

The calculation of the acoustic parameters of a curved cavity are shown in Figure 4. A mesh consisting of 30 elements across g₀ and 50 divisions across L was configured for the finite-element model (simulating only one-half of the structure in the positive x-axis due to symmetry). The profiles of the distributed force calculated by the LRF method (case I) correspond nearly exactly to those obtained by the FLNS simulation. The assertion that only the zeroth and first propagation modes are relevant for the analysis of the cavity can be verified qualitatively: the profile of the distributed force resembles a cosine function shifted by a mean value. The calculation of the friction coefficient by LRF (case I) and FEM also yield the same results. Application of the approximation for a thin actuator (case II) resulted in a slight overestimation of the real part of the friction coefficient. A further simulation possibility was explored for this case, consisting of a two-dimensional finite-element model whose gap is fitted to approximate the friction coefficient of the full structure. Such a reduction from a three-dimensional to a two-dimensional geometry would only be meaningful for the case of negligible lateral flow, for which equation (74) was derived. The pressure profile calculated for this “effective slice” should correspond to the average value, when integrating across the beam’s length, so that the projected force can be calculated as F_b = w_avgχL. Likewise, the velocity with which this slice is excited should correspond to the average velocity of the beam, so that the equivalent velocity at the midpoint can be calculated as v_b = v_avg/χ. Noticing that the ratio between distributed load and wall velocity for a two-dimensional slice in the Reynolds regime is given by equation (73), we can propose an expression for the effective gap (g_eff) of such a slice (Eq. (86)). Therefore, the question is whether the simulation of this equivalent two-dimensional structure, not with the Reynolds equations but with FLNS, results in approximately the same friction coefficient as the three-dimensional FLNS simulation. According to the results in Figure 4B, the effective gap approximation is very close to the outcome of the three-dimensional model, but requiring a significantly lower computational cost. For designs where L ≫ h, this two-dimensional “effective gap” model is an interesting alternative for a quick validation. While the three-dimensional model requires tens of minutes of computation time, the two-dimensional model may just require a few minutes and the analytic model merely the fraction of a second – hence the benefit of approximate models.

Figure 4. Calculation of (A) the distributed force and (B) the effective friction coefficient of a curved cavity (g 0 = 5 μm, s 0 = 1 μm, a 0 = 1 μm, h = 75 μm, L = 900 μm) according to the LRF model, the FLNS finite-element model, and the approximate two-dimensional model with an “effective gap”.

$$\begin{array}{cc}\hfill \frac{w_{\mathrm{avg}}\chi L}{v_{\mathrm{avg}}/\chi }& =4\mu \frac{h^3}{g_{{\text{eff}}^3}}{\chi }^{2}L=\frac{F_b}{v_b}\hfill \\ \hfill & =2\mu L{\left(\frac{h}{g_0}\right)}^3\mathrm{\Theta }(z_{\mathrm{DC}}-{\sigma }_s)\hfill \end{array}$$(79)