Increasing the sensitivity of electrical impedance to piezoelectric material parameters with non-uniform electrical excitation

To increase the robustness and functionality of piezoceramic ultrasonic sensors, e.g. for flow, material concentration or non-destructive testing, their development is often supported by computer simulations. The results of such finite-element-based simulations are dependent on correct simulation parameters, especially the material data set of the modelled piezoceramic. In recent years several well-known methods for estimation of such parameters have been developed that require knowledge of the sensitivity of a measured behaviour of the material with respect to the parameter set. One such measurable quantity is the electrical impedance of the ceramic. Previous studies for radially symmetric sensors with holohedral electrode setups have shown that the impedance shows little or no sensitivity to certain parameters and simulations reflect this behaviour making parameter estimation difficult. In this paper we have used simulations with special ring-shaped electrode geometry and non-uniform electrical excitation in order to find electrode geometries, with which the computed impedance displays a higher sensitivity to the changes in the parameter set. We find that many such electrode geometries exist in simulations and formulate an optimisation problem to find the local maxima of the sensitivities. Such configurations can be used to conduct experiments and solve the parameter estimation problem more efficiently.


Motivation
Piezoelectric effect is the physical phenomenon discovered by Pierre and Jacques Curie in 1880 that is exhibited by several crystalline and synthetic ceramic materials.When voltage is applied across certain surfaces of the solid, it exhibits mechanical strain; conversely, when mechanical stresses are applied, voltage is produced between its surfaces.
Acoustic transducers are required to construct an ultrasonic measurement device of any kind.Transducers based on circular piezoelectric ceramic disks have been manufactured and used in a wide variety of applications.
Ceramics are composed of a large number of randomly oriented crystals.The cumulative properties of each of these crystals together determine the properties of the whole solid.This means that each batch of the ceramic produced will have somewhat different material properties.Additionally, the ma-terial coefficients are dependent on geometry.If a such a ceramic disk with a specific thickness and diameter is sintered and polarised, the resulting material coefficients are different than a disk of the same material and of different thickness and diameter.Hence, post processing and treatment of the material also influences the material coefficients.Pérez et al. (2010) found the elastic coefficients to vary up to 5 % and the piezoelectric and dielectric coefficients to vary up to 20 %.These variations include measurement uncertainties concerning the determination of the coefficients and also the batch-to-batch and geometry dependencies.
Depending on the application at hand, the full knowledge of the material data set can be crucial to the whole process: for a circular disk with full width electrodes on the top and bottom (see Fig. 1) there is no shear movement; hence, e 15 is of no importance.Unlike such a disk, in the case of a shear wave transducer the knowledge of e 15 is however crucial.(Lahmer, 2008).
Analytical and numerical modelling of piezoelectricity and the resulting computer simulations (Kocbach, 2000;Unverzagt et al., 2013) have enabled more efficient design and construction of transducers in recent years.
Such numerical modelling requires the precise knowledge of the parameters of the material under analysis.To this end, new methods that utilise the mathematical concept of inverse problems to estimate the material parameters have been designed (Kaltenbacher et al., 2006(Kaltenbacher et al., , 2008)).
Previous studies (Rautenberg et al., 2011) have shown that configurations inducing resonance in radial or thickness directions (Lahmer, 2008) are insufficient for parameter estimation, when only the electrical impedance is measured, because the impedance characteristic in the frequency space shows little or no sensitivity to certain material parameters.Especially critical are c 44 , ε 11 and e 15 since in the experiments of Rautenberg et al. (2011) they show no sensitivity to material parameters at all.Low or no sensitivity means that these parameters cannot be estimated using inverse problems from the measurement of impedance characteristics only.Other, more involved measurements have been proposed in Rupitsch et al. (2009); however, they are also much more cost-intensive.
In order to improve the sensitivity of the measured impedance characteristic to the material parameters, we designed a three-electrode configuration with electrodes of various radii as shown in Fig. 2a and b and use simulations in which a non-uniform electrical excitation is applied to the ceramic in order to compute the impedance characteristic and its sensitivity to the material parameters.
In the following sections we first give a short overview of the equations governing our simulation and the excitation of the ceramic.We then define and analyse the sensitivity of the computed impedance to parameters.Finally, we formulate an optimisation problem to find a locally optimal electrode geometry that maximises the sensitivity and show that many such local optima occur.

Modelling piezoelectricity
In this paper we restrict ourselves to linear piezoelectric effect IEEE Std 176-1987.The equations of linear piezoelectricity in tensor form are given below and form the basis for a finite-element formulation.
where -D is the electrical flux density vector -E is the electrical field vector σ is the mechanical stress tensor -S is the mechanical strain tensor -c is elastic modulus tensor -e is the piezoelectric coupling tensor ε is the electrical permittivity tensor.
We shall also restrict ourselves to thin circular disks.Therefore, the use of a cylindrical coordinate system is appropriate.We shall also assume that the configuration is rotationally symmetric.Using the notation of Helnwein (2001) and Lahmer (2008)  This results in the following system of material equations from Eqs. ( 1) and ( 2): The above material equation can be extended to a full set of partial differential equations in time and space using Newton's, Gauss' and Faraday's Laws (Meschede and Gerthsen, 2010, Chap. 4 and Chap. 7).
Newton's law of motion (Slaughter, 2002, Cauchy-Navier equation) for the mechanical behaviour is where B is the differential operator relating mechanical strain to mechanical displacement.
As ceramic materials are insulators, there is no free charge and Gauss' (flux) law states Neglecting the insignificant changes in magnetic field in this case, Faraday's law states that the electric field is the negative gradient of the electric potential.
Additionally, we consider a Rayleigh damping model with positive constant α and β for the energy dissipation and arrive at a system of four partial differential equations in time and space from Eqs. (3), (4), ( 5) and ( 6): These equations are sufficient for transient analysis.The electrodes and their excitation occur as boundary conditions, in terms of free charge at the part of the boundary containing the electrode L with n being the normal vector n = (n r , n z ) at the boundary, or the applied potential φ L (t) at the electrode L. We call the remaining boundary r = L .The above equations may also be transformed for harmonic analysis into the frequency domain.Application of a Fourier transform changes the unknowns u and φ in Eqs. ( 7) and (8) to time harmonic complex variables û and φ, and the Rayleigh coefficients are scaled according to frequency.
The resulting time harmonic partial differential equations are given as follows: with ı being the imaginary unit.This transformation however contains a systematic error as α 0 and β 0 are not constants in practice but change with the central frequency of the Fourier transform.
Weak formulations and finite-element methods for the solution of the above system have been studied in Kocbach (2000) and Lahmer (2008).

Electrical excitation
In order to simulate the behaviour of a piezoelectric transducer with finite elements, it is excited by a delta function of the charge or a pulse of potential.This gives the boundary conditions to solve Eqs. ( 7) and (8).
Electrical current flow at the electrode is the time derivative of the free charge and the impedance between any two electrodes is then defined as the quotient of the potential difference and the current.In the simple case of two electrodes, with a delta function pulse of charge into one electrode with amplitude q L 0 , q L (t) = q L 0 δ(t) (see Fig. 1), this leads to the straightforward calculation where ω k is some frequency of interest, qL 0 denotes the frequency domain charge peak and φL (ω k ) is the calculated response of the potential difference using Eqs.( 7) and ( 8).
In our setup two different electric potentials are applied on two electrodes and the third is grounded.An external circuit (Fig. 2c) is required to achieve this potential difference at the electrodes.The values of the resistors R 0 and R 2 and capacitor C 2 are chosen to maintain V 1 ≥ 2V 2 .This is done by solving the equations of Kirchhoff's laws for the circuit (Meschede and Gerthsen, 2010, Chap. 7).As the circuit representation of the ceramic includes three unknown impedances, the solution of the Kirchhoff laws equations requires three finite-element solutions for Eqs. ( 7) and ( 8) using an initial guess for the material parameters and a given electrode geometry (Unverzagt et al., 2015).The chosen values for the external circuit are then kept constant for that particular electrode geometry, while we compute the sensitivity of the total impedance to material parameters and optimise it as discussed in Sect. 5.

Sensitivity
For each frequency of interest in the domain, the solution of Eqs. ( 7) and ( 8) along with the external circuit results in a total impedance, which can be measured in experiments and compared with the numerical computations in order to estimate the material parameters (Rupitsch and Lerch, 2009;Kaltenbacher et al., 2008).The first step to solve this inverse problem requires the determination of a sensitivity of the numerical solution to variations in the material parameters.This is in fact the derivative of the impedance with respect to the material parameters.
For the numerical solution we use the commercial software package CAPA, which computes the impedance for a given configuration.CAPA is a simulation tool for the numerical solution of electromechanical, coupled field problems and is therefore suitable for the analysis of most mechatronic sensors and actuators such as electromagnetic loudspeakers or piezoelectric transducers.
Besides transient analysis, which is used in this contribution, the harmonic behaviour of the piezoceramic disc can also be calculated.One disadvantage of the transient simulation method is the sole use of the Rayleigh damping model for the energy dissipation processes as mentioned in Sect. 2. This is a rather simple approximation of the damping behaviour in practice.Another disadvantage of precompiled solver packages in general is inflexibility; it is impossible to modify the computation in any way in order to be able to compute more information, like sensitivities.
In contrast to previous studies we look at the real and imaginary parts of the complex impedance separately instead of looking at the magnitude.The formulation in terms of magnitude and phase, although computationally efficient, hides some geometrical structure, which is apparent when looking at the real and imaginary components separately.This is analogous to a polar coordinate system versus a cartesian coordinate system.Figure 3a shows the complex impedance as a function of the frequency in the complex plane, whereas Fig. 3b shows the magnitude of the impedance.The 3-D plot shows more structure and we use the same representation for the sensitivities.
We used finite-difference approximations for the configuration in Fig. 8 for the sensitivity of the impedance w.r.t. the material parameters.Figure 4 shows the sensitivity of the impedance in the frequency domain w.r.t.c 44 , e 15 and ε 11 (occurring in Eq. 3) that are known to show low sensitivity in simpler configurations (Rautenberg et al., 2011).
Using a slightly different geometry (see resulting geometry in Table 1) for the electrodes but the same excitation, the sensitivity curves change both qualitatively and quantitatively (see Fig. 6 and compare with Fig. 4).However, the change is not uniform across the parameters.Hence, there is a need to find the optimal electrode configuration to ensure maximal sensitivity of the impedance to the material parameters.

Optimisation
We start by formulating a constrained minimisation problem with the electrode radii as the variables which will then be solved.
The electrode configuration is parametrised using four ring radii r = (r 1 , . .., r N r = r 4 ) and considering the outer radius as constant (see Fig. 5).Let x = (x 1 , . .., x N x ) denote the various material parameters occurring in Eq. ( 3) and let F ⊆ R denote the frequency domain.The impedance Z(f ; x, r) can be considered as a function Z : F × R N x × R N r → C. For fixed radii parametrisation r and material parameters x the impedance is a function Z(f ; x, r) : F → C which maps the argument f −→ Z(f ; x, r).This complex-valued function Z is then transformed into a two-dimensional real function z(•; x, r) : F → R 2 via the usual Euclidean mapping.We approximate the partial derivative of z towards a change of material parameter x i using a finite-differences scheme ) and for small h > 0. Hence, the sensitivity of z towards one specific parameter x i while remaining in a constant electrode configuration r can be considered as ∇ x i z(•; x, r) L 2 (F ) .We consider the norm of the corresponding function space L 2 (F) to measure the sensitivity    However, the change is not uniform across the parameters.Hence there is a need to find the optimal electrode configuration to ensure maximal sensitivity of the impedance to the material parameters. .
Discretising the frequency domain F equidistantly into intervals . Now, we wish to maximise the sensitivities with respect to all material parameters x i .
We have noticed that some parameters are very reactive towards change in geometry and others are not.Due to unevenly distributed sensitivity of the material parameters x towards optimisation, their different orders of magnitudes and the circumstance that the optimisation method used may only minimise an objective function, we introduce a weight matrix W := diag(w 1 , . .., w N x ) ∈ R N x ×N x + and a scaling matrix S := diag(s 1 , . .., s N x ) ∈ R N x ×N x depending on the orders of magnitude of the initial sensitivity evaluation and reformulate the sensitivity approximation as the minimisation problem: The scaling is only used inside the optimiser so that all components of the gradient have a similar scale and has no meaning outside the optimisation routine.For comparison of the results one needs to compare the unscaled objective function with S = I , the identity matrix.

Constraints
In the production process of these piezoceramics, a laser cuts out the electrode rings from a metal plate lying on top of the ceramic.The laser has a width of 0.3 mm; hence, some restrictions to the geometry of the electrodes apply.The minimal distance of two adjacent electrode rings is at least 0.3 mm; also, all radii must be positive and smaller than the constant outer ring radius.These constraints are all linear in the radii and can be reformulated into a vector inequality Ar ≤ b with   2008) can be used with only the measurements of the impedance required.This ces the cost of such investigations as the equipment required is comparatively cheap.order to systematically search for the maximal sensitivity in such a configuration we need to an optimisation problem with the configuration radii as variables.Due to the inflexibility of the ompiled solver we were forced to use a derivative free optimisation algorithm.More efficient rithms may be used if a more flexible finite element code with the possibility of influencing 14 .
However the results show clearly that there are many locally optimal configurations.

295
We are investigating the possibility of better sensitivity analysis by utilising the simulation software CFS++ (Kaltenbacher, 2010), being developed at the TU Vienna.Modifications in the software for this purpose are ongoing and future work.Besides sensitivities these changes can be used to compute adjoints and thus solve optimisation problems, including parameter estimation problems.
Simulations with CFS++ will also shed light on the dampening influence of the external circuit if 300 we increase the number of electrodes, thereby increasing the complexity of the external circuit, in particular the number of external impedances.
Apart from discrete sensitivities and adjoints it is also possible to formulate the sensitivity and adjoint equations for the model in function spaces and discretise these along with the primal equations and solving them.Another avenue for future development is to formulate the sensitivity maximi-

Optimisation method
Since the sensitivity itself is a finite-difference approximation, it makes little sense to use an optimisation method that requires further derivatives.For the optimisation we used Powell's latest derivative-free trust region optimiser LIN-COA (LINearly Constrained Optimization Algorithm) (Powell, 2014a, b) for linearly constrained problems.
The LINCOA method is a derivative-free optimisation algorithm for linearly constrained problems written in Fortran by M.J.D. Powell.It is based on Powell's other derivativefree optimisation algorithms with the distinction of incorporating general linear constraints.However, a detailed description of the LINCOA software has not been published yet.
According to Powell (2014a), the LINCOA method uses a quadratic model www.j-sens-sens-syst.net/4/217/2015/sation problem using shape and topology calculus instead of parametrised rings.This would help generalise the configuration of electrodes to ceramic geometries that are not radially symmetric.
available for future development.

16
as an approximation to the objective function J (x) ∈ R .However, as no derivatives are available, the quadratic model Q has to be iteratively constructed from function evaluations of J .At the beginning of the optimisation, the user chooses the amount of points m which are further used to interpolate the model Q with m ∈ {n + 2, . .., 1 2 (n + 1)(n + 2)}, a typical choice for m being 2n + 1.This leaves 1 2 (n + 1)(n + 2) − m degrees of freedom for the choice of Q which are fixed by applying a symmetric Broyden updating method to the model Q.A trust region and an active-set approach incorporating the linear constraints then determine the new points for updating the model function Q, which is iteratively minimised by the process.The trust region size in the process is decreased when certain conditions are fulfilled until ultimately the algorithm halts as the trust region size has reached a userprescribed lower boundary.

Results
In the following section we will present results of the sensitivity optimisation.The section is divided into two parts.In the first part we will concentrate on the influence of the weight matrix W on the optimisation.We have tested different weight assignments (uniform, binary and mixed) for W

Influence of the weighting matrix
As an initial point for a first optimisation we used the electrode configuration of a piezoceramic (see Fig. 8) we physically possess for measurement purposes.This ceramic configuration has been developed in Unverzagt et al. (2015) using statistical methods.We have made the experience that this configuration has exceptionally good properties with regard to optimisation as opposed to other geometries tested.We shall call this configuration the reference initial point.For the optimisation process we used the uniform weights W = diag(1, . .., 1).A summary of the results can be found in Table 1 along with the final geometry.This shows an overall in sensitivity ∇Z 2 by more than 41 % and some partial sensitivities were increased by up to 307 %.However, not all partial sensitivities are very reactive towards changes in the electrode configuration, i.e. the sensitivity gain regarding e 15 is −6.1 %.
As a reaction to this slight partial decrease in sensitivity regarding e 15 we chose to neglect the partial sensitivities of those parameters which have a positive gain and focus on e 15 by setting the corresponding weights to 0 and 1, respectively.Through this binary weight-matrix setting, the sensitivity with regard to e 15 was increased by 31 %.(see Table 2) However, following our expectations, the sensitivities regarding the other parameters have mainly decreased.The evaluation histories can be seen in Fig. 7.
These two resulting geometries combined demonstrate the feasibility of optimising the sensitivity with regard to all parameters, i.e. the two geometries show a combined increase in sensitivity for all parameters.
We have experimented with different mixed weights (see Table 3) with the aim to find a single globally optimal electrode configuration, however, it has not been possible to find such an electrode configuration.Although it is not globally optimal in the sense discussed above, it is still possible to use the configuration from Table 3 for the purpose of solving an inverse problem in the material parameters since each parameter shows some (non-zero) sensitivity.Due to long computational time for each optimisation process we did not try to fine-tune the selection of the weight matrix.

Multiple starting points
To examine the influence of different initial electrode configurations to the optimisation process we have chosen a wide range of barely feasible configurations (see Tables 4a-c).All these are optimised using uniform weights and the resulting optimal is compared to the resulting optimal when starting from the reference initial point above with uniform weights.The resulting geometries are all further away from infeasibility than their initial configuration and locally optimal.These three cases demonstrate that the problem of identifying a single globally optimal geometry is hard, since there are many locally optimal configurations.

Conclusions and future work
In contrast to Rautenberg et al. (2011) with the two electrode approach we have shown that the sensitivity of the impedance to various critical material parameters is non-zero in all our ring electrode configurations with non-uniform excitation.Therefore, parameter estimation techniques as in Kaltenbacher et al. (2008) can be used with only the measurements of the impedance required.This reduces the cost of such investigations as the equipment required is comparatively cheap.
In order to systematically search for the maximal sensitivity in such a configuration we need to solve an optimisation problem with the configuration radii as variables.Due to the inflexibility of the precompiled solver we were forced to use a derivative-free optimisation algorithm.More efficient algorithms may be used if a more flexible finite-element code with the possibility of influencing the internal computations were available, so that one could compute derivatives simultaneously.However the results show clearly that there are many locally optimal configurations.
We are investigating the possibility of better sensitivity analysis by utilising the simulation software CFS++ (Kaltenbacher, 2010) being developed at the TU Vienna.Modifications in the software for this purpose are ongoing.Besides sensitivities, these changes can be used to compute adjoints and thus solve optimisation problems, including parameter estimation problems.Simulations with CFS++ will also shed light on the dampening influence of the external circuit if we increase the number of electrodes, thereby increasing the complexity of the external circuit, in particular the number of external impedances.
Apart from discrete sensitivities and adjoints, it is also possible to formulate the sensitivity and adjoint equations for the model in function spaces and discretise these along with the primal equations and solving them.Another avenue for future development is to formulate the sensitivity maximisation problem using shape and topology calculus instead of parametrised rings.This would help generalise the configuration of electrodes to ceramic geometries that are not radially symmetric.

Figure 2 .
Figure 2. (a) Arrangement of electrodes in thin rings; (b) arrangement of electrodes in wide rings; (c) schematic diagram of the whole circuit with the piezoceramic shown as a T-Network.
(Voigt notation)  in the rotationally symmetric case we can rewrite the tensor notation into matrix formulation:

Figure 4 .Figure 5 .
Figure 4. Sensitivity of the impedance for the geometry of the ring electrodes shown in Fig. 8 to various material parameters against frequency w.r.t.(a) c 44 , (b) e 15 , and (c) ε 11 .

Figure 6 .
Figure 6.Sensitivity of the impedance for the resulting geometry of the ring electrodes shown in Table 1 to various material parameters against frequency w.r.t.(a) c 44 , (b) e 15 , and (c) ε 11 .
Uniform weights: Improvement of sensitivity of impedance w.r.t material parameters using W = (1, . . ., 1): The sensitivity has increased for all the parameters with exception of e15.These are the best onent-wise results we have computed so far.Objective function values are shown unscaled.Figure shows ting geometry utilising uniform weight matrix in mm r1 = 4.01; r2 = 4.33; r3 = 2.20; r4 = 3.82; onclusions and future work ntrast to Rautenberg et al. (2011) with the two electrode approach we have shown that the itivity of the impedance to various critical material parameters is non-zero in all our ring elecconfigurations with non uniform excitation.Therefore parameter estimation techniques as in enbacher et al. ( Binary weights: Improvement of sensitivity of impedance w.r.t.material parameters using W = diag(0, 0, 0, 0, 0, 0, 0, 0, 0, 1): The sensitivity w.r.t e15 has increased.However, the sensitivity toward other parameters has mainly decreased.Objective function values are shown unscaled.Figure shows resulting geometry utilising binary weight matrix in mm r1 = 3.4; r2 = 3.7; r3 = 1.23; r4 = 3.54; r) 2 w.r.t.Ar ≤ b.

Table 1 .
Uniform weights: improvement of sensitivity of impedance w.r.t material parameters using W = diag(1, . .., 1).The sensitivity has increased for all the parameters with exception of e 15 .These are the best component-wise results we have computed so far.Objective function values are shown unscaled.Figure shows resulting geometry utilising uniform weight matrix in millimetres; r 1 = 4.01; r 2 = 4.33; r 3 = 2.20; r 4 = 3.82.
The sensitivity w.r.t e 15 has increased.However, the sensitivity toward other parameters has mainly decreased.Objective function values are shown unscaled.Figure shows resulting geometry utilising binary weight matrix in millimetres; r 1 = 3.4; r 2 = 3.7; r 3 = 1.23; r 4 = 3.54.

Table 3 .
Mixed weights: improvement of sensitivity of impedance w.r.t.material parameters using W = 1 With exception of c 44 and e 15 all partial sensitivities are increased.Objective function values are shown unscaled.Figure shows resulting geometry utilising mixed weight matrix in millimetres; r 1 = 4.01; r 2 = 4.34; r 3 = 2.17; r 4 = 3.93.