Numerical Shape Optimization of Photoacoustic Sample Cells: FirstResultsBernd Baumann∗,1, Bernd Kost1, Marcus Wolff1,2, Hinrich Groninga2, Tanja Bloß1 and StefanKnickrehm1

1Hamburg University of Applied Sciences, 2PAS-Tech GmbH∗Corresponding author: Prof. Dr. B. Baumann, HAW Hamburg, Faculty of Engineering and Computer Science,Berliner Tor 21, 20099 Hamburg, Germany, bernd.baumann@haw-hamburg.de

Abstract: First results of an automatic shape optimiza-tion of a photoacoustic sample cell are described. The aimis to maximize the sensor’s signal strength. The approachconsiders all shapes that can be represented by a numberof axis-symmetrical truncated cones which are connectedin a continuous way. In addition, the cell is subjected tocertain constraints (e.g. laser beam is not blocked duringpassage through cell). The acoustic pressure at the micro-phone represents the objective function and is calculatedusing an eigenmode expansion combined with a finite ele-ment analysis. The performance of different optimizationmethods is investigated.Keywords: Photoacoustics, gas sensor, acoustic res-onator, shape optimization

1 IntroductionThe monitoring of trace gases has become an importanttask in many fields. In situ environmental gas detectionand the diagnosis of diseases by breath tests are promi-nent examples. Photoacoustic sensors have a long andsuccessful history concerning the detection of very lowgas concentrations.

The mode of operation of a photoacoustic sensor is sim-ple and can be understood with the aid of Figure 1. Thephotoacoustic cell contains the gas sample. It is passed bya laser beam. The wavelength of the laser is adjusted toone of the absorption lines of the molecules one intendsto detect. If molecules of this type are contained in thesample absorption takes place and the gas heats up due torelaxation through molecular collisions.

If the laser beam is modulated by a mechanical or elec-tronical chopper, the resulting local and transient heatingresults in the generation of thermal and sound waves. Thelatter can be detected by a microphone. The chopper fre-quency and the frequency of the waves are identical. It is a

natural idea to amplify the microphone signal by adjustingthis frequency to an acoustic resonance of the measuringcell: The photoacoustic sensor is said to be operated inthe ’resonant mode’.

For the detection of low molecule concentrations it isimportant to obtain a large signal amplification in theacoustic resonator. It is evident, that the volume andshape of the resonator influence the signal strength. Inthis paper we present first results of an automated opti-mal shape search with respect to the signal strength ofphotoacoustic cells of axial symmetry. The search spaceconsists of all shapes that can be represented by 6 axis-symmetrical truncated cones which are connected in acontinuous way. A typical configuration is depicted inFigure 2. To parametrize these cell shapes 7 radii areneeded. These radii Ri, i = 1, . . . , 7 represent the de-sign variables of this optimization problem.

For the numerical optimization three different searchstrategies have been applied: Sequential Quadratic Pro-gramming (SQP) [1, 2], Genetic Algorithm (GA) [3] andEvolution Strategy (ES) [4, 5].

2 Calculating the Objective Func-tion

Various types of photoacoustic cells have been investi-gated and applied to the detection of trace gases [6]. One

Figure 1: Design of a photoacoustic sensor with cylindrical res-onator cell.

1

Figure 2: Exemplary resonator cell in the search space. Each ofthe truncated cones is of length 12.5 mm.

of the most important cell types are cylinder shaped cells.Even for this simple geometry it is not possible to ana-lytically calculate all quantities which are required for thedetermination of the photoacoustic signal. For cells ofmore complex shape, an analytical calculation is impossi-ble and numerical methods have to be applied. Recently,a complete numerical modeling of the photoacoustic sig-nal generation process as well as the calculation of signalstrength has been achieved [7, 8, 9, 10]. The calculationof the signal strength is a prerequisite for an automatedoptimization, since enhancing the signal strength is theobvious objective for photoacoustic sensors. In this sec-tion we give a short overview of the method used for thecalculation of the signal strength. Details can be found inthe preceding references.

The signal strength at point ~r is determined through thestrength of the excitation and the amount of loss. It can beexpressed as a sum over all normalized resonator modespj(~r) [11]:

p(~r, ω) =∞∑

j=0

Aj(ω)pj(~r). (1)

Sound hard walls of the cell are assumed. The amplitudescan be calculated from

Aj(ω) = iAjω

ω2 − ω2j + iωωj/Qj

(2)

with the eigenfrequencies ωj and the corresponding qual-ity factors Qj . The latter account for various loss mecha-nisms (see below).

The excitation of the sound waves is described by

Aj =α(γ − 1)

VC

∫VC

p∗j · I dV . (3)

α denotes the absorption coefficient and I = I(~r) the spa-tial distribution of the laser intensity within the resonator.γ denotes the ratio of the specific heat at constant pres-sure cp to the specific heat at constant volume cV . VC isthe cell volume.

Various mechanisms contribute to loss (for details con-sult [9]): The Stokes-Kirchhoff loss due to viscosity andthermal conduction in the gas volume, thermal conduc-tivity surface loss and surface loss due to viscosity. Thecombined effect of all loss mechanisms leads to a result-ing quality factor which is calculated from

1Qj

=∑

i

1Qi

j

. (4)

The objective function is determined by searching forthe maximum value of |p(~rM, ω)| over the frequency rangefmin < ω/2π < fmax (step size 1 Hz). The vector ~rM de-notes the position of the microphone. The investigationsdescribed here have been performed for a microphone po-sition located on the cell surface at half cylinder length(see Figure 1).

Certain assumptions have to be made concerning theintensity of the laser beam (see Equation (3)). The beamprofile has been modeled according to

I(~r) = I0 exp

[−2

(r⊥rB

)2]

. (5)

We assume that the x-axis coincides with the symmetryaxis of the resonator and r⊥ =

√y2 + z2 is the distance

perpendicular to this axis. rB is the beam radius which isassumed to be approximately 2 mm. In trace gas analysisthe absorption of the radiation field is very small. There-fore, the x-dependency of I(~r) has been neglected. Ab-solute values of the signal are not important in the presentcontext and the product αI0 can be chosen arbitrarily.

3 Shape Optimization of Photo-acoustic Resonator

For the simulations described in this paper the gas pa-rameters of n-butane under atmospheric conditions have

2

been used (Table 1). As has been reported in a previouspaper it is necessary to use a maximum element size of0.05c/fmax instead of the 0.2c/fmax-rule usually appliedin acoustic simulations [9]. The reason is that not only thesound waves, but also the laser beam has to be resolvedproperly by the finite element mesh.

Typical cylinder-shaped photoacoustic cells possess adiameter to length ratio much smaller than one. It iscertainly a reasonable idea to start the shape optimiza-tion from a cylinder cell with this property. The eigen-modes and the eigenfrequencies of cylinder cells can beobtained by analytical means [12]. The modes can beclassified as longitudinal, azimuthal, and radial. In ourcase, where a laser beam travels along the symmetry axisof the resonator, longitudinal and azimuthal modes cannotbe excited, since I and pj in Equation (3) are orthogonal.Therefore, longitudinal and azimuthal modes do not con-tribute to the sum in Equation (1). Due to the small di-ameter of the resonator standing waves in radial directionhave small wavelengths and large frequencies. Since thephotoacoustic signal is inversely proportional to the mod-ulation frequency the contribution of the radial modes tothe signal is small . Therefore, the only mode, which canbe excited efficiently is the trivial mode correspondingto eigenfrequency 0. This mode corresponds to a space-independent pressure variation, i. e. p0(~r, ω) = p0(ω). Ifthe cell is operated at low frequencies one speaks of the’nonresonant mode of operation’. In this mode large sig-nals can be obtained, however its usefulness is restricted,because noise also increases approximately according to1/f . These facts have to be kept in mind for the interpre-tation of the result of our first attempt of a shape optimiza-tion.

The code used for the calculation of the objective func-tion is a Matlab-script, which calls COMSOL Multi-physics commands. The easiest way to implement an op-timization is to use the Matlab standard optimization toolfmincon [13]. This tool is based on the SQP-optimizationmethod. The radii have been restricted to 2 mm ≤ Ri ≤13 mm. The lower bound has been chosen equal to thelaser beam radius rB. This guarantees that the main partof the laser intensity travels through the resonator cavity.The upper bound has been introduced, since it turned outthat without an upper bound on the Ri very large cellsare generated in the course of the simulation and the cor-responding large number of degrees of freedom is very

demanding in respect of computer time and memory.For the evaluation of Eq. (1) it has to be examined,

how many eigenmodes contribute significantly to the re-sponse function. The simulations have been carried outwith four modes. In order to check the consistency ofthe achieved results some calculations have been repeatedwith 8 modes without relevant differences.

Furthermore, the calculation of the response functionEq. (1) is based on the frequency step size ∆f = 1 Hz inthe frequency range fmin ≤ f ≤ fmax, with fmax =3500 Hz. For the same reason as described abovesome simulations have been checked using a maximal fre-quency of 7000 Hz to ensure that modes with frequencies> fmax have no significant influence on the signal.

When starting the optimization from a cylindrical cellwith a radius in the allowed range and fmin near zerothe resulting optimal cell is always a cylinder of radius2 mm. In this case the frequency of maximal signalstrength is fmin (nonresonant mode of operation, seeabove). This is understandable because the signal strengthis inversely proportional to the cell volume and the opti-mization shrinks the cylinder to its minimal allowed size.So the optimization works properly but the result is notsatisfying because of the above mentioned noise problem.

In order to circumvent this problem, we performed sim-ulations where we restricted the search for the signal max-imum to frequencies above the influence of the trivialmode resonance, i. e. fmin has been set to 0.8× the low-est nontrivial eigenfrequency. Starting again from cylin-der shaped cells of various radii the optimal cell found isdepicted in Figure 3.

The important question now is if the found cell shapeconstitutes the global optimum or a local optimum only.To pursuit this question simulations using different initialcell shapes have been performed. An fmincon optimumsearch starting from the cell depicted in Figure 4 resultsin the cell, which is shown in Figure 5 and starting fromone of the two cells depicted in Figure 6 results in the so-lution of Figure 7. The shape of Figure 7 constitutes theresonator producing the largest signal found during thecourse of the investigation presented in this paper. With-out going into detail we would like to mention, that theoptimal shapes have been found in less than 200 functionevaluations. The quality of the obtained solutions can bejudged from Table 2. It is worth mentioning, that this isan improvement by a factor of about 1000 with respect to

3

density ρ 2.376 kg/m3

sound velocity c 210.59 m/sviscosity η 7.6 · 10−6 Pa scoefficient of heat conduction κ 1.62 · 10−2 W/m Kspecific heat capacity at constant volume cV 1.598 · 103 J/kg Kspecific heat capacity at constant pressure cp 1.735 · 103 J/kg K

Table 1: Gas parameters of n-butane under atmospheric conditions (T = 300 K, p = 1013 hPa) [14].

Figure 3: Best cell found by fmincon when starting from acylinder shaped cell (Ri = 2.00, 2.48, 2.76, 5.37, 2.68, 2.47,2.00 mm).

cylindrical resonators. The exact value depends naturallyon the radius of the cylinder.

Figure 4: Alternative initial cell shape. The radii of the trun-cated cone vary linearly from 2 to 12 mm.

As an SQP-method, fmincon is a derivative basedmethod and in danger to get stuck in local optima. Thepreceding paragraphs show that this is what happens inthe shape optimization of axis-symmetrical photoacousticresonators. This problem has been circumvented by the

Figure 5: Best cell found by fmincon when starting from theinitial configuration of Figure 4 (Ri = 2.00, 2.00, 2.65, 2.73,2.00, 13.0, 13.0 mm).

use of initial configurations of various shapes. However,as always, if the resonator of Figure 7 is the global opti-mum is uncertain.

What one would like to have is an optimization methodwhich does not need such heuristic user support. Suchmethods are available and are known as global optimiza-tion methods. The evolutionary algorithms like GA andES are of this type. Here results of preliminary investi-gations of the shape optimization problem with GA andES are presented. To give a first rough impression of the

0.192 0.598 1.000

Table 2: The relative performance of the found optima can bejudged from the signal ratio, i. e. the signal of the respective res-onator divided by the signal of the best resonator (the rightmostone).

4

Figure 6: Second alternative initial cell shapes. The radii vary alternatively between 2 and 12 mm.

Figure 7: Best cell found by fmincon when starting from oneof the initial configurations of Figure 6 (Ri = 13.0, 3.94, 2.00,2.00, 2.00, 3.95, 13.0 mm).

performance of these methods the results of four runs forboth methods using different random number sequencesand each comprising 100 function evaluations are pre-sented in Figures 8 and 9.

For the GA results (Figure 8) the Matlab GA-toolbox[13] with standard parameters has been used. All runshave been started from a cylinder with the maximal al-lowed radius 13 mm.

The ES results (Figure 9) have been obtained with asimple (3 + 5)-ES with general step size adaptation usinga modified barrier function approach to include the con-straints [2]. As usual in ES, randomly generated vectorsin the allowed range are used as start configurations.

Figure 8 shows, that despite the cylindrical initial con-figuration and the low number of function evaluations GAis en route to the best solution found by the SQP method.ES as well shows a clear tendency towards the best SQP

solution. However, in the second run it seems, that ESis in the vicinity of the local optimum depicted in Figure3. It should be mentioned, that in contrast to the commer-cial GA-tool with its tuned strategy parameters, for the ESmethod an investigation of suitable search parameters hasstill to be done.

We expect, that a larger number of function evaluationswould lead to the solution of Figure 7. Alternatively, apromising optimization method would emerge by switch-ing from GA or ES to SQP after a certain number of func-tion evaluations (hybrid method).

4 Summary and Conclusions

In the present article we have described the first steps to-wards an automatic shape optimization of the resonatorcell of axially symmetric photoacoustic sensors. The SQPmethod is not able to find the global optimum when theinitial cell is of cylindrical shape. In contrast, GA andES are more flexible with respect to the initial configura-tion. We propose, that it might be an efficient strategy forthe practical design of photoacoustic resonators to startthe optimization with GA or ES and later switch to SQP,which is supposed to have its strength in the fine tuning.

The results of this investigation give hints, how efficientaxially symmetric resonators should be shaped (Figure 7).It should, however, be kept in mind, that real photoacous-tic cells have gas inlets, gas outlets, boreholes for the mi-crophone etc. which affect the signal. Such effects havenot been included in the investigations presented here.

5

Figure 8: Variation of radius along symmetry axis of the bestshapes found by GA within 100 function evaluations. For a goodcomparability the best shape found by the SQP method has beenincluded.

Figure 9: Variation of radius along symmetry axis of the bestshapes found by ES within 100 function evaluations. For a goodcomparability the best shape found by the SQP method has beenincluded.

References

[1] R. Fletcher, Practical Methods of Optimization,John Wiley and Sons (1987)

[2] J. Nocedal, S. Wright, Numerical Optimization,Springer, Berlin (2006)

[3] J. H. Holland, Adaptation in Natural and Artifi-cial Systems, University of Michigan Press, AnnArbor (1975)

[4] I. Rechenberg, Evolutionsstrategie - Opti-mierung technischer Systeme nach Prinzip-ien der biologischen Evolution, Frommann-Holzboog, Stuttgart (1973)

[5] H.-P. Schwefel, Evolution and Optimum Seek-ing, Wiley and Sons, New York (1995)

[6] A. Miklos, P. Hess and Z. Bozoki, Applicationof acoustic resonators in photoacoustic trace gasanalysis and metrology, Rev. Sci. Instr., Vol. 2,No. 4, 1937-1955 (2001)

[7] B. Baumann, M. Wolff, B. Kost, and H.Groninga, Calculation of Quality Factors andAmplitudes of Photoacoustic Resonators, Proc.Comsol Conf. Frankfurt 2006, 134-138 (2006)

[8] B. Baumann, B. Kost, H. G. Groninga, M. Wolff,Eigenmodes of Photoacoustic Sensors via Fi-nite Element Method, Review of Scientific Instru-ments, 77, 044901, (2006)

[9] B. Baumann, M. Wolff, B. Kost, and H.Groninga, Finite element calculation of photoa-coustic signals, Appl. Opt. 46, 1120-1125 (2007)

[10] B. Baumann, M. Wolff, B. Kost, and H.Groninga, Solving a coupled field problem byeigenmode expansion and finite element method,Int. Jnl. of Multiphysics, Vol. 1, Nr. 3, 303-315(2007)

[11] L. B. Kreuzer, The Physics of Signal Generationand Detection, Optoacoustic Spectroscopy andDetection, Pao, Y.-H. (Ed.), Academic, London,1-25 (1977)

[12] S. Temkin, Elements of Acoustic, John Wiley andSons Inc., New York (1981)

[13] www.mathworks.com

[14] VDI-Warmeatlas, 9. Auflage, Springer Verlag,Berlin (2002)

6