Mixed finite elements applied to acoustic wave problems in compressible viscous fluids under piezoelectric actuation

In the present contribution, we develop a mixed finite element method capable of the coupled multi-field simulation of a viscous fluid actuated by a piezoelectric resonator. Several challenges are met with in this setting, among which are the necessity of correct interface coupling, near incompressibility of the fluid, adverse geometric dimensions of flat piezoelectric transducers and different length scales of shear and pressure wave. Assuming small deformations and velocities, we present a mixed variational formulation with consistent interface coupling conditions in (mechanic) frequency domain. Both fluid and piezoelectric solid domain are discretized using Tangential-Displacement Normal-Normal-Stress elements. These elements model not only the deformation, but add an independent tensor-valued stress approximation. The method has been rigorously proven to be free from shear locking for flat prismatic or hexahedral elements. Thus, modeling of the flat geometry of piezoelectric resonators as well as resolution of the fastly decaying shear wave are facilitated. To circumvent the problem of volume locking due to the near incompressibility of the fluid, an additional independent pressure field is introduced. We present computational results indicating the capability of the method.

the excitation of attenuated shear waves in the fluid, which can be used to sense the viscosity of the liquid as it affects the resonance frequency as well as the Q-factor of the resonance. The effect can be evaluated using a simple 1D-model [16], which is justified by the fact that the diameter of the disk is comparatively large to the characteristic length dimensions involved with the viscous entrainment, particularly the penetration depth elements proposed for the two domains. Hybridization and other implementational issues are discussed. Finally, we present computational results for a viscous fluid under actuation from a circular piezoelectric transducer in Sect. 4.

Theory
The focus of this work is to develop a mixed finite element formulation for viscous fluids and piezoelectric solids in the acoustic regime, such that subsequent coupling of the two physical domains through their common interface is facilitated. In the sequel, the partial differential equations governing the (acoustic) flow in viscous fluids as well as equations for piezoelectric solids in the linear regime are summarized. Furthermore, we devise interface conditions for their mutual coupling.
In the sequel, we assume Ω ⊂ R 3 to be the domain of interest, which contains in part the piezoelectric solid, in part the viscous fluid. We denote this spatial splitting by Ω = Ω p ∪ Ω v with Ω p ∩ Ω v = ∅, using Ω p for the piezoelectric solid and Ω v for the domain of the viscous fluid. We define their common interface by Γ i = ∂Ω p ∩ ∂Ω v .

Balance equations for viscous fluids
We are interested in the velocity field v in a general compressible and viscous fluid in the domain Ω v with boundaries Γ v = ∂Ω v . The conservation of mass in fluids is given by with ρ denoting the density of mass. The conservation of momentum reads with the body force per unit mass k and the Cauchy stress tensor T. As we focus on acoustic effects only, it is justified to assume the velocity v to be a small perturbation around a vanishing hydrodynamic velocity. As a consequence, we identify spatial and material points, which enables the definition of a displacement vector u. Velocity and displacement are related via derivation by time, v = ∂u ∂t .
Under this assumption of small velocities, the convective terms in (2) can be neglected, which leads to the linearized mechanical balance equation with acting forces per unit volume f = ρk. According to Filippi et al. [6], we introduce the static pressure p as well as the viscosity stress tensor σ and rewrite the Cauchy stress tensor with I denoting the unit tensor. The constitutive equations for Newtonian fluids connect the strain rate tensor D v = 1 2 (∇v + ∇v T ) with the viscous stress tensor with μ the (first) viscosity coefficient and λ B the second viscosity coefficient. In analogy to the constitutive equations for linear elastic solids, we implicitly defined the tensor of viscosity C v . Often, the constitutive equations are cited using the bulk viscosity μ B = λ B + 2/3μ as second material parameter, Note, that as tr(σ ) = 0 the viscosity stress tensor is not introduced as a deviatoric tensor, which furthermore implies tr(T) = −3 p . The constitutive equation for the pressure is introduced via the speed of sound c [32,39], where the subscript ( ) s denotes adiabatic changes of state. The compressibility modulus K s is defined via the Newton-Laplace relation In acoustic flows, pressure and density may be assumed as infinitesimal perturbations p e and ρ e of constant, static values of state p 0 and ρ 0 , In the examples, we have T · n = 0 on some parts of the boundary. This necessarily implies p 0 = 0. Inserting these assumptions into (8) leads in combination with (9) to the relation We substitute the assumption of small-density perturbations (10) and small velocities in the equation of mass balance (1) and find Integration by time on both sides together with (12) leads to the constitutive equation for the pressure At this point, we mention that a linear relation between dilatation div u and pressure p may also be assumed a priority for barometric fluids as it is done, e.g., in [42]. Summing up, the constitutive equations for the Cauchy stress tensor now read

Balance equations for piezoelectric solids
Below, we summarize balance and constitutive equations for piezoelectric solids. For a detailed introduction, we refer the interested reader to the monograph [40]. Let Ω p denote the domain associated with the piezoelectric solid. In addition to the displacement field u, we are interested in computing the electric field E. We assume to stay in the linear regime of small deformations and electric fields, such that we may use Voigt's theory of linear piezoelasticity [38]. If each part of Ω p is simply connected, Faraday's law in the electrostatic regime ensures that the electric field is a gradient field, connected to the electric potential φ via E = −∇φ. Its work conjugate is the dielectric displacement D. Moreover, due to the assumption of small deformations, it is justified to use the linearized strain tensor ε given by as work conjugate of the total stress tensor T. Voigt's theory of linear piezoelasticity postulates a linear relationship between mechanical stresses T and electric field E as well as strain ε and dielectric displacement D, which is usually cited in the form Above, S E denotes the mechanical compliance tensor measured under constant electric field, T is the dielectric permittivity tensor at constant stress, and d is the piezoelectric tensor which describes the electromechanical coupling.
In correspondence with our notations in the fluid domain, we denote external body loads by f. Concerning electric loads, we assume that the interior of the non-conducting piezoelectric body is devoid of free charges. Under these assumptions, the mechanical balance equation in solids and Gauss' law of zero charges read in differential form Together, (17) and (18) yield a system of partial differential equations to be satisfied within the piezoelectric domain Ω p .

Harmonic excitation
We consider time-harmonic processes with excitation signals of general form, Here,ḡ depends only on the spatial coordinate x, ω is the angular frequency, and j = √ −1 is used to denote the imaginary unit. A linear partial differential equation thus allows for a time-harmonic displacement ansatz Having simplicity of notation in mind, we will omit this exact notation in the following, writing, e.g., u instead ofū for the time-harmonic displacement field. The velocity is given by v = jωu. The strain rate tensor D v is connected to the linearized strain tensor in the same manner, as D v =ε = jωε. We are interested in computing the displacement, respectively, velocity field in the fluid as well as the piezoelectric solid part. The electric potential is considered as unknown field in the piezoelectric solid Ω p only. Thereby we neglect all effects of the electric field within the viscous fluid or the surrounding air. Moreover, within the fluid domain Ω v pressure perturbations p e are of interest. Given the static density ρ 0 , density perturbations ρ e can then be evaluated by (12). To simplify notation, we will omit all indices and use p = p e to denote the pressure perturbation and ρ = ρ 0 for the static density.
Inserting the time-harmonic ansatz (20) into (4) and (15) the balance and constitutive equation in the fluid A more general derivation of a formulation for harmonically excited viscoelastic materials can be found in [8, p.188ff]. In the sequel, the pressure p and the viscous stresses σ will be used as independent variables in the numerical discretization.
For the piezoelectric solid domain Ω p , the following system of partial differential equations for unknown displacement u, electric potential φ, and total stresses T is obtained,

Boundary and interface conditions
In addition to the partial differential equations (21), boundary conditions on Γ = ∂Ω as well as interface conditions on the common interface Γ i of piezoelectric solid and fluid have to be defined. For the sake of immediate readability, we introduce only the most basic cases below, such as zero displacement conditions or external forces. More advanced settings including absorbing boundary conditions are discussed in Sect. 4.1.2. We first consider the boundary of the viscous fluid Γ ∩ ∂Ω v . To model zero surface velocities v = 0 on the boundary part Γ v,u , the displacement u is prescribed as zero. Then, this condition coincides with the standard condition on clamped solid surfaces Γ p,u ⊂ Γ ∩ ∂Ω, such that Surface tractions and pressures t 0 on piezoelectric as well as fluid domain boundaries are prescribed via Note that this includes free boundaries by setting t 0 = 0. Last, on the common interface Γ i , continuity of displacement, respectively, velocity and surface stresses has to be assumed, Here we used n v and n p as the outer unit normals of Ω v and Ω p , respectively. On the common interface, we have n v = −n p . On the boundary of piezoelectric domains, ∂Ω p , additionally electrical boundary conditions are given. On a non-vanishing part of the boundary Γ φ , the electric potential is prescribed, This condition models an electrode with given applied potential. On the remaining boundary, Γ D = ∂Ω p \Γ φ , surface charges ρ e 0 are given This includes the case of insulated boundaries when setting ρ e 0 = 0 and is used as a boundary condition toward the surrounding air.

Discretization with mixed finite elements
In this Section, we briefly introduce the Tangential-displacement normal-normal-stress (TDNNS) finite element method for linear elasticity, to follow up with a formulation for nearly incompressible materials introducing an independent variable for the pressure. Finally, the extension of TDNNS method to linear piezoelectric materials is summarized. For more detailed descriptions and analyses, we refer to the works of Pechstein et al. [25,26,28] for the linear elastic case. Piezoelectric problems are treated in [17,27].
As a prerequisite, we introduce some notation on normal and tangential components of vector and tensor fields. On a (boundary) surface with (outer) unit normal vector n, a vector field a can be split into a normal component a n = a · n and a tangential component a t = a − a n n. Tensor fields B such as the stress field have a well-defined surface vector b n on any (boundary) surface, which is obtained by multiplication with the normal vector, b n = B · n. This surface vector can again be decomposed into its normal component b nn = (n · B) · n and its tangential component b nt = b n − b nn n.
Let T = {T } be a finite element mesh subdividing the domain of interest Ω into finite elements of tetrahedral, prismatic or hexahedral form. Furthermore, the finite element mesh T = T v ∪ T p shall respect the subdivision into a viscous fluid and piezoelectric solid part. The outer unit normal on element or domain boundaries shall be denoted by n without any index if the domain is clear from the context.

The TDNNS method for linear elastic solids
The TDNNS finite element method for linear elastic solids is a mixed finite element method with independent approximations for displacements u and stresses T. It is based on a generalized form of Reissner's principle (see (29)), but uses the tangential displacement u t and the normal component of the normal stress vector t nn as degrees of freedom. Accordingly these quantities are continuous across element interfaces. Note, that the normal displacement u n is discontinuous, and gaps may arise between the elements.
Tangentially continuous elements were first introduced by Nédélec [19,20] in the context of finite element methods for Maxwell's equations. Within the TDNNS method, it is proposed to use these elements originally designed for the vector potential of the magnetic flux density to represent the displacement. For the stresses, tensor-valued elements with normal-normal continuity were presented in [17,25,26]. To follow the deductions below, a detailed knowledge of implementational issues surrounding these choices is not mandatory, we refer the interested reader to the contributions cited above and also [28] for further stability results.
The mixed finite element method is based on the following Reissner-type variational formulation: find u and T piecewise polynomial on the finite element mesh T satisfying continuity conditions: u t and t nn are continuous at element interfaces, (28) and the essential boundary conditions u t = 0 on Γ u and t nn = t 0,n on Γ T , respectively, such that for all admissible virtual displacements δu and virtual stresses δT. Note that, in the above context, admissible virtual displacements and stresses satisfy the respective continuity conditions (28) and also the homogeneous essential boundary conditions, i.e., δu t = 0 on Γ u and δt nn = 0 on Γ T . For the exact choices of polynomial ansatz functions, we refer to [17,25,26] for the various element types. As a consequence of the discontinuous displacement field, the strain only exists in distributional sense. While the elastic work pair is usually defined in the sense of the Lebesgue integral Ω ε(u) : T dΩ for weakly differential displacements, the duality product ε(u), T has to be considered. This duality product is well-defined in infinite-dimensional Sobolev spaces as well as for finite element functions [28]. For piecewise smooth finite element functions, T and u satisfying (28) the duality product can be evaluated by The equivalence of (30) and (31) can be shown by integration by parts on each element. For more detailed explanations as well as for inhomogeneous boundary conditions, we refer to [25,26,28]. In [22,33], a hybridization technique is proposed and analyzed for the geometrically linear and nonlinear case, respectively. It has been shown that the proposed hybridization improves the condition of the assembled stiffness matrix, and (for linear elastic problems after static condensation of the stresses) leads to a symmetric positive definitive system matrix. In this hybridization approach, the normal-normal continuity of the stress tensor is broken and imposed by Lagrangian multipliers defined on the element interfaces. These Lagrangian multipliers resemble the normal displacement u n as is detailed below. If the Lagrangian multipliers are chosen of the same polynomial order as the normal-normal component of the stress tensor, the original and the hybridized system are equivalent. As the stresses are then purely local not explicitly satisfying the inter-element continuity constraint (28), these unknowns can be eliminated element by element in a static condensation procedure.
The Lagrangian multiplier λ is a vector valued finite element function defined uniquely on the element (inter-)faces pointing into the normal direction (λ × n = 0). The according finite element space can be implemented by using a facet space equipped with a normal vector [4,21,22]. The hybridized variational problem reads: find u, T, and λ, with u and λ satisfying the essential boundary conditions u t = 0 and λ n = 0 on Γ u such that for all admissible virtual displacements δu, admissible δλ, and all δT. The continuity constraint (28) on δT is dropped in this formulation. Also all stress boundary conditions on Γ T are now natural boundary conditions appearing as virtual work of external forces. Note that only the normal component of λ n = λ · n enters into the variational formulation.
We provide a short justification of the above variational Eq. (32). Considering the virtual work of the stresses δT, after reordering the element-wise integrals of the duality product ε(u), δT given in (30), we obtain We thereby recover the constitutive equation T : S = ε in weak form. Normal displacement u n and the Lagrangian multiplier component λ n are identified.
Reordering and summing up the remaining terms in the variational formulation (32), this time using the second line (31) to express the virtual work of elastic forces, leads to Thus, equilibrium is recovered in weak sense. The (continuous) virtual tangential displacement δu t acts as a Lagrangian multiplier for continuity or boundary values of t nt , while δλ n ensures continuity or boundary conditions of t nn exactly.
The discussed hybridization technique will be crucial for the following deductions. Finally we want to mention that Pian's method of assumed stresses [29] can be seen as a similar approach as the TDNNS method in combination with hybridization.

TDNNS for nearly incompressible fluids
In this Section, a finite element formulation for harmonically excited nearly incompressible fluids is derived, where the pressure and the viscosity stress tensor are treated as independent variables. In order to get a variational formulation, the set of Eq. (21) for viscous fluids is rewritten in the form with S v = (C v ) −1 . In analogy to the discussed hybridization technique, the normal continuity of the total stress vector t nn = σ nn − p and the boundary condition t nn = σ nn − p = t 0,n on Γ T are imposed by Lagrangian multipliers, Multiplying the constitutive equation for the pressure (36) by a piecewise smooth (polynomial) but discontinuous virtual pressure and integrating over Ω v yields Element-wise integration by parts and replacing u n by λ n leads to Next, we pay attention to the balance Eq. (35). Multiplication by an admissible virtual displacement δu satisfying (28) and using the algebraic identity div T = − grad p + div σ as well as the duality product identity (30)- (31) gives Last, the constitutive equation for the viscosity stress tensor (37) is multiplied by a virtual viscosity stress δσ which is piecewise smooth (polynomial) but not satisfying any continuity assumptions or boundary conditions. In analogy to (33), the following variational equation is obtained: Summing up Eqs. (38), (40), (41) and (42) and reordering terms, we arrive at the following structurally symmetric variational problem: find u, σ , p, and λ satisfying the essential boundary conditions u t = 0 and λ n = 0 on Γ u , respectively, such that for all admissible virtual quantities δu, δσ , δp, and δλ. By admissible, we mean that δu is tangentially continuous as in (28), δλ is unique and pointing in normal direction on each element interface, and δu t = 0, δλ n = 0 on Γ u . The virtual stress quantities δσ and δp do not satisfy any continuity or boundary conditions.

TDNNS for linear piezoelectric materials
In this Section, we briefly summarize the derivation of a mixed variational formulation for piezoelectric materials, which is shown in detail in [17,27]. Multiplying Eq. (22), with a virtual displacement δu, Eq. (22) with a virtual stress δT and Eq. (22) with a virtual electric potential δφ, and integrating over Ω p , a variational formulation is obtained. Standard continuous elements (e.g., nodal elements) are chosen for the electric potential. After stress hybridization, the following problem is obtained: find u, T, λ, and φ satisfying the essential boundary conditions u t = 0, λ n = 0 on Γ u and φ = φ 0 on Γ φ such that for all admissible virtual quantities δu, δλ, δT, and δφ. Here, we mean by admissible that δu satisfies the tangential continuity assumption (28), δλ is unique on element interfaces pointing in normal direction, on the clamped boundary δu t = 0, δλ n = 0 on Γ u , and the δφ is continuous with δφ = 0 on Γ φ . There are no boundary or inter-element continuity constraints assumed on δT.

Examples
This Section is dedicated to the presentation of computational results. We present results from two different computations concerning the problem of viscosity measurements by piezoelectric actuation. To this end, a piezoelectric circular AT-quartz resonator is applied to the surface of the fluid of interest and actuated in shear mode. The resonance frequency of these quartz actuators is slightly reduced due to liquid loading. The decline of the resonance frequency is a measure for the viscosity coefficient. In our numerical results, we aim at computing the frequency response of the resonator when applied to fluids of different viscosity. Within the first example, a simplified setting is considered, where the quartz resonator is modeled via the boundary condition. For this case, analytic solutions obtained by Beigelbeck and Jakoby [3] are available. We compare our results to their solution. In the second example, the piezoelectric resonator is modeled by a finite element discretization, and the two computational domains (piezoelectric solid and viscous fluid) are coupled as described above. The frequency response for fluids of different viscosity is computed. All computations are carried out in the framework of the open-source software package Netgen/NGSolve [1]. This package provides all required finite elements. By a python interface, different finite elements can be linked, and variational equations can be entered symbolically.

Comparison with analytic results
In this first example, the fluid domain Ω v ∞ = {(x, y, z) ∈ R 3 : z < 0} is assumed the infinite half-space. To model the behavior of a piezoelectric resonator with circular electrodes of radius R, the tangential displacement is prescribed on the fluid's top surface. The displacement in x-direction is assumed to be Gaussian, while there is no displacement in y-direction. Due to the non-uniform displacement at the fluid boundary both shear and pressure waves are excited.
Beigelbeck and Jakoby [3] presented an analytical solution for the according two-dimensional problem, solving for the amplitudes of the acoustic pressure wave and the shear wave. While the pressure wave, and accordingly the deformation component u z travel the fluid nearly undamped, the shear wave is highly damped. The penetration depth κ of the shear deformation is related to the viscosity, the angular frequency (ω = 2π f ), and the density via Along the z-axis (x = 0, y = 0. z ≤ 0), the shear displacement is analytically given by In [3], a penetration depth of κ = 0.23 µm is reported in water at a frequency of f = 6 MHz. We aim at reproducing these results in a 3D finite element computation.

Geometric setup
A sketch of the geometric setup is shown in Fig. 1. At the upper boundary of the fluid (z = 0), the displacement is prescribed as pointing in x-direction with its absolute value following the two-dimensional Gaussian, The displacement along the x-axis (y = 0) is visualized in Fig. 1. The radius of the resonator is set to R = 3 mm, which is a typical value for AT-quartz patches. For our computations u 0 = 1 µm is used. Following the reference, in the sequel all results will be presented dimensionless.
In contrast to the reference solution, the finite element computations are carried out in a bounded domain Ω v . We choose Ω v = {(x, y, z) : x 2 + y 2 < r ∞ , −h < z < 0}, setting the width r ∞ = 7.5 mm and height h = 2.471 mm = 10λ corresponding to ten wavelengths of the excited pressure wave. The following material parameters of water are taken from [3]: density ρ = 998.2 kg m −3 , speed of sound c = 1483 ms −1 ; and dynamic viscosity μ = 1 mPas. The compressibility modulus is evaluated via (9) K s = 2195 Nmm −2 . The symmetric nature of the problem allows for reduction of the geometry by half, and accordingly a reduction of degrees of freedom. On the plane of symmetry y = 0, the symmetric boundary condition u y = 0, t x = t z = 0 is imposed. The vertical outer face of the cylinder (r = r ∞ ) as well as the bottom surface are assumed to satisfy an absorbing boundary condition for p with vanishing viscous stress σ · n = 0, see al Sect. 4.1.2. We use a layered mesh of prismatic elements for the finite element analysis. To resolve the onset of the displacement wave and the damped behavior of the shear displacement near the actuated surface (z = 0), three layers of thickness 0.2 µm, 0.4 µm, and 0.6 µm are introduced at the top of the geometry. The remaining fluid body is discretized by 30 equal-sized slabs of prismatic elements. Thus, we use approximately three elements per pressure wave length. On a half-circle with radius R located at the significant support of the actuating boundary condition, the in-plane mesh size is set to 1.2 mm, whereas toward the boundaries the mesh size increases up to 3 mm. The mesh is provided in Fig. 2. Note that for the used mesh parameters, the length-to-height ratio of the elements is up to 1000:1 near the top surface. The final mesh consists of 2046 prismatic elements. We use finite elements of polynomial order k = 2 on the described finite element mesh. For this choice, we count 122,535 coupling degrees of freedom. In our problem, we assume the viscous stress to vanish far from the upper boundary. Substituting the displacement into the Eq. (48.2) and using the relation (9), the wave equation for the (acoustic) pressure is found For this equation, absorbing boundary conditions are given by [5,30] jω with ∂ p ∂n = (∇ p) · n. Using the relation ∇ p = ρω 2 u from (48.1,2) and the Lagrangian multiplier λ resembling the normal displacement, we reformulate This absorbing boundary condition (50) can be embedded in the standard way into our variational formulation (43). It results in adding the following boundary integral to the left-hand side of (43): Γ abc jωcρλ n δλ n dΓ.

Results
The results obtained for the finite element problem are presented in Figs. 3, 4, 5, 6, 7. Pressure p as well as shear deformation u x and transverse displacement u z are evaluated at the plane of symmetry (y = 0). In Figs. 3, 4, 5, near-field results for z ≥ −1 µm are shown: recall that at the boundary (z = 0) the shear deformation u x is prescribed as a Gaussian with maximum value u 0 . As shown in Fig. 3, it decreases rapidly and nearly vanishes within a distance of 1 µm. These results correspond to the reported penetration depth of 0.23 µm.
The transverse displacement u z is-according to the boundary condition-zero at the top surface (z = 0), but shows a steep gradient toward the nearest local maximum, see Fig. 4. The near field of the pressure p is virtually constant, compare Fig. 5. Far field results for u z and p are shown on a different length scale for z ≥ −3λ with λ = c f = 0.247 mm the wave length. For the pressure wave, a decay length of 2.37 m is reported in [3]. For the given geometry, this results in nearly undamped waves, as it can be seen in Fig. 6 for the pressure and Fig. 7 for the transversal displacement u z . While here the plots are shown for a distance of three wavelengths (3λ), the calculation is carried out on a thickness of the fluid of ten wavelengths (10λ). We observe that, as expected, maximal     The next example is an extension of the previous one. A circular piezoelectric AT-quartz patch resonator is submerged in viscous fluid. The patch radius is r ∞ = 7.5 mm, and two circular electrodes of radius R = 3 mm are positioned at the patch center. The patch thickness is assumed as h p = 0.2 mm. The patch is clamped around its circumference. Due to the smallness of the electrodes in comparison with the patch radius, it is sufficient to model the fluid within the (infinite) cylinder of radius r ∞ . Moreover, we use the symmetric nature of the problem, such that only half of the patch in thickness direction and the fluid below the patch are modeled. Thereby, we introduce the symmetric boundary Γ sym . Last, the infinite cylinder is truncated at h = 1 mm, and an absorbing boundary condition is prescribed on the according artificial boundary Γ abs . A schematic view is provided in Fig. 9. We evaluate the frequency response for fluids of different viscosity. The material data taken from [40] are summarized in Table 1. An electric potential difference of 2V is applied to the electrodes. In the simulation, we assume that the potential is gauged such that it vanishes on the plane of symmetry, i.e., φ = 0 on Γ sym . The remaining, nonelectroded boundaries of the patch actuator are assumed free of charges (D · n = 0). The actuator is assumed to be clamped around its circumference. For a shear actuator, the symmetry condition implies that u = 0 on Γ sym as well. Additionally, we assume that there are no displacements at radius r ∞ for the fluid, while absorbing boundary conditions are imposed on Γ abs as described in Sect. 4.1.2.
We discretize the whole computational domain by a layered finite element mesh stemming from an unstructured triangular in-plane mesh. The maximum in-plane mesh size within the radius of actuation is chosen as 1.2 mm and grows up to 3 mm toward the fluid domain boundary. With respect to the thickness direction, the piezoelectric patch is discretized using two equal-sized layers of prismatic elements. Two layers of elements are introduced in the uppermost 1.5 µm of the fluid domain, whereas the remaining domain is divided into 20 slabs of equal thickness. This results in a finite element mesh consisting of 2162 elements. We use second-order displacement, stress and pressure elements, as well as third-order elements for the electric potential. In total, we count 130,633 coupling degrees of freedom.
In Table 2, we provide the eigenfrequencies of the unloaded, clamped actuator. Open-circuit and closedcircuit simulations are compared. We observe several additional eigenfrequencies near the intended shear mode at f 1 . Next, we carry out frequency sweeps for the resonator both when unloaded and submerged in fluids of different viscosity. As results, the average of the absolute value of the shear displacement at the bottom electrode and the admittance of the piezoelectric patch are evaluated in a frequency range around the resonance frequency of the unloaded actuator. The calculations are carried out for four different values of viscosity (0.05 mPas, 0.5 mPas, 1 mPas, 3 mPas). The density as well as the speed of sound are kept constant, using the values from our first example (ρ = 998.2 kgm −3 , c = 1483 ms −1 ). The average displacement u el x on the electrode is evaluated via with A el = R 2 π denoting the area of the electrode. The admittance Y is evaluated via the electric current and the applied potential difference φ = 2V, The results are shown in Fig. 10 for the displacement u el x and in Fig. 11 for the admittance Y . Both curves show a decline of the resonance frequency for higher values of the viscosity coefficient and higher damping. The change of the resonance frequency is used in sensing applications as a measure for the viscosity coefficient.

Conclusions
We presented a mixed finite element method suitable for the coupled discretization of nearly incompressible viscous fluids and piezoelectric actuators. Due to geometric dimensions of the piezoelectric patch, as well as the rapidly decaying nature of the shear waves, a discretization using elements with high aspect ratio is advisable in solving the problem. TDNNS elements have been known to be free from shear locking on flat prismatic or hexahedral elements. To make elements of high aspect ratio robust when approaching the incompressible limit, we proposed an enrichment by additional degrees of freedom for the pressure. A variational formulation for the hybridized system including independent viscosity stress tensor and pressure was derived. The performance of the presented method was shown in two examples, including the frequency response of a liquid-loaded AT-quartz resonator.