Analysis of the sensitivity to pressure and temperature of a membrane based SAW sensor

ABSTRACT This paper presents a FEM analysis of a membrane-based Surface Acoustic Wave (SAW) sensor. The sensor is a 2.45GHz Reflective Delay Line (R-DL) based on Lithium Niobate (LiNbO3). As the wave propagation time is much smaller than the typical time constant of the phenomena to be monitored (deformation, temperature change etc.), the analysis can be performed in three successive steps. First, a static FEM study of the complete sensor (housing included) is carried out, to compute the temperature, stress and strain fields generated in the sensitive area by the measured parameters (pressure, temperature, etc.). Then, a dynamic electro-mechanical study of the R-DL is performed. The simulation takes the previously computed fields into account, which makes it possible to compute the sensor sensitivity to the measured parameters. The model takes advantage of the periodicity of the components of the R-DL to compute phenomenological parameters (Coupling-of-Mode parameters), which can later on be used to compute the electrical response of the sensor (step 3). In this paper, we focus on the first two steps. The COM parameters are extracted, under simultaneous thermal and mechanical stresses. Especially, the sensor sensitivity is obtained from the evolution of the velocity, under various stress configurations.


Introduction
Surface acoustic wave (SAW) devices use interdigital transducers (IDT, comb shaped electrodes) to generate waves at the surface of a piezoelectric substrate. The velocity of the generated waves changes with the surface state (stress and strain, mass deposition or temperature), which makes it possible to use these devices as sensors by arranging the IDTs in delay lines or resonator structures. Due to their simplicity, it is possible to build robust SAW sensors. In addition, their passiveness and capability to be interrogated wirelessly make them good candidates for harsh environment sensing applications. Besides being of potential interest for structural health monitoring [1], pressure measurement is a promising application for the SAW devices. Authors have already proven the feasibility of such sensors using both resonator [2,3] and reflective delay line [4][5][6][7] designs.
An accurate finite element simulation of such devices is of interest when it comes to reducing their development time and improving their sensitivity. However, due to the high aspect ratio of the device (the whole sensor is usually in the millimeter range while a single finger of an IDT is in the micrometer range), the simulation of a whole device using the finite element method (FEM) would be too time consuming. That is why the simulation of a whole SAW sensor should be carried out in three steps which is made possible because the effect to measure is slow compared to the wave propagation phenomenon. In a first step, the effect of the measurand (in terms of temperature, strain and stress fields) is computed along the propagative path of the sensor. It is possible here to take the influence of the sensor housing into account. In a second step, the socalled coupling of mode (COM) parameters are extracted from the simulation of a simplified infinite grating of electrodes [8,9]. These phenomenological parameters are computed for all the surface states present in the acoustic path (i.e. IDT, reflector gratings, free surface. . .). They can already be used to compute the sensitivity of the sensor with respect to temperature, strain and stress. Finally, in a third step, the frequency response of the whole device can be computed using the obtained sets of COM parameters and the device geometry (position and number of electrodes, aperture. . .). The obtained complete set of COM parameters allows a fine tuning of the RF response of the sensor [4]. A similar approach has been developed in [10] to optimize a mass sensitive sensor based on a resonator design. This paper focuses on the two first steps since they are sufficient to compute the sensor sensitivity. They are illustrated in this paper using a CTR sensor [11] as an example ( Figure 1). The latter is aimed at detecting changes in intra-cranial pressure. It is made of a thin membrane that is deformed under pressure. This deformation generates stress and strain fields that affect the wave velocity. The sensitivity of the sensor to pressure and temperature has been computed. It is compared to experimental results, at the end of Section 4. In our case, the computation of the 'free surface' wave velocity is sufficient to calculate the sensor sensitivity. However, the whole set of COM parameters is computed to illustrate the general computation procedure, as described above.

Coupling of mode theory
The COM model makes it possible to simulate SAW devices in an efficient way. Nevertheless, the accuracy of the model depends on its parameters (the so-called COM parameters). COM parameters can be determined experimentally. It is, however, a tedious procedure. For this reason, it is advantageous to extract them from a numerical model.
The COM parameters are determined from an analysis of the behavior of waves propagating in an infinite grating. The main parameters of such a grating are illustrated in Figure 2 (p is the pitch, W is the aperture, h and a are respectively the electrode thickness and width).
The COM theory considers two mechanical waves u þ and u À propagating in the positive and negative directions. When traveling through the grating, these waves can be reflected and contribute to the wave propagating in the opposite direction. The COM model considers only the reflection of the first harmonic. It also assumes that the amplitudes may decrease due to attenuation and that there is a linear transduction coefficient due to piezoelectricity. This leads to the following equations [8,9]: Where v SAW (m=s) is the SAW velocity, γ(m À1 ) is the attenuation, κ 12 (%=m) is the reflection coefficient, ζ (Ω À1=2 ) is the transduction coefficient and C (F=period) is the capacitance.
Assuming that the waves are 2p-periodic: In this matrix, θ u ¼ ω v SAW À π p À jγ (rad=m) is the detuning factor. The behavior of the system can be fully described using the five COM parameters (v SAW , κ 12 , ζ, C, γ). These parameters are determined from the simulation of an infinite periodic grating, using different electrical boundary conditions, as described in details in [8,9].

Short circuit behavior
When considering no attenuation (γ ¼ 0), the equations system (6) under short-circuit (SC) conditions (i.e. V ¼ 0) corresponds to an homogeneous system. The grating behaves like if there are no electric field and transduction (ζ ¼ 0): We look for solutions of the form: Substituting the solution in (7), we get this system: Which has non-zero solutions only if the matrix determinant is zero: The general solutions are: Where A and B are constants determined by the boundary conditions and q is the wave number of the existing mode, given by: Substituting the detuning factor (θ u ) in the dispersion relation (12) we find a quadratic polynomial equation in terms of frequency. The two solutions of this polynome are the two frequencies that define the stop band edges of the dispersion relation for the short circuit condition: The average of these two frequencies and their difference give the SAW velocity and the reflection coefficient:

Open grating behavior
In open grating (OG) condition, the current is zero. Thus, the third line in (6) gives this expression of the voltage [8,9]: Substituting this relation in the first two lines of Equation (6) yields: A non-trivial solution exists if the matrix determinant is zero: This gives the two following solutions: These two frequencies define the stop band edges in open grating condition. In addition, we can assume that the coupling and transduction coefficients are real [8].
Then, adding the two stop band frequencies for the open grating condition yields: Therefore, from the latter and the first of (14), we get the relation between ζ and C: We need another relation between ζ and C to determine them. This can be obtained from the harmonic admittance of the IDT.

Harmonic admittance
To determine the harmonic admittance in the case of no attenuation (γ ¼ 0), we apply a harmonic driving voltage [8,9]. Since the potential on the electrode fingers does not depend on the position x, a particular solution of the mechanical equations (first two lines) in (6) has the form:R Substituting this in the first two lines of (6) provides: Since dV dx ¼ 0 and after simplification by V we get: Inverting this system, yields fA þ A À g. Thus, the solutions (21) become: Substituting further this solution in the third line of (6) yields the harmonic admittance: This harmonic admittance shows a resonance and an anti-resonance at two different frequencies. After substituting the already known COM parameters (v SAW , κ 12 , ζðCÞ), it becomes possible to extract the capacitance (C) by fitting the simulated harmonic admittance using Equation (26).

FEM simulation of an infinite periodic grating on LiNbO 3
Simulation of the infinite periodic grating can be done using COMSOL Multiphysics [12]. The model is implemented in 3D. From this model, COM parameters can be extracted using specific boundary conditions reproducing the three configurations (short circuit, open grating, harmonic admittance), described in the previous section.

Geometry
It is possible to reduce the model geometry to a pair of fingers, using periodic boundary conditions on potential and displacement. In addition, as the aperture of the device is considered infinite, it can be (in the model) reduced to a fraction of the wavelength by using periodic boundary conditions. The substrate thickness is also reduced to a few wavelengths, using a perfectly matched layer (PML) placed at the bottom of the substrate [13]. This PML avoids reflections of acoustics modes at the bottom of the substrate. An air layer is also taken into account above the substrate. Nevertheless, LiNbO 3 has a high permittivity. Thus, air is expected to have a relatively weak influence on the results. The resulting model is often called 'unit cell model'. It is shown in Figure 3.

Governing equations
As the electrodes have a high conductivity compared to the piezoelectric substrate or the air layer, their conductivity is assumed to be infinite. Consequently, the electrodes are not simulated from an electrostatic point of view. Also, the air layer is supposed to have no stiffness or mass. As a result, the mechanical equation is not solved in the air layer. The electromechanical behavior of the piezoelectric material is described by the following constitutive equations that couple the stress σ and strain S to the electrical field E and displacement D: Where c E is the short circuit stiffness matrix, e is the stress piezoelectric matrix, and ε S is the blocked permittivity matrix. The material constants used in this study are presented in Appendix. In addition, assuming small displacements, the strains are S kl ¼ 1 2 ðu k;l þ u l;k Þ, and E k ¼ ÀV ;k . Thus the mechanical and electrostatic equations are (in the absence of body and surface charges and body force): The boundary conditions used for the electrostatic and mechanical equations are detailed in Table 1.
It is also possible to compute the COM parameters for an initially pre-stressed structure. Under the assumption that the displacement causing the waves is small compared to the initial displacement responsible for the pre-stress field (σ 0 ), Equation (28) becomes [14]: In addition, to consider the effect of the initial strain on the material stiffness and density, the stiffness tensor and the density must be modified according to their respective sensitivity to the applied strain [15]. This is discussed hereafter.

Material data perturbations
The effect of initial strain on the stiffness constants is described as follows (linear dependency) [15]:c E pq ¼ c E pq þ c pqr S r p; q; r ¼ 1; :::; 6 (31) The constants c pqr are called third-order stiffness coefficients. They are determined experimentally [16] and generally presented in a 6 Â 6 Â 6 matrix. However, due to crystal symmetries, only few coefficients are needed to depict the entire matrix. In the case of LiNbO 3 , 14 coefficients are sufficient. Besides, the strain field affects the density. The change in density is directly proportional to the change in volume [15]: The effect of temperature is easier to consider since it only affects material parameters. Assuming that the material properties are linearly dependent on the temperature variation (ΔT), we get: The temperature dependent coefficients (TC) are determined experimentally and are well known for LiNbO 3 [17]. All the constants used to carry out this analysis have been reported in Appendix. As experimental values are not available in the literature [16,17], temperature variations similar to (33)-(35) are not used for the third-order nonlinear elastic constants of (31). In general, the temperature derivatives of the third-order elastic constants are not known [18].

Mesh
The structure is meshed using prismatic elements having quadratic interpolations. The mesh has been refined close to the electrode edges to better consider the charge singularity that occurs there. In addition, the mesh is made coarser along the depth of the substrate to follow the attenuation of the Rayleigh wave. An example of the mesh used for this study can be seen in Figure 4. This particular mesh has 2095 elements and 34602 nodes. It has been refined iteratively, thus increasing the number of degrees of freedom (nDOF) to ensure a proper convergence of the solution. The convergence is monitored taking the COM parameters as figures of merit (see Figure 4).

Simulation results
The computational procedure described above has been tested on a sensor that is currently in development at CTR (see Figure 1). This implantable device can detect very low pressure variation in the liquid surrounding the brain. It can also sense temperature. It is made of a thin LiNbO 3 circular membrane deposited on a silicon substrate. The acoustic path (the direction along which the acoustic waves travel) is oriented along the crystallographic Z axis. The sensitivity of the sensor to pressure and temperature can be computed and compared to the experimental results.

Simulation of the SAW sensor sensitivity to pressure
To determine the SAW sensor sensitivity to pressure only, the simulation is carried out at the reference temperature T 0 (in our case, T 0 ¼ 20 C). When the membrane is deformed under pressure, the generated strain and stress fields affect the wave propagation. As a first step, the whole sensor is simulated in static state, without considering the acoustic phenomena. Taking advantage of symmetries in material, geometry and load, only one forth of the SAW sensor is modeled. The boundary conditions are illustrated in Figure 5. From this simulation, the stress and strain fields are extracted along the acoustic path. The latter is divided into two sections that are designed to optimize the sensor sensitivity (see Figure 5). Sections 1 and 2 correspond to two different pressure and temperature sensors, which are designed to work in a so-called dual mode. They are identically sensitive to temperature variations, but are subjected to strain fields of opposite signs. By subtracting the signals of both sensors, it is possible to minimize the temperature sensitivity, while increasing the pressure sensitivity. More design details are given in [11]. Since we consider only a linear effect of the perturbation field on the material properties, only the average values over the section lengths of the fields are taken as inputs for the computation of the COM parameters.
We summarize in Table 2 the obtained COM parameters, for the two sections, when a 3000Pa pressure is applied on the membrane. For comparison purposes, we also computed the initial state, where no pressure is applied.
It can be noticed from these results that loaded section 1 COM-parameters values are bounded from below by values for loaded section 2 and from above by those of the unloaded section 1. This is true for the first three parameters while the reverse bounding happens for the fourth one. It can be explained by the fact that the COM parameters are more influenced in section 2 than in section 1, although section 1 is locally subjected to  a higher level of stress (notably near the edges). This is because section 1 extends over the substrate where there is almost no stress/strain at all. As already mentioned above, the strain and stress fields are first averaged over the whole section before computing the COM parameters, which leads to section 1 COM parameters being less affected by the stress/strain fields than section 2 COM parameters, although section 1 comprises the point of maximum stress/strain. As mentioned in the introduction, the sensitivity (the change in phase due to the load) is computed from the velocity by using the equation: Where ϕ 0 and v 0 are the initial phase and SAW velocity in the section (when no load is applied), S and v are respectively the strain in the propagation direction and velocity in the loaded section. The sensitivity is compared to the experimental results (see Figure 6, left).
It is worthy to mention that, although pressure and temperature variations were considered here individually, the present analysis can be used for their simultaneous variations. In this case, the strain S, in (36), should be replaced by S tot ¼ S þ S th , where S th is the thermally induced strain while S is due to the pressure variation.

Simulation of the SAW sensor sensitivity to temperature only
The sensitivity to temperature is computed without any pressure applied. As the sensor is made of two different materials, strain and stress fields are generated when the sensor is heated up (bimorph effect). This affects the wave propagation conditions, in addition to the pure temperature effect itself.
A static thermomechanical study can be carried out to determine the initial stress and strain fields due to the bimorph effect, using the same model as previously. The resulting fields are significant. They have the same order of magnitude as that computed when applying pressure only (see Figure 7).
Then, the acoustic behavior can be computed using the unit cell model. The COM parameters were extracted for a temperature T 1 ¼ 30 C. The parameters are given in Table 3 and compared to the COM parameters of a reference pure LiNbO 3 at the same temperature (T 1 ).  Figure 6. Comparison of the experimental and simulated sensitivities to pressure (left) and temperature (right).
These parameters (Tables 2 and 3) allow to compute the sensor sensitivity to pressure and temperature taking the bimorph effect (thermal stress) into account. First of all, we can see that the pressure sensitivity is in almost perfect agreement with the experimental data. Concerning the temperature, the results demonstrate that taking the thermal stress into account improves the accuracy of the simulation. However, there is still a small discrepancy between computed and experimental data which could be explained by the unknown pre-stressed state of the membrane due to manufacturing process or to the inaccuracy of the third order material coefficients. This small deviation may be also due to the incompleteness of the changed elastic constants in the strained material (31). Indeed, differently from [15] and [16] where expanded internal energy, thus thermodynamic tension, were used, it was shown in [18][19][20], using proper (rotationally invariant) nonlinear electro-elastic equations and a special thermodynamic functional, that the effective nonlinear elastic constants present two additional terms that are proportional to the displacement gradients. Nevertheless, the latter are not direct outputs of commercial finite element codes. Thus, they need non-worthy (compared to the expected gained accuracy versus experimental results) heavy numerical and non-accurate (due to the required numerical differentiations) post-treatments. For these reasons, it was preferred here to use the simpler relation (31). This practice is admitted by the ANS/IEEE standard on piezoelectricity [21] (see page 42 therein).

Conclusion
This paper reports an efficient and accurate computational procedure for designing SAW sensors using the finite element method; it is demonstrated on a sensor sensitive to both  Figure 7. Stress in the principal directions over the acoustic path generated by a 10 C increase in temperature (σ zz ¼ 0). pressure and temperature. A static simulation of the complete sensor is carried out to consider the effect of the measurand. Then, the COM parameters are extracted for a substrate stressed by temperature and/or pressure. Among these parameters, the velocity is used to compute the sensitivity of the sensor and it is compared to experimental measurements. The FEM analysis shows a 0:08rad=mbar sensitivity to pressure and a 3:75rad= , that is in good agreement with the experimental results within the respective range, 0 À 30mbar and 20 À 30 C. The procedure is accurate enough to enable further improvement of the sensor design. In particular, it would be interesting to study the effect of the certainly occurring buckling of the membrane when the temperature increases.

Appendix. Material constants
The material constants used in the study are taken from [22]; they are given in Table A1. The Table A2 summarizes the 14 third order coefficients needed to determine the dependency of LiNbO 3 to a strain field. The complete set of constants can be calculated using the crystal symmetries; the relations can be found in [23].
The Table A3 gives the temperature dependent coefficients of LiNbO 3 together with the expansion coefficients [17].