Modelling of open photoacoustic resonators.

Photoacoustic spectroscopy employs acoustic resonators for signal amplification. Resonators are usually closed, however, in some applications, open resonators are preferred. The opening deteriorates the photoacoustic signal hence reducing the sensitivity of the photoacoustic measurement. We present two new approaches for simulating the photoacoustic signal in open resonators using finite element modelling. The approaches are based on the amplitude mode expansion model and the viscothermal model with the opening modelled using perfectly matched layers and the boundary element method respectively. Additionally, the performance of the viscothermal model using perfectly matched layers for simulating open resonators is extended to the ultrasound region. The simulation results are verified by comparing them to photoacoustic measurements. The approaches provide an accurate basis for designing and optimizing open resonators with high sensitivity.

Introduction

Photoacoustic spectroscopy (PAS) is a technique that detects acoustic waves generated from thermal de-excitation of molecules upon absorption of modulated electromagnetic radiation. The technique provides high sensitivity, high selectivity, large dynamic detection range making it suitable for trace gas applications [[1], [2], [3], [4]] and, unlike other optical techniques, can measure optically opaque samples [[5], [6], [7]].

Acoustic resonators are used to further increase the sensitivity of the photoacoustic (PA) signal. The modulation frequency of the electromagnetic radiation is matched to an acoustic eigenfrequency of the resonator thus exciting the corresponding acoustic mode. Conventional photoacoustic spectroscopy usually applies closed resonators (using windows): gaseous samples are filled into a tightly sealed resonator while solid samples are placed inside a closed resonator filled with a non-absorbing gas. However, the sensitivity of measurements from closed resonators is limited by window absorptions which produce a coherent background signal. This background signal can vastly affect measurements of samples of low concentration where the PA signal can be weaker than the background signal. To avoid the problem, authors have proposed reducing the background signal by placing the windows at the acoustical nodal points [8], attenuating the background signal by use of buffer volumes next to the window, making 2- measurements and using windowless resonators which eliminates the problem.

Windowless or open resonators have the additional advantage of enabling continuous exchange of fluids between resonator and environment. This is particularly useful for the continuous PA measurement of blood glucose levels since closed resonators result in the build-up of humidity inside the cell [9]. Open resonators are instead used to allow the humidity to evaporate, minimize the effects of pressure and temperature fluctuations during measurement hence increasing the long-term stability of the measurements [10]. In fluid applications, open resonators have been reported to eliminate the need of pumps and valves due to natural flow of fluids thus simplifying the PA set up [[11], [12], [13]]. The consequence of having an open resonator is the deterioration of the PA signal and penetration of external acoustic signal into the resonator.

A strong PA signal is crucial for the development of PA sensors. The strength of the PA signal is mainly influenced by resonance amplification. Accurate modelling of the photoacoustic signal of an open resonator would provide a good framework for optimizing its geometry with reference to the signal amplification, hence increasing the sensitivity. To the best knowledge of the authors, only Lang etal. have simulated the photoacoustic signal of an open cell [14] employing the viscothermal (VT) model with perfectly matched layers (PML). Their model estimated the signal of nitrogen dioxide in a simple cylindrical resonator. A Gaussian shaped heat source was implemented diagonally across the resonator to mimic the absorption of the laser beam by the nitrogen dioxide molecules. Their measurement was restricted to a single resonance peak in the audible region. Miklós et al. [11] used an analytic solution in combination with the transmission line model to calculate the photoacoustic signal of a simple open cylindrical resonator. However, the method is limited to resonators of simple geometry and is not appropriate for resonators with complex geometry.

In this paper, we simulate the photoacoustic signal of solid samples in a complex open T-cell resonator which has lately been used for non-invasive photoacoustic measurements of blood glucose levels [10,15]. The resonator is designed for ultrasound operation. We have previously described and demonstrated the amplitude mode expansion model (AME) as an accurate method for simulating the photoacoustic signal in closed cell macro-resonators [16]. Here, the method is extended to open resonators by using PML at the opening. Additionally, we extend the study of the viscothermal method with a PML into the ultrasound region and present a new method for modelling the surrounding of the resonator with the boundary element method (BEM). The simulations are done over a wide frequency range between 10 kHz and 60 kHz and are compared to the experimental measurements of the resonator’s PA signal.

Method

Simulation

This section will briefly describe how the models are realised and implemented. The investigated T-cell resonator consists of 3 interconnected cylinders as seen in Fig. 1.

Fig. 1

Left: The T-cell resonator with the openings highlighted in blue. Right: the schematic of its cross-section (light grey) showing the size of each cylinder in mm.

The description of how the cylinders are interconnected and their dimensions are published elsewhere [16]. The resonator has three openings: A small opening where a person would place their hand during blood glucose measurements (Fig. 1 right: absorption cylinder) and a small opening at the end of the resonance cylinder where the microphone is located. Due to the conjuncture that these openings are sealed by the hand of the diabetic patient and by the microphone respectively, these openings can be considered as effectively closed. This leaves the big opening at the opposite end of the cavity cylinder for ventilation.

Viscothermal model

The model is set up using the thermoviscous acoustic module of COMSOL Multiphysics software which solves the linearized Navier-Stokes equation supplemented by the continuity equation and a balance equation for the energy:

, , , , , , and are density, velocity, pressure, dynamic viscosity, bulk viscosity, specific heat at constant pressure, temperature and the coefficient of thermal expansion of the fluid, respectively. The field describes the spatial distribution of the heat source. It is defined at the end of the absorption cylinder with its size restricted within the absorption cylinder. The spatial distribution can be found in our article on the simulation of closed PA resonators [16].

The effects of thermal conductivity and viscosity are explicitly included in the equations. This is important since the resonator is of small size and surface thermal and viscous losses are the main loss effects in small resonators. The derivation of these governing equations has been discussed in depth elsewhere [17].

Air is selected as the propagating fluid from the COMSOL Multiphysics database. The temperature is set to 20 °C and the static pressure to 1013 hPa. The solid sample in the model is simulated with a heat source term to mimic the thermal de-excitation of the sample’s glucose molecules after laser excitation. The heat source is defined at the small opening of the resonator and limited within the cross-sectional area of the opening.

The big ventilation opening of the resonator is modelled using two different approaches: boundary element method and the perfectly matched layers.

The viscothermal simulations were done using 2 Intel Xeon CPUs with 14 cores each, running at 2.6 GHz with 384 GB of RAM.

Perfectly matched layers

The PML is an absorbing domain that is added to the exterior of the model. The ventilation opening of the resonator is modelled by creating an additional domain using a hemisphere to represent the free space. The PML form the end of the hemisphere that supresses reflection of the radiated sound waves into the resonator and models the extension of the free space to infinity (Fig. 2).

Fig. 2

Grey: The resonator and the hemisphere representing the open end. Blue: The PML domain.

The hemisphere’s radius is four times the radius of the opening while the thickness of the PML is twice the radius of the opening. The walls of the resonator (including the small openings sealed by the skin of the hand and the microphone) and the flanged edge are set with no slip and isothermal boundary conditions. The method is used to simulate the PA signal of the resonator between 10 kHz and 60 kHz with 50 Hz increments.

A swept mesh is applied to the PML while the resonator is meshed using prism elements. Boundary layers are generated along the interior resonator walls to accurately map the thermal and viscous surface loss effects. A convergence study is performed using two different meshes as seen in Table 1.

Table 1

Properties of the meshes generated for the PML simulations.

	Fine mesh	Finest mesh
Max. element size (mm)	1.429	1.2
Min. element size (mm)	0.481	0.481
Growth rate	1.4	1.1
DOF	654,416	1,176,891

Additionally, the convergence of the PML is verified by increasing and decreasing the thickness of the layer by 10 % and simulating the PA signal of a resonance between 46 kHz and 51 kHz with 50 Hz increments. The non-invasive blood glucose measurements are intended to utilize the resonance amplification from this specific resonance hence its selection.

The boundary layer used in this study consists of 4 layers with the layer closest to the resonator wall having a thickness of 0.02 mm. This has been found to be sufficient even at 60 kHz. The confirmation is done by simulating the high frequency range of our study using 8 layers with the thickness of the layer closest to the wall being a third of the thermal penetration depth. Such a refined boundary layer requires long simulation time and we did not simulate the entire frequency range with this mesh.

Boundary element method

The acoustic boundary element method in COMSOL solves the Helmholtz equation for pressure variations in fluids. The method is suitable for modelling large and infinite domains such as the free space exterior to the resonator domain. At the open boundary, a multiphysics coupling between the viscothermal module and the BEM method is established. The method reduces the number of degrees of freedom since only the open boundary needs to be meshed unlike the PML where the hemisphere is also meshed.

The method is used to simulate the PA signal of the resonator between 10 kHz and 60 kHz with 100 Hz intervals. A convergence study of the method is performed using two different meshes as seen in Table 2. Details of the convergence studies are discussed later.

Table 2

Properties of the meshes generated for the BEM simulations.

	Fine mesh	Finer mesh
Max. element size (mm)	1.429	1.2
Min. element size (mm)	0.481	0.481
Growth rate	1.4	1.1
DOF	205,272	407,242

Amplitude mode expansion model

The AME method calculates the acoustic pressure at the location of the microphone and the modulation frequency as a superposition of the resonator’s acoustical eigenmodes Here, is set to 50 modes. It has been confirmed that increasing the value of does not affect the results. The eigenmodes and the eigenfrequencies of the resonator are calculated by solving the Helmholtz equation with a sound hard boundary condition at the resonator’s interior walls. As the AME method is usually employed for closed cavities, we have extended the method by adding a PML. The ventilation opening is modelled using a hemisphere which is bounded by a PML. So-called spurious modes that exist mainly in the hemisphere domain are removed and only the real resonator modes are selected during modal summation of Eq. (4). This is achieved by comparing the integrated absolute pressure within the resonator and the hemisphere. The modes that have higher integrated absolute pressure in the hemisphere compared to the resonator are defined as spurious modes. From the 50 modes, only 29 modes qualify as real modes and used in the summation. Some examples of spurious modes are shown in Fig. 3. The BEM approach was not implemented in this model since COMSOL Multiphysics does not support eigenfrequency studies with the BEM module.

Fig. 3

Distribution of the absolute value of the pressure for some spurious mode (from left to right: 18.5 kHz, 22.8 kHz, 31.8 kHz). Green represents nodes, red antinodes. Absolute values are of no importance, since in the AME model the modes need to be normalized.

The amplitudes are calculated usingwhere the describe the excitation of the sound waves and are obtained from

is the volume of the resonator, is the ratio of isobaric and isochoric heat capacity and is the absorption coefficient of the sample. The laser intensity profile in the sample , like the heat source term in the viscothermal model, represents the thermal de-excitation of the sample’s molecules and is similarly defined in the resonator’s absorption cylinder. Loss is introduced by loss factors in Eq. (2), which account for the thermal and viscous losses at the bulk of the fluid and at the resonator’s surface. Detailed descriptions of the formulas can be found elsewhere [18].

Air is selected as the propagating fluid with the parameters in Table 3.

Table 3

Air parameters at 20 °C and a static pressure of 1013 hPa [19].

Density ρ	1.2044 kg/m^3
Sound velocity c	343.2 m/s
Viscosity μ	1.814 × 10−5 Pa s
Coefficient of heat conduction α	2.58 × 10−2 W/m K
Specific heat capacity at constant volume C_v	7.1816 × 102 J/kg K
Specific heat capacity at constant pressure C_p	1.0054 × 103 J/kg K

A swept mesh is applied on PML while prism mesh is generated on the resonator. Since the finite element software COMSOL is only used to calculate the eigenfrequencies and modes, boundary layers are not generated for this model. To reduce the number of spurious modes, the size of the hemisphere is reduced in comparison to the hemisphere in the VT model described above. The simulations are performed with hemispheres of different sizes to check for convergence of the model. The radii of the hemispheres are 1.610, 1.785 and 1.800 times the radius of the opening as seen in Fig. 4.

Fig. 4

Mesh of the resonator with three different hemisphere sizes. From left to right: HS 1, HS 2 and HS 3 are hemispheres with radius 1.610, 1.785 and 1.800 times the radius of the opening respectively.

The simulations of the AME model were done on a 64-bit computer with a processor speed of 2.5 GHz and 32 GB RAM.

Experiment

To evaluate the results from the simulation methods described above, we measured the resonator’s frequency response using a set up described in our previous work [16]. The window sealing the cavity cylinder is removed for the measurements as described in section 2.1. An average of 10 acquisitions were recorded when the laser modulation frequency is swept between 10 kHz and 60 kHz with an increment of 10 Hz.. The frequency response of the microphone is fairly flat between 35 kHz and 60 kHz while between 10 kHz and 35 kHz it has a resonance with a peak at around 26 kHz as can be seen in Fig. 5.

Fig. 5

Frequency response of the microphone used for the measurements.

Results and discussion

Convergence study

The results from the convergence study of simulations with different mesh parameters for the VT model (see Table 1, Table 2) and different hemisphere size for the AME model (see Fig. 4) are presented in Fig. 6. The plots have been rescaled to avoid overlapping of the curves of the different approaches.

Fig. 6

Frequency response of the different models. The upper solid lines represent the VT-PML results, the middle-dashed lines represent the results from the VT-BEM approach (blue: mesh 1, red: mesh 2). The lower full lines represent the results from AME-PML model with different hemispheres (grey, red and green refer to HS 1, HS 2 and HS 3 respectively).

The plots from the VT simulations with the PML and BEM approach show that the model is converging towards the solution as the different meshes produce similar spectral features. The curves slightly deviate at frequencies above ca. 55 kHz, which indicates that the mesh resolution is no more sufficient in this region. The plots from the AME model with different hemisphere sizes have also similar spectral features indicating a convergence of the model towards the solution. The different hemispheres have the same resonant frequencies albeit with small differences in the amplitudes.

The convergence of the PML is also checked by changing the thickness of the PML layer and simulating the PA signal between 46 kHz and 52 kHz with the results shown in Fig. 7. When the thickness of the layer is reduced by 10 %, the size of the hemisphere’s radius is increased by 10 %.

Fig. 7

Frequency response of the VT model with different PML thickness against the experiment.

Changes in the PML thickness do not significantly change the results. The amplitude of the resonance is slightly changed, however, the resonance frequency is the same. It is noted that while the magnitude of the amplitude decreases a peak-like feature between 48 kHz and 49 kHz becomes more prominent. The result obtained with the larger hemisphere (larger air volume) resembles the measurement more than the other two configurations as observed in Fig. 7 where the peak-like feature of the resonance between 48 kHz and 49 kHz is also not very prominent.

Comparison of the simulation models

Fig. 8 shows the frequency response of the measurement and the simulation of the open T-cell resonator using the three different approaches described. The experimental data has been smoothed using the Savitzky Golay Matlab function [20].

Fig. 8

Frequency response of the experiment against the simulations. The blue and red plots represent the VT model using boundary element method and perfectly matched layers respectively. The green curve is the result of the AME model and the blue dashed line shows the experimental measurements.

The VT and the AME models predict the same resonances with a difference in the resonant frequencies of less than 0.9 %. The resonant frequencies of the AME model are slightly shifted to higher values compared to the VT model. This behaviour is expected since the thermal and viscous losses in the AME model are known to be underestimated when using quality factors, whereas the VT model accurately captures the losses using boundary layers [18]. This is further reflected by the narrower resonance widths of the AME model.

The resonances of the two different VT approaches are in very good agreement. The resonant frequencies differ by less than 0.2 %. The amplitudes of the BEM approach are higher compared to the PML indicating that the loss predicted by this approach is slightly lower than the overall losses of the PML. Consequently, the resonant frequencies are also slightly shifted to higher frequencies compared to the frequencies of the model with the PML.

The wide resonance between 45 kHz and 52 kHz is due to superposition of two resonances. The corresponding modes are depicted in Fig. 9. The BEM approach and the AME method predict a small peak like feature on the wing of the resonance at around 48.9 kHz. However, the second peak is diminished with the PML approach. The difference is attributed to the fact that the BEM approach and the AME model underestimate the loss effects since the additional peak is also rather minor in the experimental measurement. The VT model with the PML is therefore rated as the most reliable of the three different simulation approaches since it reflects the experimental measurements most accurately.

Fig. 9

The acoustic modes showing the absolute value of the pressure at 49 kHz and 49.9 kHz. The plots are normalized by the highest value of the individual plot.

The AME model has already been established as a very fast method, taking only a few minutes to obtain a complete response curve [18]. The VT model with BEM takes ca. 21 h to simulate a 10 kHz range with 10 Hz increments, while the PML approach takes ca. 102 h to simulate the same range. The difference is mainly attributed to the extra number of degree of freedoms that the PML approach generates due to the hemisphere created at the opening.

Comparison of the models with measurements

The simulation results show good accordance with the measurements. All the major resonances that the models predict are experimentally measured. The relative difference in the peak resonance frequency between experiment and simulation is not more than 1.04 %. In this section, the measured results are compared against the VT approach with the PML since this model shows the best agreement with the measurement. The base noise level of the measurements is high and as a result, the measured resonances with small amplitudes are difficult to distinguish.

The signal to noise ratio (SNR) of the measurement at the resonance with the highest amplitude is 16.46. It is defined as the ratio of the mean to the standard deviation of the measurements. The amplitude of the simulated resonance at 28 kHz is three times smaller compared to the highest amplitude at 49.8 kHz. However, the measured amplitude at 28 kHz is more than half the amplitude of the resonance at 49.8 kHz. The difference is attributed to additional amplification by the microphone. As previously mentioned in Section 2.2, the microphone exhibits a resonance peak in this frequency range hence the higher amplitude is as a result of synergetic amplification by the resonator and the microphone.

The widths of the simulation resonances are narrower than those of the measurements. This indicates that the damping in the measurement is higher than what the model predicts. The higher damping is suspected to result from signal leakage from the microphone mount which does not seal the resonator tightly. The higher damping in our measurement is further indicated when the factor of the two prominent resonances at 28 kHz and 49.8 kHz are compared against the simulation. The factor is defined as the ratio between the resonance frequency and the width of the resonance peak at half its resonance amplitude (FWHM). It is estimated usingwhere denotes the resonance frequency while and are the frequencies at which the value of the pressure amplitude has decreased to half of its value. The factors of the measurement are 9.74 and 22.46 while the VT model yields 51 and 30 with PML and 54.88 and 29.77 with BEM, respectively.

Conclusion

The photoacoustic signal of an open resonator of complex geometry has been accurately modelled using three different approaches. The AME method especially stands out due to its significantly faster computational time and reasonable precision. The VT model is more accurate than the AME model however it is computationally much more demanding. To reduce the computational time of the VT model, it is reasonable to use the BEM approach. The methods described are not limited to simulating the PA signal of T-shaped resonators only. They can be used to simulate open resonators of arbitrary shape. The ability and accuracy of the methods have been demonstrated in the audio and ultrasound region. With small changes the methods can be adapted to gaseous samples.

The methods described can be used in the design and optimization of the photoacoustic signal of open resonators. Future studies are planned where the methods will be used to optimize the geometry of the T-cell resonator to obtain maximum signal amplification. We intend to improve the sensitivity of the PA measurements employing the T-cell resonator leading to the development of a continuous non-invasive blood glucose sensor.

Funding

We would like to thank the Hamburg University of Applied Science Promotionsförderung for funding.

Declaration of Competing Interest

The authors declare no conflict of interest.