Acoustics on small scales

We present a modelling strategy based on the finite element method to describe flexible, piezoelectric structures surrounded by a compressible fluid, including viscosity. Non-conforming interfaces based on the Mortar method are used to couple the different physical domains. Finally, we present an application example of a piezoelectrically actuated MEMS structure to illustrate the modeling procedure and the impact of viscous effects.


Introduction
In classical acoustics, one uses the wave equation to describe the physical phenomena. In its derivation, the fluid is assumed as compressible but inviscid. Thus, viscosity, which is present in any real fluid, is usually neglected. This is a reasonable simplification if the viscous boundary layer is thin compared to both the acoustic wave length and relevant feature dimensions of the problem under investigation. The thickness of the viscous boundary layer, which will develop in any real fluid, e.g. when particles move tangential to a wall, can be estimated based on the Stokes boundary layer thickness δ μ = μ/(π fρ) (see e.g. [1]). Given usual parameters (dynamic viscosity μ and density ρ of air), we obtain a boundary layer thickness of approximately 70 µm at a frequency f of 1 kHz. Thus, viscous effects become important for problems in the micro-scale like MEMS devices.
Several methods exist to model viscous effects e.g. impedance-like boundary condition models [2][3][4] and low reduced frequency models [5][6][7]. Despite being computationally reasonable, these models have a limited range of applicability. The impedance-like boundary condition model has geometric constraints, and its use is recommended for geometries that are large compared to the viscous boundary layer thickness. The low reduced frequency model is only suggested for cases where the acoustic wavelength is considerably larger than the length scale and the boundary layer thickness. In modelling the fluid-structure interactions, these models do not consider the coupling of the shear (tangential) velocities with structural mechanics. Using a full linearised Navier-Stokes (FLNS) formulation is computationally more challenging. On the other hand, the FLNS formulation meets the full fluid-structure coupling conditions and can simulate a compressible viscous fluid [7][8][9][10][11]. For typical MEMS devices such as microphones [12,13], viscosity and density sensors [14][15][16], the interaction between flexible structure, e.g., elastic cantilevers, surrounded by a compressible and viscous fluid, and the electric excitation mechanism need to be considered. A qualitative sketch of such a typical geometry is given in Fig. 1.
We propose an efficient modelling strategy based on the finite element method, using the non-conforming interface technique for coupling solid and fluid mechanic domains with flexible discretization. The modelling strategy is illustrated in an example problem.

Physical modelling
The solid mechanics problem comprises the domain of the flexible structure Ω s and the active piezoelectric material Ω p . We need to satisfy the balance of momentum where u denotes the displacement vector, ρ i the density of the solid material, and σ i the stress tensor. The index i ∈ {s, p} is used to distinguish between flexible structure and active piezoelectric material. We assume a linear strain-displacement relationship, i.e. the strain tensor is computed by and consider the material behavior in the structure domain as linear elastic but anisotropic, relating stress and strain tensor through the material stiffness tensor C by In the piezoelectric domain we use the linearised piezoelectric constitutive law (Voight piezoelectricity), i.e.
where E and D denote the electric field and flux vectors, respectively, e is the piezoelectric coupling tensor and the electric permittivity tensor. The assumption of a linear constitutive law is valid for small perturbations from a reference state, which is appropriate for the harmonic problems considered in this work. Furthermore, we use Gauss' law to describe the electric flux density by and describe the electric field E = −∇φ by the electric scalar potential φ, thus identically fulfilling Faraday's law for the electrostatic case, ∇×E = 0. The final linear piezoelectricity formulations are obtained by inserting (4) into (1) and (5). The viscous compressible fluid is modeled by the linearised balance of mass and momentum respectively. The variables p v and v denote the fluid pressure and fluid velocity vector, respectively, which are perturbation quantities from a reference state. The speed of sound c = √ K/ρ 0 is defined via the adiabatic compression modulus K, and the unperturbed fluid density ρ 0 . Note that we have assumed zero background velocity in the reference state. The viscous effects enter by assuming an isotropic Newtonian behavior for the fluid, i.e. we describe the fluid stress tensor by where λ and μ are the fluids bulk and dynamic (shear) viscosity, respectively, and I is the unit tensor. Neglecting viscous effects is usually appropriate in regions outside the viscous boundary layer. In this acoustic far-field region Ω a we assume μ = λ = 0, which allows to eliminate the velocity from (6) obtaining the acoustic wave equation for the acoustic pressure where the variable p a denotes the acoustic pressure perturbation. Finally, we need to enforce boundary and coupling conditions. On the wall boundaries of fluid and solid we require zero velocity and displacement, respectively. This corresponds to homogeneous Dirichlet conditions. At the boundary of the acoustic fluid region we assume homogeneous Neumann boundary conditions. This physically correspond to zero normal acceleration of the boundary, thereby modelling a stationary wall as well as symmetry.
The coupling conditions at the interface between viscous fluid and solid or piezoelectric domain are the dynamic and kinematic condition respectively, where n is the outer normal vector of the fluid domain and reverse direction of outer normal vector of solid or piezoelectric domains (n = n v = −n i ). These transmission conditions ensure the traction and velocity continuity at fluid-structure interfaces. Similarly, the coupling conditions at the fluid-fluid interface between viscous and acoustic formulation are formulated as The outer normal n is defined as n = n v = −n a , with the outer normal of acoustic n a .

Finite element formulation
To obtain the weak form for the finite element solution, we multiply the governing PDEs by appropriate test functions, denoted by () in the following, and integrate them over the whole computational domain.
In the solid mechanic domain Ω s we consider the balance of momentum (1) and in the piezoelectric domain Ω p we add Gauss' law (5). After integration by parts and insertion of the constitutive relations (3) and (4) as well as the boundary conditions one obtains where the surface terms at the interfaces have been retained to illustrate the surface coupling. In the viscous fluid region we consider balance of mass and momentum (6), where the latter was integrated by parts yielding Similarly we have incorporated boundary conditions but retained the interface terms. In the acoustic fluid region we obtain the weak form of the wave equation (8) Ωa where we have ready incorporated the homogeneous Neumann boundary conditions, but retained the interface term.
Coupling of the domains is done via non-conforming interfaces using the Mortar method, e.g. [17,18]. At the interface between viscous fluid and acoustic region Γ va we can directly insert the interface conditions (10). For the interfaces between structure or piezoelectric domain and viscous fluid, Γ sv and Γ pv , respectively, we use the Lagrange multiplier method. We introduce an additional unknown t for the traction at the interface, i.e. we have Inserting the interface conditions into (11), (14) and (15) yields Continuity of the velocities, i.e. the kinematic condition in (9) is enforced in a weak sense at the interface, i.e. we require which also provides the additional equations necessary for the determination of the introduced unknowns at the interface. The final weak form is given by equations (17), (12), (18), (13), (19) and (20) which form the basis for the finite element implementation, which was done in the open source finite element framework openCFS [19]. One can use standard nodal finite elements with ansatz functions of variable order. In order to fulfill the inf-sup condition in the viscous fluid region, it is recommended to reduce the order for the pressure degree of freedom, e.g., use quadratic ansatz functions for velocity and linear ansatz functions for the pressure degrees of freedom. The evaluation of the interface integrals is done using an intersection mesh between the two non-conforming interfaces.

Example problem
As an application example for the described modelling strategy, we consider a cantilever structure. For the sake of simplicity it is treated as a 2D model, which is accurate for structures that are considerably wider than long. The model is analyzed in the frequency domain, i.e. we compute the steady state response to periodic forcing.

Model description
The cantilever consists of a silicon substrate and a piezoelectric layer, which is placed between two electrodes. The structure is surrounded by air. As discussed in Sect. 1, the effect of viscosity is especially significant in viscous boundary layer regions. Therefore, the fluid close to the structure (viscous fluid region) is modelled using the fluid viscous formulation. Further away from the speaker it is possible to  neglect the effect of viscosity. In this region (acoustic region), the air is represented by the acoustic wave equation. An overview of the model and some of the most important dimensions can be found in Fig. 2 and Table 1. The used material properties are given in Table 2.
The piezoelectric cantilever is excited with a harmonic voltage signal of 1 V across the electrodes. The harmonic solution was calculated for excitation frequencies from 200 Hz to 200 kHz. Within this frequency range 500 frequency steps (logarithmic sampling) were calculated. Apart from the excitation, the cantilever is also fixed at the nodes along the left boundary.
The velocity at the walls (left boundary below cantilever and floor) is set to zero in the viscous fluid region. At the left boundary above the cantilever, the normal velocity is zero to satisfy the symmetry condition.
In the acoustic region, we apply homogeneous Neumann boundary conditions for the symmetry (left boundary) and walls (floor). There should be no reflection at the right and the top boundary as the device is considered in an open domain. This is accomplished by adding a perfectly matched layer (PML) region, in which the acoustic waves are damped (i.e. absorbed).
The model was calculated with the quadrilateral mesh shown in Fig. 3. When meshing the viscous fluid region, it is important to consider the boundary layer thickness δ μ . Within this region the velocity gradients are especially high. In order to resolve these gradients, a boundary layer mesh was added. This boundary layer mesh has the thickness of the Stokes boundary layer at the maximum excitation  frequency of 200 kHz. It consists of 5 element layers with a mesh bias, which gradually increases the thickness of the elements in the normal direction of the mechanics-fluid interface (see Fig. 3). During tests with different meshing schemes, it was found that the results are also sensitive to the tangential resolution within the boundary layer. Thus, we aim to use approximately square cells (as shown in the top right detail of Fig. 3) and larger cells outside the boundary layer region (see top center detail of Fig. 3). Using this meshing strategy made it possible to reduce the computational effort without sacrificing accuracy in the result (see Fig. 4 and Table 3). For the acoustic region, the mesh size was determined using the wavelength λ min at the maximum excitation frequency. In particular, a mesh size of λ min /10 was applied across the whole region.
The different meshing schemes result in inconsistent mesh sizes at the boundaries between the various domains (mechanic, viscous fluid and acoustic regions). Instead of creating a transition mesh, we use non-conforming interfaces in order to reduce the complexity and distortion of the mesh (see Fig. 3).

Results
Parameter studies were conducted to investigate the influence of the gap thickness below the cantilever, and the effect of the dynamic viscosity of the fluid. The computed response for the displacement of point P (see Fig. 2) at the end of the cantilever is shown in Fig. 5. When comparing the results for different gap thicknesses (see Table 3. Computational effort for the locally and globally refined meshes in Fig. 4 5a) it is apparent that the peaks in the amplitude diagrams are less prominent for smaller gap heights. This effect is frequency dependent. At low frequencies the peaks are shifted to the left and more attenuated. As the frequency increases, the curves converge. For example, the first peak is located at ≈ 8 kHz for h g = 40 µm. For h g = 20 µm, and h g = 10 µm, the first peaks have turned into saddle points, at ≈ 1 kHz and < 200 Hz At the last peak the difference between the transfer functions is less extreme. These results suggest that there is a damping effect which is dependent on the thickness of the boundary layer relative to the gap height. By varying the dynamic viscosity, it is possible to directly affect the ratio between boundary layer thickness and structure dimension, since the thickness of the viscous boundary layer δ μ is square root proportional to dynamic viscosity μ. In Fig. 5b, a similar effect as in Fig. 5a can be observed. At low frequencies the peaks are strongly dampened. As the frequency increases and the boundary layer de- creases in thickness, the damping effect is reduced. Furthermore, the peaks at lower viscosities are less prominent, which supports our findings from the study of the gap thickness. We can clearly see that damping effects increase with relative boundary thickness compared to gap height. Figure 6 shows fluid velocity in the viscous fluid region at two characteristic frequencies: 1000 Hz and 69000 Hz (first saddle point and peak at h g = 20 µm, μ = μ air ). At both frequencies the flow underneath the cantilever is almost unidirectional. The magnitude of the fluid velocity in this region is much greater compared to the rest of the domain. As the air escapes the area underneath the cantilever it spreads out and the velocity magnitude decreases. Close to the wall the viscous boundary layer is visualised. At the low frequency the thickness of this region is greater compared to the result at the higher frequency.

Mesh
When looking at the acoustic pressure in the far field in Fig. 7 it is possible to observe spherical wave fronts, typical for point sources. Comparing the distance between the wave fronts, it appears that the wavelength is smaller for the higher frequency, as expected. It is also interesting to note that there is a high pressure region underneath the cantilever. Underneath the cantilever the boundary layers impede the flow of the air. As a result the fluid is compressed and a high pressure region is formed.

Conclusion
The presented modelling strategy allows to describe viscous effects in acoustic phenomena. By using a linearised Navier-Stokes type formulation in the viscous fluid region, one can gain detailed insight into pressure and velocity fields in the arising viscous boundary layers. Together with the ability to couple the viscous fluid formulation to solid mechanic and piezoelectric formulations, one can model typical MEMS devices. The presented application example illustrated the pronounced damping effect of the viscous fluid on the response of the piezoelectric actuator. Selected parameter studies clearly demonstrated that the damping effect is more pronounced if the ratio between the viscous boundary layer thickness and a characteristic structure dimension is larger. Several strategies have been suggested to minimize the computational effort arising in the viscous fluid formulation: Selective refinement of the mesh only in boundary layer regions is very effective compared to a uniform mesh size in the viscous fluid domain. Furthermore, the viscous fluid domain should be limited to regions where viscous effects need to be considered. In the acoustic far field, one can neglect viscous effects and use the acoustic wave equation, thereby reducing the computational effort significantly. Finally, the presented non-conforming interface formulation allows for independent meshes on different domains, which greatly simplifies the meshing procedure in many cases.