Higher-Order Models for Resonant Viscosity and Mass-Density Sensors

Advanced fluid models relating viscosity and density to resonance frequency and quality factor of vibrating structures immersed in fluids are presented. The numerous established models which are ultimately all based on the same approximation are refined, such that the measurement range for viscosity can be extended. Based on the simple case of a vibrating cylinder and dimensional analysis, general models for arbitrary order of approximation are derived. Furthermore, methods for model parameter calibration and the inversion of the models to determine viscosity and/or density from measured resonance parameters are shown. One of the two presented fluid models is a viscosity-only model, where the parameters of it can be calibrated without knowledge of the fluid density. The models are demonstrated for a tuning fork-based commercial instrument, where maximum deviations between measured and reference viscosities of approximately ±0.5% in the viscosity range from 1.3 to 243 mPas could be achieved. It is demonstrated that these results show a clear improvement over the existing models.


Introduction
When electromechanical transducers are immersed in a fluid, the mechanical properties of the fluid influence the resonance spectrum measured at the electrical terminals [1]. Based on the theory of Sauerbrey [2], quartz thickness shear mode resonators have long been used to measure smallest mass depositions on these sensors, e.g., in vapor deposition chambers under close to vacuum conditions. An equivalent to the Sauerbrey equation for fluids has been developed by Kanazawa and Gordon [3] and found wide application in measurement practice following the work of Martin [4]. In the subsequent years, the sensitivity of acoustic fluid sensors was enhanced by replacing these bulk acoustic wave devices by surface acoustic wave sensors [5]. The drawback of using such sensors featuring dominant shearing motion, is that that according to the one-dimensional model, only the density-viscosity product can be measured, but not the individual parameters. The ability to separate density ρ and viscosity η due to spurious effects (e.g., finite area of vibration and boundary effects) is very limited for shear sensors and yields large measurement errors on the individual parameters [6]. Although this problem could be solved by adding liquid traps (see e.g., [7,8]), such sensors are not widely used in measurement practice, likely due to demanding cleaning requirements.
Piezoelectric cantilevers [9][10][11][12][13] and quartz tuning fork sensors [14][15][16][17] are excellent candidates for simultaneous determination of density and viscosity, when the resonance characteristics of the immersed sensors are evaluated. Contrary to shear wave sensors, where the relationship between the resonance parameters and ρ, η is simple, the problem is profound e.g., for vibrating cantilevers in general. Even when the cantilever is assumed with infinite length, such that the problem can be reduced to two dimensions, closed form solutions are only available for cylindrical cross-sections [18,19].
A numerically determined correction factor, boundary element methods, and tabulated values for the hydrodynamic loading for rectangular cross-sections are found in [11,20,21], respectively. However, the numerical modeling approach for vibrating rectangular cross-sections [21] is already involved, and closed form expressions are only available in some limiting cases, such as for zero viscosity or for infinitely extended plates.To our best knowledge, finding closed form solutions for arbitrary geometries is a hopeless endeavor, but the significance of such forms is questionable in measurement practice, in particular since fabrication and material tolerances are encountered. In order to achieve measurement accuracies which are on par with lab grade instruments, additional model parameter adjustments using test liquids are always required. It is therefore more suitable to work with parametrized models which approximate the physical problem well and can be adjusted by calibration measurements. In the derivation of such parametrized models, the hydrodynamic function [11], which essentially describes the fluid resistance acting against a vibrating body, is central. This fluid force constitutes a boundary condition for the resonator changing its resonance frequency and Q-factor. It is therefore possible to establish a relationship between the physical fluid properties and the measured resonance parameters. The task of the sensor designer is to guarantee that the resonance parameters can be correctly measured and that a fluid model is available which can be transformed to yield unique and accurate fluid parameters from the resonance parameters. As the exact hydrodynamic function is only known for some special cases, the conventionally used fluid models derived, e.g., by [22][23][24][25][26] are all based on an approximation of the hydrodynamic function of a vibrating prismatic beam featuring a cylindrical cross-section by means of a truncated power series. A remarkable fact, shown by a many researchers in experimental work [22][23][24][25][26][27][28][29][30][31], is that the resulting reduced order models with adjustable parameters are suitable for representing various different types of fluid loaded vibrating structures, including cantilevers of different cross-sections, large steel and miniaturized quartz tuning forks, spiral springs, platelets, U-shaped wires, piezo buzzers in fluids, etc. However, the limitations of the approximation being made to obtain models of manageable complexity, narrows the accurately measurable range of fluid parameters. The unified fluid models presented in this work, follow in a straightforward manner from the work of [11,22,27]. It encompasses the previous ones but can be extended to arbitrary order of approximation. The key features of these models are that any order of approximation can be achieved, given enough calibration measurements. Furthermore, the complexity of the fluid parameter determination from resonance parameters remains manageable for any order. As a side effect, also a viscosity-only model is derived which needs no information about the fluid density for calibrating the model parameters.
The remainder of the paper is structured as follows: In Section 2, the theoretical background is established. The two new arbitrary order fluid models in Equations (16) and (24) are derived in Sections 2.4 and 2.5. Numerical experiments for the vibrating cylinder and analyzes of measured data obtained using a piezoelectric tuning fork sensor are shown in Section 3. Two fundamental, but longer treatises on dimensional analysis and resonance parameters are found in the Appendixes A and B.

Materials and Methods
The closed form solution for the vibrating cylinder is chosen as a vehicle to reveal functional dependencies which apply also for more complex cross-sections, but for which no closed forms are known.

Vibrating Cylinder Immersed in Viscous Fluid
The complex valued fluid force amplitude per length F acting against the oscillatory motion of an infinitely long prismatic cylinder is given by [19,23] with The fluid is considered incompressible and Newtonian. The radius of the cylinder, the velocity amplitude of the rigid cross-section, and the angular frequency are denoted by R, v, and ω, respectively. The relation between time domain and complex amplitude is given, e.g., for the vibration velocity by v(t) = R{v exp(jωt)}. K 0 and K 1 denote the modified Bessel functions of second kind and j is the imaginary unit √ −1. The function Γ assumes the value of 1 for vanishing viscosity (i.e., β → ∞). The associated flow field in this case is that of the potential flow [32], where the mass/length m is periodically displaced, but no viscous losses occur. For non-zero viscosity, the real part of Γ encompasses how viscous effects add additional mass drag and the imaginary part represents viscous losses. The function is therefore also termed hydrodynamic function. We demand that the vibration amplitudes are small enough that the relation between force and velocity amplitudes remain in the linear range. This means that the convective part in the Navier-Stokes equations causing non-linearity and fluid instability is neglected. Although not generally permissible, this assumption holds for many vibrators, especially for QCMs [33], QTFs and micro cantilevers [27]. It is then convenient to eliminate the driving velocity and to consider the fluid impedance per unit length henceforth. It is essential to note that Γ and β are both dimensionless variables, where the latter is termed Reynolds number [11], or non-dimensional frequency [18]. Furthermore, it is a coincidence that the equivalent moved mass per length m equals the area of the cylinder cross-section times the fluid density. This is generally not the case for other geometries. The derivation of Equation (2) for this very basic geometry is already involved with details shown, for example, in [19]. The determination of the fluid parameters from the fluid loaded resonator is difficult, as ρ and η appear also in the arguments of two Bessel functions. The fluid models presented in this work, are based on an expansion the hydrodynamic function in Equation (2) into a power series around β → ∞, yielding For the reduced order models [22][23][24][25][26], the approximation of the fluid force in Equation (7) is added to the equation of motion of the resonator (see Section 2.3), yielding models for the resonance parameters of the fluid loaded resonator shown in Equation (12). A consequence of the agreement between Equation (7) and Equation (2) being best for large β is that the model approximation is more accurate for thicker cylinders vibrating at higher frequencies. By introducing the characteristic decay length of plane shear waves δ it also becomes clear that the decay length should be small compared to the dimension of the vibrating cylinder. Consequently, the reduced fluid models are only applicable up to a certain maximum viscosity for a given resonator.

Arbitrary Cross-Sections Vibrating in Fluid
In general, hydrodynamic functions for different cross-sections must be determined numerically. For the numerical modeling of complex geometries the characteristic variables are not easily recognized in the modeling approach. By means of a dimensional analysis, with details given in the Appendix A, however, it can be shown that the dependence of the hydrodynamic function on the fluid properties is in all cases determined exclusively by a single Reynolds number β. The additional characteristic parameters required to describe the problem completely can be expressed by additional aspect ratios α i as will be illustrated in the following. For the simple cylinder versus a rectangular cantilever vibrating in a tube (see Figure 1), for instance, follows cylinder: cantilever in tube: As the additional aspect ratios α 1 and α 2 are fixed for a given sensor geometry, the hydrodynamic function can be considered to be depending only on β for any cross-section and configuration. The method to derive a reduced order model from the simple cylindrical geometry can therefore be generally applied, but it must be kept in mind that the shape of the hydrodynamic functions is unknown and may be more complex, requiring higher orders of approximation. An indicator for requiring higher orders is, if any involved dimension is comparably small to the decay length δ given in Equation (8).

Resonator Model
The fluid forces F represent a boundary condition for the immersed resonator changing its resonance parameters. These resonance parameters are defined as characteristics of the mechanical admittance spectrum and the dependence of these parameters on the hydrodynamic fluid loading is outlined in the following.
In the vicinity of a resonance, a vibrating structure can be approximated by an equivalent spring-mass-damper system with the respective parametersk,m, andc and the additional fluid mass loading m f and damping c f . The acoustic admittances Y = v/F of such a system is given by The Nyquist plots of such admittances resemble circles as is shown in Figure 2. The resonance frequency is defined at the point where the function is real and maximum. The bandwidth ∆ω used to calculate the Q-factor by Q = ω r /∆ω lies between the ±45 • phase angle points. For driven resonances, it is therefore favorable to use the admittance Y. (Contrary to when force-displacement relations of driven resonators, or free vibrations, e.g., measured in ring-down mode, are considered.) The resonance parameters are given by A relation between resonance parameters and fluid parameters can be obtained by splitting Equation (7) in real and imaginary part and attributing them to added fluid mass m f and damping c f (see e.g., Appendix B for details). Some of the established models, which are all based on similar derivations and that can in principle be converted into each other (In case of free instead of forced oscillations, expressions for the resonance parameters have to be adapted.) are: Heinisch et al.: Youssry et al. : Toledo et al., Dufour et al.: Zhang et al.: Here, Q r , ω r = 2π f r and ω 0 denote Q-factor and angular resonance frequency in the fluid and in vacuum, respectively. g 1 and g 2 are functions of the resonance parameters which correspond to the real and imaginary parts of the hydrodynamic function. It is apparent that the number of parameters of the models m 0k , m ρk , m ηρk , c 0k , c ηk , c ηρk for [23], C 1 , C 2 , C 3 , C 4 for [22,25], M L , a 1 , a 2 , b 1 , b 2 for [24] or A, B, C for [26] differ. For instance, the model from Heinisch et al [23] where the vacuum resonance parameters ω 0 and Q 0 are fitted in the calibration procedure (m 0k , c 0k ), are considered known in the model of Zhang et al. [26]. Furthermore, c ηρk and m ηρk are allowed to differ, while according to Equation (7), the same 1/ β dependency in real and imaginary part of the hydrodynamic function suggests that c ηρk and m ηρk should be equal as in the model of Zhang et al. [26] (parameter A). Consequently, the model from Heinisch et al., yields lower deviations between the fluid reference values and the measured values, as is also outlined in [26]. For the above mentioned models, the closed form solutions for the fluid force of vibrating cylinders have been approximated by a power series in the Reynolds number, where a low order of approximation was chosen in order to obtain models which can be rearranged yielding simple expressions for η and ρ [29].
Nyquist plots of the mechanical admittance of three resonances.

Higher-Order Fluid Models
The extension to higher-order approximation is motivated by using the expression for the resonance parameters derived in Appendix B, where ρ s represents the resonator mass-density These relations can also be found in [11], for instance. However, in this work, alternative scales are used which render the equations more convenient for our approach. Γ R and Γ I represent the real and the negative imaginary parts of the hydrodynamic function i.e., Γ = Γ R − jΓ I . Following the usual derivation of the fluid model for the cylindrical cross-section, but without the early truncation of Equation (6), yields From Equation (8) follows the proportionality 1/ β ∝ √ ν/ω r , with ν denoting the kinematic viscosity ν = η/ρ. It can be expected that for vibrating cross-sections of alternative shape, the constants of the polynomial are different. By introducing the new variable ξ = √ ν/ω r , attributing ρ s to the constants (a i and b i ) and neglecting higher orders than N a and N b , the general model follows Although g a and g b are equal to Γ R and Γ I , i.e., a distinction is made, to emphasize that g a and g b are functions of the resonance parameters that are measured. It is assumed that the constants a i and b i are known from a model calibration procedure and g a and g b are calculated from the resonance parameters provided by a resonance estimation algorithm.
ξ and density ρ are therefore the two unknowns to be determined. The density ρ can be eliminated from Equation (16) yielding a polynomial in ξ alone Additionally, Q r = Q 0 at ξ = 0, and therefore b 0 = 0. The constants of the polynomial are all known such that an estimate for the kinematic viscosity is obtained byν = ξ * 2 ω r , where ξ * denotes the correct real root of Equation (18). Subsequently, an estimate of densityρ is obtained by substituting ξ * in Equation (16) and dynamic viscosity follows fromη =νρ.
A comparison with the models in Equation (12) shows that they all can be reduced to respectively.

Model Calibration
Given M calibration measurements, Equation (16) can be transformed into matrix form The ranges of viscosity and density of the test fluids define the valid calibration range. A least square fit for the parameter vector a and b can be obtained using the Moore-Penrose inverse (where W is an identity matrix)â However, the least squares optimization favors larger relative deviations in the lower viscosity range, so that W should be adjusted in order to achieve the wanted distribution of the unavoidable model deviations for M > N a + 1 and M > N b . With a sufficiently high number of calibration measurements available, the order of the approximation can be expanded as desired.

Viscosity-Only Model
For the calibration of the previous model, the densities and viscosities of various test liquids where required. However, it is possible to perform model calibration and measurement of the kinematic viscosity without knowledge of fluid density. Rearranging Equation (13) for the fluid loss factor g = Γ I /Γ R yields a model which depends only on the kinematic viscosity ν but not on density ρ The rational function of the fluid loss factor can itself be approximated by a polynomial in ξ.
The unknown ξ can be calculated for given constants c i and a measurement point g by finding the correct root ξ * of Equation (24). Subsequently, the kinematic viscosity results withν = ξ * 2 ω r . The dynamic viscosity η, however, cannot be determined using this particular model.
The constants are again obtained by a weighted least squares fit It is noteworthy that although the resonance parameters depend on both density and viscosity, a fluid model can be found which depends solely on the kinematic viscosity ν. The question of whether dynamic (η) or kinematic viscosity should be regarded as the primary viscosity parameter can be decided on the basis of this model in favor for the kinematic viscosity ν.
In this section, it was shown that extending the usual fluid models to higher polynomial order is straightforward. Given Equations (13) and (A17), it is apparent that also the extended models rely on power series approximations of the real and imaginary parts of the hydrodynamic function. Although polynomials in ξ were used for the presented models, it should be kept in mind that any function suitably approximating the hydrodynamic function could be used instead.

Results and Discussion
Synthetic data, generated using the closed form solution for the cylinder, as well as experimental data obtained with quartz tuning fork sensors, are analyzed. For the latter, the hydrodynamic function is estimated from measured data and compared to theoretical predictions. The deviations between the reference values of certified standard fluids and estimates using different model orders, are analyzed in detail. It will be shown that with increasing order, the systematic model deviations reduce to a point where measurement errors and presumably deviations between actual and nominal fluid parameters or unnoticed fluid degradations become dominant.

Numerical Analysis for the Vibrating Cylinder
Equations (2) and (13) were used to generate resonance parameters ω r and Q r for a resonator with vacuum resonance frequency 32.768 kHz (i.e., ω 0 = 205.89 · 10 3 s −1 ) and unloaded Q-factor of Q 0 = 10 4 . The fluid properties correspond to that of a N140 fluid standard at various temperatures, as shown in Figure 3. The resonator radius and density are defined with R = 0.1 mm and ρ s = 2800 kg/m 3 , respectively. The correct hydrodynamic function expression from Equation (2) (solid line in the left Figure 4) was used to generate the resonance parameters ω r and Q r for the fluid parameters of N140 in the temperature range from 10 • C to 100 • C. The dashed lines in Figure 4a represent the power series approximation of different order in Equations (14) and (15) ( , , •, * ). In (b), the deviations between reference and calculated viscosities and densities are shown. It is apparent that the first order power series model is accurate only around β → ∞ which corresponds to δ/R → 0 (see Equation (8)). The approximation ( ) [0,1], [1,2] is of the same order as the power series model ( ), but the parameters a 0 , a 1 , b 1 and b 2 were adjusted such that the relative deviations are minimized on average, yielding deviations which are approximately 5 times smaller at the higher viscosities. The deviations due to model approximation can be made arbitrarily small by increasing the order of the approximation, which is a confirmation for the validity of the approach.  [34]. The numbers of the standards coarsely agree with the kinematic viscosity in cSt at 40 • C, i.e., ν of N14 is approx. 14 cSt at 40 • C.

Quartz Tuning Fork Measurements
The measurement results were obtained with the VDC100 measurement cell and the MFA200 resonance analyzer, both provided by MicroResonant [35]. Details on the instrument and the dimensions of the QTF can be found e.g., in [36,37]. The parameters of the reference fluids and the measured resonance parameters are shown in Table 1. The standard fluids were measured at different temperatures to enhance the use of the fluids.

Extended Viscosity and Density Model
The set of four fluids # 2, 13, 19, and 23 in Table 1 were used for calibrating the parameters a 0 . . . a 3 and b 1 . . . b 4 of the model g a = ρ a 0 + a 1ξ + a 2ξ 2 + a 3ξ Scalingξ = ξ/max(ξ) was introduced so that the constants have the same physical units and the orders of magnitude are more similar. Using the estimated vacuum resonance parameters ω 0 = 205.818 s −1 and Q 0 = 14100, the determined constants are (in m 3 /kg) a 0 = 2.9983 × 10 −4 , a 1 = 2.2803 × 10 −4 , a 2 = 5.1036 × 10 −6 , a 3 = 6.0255 × 10 −8 , Figure 5 shows the deviations between reference and measured values of dynamic viscosity and density for different orders of the polynomials. The order which corresponds to the models in Equation (12) ., a 0 , a 1 , b 1 , b 2 ) showing deviations in the range of ±1.5%.Order extension reduces the deviations gradually, such that deviations of −0.57% to 0.22% are achieved over the full viscosity range from 1.27 mPas to 242.9 mPas using the parameters in Equation (28). Adding a further point to the calibration (#22) and increasing the orders by one does not enhance the agreement significantly. Interestingly, the deviations on density reduce only slightly with increasing order. In Figure 6 the various models in Equation (12) are compared with an extended model presented in this work. The models were calibrated using the same fluid as before with no weighting (i.e., W is an identity matrix, see Equation (23)) applied. The deviations are differently distributed for the various models, but which is mostly coincidentally. Although all models yield reasonable results, it apparent that the models using more parameters are superior for minimizing deviations.

Viscosity-Only Model
The function g is calculated from the resonance parameters in Table 1 by and is plotted over ξ in the left Figure 7. The fact that all measured fluids lie closely on a smooth line confirms the suitability of this model. The filled dots in Figure 7a indicate the fluids selected for model calibration. In (b), the deviations between calculatedξ and that of the reference values from Table 1 are shown for a third order polynomial. The deviations vanish therefore for these selected fluids.
The blue results represent the deviations when all fluids were used to calibrate the four parameters of the polynomial. It can be observed that the agreement is not improved drastically, which is an indicator that these deviations are not model errors, but measurement errors and alterations of the test fluids from the certificate values. Figure 8 shows the consequence of reducing the order. The dashed lines on the left are the measured values g and the solid lines are the models for polynomial orders ranging from 1 to 5. (All fluids were used for these parameter fits.) The intersections of dashed and solid lines represent ξ * . Although the number of roots of the used polynomials is equal to the order of the polynomial, it is observed that the intersections are unique, i.e., there is only one real root in the considered range, also for higher orders. In Figure 8b, the deviations between reference and estimated values of the kinematic viscosity are shown. Orders higher than four do not lead to further improvements in this case.  The relative deviations between the reference (ν) and estimated (ν) kinematic viscosities decrease with increased order, but do not improve noticeably above an order of 4. The cut-off extreme deviations for first order are +10% and −22%.

Measurement of the Hydrodynamic Function of the QTF
As was mentioned in Section 2, the approximation to the hydrodynamic function may be any suitable function and not necessarily a polynomial as used in this work. If for a given resonator a sufficient amount of measurements in various test fluids is available, the shape of the sampled hydrodynamic function can be obtained from resonance parameters by evaluating i.e., Equation (A17) derived in Appendix B. The results are plotted in Figure 9 versus the Reynolds number β defined as The density of the quartz material is ρ s = 2649 kg/m 3 and the width W and thickness H of one prong of the QTF are 610 µm and 350 µm, respectively (i.e., α = 1.74). The result is compared the a 2D finite element simulation which was performed for a cantilever with the same cross-section using the solid mechanics module of COMSOL 5.3a ( The fluid properties were implemented by equivalent complex elastic moduli [33,38].). The ( * ) represent the tabulated values from Brumley et al. [21] for an aspect ratio of α = 2, where for the alternative definitions of the Reynolds number and the hydrodynamic function is accounted for by Γ = απΓ * Brumley /4 and β = 4β Brumley , with * denoting the complex conjugate. The deviations between 2D FEM simulation of the cantilever and the measurement results for the QTF may be attributed to the interaction and finite length of the QTF prongs. Furthermore, the not considered contribution of the fluid interaction at the end faces of the QTF is particularly strong, as vibration amplitudes are largest at the ends [39]. In addition, the resonator model was considered purely mechanical, but the neglected electrical interaction also affects the resonance parameters [40]. It can therefore be concluded that in most cases using the theoretical hydrodynamic functions without model calibration will not allow satisfactory accuracies for the determination of the fluid properties. Figure 9. Estimated hydrodynamic function using Equation (A17) ( • ) and comparison with 2D simulation of rectangular cross-section (-·-) and tabulated values from [21] for aspect ratio α = 2, which is the closet match to the actual α = 1.74.

Conclusions
Fluid models providing a higher degree of approximation of the hydrodynamic function than the conventionally used ones were derived and applied to measurement results obtained with a quartz tuning fork sensor setup. These models demonstrate a clear advantage in terms of reduced model deviations for this particular problem, but it is assumed that they can be applied to many other types of vibrating sensors as well. Increased orders, however, do not significantly increase the complexity of the model calibration method and the viscosity and density calculation from measured resonance parameters.

Conflicts of Interest:
The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

Appendix A. Dimensional Analysis
Dimensional analysis allows insights into the nature of a problem solely on basis of the involved variables and their physical units [41]. For this, all involved physical properties are collected and their scaling with length, time, and mass are noted. (A1) The square bracket notation, e.g., [R] = L reads as "the dimension of R is L". In this case L represents a dimensionless scaling factor that accounts for a change of units e.g., from meters to inch. By taking the logarithm of the dimensions of the physical quantities, a matrix H can be established According to the Buckingham Π theorem [41], there are two dimensionless numbers in our case, forming the so-called π-group. The dimension of two is a consequence of having five parameters involved (the dependent variables) and only three basic physical scales (length, time and mass), which represent the independent variables. One, but not unique, representation of the π-group is obtained by calculating the kernel (null space) of the transpose of H The number of dimensionless variables is therefore equal to the nullity of H T . Each column vector contains the exponents of the dependent variables, such that their product is dimensionless, i.e., As is outlined in [41], a function f 1 can be constructed such that which implies that one dimensionless variable is a function of the (all) other, i.e., Solving Equation (A4) for Z yields therefore By comparing this result with the solution in Equation (1), π 1 is identified as jπΓ and π 2 = β as the Reynolds number. Therefore, the basic structure of the fluid impedance could be revealed by dimensional analysis alone. The idea is now applied to the more complex situation of a cantilever of rectangular cross-section mounted in a tube of radius R as shown in Figure 1. Performing a dimensional analysis with these parameters adds two aspect ratios as dimensionless variables e.g., π 3 = W/H = α 1 , and π 4 = W/R = α 2 . (A8) In this case, the fluid impedance per length, is Z = jωAρ Γ(β, α 1 , α 2 ).
Consequently, the hydrodynamic function Γ is, no matter how complex the geometry actually is, always a function of the Reynolds number β and all the aspect ratios which appear in the considered problem. For fluid sensing applications, the geometry of the sensor is fixed, but the fluid parameters ρ and ν appearing in the fore factor and in the hydrodynamic function are variable. Therefore, for any vibrating cross-section and configuration the form can be expected, where A is a parameter of dimension L 2 which is chosen to coincide with the area of the vibrating cross-section, as it will turn out to be convenient in Appendix B. The fluid dependent factor in the Reynolds number β defines a new variable ξ = √ ν/ω, which is used henceforth. It is mentioned for the sake of completeness that properties which are associated with the finite compressibility of the fluid, such as the bulk modulus and the longitudinal viscosity [42] were not considered, but would extend the π-group. Fluid sensors for density and viscosity sensing, are typically designed such that these influences need not to be considered.

Appendix B. Resonance Parameters
The deflection amplitude u(x, ω) of a single or double clamped prismatic beam of length L can be expressed by superposition of modeshapes φ i (x) times generalized coordinates q i (ω) The equivalent mass and the generalized force due to fluid loading for one particular mode of vibration can be determined from considering energy balances, yieldinḡ When the area A in the fore factor in Equation (A9) is chosen to coincide with the area of the resonator cross-section, they cancel out. Splitting the complex hydrodynamic function in Γ = Γ R − jΓ I and adding the fluid force to the spring-mass damper system driven with the generalized force F 0 , yields As a consequence, the resonance parameters of the fluid loaded resonator can be given solely in terms of vacuum resonance parameters ω 0 and Q 0 and the hydrodynamic function scaled with the ratio of fluid density ρ to resonator density ρ s Given Equation (A16), the hydrodynamic function can be determined from measurements in known test liquids by For unknown sensor parameters, the scaling of Γ and β remain unknown, but the complexity of Γ and the required order of approximation can readily be assessed from measurements. These scaling factors are implicitly accounted for by the model parameter calibration procedures in Sections 2.4 and 2.5.