Solving a coupled field problem by eigenmode expansion and finite element method June 18 , 2007

The propagation of sound in fluids is governed by a set of coupled partial differential equations supplemented by an appropriate equation of state. In many cases of practical importance one restricts attention to the case of ideal fluids with vanishing transport coefficients. Then, the differential equations decouple and sound propagation can be described by the wave equation. However, when loss mechanisms are important, this is in general not possible and the full set of equations has to be considered. For photoacoustic cells, an alternative procedure has been used for the calculation of the photoacoustic signal of cylinder shaped cells. The method is based on an expansion of the sound pressure in terms of eigenmodes and the incorporation of loss through quality factors of various physical origins. In this paper, we demonstrate that the method can successfully be applied to photoacoustic cells of unconventional geometry using finite element analysis.


INTRODUCTION
Photoacoustic spectroscopy finds many applications in the field of concentration measurements of gaseous compounds.The photoacoustic effect utilizes the fact that the excitation energy of light absorbing molecules is transfered into a local pressure increase in the absorbing gas.If the excitation source is modulated, a sound wave is generated that can be detected by a microphone (Figure 1).The signal detection sensitivity of photoacoustic sensors can be considerably improved by taking advantage of acoustical cell resonances, i.e., the radiation is modulated at a frequency equivalent to an acoustical eigenmode of the measuring chamber [1].
During the last decades, considerable effort has been made to develop photoacoustic sensors.One of several fields of interest in this context is the investigation of the propagation of the sound waves in the photoacoustic cell [2].While for cylinder shaped cells some acoustical quantities of interest can be calculated analytically other quantities are not amenable to an analytical treatment [3].In recent years, photoacoustic cells of unconventional cell geometry have come into focus [4,5].In this case, analytical calculations are generally not possible.
The finite element method on the other hand offers the possibility to perform calculations where analytical results cannot be obtained.First attempts to utilize the finite element method (FEM) for the investigation of the photoacoustic signal generation have been described in [6,7,8].However, only modern software tools and hardware allow an efficient application of this method.Recently, first results have been published [9,10].
In the next section the dominant loss mechanisms relevant in photoacoustic cells are presented.Thereafter, the differential equations describing the photoacoustic signal propagation are formed under consideration of the loss mechanisms.Some approaches for the solution as well as their limitations are briefly discussed [11].In section 4 a more appropriate solution method is described.This technique, combined with the finite element method, is used for the calculation of acoustical quality factors and amplitudes.

LOSS MECHANISMS
The resonance amplitude of the photoacoustic signal is determined by the strength of the excitation and by the amount of loss.Therefore, the wave equation is not sufficient for the calculation of the sound pressure at the diaphragm of the microphone and a more comprehensive treatment of the fluid properties is required.The following discussion is based on Morse, Ingard and Kreuzer [11,2].
The wave equation for sound waves is derived under the assumption that the fluid behaves adiabatically, i.e. thermal conductivity of the fluid is negligible.Then, no heat is exchanged between neighboring pressure maxima and minima.If on the other hand this heat exchange cannot be neglected, the energy density fluctuations will die out and the sound wave dissipates (Figure 2).A second loss mechanism is due to viscosity.Neighboring layers of the fluid move at different velocities when the sound wave passes (Figure 3).This generates shear stress and therefore, viscous friction.
The characteristic lengths of these loss mechanisms can be estimated by and with η being the coefficient of viscosity and κ the thermal conductivity.ρ represents the gas density, c the sound velocity and c p the heat capacity at constant pressure.These lengths are very small for many fluids of practical interest, which means that a substantial attenuation occurs only after the wave has traveled over large distances.It is well known that for small cavities like photoacoustic measuring cells volume loss is not the dominant source of loss.The dominant loss effects happen at the wall of the photoacoustic cell.They are due to heat conduction and viscosity as well.Typical cells are manufactured out of metal which has a high coefficient of thermal conductivity compared to the fluids.Therefore, temperature waves occur far away from the cell wall whereas the wall itself constitutes a region of isothermal behaviour.This results in a heat exchange as is indicated in Figure 4 (left part).The transition from the adiabatic to the isothermal behaviour takes place in a boundary layer of thickness (2) (ω is the angular frequency).
If viscosity is taken into account, one has a "no slip" condition at the wall, whereas the fluid moves due to the wave propagation in the interior of the container.Again, neighboring layers of the fluid have different velocities which results in viscous friction (Figure 4, right part).The thickness of Prandtl's boundary layer is There exist further loss mechanisms which are not discussed in this article [2,11].

DIFFERENTIAL EQUATIONS
For the calculation of the pressure distribution inside a cavity and at resonance, it is required to include loss.Therefore, the lossless wave equation is not an appropriate description of the fluid dynamics.Balance equations, however, represent an adequate starting point.In this case, the continuity equation (mass balance) (D/Dt substantive derivative, velocity) the Navier-Stokes equation (momentum balance) ( p pressure, µ expansion coefficient of viscosity) and the equation that expresses the energy balance (1.law of thermodynamics) are required.In the present context the first law is most useful in the form [12] 306 Solving a coupled field problem by eigenmode expansion and finite element method

Adiabatic region Gas
Prandtl's boundary layer Metal (T absolute temperature, β coefficient of thermal expansion, strain tensor).These partial differential equations have to be supplemented by an equation of state 1 , i. e.
This results in six equations for the six field quantities , T, ρ and p.
A first step towards the solution of these equations is linearization.All field quantities except are split into a large static and a small varying part: , , ( 5) is assumed to be small.Insertion into the differential equations and dropping all terms, that are proportional to the small quantities, lead to the linearized counterparts of the balance and state equations.When combined, the following two equations can be derived [11]: (7) (modified wave equation) and (8) (modified heat equation).γ denotes, as usual, the ratio of the specific heat at constant pressure c p to the specific heat at constant volume c V and α V = . The quantity of greatest interest, decouples from all field quantities but from T ~.
In order to avoid solution of the coupled system (Equations 7 and 8), different possibilities of simplification are possible: 1. Disregard heat conduction (κ = l κ = 0).In this case, temperature and acoustic pressure are proportional to each other (T ~= (γ -1)p ~/α V γ ) and,γ (p ~-α V T ~) equals p ~, which implies The problem decouples and the resulting wave equation includes a dissipative term.2. Maintaining heat conduction and assuming l κ ≈ l η leads to the following modified wave equation (10) where l is the mean of l κ and l η (see [3]).The mathematical structure of the two differential equations is identical.Performing a Fourier Transformation, we obtain (11) where Λ = l η in the first case and Λ = γ l in the second case, p(r r , ω) is the Fourier transform of the acoustic pressure is the wave number.In the case of photoacoustic cells, the sound waves are a result of the interaction of the laser beam and the molecules.The differential equations above have to be supplemented by a source term which accounts for the sound generation.

THE QUALITY FACTOR APPROACH
In this paper, we discuss a third approach that is particularly appropriate if one wants to include loss effects at the boundaries.It has been used in the theoretical treatment of photoacoustic cells of cylinder shape [2] and we show that it can be combined with the finite element method and applied to cells of arbitrary geometry.
Starting point is the loss free Helmholtz equation which describes the propagation of sound waves: H (r r , ω) is the Fourier transform of the power density (r r , t).Assuming that the absorbing transition is not saturated and that the modulation frequency of the light source is considerably smaller than the relaxation rate of the molecular transition, the relation H (r r , ω) = α I (r r , ω) applies, where I(r r , ω) is the Fourier transformed intensity of the electromagnetic field and α is the absorption coefficient [2].It is assumed that the walls of the photoacoustic cell are sound hard which leads to the boundary condition (13) i.e., the normal derivative of the pressure is zero at the boundary.
It is well known that the solution of the inhomogeneous wave equation can be expressed as a superposition of the acoustical modes of the photoacoustic cell: (14) The modes p j ( ) and the according eigenfrequencies, ω j = ck j , are obtained by solving the homogeneous Helmholtz Equation (15) under consideration of the boundary condition (13).In order to use Equation ( 14), it is necessary to normalize the modes according to (16) where V C denotes the volume of the photoacoustic cell and p * i the conjugate complex of p i .The amplitudes A j (ω) of the photoacoustic signal are determined from (17) with (18) The inhomogeneous Helmholtz Equation (12) does not contain terms that account for loss.Therefore, loss effects are included via the introduction of quality factors Q j in the amplitude Formula (17): The Stokes-Kirchhoff loss due to viscosity and thermal conduction in the gas volume is described by (20) Thermal conductivity surface loss results in (21) and surface loss due to viscosity in p j is the component of the pressure gradient tangential to the cell wall.The combined effect of all loss mechanisms lead to the following resulting quality factor (23)

RESULTS
Two types of photoacoustic cells are investigated with the approach described in the previous section: The first type is of cylindrical shape.Cylinder cells are widely used in photoacoustic spectroscopy.Its quality factors can be calculated analytically and we can use these cells to test our algorithm.The second type, the so called T cell, is depicted in Figure 5 and the geometrical dimensions are compiled in Table 1.The advantage of the T cell compared to the cylinder cell is that this design allows independent optimization of the key parameters affecting the signal strength.An analytical treatment of T cells does not appear to be feasible.The numerical calculations presented in this section have been performed with the finite element tool COMSOL Multiphysics2 [13].
The discussion of Section 4 shows, that the acoustic eigenmodes of the photoacoustic cell are needed.An investigation of the modes of the T cell has been published elsewhere [9].Here, we include only Figure 6 for demonstration purposes.Figure 5 Interior of a T cell consisting of an absorption cylinder (large diameter) and the perpendicularly mounted resonance cylinder (small diameter).The resonance cylinder carries the microphone at its upper end.For some technical reason there is a diminution at the connection of the two cylinders.The laser beam is along the axis of the absorption cylinder.Upon determination and normalization of the eigenmodes, it is possible to calculate the different components of the quality factor according to section 4. For cylinder cells the associated surface integrals can be evaluated analytically.The agreement between the exact and the finite element results was found to be excellent [10].
For the numerical modeling of sound waves with FEM the element diameter is required to be smaller than h = 0.2c/ f, where c is the speed of sound and f the frequency of the sound wave.In our case the largest frequency of relevance is 3500 Hz.When calculating amplitudes, not only the sound waves but also the heat generating laser beam has to be resolved in the finite element mesh.Here, element diameters of h = 0.05c/ f result in a satisfactory numerical accuracy.Orders higher than eight have no considerable contribution to the simulation.Therefore, we truncated the sum ( 14) after the 8th term.
We have extracted resonance frequencies, quality factors and amplitudes by fitting Lorentzian profiles (24) to the data.The quality factor is related to the half width (FWHM) of the resonance by ∆ f j = f j /Q j .Two sets of data, together with the fit functions are depicted in Figure 7.
In Figure 8, values of the quality factors for resonance cylinders of length 1 cm to 14 cm are depicted.The obvious parallelism of the two curves demonstrates the suitability of the model for the description of sound propagation in photoacoustic measuring cells.For very short resonance cylinders only a deviation of the parallelism is observable.The fact, that experimental values are consistently smaller than finite element values, is understandable because the FE simulation does not include all effects that contribute to loss of the photoacoustic signal.For example, the microphone and its surrounding are not modeled in detail.The sharp edges and non-sound hard parts are not considered.Also, we did not model the gas in-and outlet of the cell, which also add sharp edges.There might be impurities in the gas, that we do not have controll of.Also, non negligible coupling of the sound wave with the gas container might occur.These and other origins of loss are not included in the finite element model.Therefore, we expect higher loss, i.e. smaller quality factors, in the experiment than in the simulation.This is exactly what we observe.
Furthermore, the amplitudes according to Equation (18) have been calculated.In order to evaluate the volume integral, we use the following description of the radiation intensity of the laser beam: (25) We assume that the x-axis coincides with the symmetry axis of the absorption cylinder and is the distance perpendicular to this axis.w is the beam radius of approximately 2 mm.The analytical calculation of the amplitudes according to Equation (18) with the Gaussian laser beam profile (25) is not even possible for the simple case of cylinder cells and numerical methods are necessary.
Since we compare our finite element results to the phase-sensitively amplified photoacoustic signals of T cells, absolute pressure values are not relevant.Therefore, the value for the product αI 0 can be used to rescale the numerical data.We performed a least square fit to adjust the parameter αI 0 .Results are displayed in Figure 9.

CONCLUSION
We have performed a calculation of the signal of a photoacoustic sensor of unconventional shape.The associated coupled field problem could be avoided by taking advantage of an approach based on the Helmholtz equation and the incorporation of loss through quality factors.The calculations have been carried out with a commercial finite element tool.The achieved results compare very well with experiment.

Figure 2 Figure 3
Figure 2 Sound waves are accompanied by temperature waves.Heat flows from regions of elevated temperature to regions of lower temperature.

Figure 4
Figure 4 Boundary loss of temperature wave due to heat conduction (left) and velocity profile in the vicinity of Prandtl's boundary layer (right).

Table 1 Figure 6
Figure 6 Laser beam and two excited modes.

Figure 7 Figure 8 QFigure 9
Figure 7Experimental response function (dashed line) and corresponding fit (solid line) for the L R = 40 mm (upper) and L R = 80 mm (lower) T cell.