Multi-Phase Multi-Physics Finite Element Model Updating of Piezoelectric Transducer

The overall objective of this dissertation is to propose a quick and accurate method of updating transducer models. When one makes measurements on a piezoelectric transducer, oftentimes only the impedance function is measured as the experimental data. Thus, to perform model updating for transducers, two major tasks will be covered: (i) developing and verifying an efficient method for estimating the electric impedance function of a transducer, and (ii) developing and testing a FE model updating method for piezoelectric transducers. The proposed method to estimate the impedance function of a transducer is a Laplace domain method. It expresses both the voltage and current in their partial-fraction forms in the Laplace domain, and obtains the impedance function of the transducer from the ratio of the voltage and current. The Prony-SS method is employed to extract the poles and residues of the voltage and current signals. Compared with traditional methods, the proposed method uses the transient signals, and will not suffer any leakage problems or resolution issues. In addition, this method requires only very short signals to obtain the impedance function, and is excellent for rejecting noise. This proposed model-updating method is a multi-physics FE model-updating method, including the correction of the elastic material properties based on a short-circuit model, and the correction of dielectric and piezoelectric parameters based on an open-circuit model. The fundamental updating algorithm employed in both steps is the cross-model cross-mode (CMCM) method. In addition to its accuracy and efficiency, this method has the advantages of both the direct matrix methods and indirect physical property-adjustment methods. Implementing the CMCM algorithm requires a knowledge of both the measured modal frequencies and the corresponding mode shapes, but a measured impedance function could provide only modal frequencies for shortand open-circuit transducers. When dealing with the incomplete modal information, an iterative procedure is taken. In each iteration, the “measured” mode shapes are approximated by the mode shapes obtained from the previous iteration’s updated FE model. In this study, we employed a tube transducer, which is made of piezoceramic material, to develop and test new methods of estimating the impedance function and updating piezoelectric constitutive properties. Both computer simulations and lab experiments have been conducted to verify the accuracy and efficiency of the proposed methods.


PREFACE
This dissertation follows the University of Rhode Island Graduate School guidelines for the preparation of a dissertation in standard format. The material presented in this thesis is divided into eight chapters.
• Chapter 1 serves to provide an introduction to the studies within this dissertation, as well as a review of the relevant research in FE model updating.
• Chapter 2 is an introduction to piezoelectric transducers, especially tube transducers. We begin with a brief overview of transducers and piezoelectric material. Both full and reduced constitutive equations are introduced to explain the performance of thin-tube transducers, and elaborate the relationship between piezoelectric constitutive properties. In addition, we address the analytical solution for ring (single DOF vibrators) and the FE solution for tube (multi DOF vibrators), and their electrical representations, i.e., equivalent circuits, are extended from those solutions.
• Chapter 3 presents the electromechanically coupled FE formulations and the procedures for modeling tube transducers with 2D and 3D elements using a commercial package. Furthermore, we present numerical verifications of the FE model, including comparisons with an approximate analytical solution, among two commercial FE packages, and with physical measurements.
• Chapters 4 and 5 show techniques and numerical studies for the estimation of impedance functions. In addition to a review of the traditional methods for estimating the impedance function of a piezoelectric transducer, in Chapter 4, we also newly propose an improved method using the Laplace domain. The newly proposed method decomposes the signals into their complex exponential components using the Prony-SS method, and it is particularly effective v for handling measured voltage and current signals that have been contaminated by noise. In Chapter 5, we discuss a numerical study of the estimation of the impedance function. To test the effectiveness of the proposed Laplace domain method, both numerical simulations and lab experiment are included. The impedance function that was estimated from the proposed method is compared with those obtained using traditional methods. The effectiveness of noise cancellation for both the proposed and traditional Fourier-based method is also investigated.
• In Chapter 6, we describe procedures for updating piezoelectric transducers.
Besides a brief review of traditional model updating methods, a new approach which is based on the extended cross-model cross-mode method is proposed.
The proposed model-updating method is a two-step method, including the correction of the elastic material properties based on a short-circuit model in the first step, and the correction of dielectric and piezoelectric parame-   The comparison of modal frequencies between 2D and 3D models 41

Statement of the Problem
The accuracy of the finite-element (FE) model of a piezoelectric transducer depends on the accuracy of the material's constitutive properties, but the values provided in the manufacturer's specification sheets for piezoelectric materials are often incomplete or inaccurate (Abboud et al., 1998). Because of the polycrystalline nature of piezoceramics and the statistical variations in their composition and production-related influences, the standard tolerance range of piezoelectric material is ±20% (Ltd, 2011). In addition, depolarization influences the accuracy of the material properties, which may result from operating temperatures that are too near the Curie temperature, high static-pressure cycling in deep-water applications, high alternating electric fields, or a slight extent from the passage of time (Sherman and Butler, 2007). Therefore, the properties of different transducers made with the same piezoelectric material may vary, and it is necessary to obtain a quick and accurate method of updating transducer models.
FE model updating has been defined as the systematic adjustment of FE models using the experimental data. When one makes measurements on a piezoelectric transducer, oftentimes only the impedance function is measured as the experimental data. Thus, to perform model updating for transducers, two major tasks will be covered: (i) developing and verifying an efficient method for estimating the electric impedance function of a transducer, and (ii) developing and testing a FE model updating method for piezoelectric transducers.
The electrical impedance function is a function of frequency, and it is generally defined as the total resistance offered by a device or circuit to the flow of an alternating current at specific frequencies. As an electrical feature, it is relatively easy to be estimated, and includes many important transducer properties such as the resonance frequency, capacitance, bandwidth, and damping. Traditionally, there are two methods for estimating the impedance function: (i) stepped harmonic analysis (i.e., single frequency) method, and (ii) Fourier-based analysis method using transient signals. The stepped harmonic analysis involves repeatedly applying a harmonic input signal to obtain the corresponding steady-state response, then computing the impedance of the given frequency. The Fourier-based method computes the impedance function from the complex ratio of the Fast Fourier Transform (FFT) of the input and output time-domain signals. However, both methods above have drawbacks. The stepped harmonic analysis method is time consuming and requires specialized equipment, while the Fourier-based analysis method suffers from problems related to frequency leakage and resolution (Lis and Schmidt, 2004;Lis et al., 2005;Su et al., 2014).
In terms of model-updating technologies, the traditional methods have been classified into two major types: direct matrix methods and indirect physical property adjustment methods (Friswell and Mottershead, 1995). The first of these two groups generally involves non-iterative methods, all of which are based on computing changes that are made directly to the mass and stiffness matrices. Such changes may have succeeded in generating modified models that have properties close to those measured in the tests, but the resulting models become abstract "representation" models, and cannot be interpreted in a physical way. The second group comprises methods that are in many ways more acceptable in that the parameters that they adjust are much closer to physically realizable quantities.
The methods in this second group, which seek to find correction factors for each individual FE or for each design parameter related to each FE, have generally been viewed as the main hope for the development of updating technology. The pub-lished model-updating methods that were implemented for ultrasound transducers belong to the second group of methods (Piranda et al., 1998;Piranda et al., 2001).
Because these methods are all iterative, a much greater computation effort is required, and the solutions for the physical properties are often prone to diverge (Hu et al., 2007).

Methodology and Procedures
In this study, we employed a tube transducer, which is made of piezoceramic material, to develop and test new methods of estimating the impedance function and updating piezoelectric constitutive properties. The effect of the backing layer and the transducer housing is not involved in this study, and the general procedure is also applicable to other types of piezoelectric transducers. method is employed to extract the poles and residues of the voltage and current signals (Hu et al., 2013). Compared with traditional methods, the proposed Pronybased method uses the transient signals, and will not suffer any leakage problems or resolution issues. In addition, this method requires only very short signals to obtain the impedance function, and is excellent for rejecting noise. Both computer simulations and lab experiment have been conducted to verify the accuracy and efficiency of this Laplace domain method.
The proposed model-updating method is the multi-physics FE model-updating method. This method is a two-step method, including the correction of the elastic material properties based on a short-circuit model in the first step, and the correction of dielectric and piezoelectric parameters based on an open-circuit model in the second step. The fundamental updating algorithm employed in both steps is the CMCM method (Hu et al., 2007;Hu and Li, 2008). In addition to its accuracy and efficiency, this method has the advantages of both the direct matrix methods and indirect physical property-adjustment methods. Implementing the CMCM algorithm requires a knowledge of both the measured modal frequencies and the corresponding mode shapes, but a measured impedance function could provide only modal frequencies for short-and open-circuit transducers. When dealing with the incomplete modal information, an iterative procedure is taken. In each iteration, the "measured" mode shapes are approximated by the mode shapes obtained from the previous iteration's updated FE model.
Specific to piezoelectric tube transducers poling in the radial direction, which obeys the thin-tube assumption, elastic material parameters associated with a tube's thickness direction would not affect the functionality of the piezoelectric tube. In turn, the corresponding impedance function depends only on four independent material parameters: two elastic parameters Y E 11 and ν 12 (or s E 11 and s E 12 ), the permittivity under constant stress ε T 33 , and the piezoelectric parameter However, with respect to References. (Hu et al., 2007;Piefort, 2001)

Introduction
Transducers refer to devices or agencies that convert energy from one form to another. To convert electrical energy to acoustical energy, or vice versa, we use an electroacoustic transducer (Burdic, 1991). The electroacoustic transducer may be used operated in either transmitting (projector) or receiving mode (hydrophone), or both. In transmitting mode, an ultrasonic wave is generated by an applied electrical voltage or current, while in receiving mode, an electrical signal is generated by an incoming acoustical wave. While loudspeakers and microphones are transducers that are used as sources and receivers of sound in air, the counterparts in water are projectors and hydrophones for sources and receivers, respectively.
A transducer may also be operated in either continuous-wave mode, where the transducer is harmonically operated at a specified frequency, or in transient mode, where the transducer is used to send out pulses (Kocbach, 2000).
Transducers were first used to obtain the direct measurement of the speed of sound in the fresh water of Lake Geneva in Switzerland in 1826. At that time, there were no electroacoustic transducer to generate sound in the water. Instead, the projector was a mechano-acoustic transducer, i.e., the striking of a bell under water. In 1915, Paul Langevin et al. began work in France using an electrostatic transducer as a projector and a waterproofed carbon microphone as a hydrophone (Sherman and Butler, 2007). They realized that the use of the piezoelectric effect in quartz had the potential for realizing improved transducers. Then, the search for new man-made transduction materials was further prompted by the potential for application during both World Wars and the Cold War. In 1944, piezoelectricity was discovered by A.R. von Hippel in permanently polarized barium titanate ceramics, and in 1954, even stronger piezoelectricity was achieved in polarized lead zirconate titanate ceramics (PZT), which are still being used in most underwater sound transducers (Jaffe et al., 1955). The development of underwater electroacoustic transducers expanded rapidly during the twentieth century, and it continued to be a growing field of knowledge (Sherman and Butler, 2007 acquire data related to a wide variety of topics such as the Acoustic Thermometry of Ocean Climate project (ATOC) and the Sound Surveillance System, which has been used to study the behavior of sperm whales, and to detect earthquakes and volcanic eruptions under the sea (Sherman and Butler, 2007).

Tube Transducers
The cylindrical shell is the basic shape of transducers that are widely used in many engineering applications. Three possible modes of cylindrical transducers are shown in Figure. 2. Figure. 2(a) shows a cylindrical shell that operates in the radial mode. The inner and outer cylindrical surfaces are coated with conducting material, and the poling field is applied between these surfaces. This geometry is characterized by moderate sensitivity and high output capacitance. Because of the large capacitance, the design is capable of driving long cables directly without suffering much loss in sensitivity. Figure. 2(b) shows the tangential mode of operation. Conducting stripes that are formed longitudinally on the ceramic cylindrical surface divide the cylinder into an even number of curved segments. Alternate stripes are electrically connected and a poling field applied in the circumferential direction. The cylindrical segments are thus electrically in parallel but mechanically in series. This configuration is characterized by high sensitivity and low capacitance, and is an excellent design where high sensitivity is required, and where suitable electronic isolation from the cable can be provided close to the hydrophone element.
The configuration shown in Figure. 2(c) operates in the longitudinal mode.
The ends of the ceramic cylinder are made to conduct electrically and the poling field is applied in the direction parallel to the cylindrical axis. This configuration has a somewhat higher sensitivity than the radial mode, and as with the tangential mode, the electrical impedance of the device is high (Burdic, 1991).
The present study focuses on the radial mode transducer. Only the active element (thin piezoelectric tube) which is sandwiched between a front layer and a backing layer is considered, and the effect of assembling and the transducer housing is not involved in the analysis. As a varying voltage is applied over the electrodes, the cylindrical transducer will begin to vibrate. This vibration is related to the geometry of this tube and the material parameters.
The present study focuses on the radial-mode transducer. We consider only the active element (thin piezoelectric tube), which is sandwiched between a front layer and a backing layer, and the effects of the assembly and the transducer housing were not considered in the analysis. As a varying voltage is applied over the electrodes, the cylindrical transducer will begin to vibrate. This vibration is related to the geometry of the tube and the material parameters.

Piezoelectric Material
At present, piezoelectric materials are the most commonly used material when developing transducers. When mechanical stresses are applied to a piezoelectric material solid, a voltage is produced between its surfaces, called a piezoelectric effect. Conversely, when a voltage is applied across certain surfaces of the solid, the solid undergoes a mechanical distortion, which is called the inverse piezoelectric effect. Most of the piezoelectric materials are crystalline solids. They can be single crystals, which are either formed naturally or by synthetic processes, or polycrystalline materials such as ferroelectric ceramics, which can be rendered piezoelectric in any chosen polar direction by the poling treatment. During the manufacture of piezoceramics, a suitable ferroelectric material is first fabricated into the desired shape. The piezoceramic element is then heated to an elevated temperature (the Curie point) while in the presence of a strong DC field.
The best known piezoelectric material is lead zirconate titanate (PZT) which is an intermetallic inorganic compound. It is a white solid that is insoluble in all solvents, and it has a perovskite-type structure. The remanent polarization gives PZT an approximately linear response to an alternating electric field. It is very stable and large enough to give a strong piezoelectric effect. However, this also causes depolarization to significantly affect the performance of PZT transducers. This depolarization may result from operating temperatures that are too near the Curie temperature, high static pressure cycling in deep-water applications, high alternating electric fields, or a slight extent from the passage of time.
Thus, while the properties of true piezoelectrics are determined by their internal crystal structure and cannot be changed, the piezoelectric properties of polarized electrostrictive materials depend on the level of remanence achieved during the polarization process, and may be changed by the operating conditions. Therefore, although PZT material has many advantages compared to other materials for underwater transducer applications, its limitations must also be considered during the design process (Piefort, 2001).

Constitutive Equations 2.2.1 Full Constitutive Equations
The linear piezoelectric effect in a piezoelectric medium can be described using a set of constitutive equations. When it is expressed in the so-called e-form, one has: where T ∈ R 6×1 is stress vector, S ∈ R 6×1 strain vector, E ∈ R 3×1 electric field vector, D ∈ R 3×1 electric displacement vector, c E ∈ R 6×6 elastic stiffness matrix, e ∈ R 3×6 piezoelectric stress matrix (superscript t stands for the transpose) and ε S ∈ R 3×3 permittivity matrix. The superscript of c E indicates that the elastic compliance coefficients matrix is obtained with E being held constant, and the superscript of ε S indicates that the permittivity coefficients matrix is obtained with S being held constant (Sherman and Butler, 2007;Piefort, 2001).
Generally, the coefficient matrices in Eqs. 1 and 2 include 45 independent coefficients. However, for polarized piezoelectric ceramic material (such as PZT) which belongs to the 4 mm (C4V) crystal class (IEE, 1987) , many of the coefficients are zero or related, leaving only 10 independent coefficients, and the corresponding equations are explicitly expressed in the e-form as: 4) where c E 66 = (c E 11 − c E 12 )/2. The subscript 3 is the direction in which ceramic element is polarized. Subscripts 1 and 2 represent directions perpendicular to poling direction, and 4, 5 and 6 refer to shear directions that follow IEEE convention, namely, 23 → 4, 13 → 5, and 12 → 6 (Wik, 2015).
A widely used alternative and equivalent representation consists in writing the constitutive equations in the so-called d-from: where s E ∈ R 6×6 is elastic compliance matrix, d ∈ R 3×6 piezoelectric strain matrix and ε T ∈ R 3×3 permittivity matrix at constant stress. Combined into one equation, the expanded form for 4 mm crystal class is explicitly expressed as (Balmes and Deraemaeker, 2013): From the constitutive equations of d-form and e-form, it is straightforward to obtain the following relationships: and Furthermore, Eq.10 can also be expressed as either or For quantifying the elastic parameters, many people prefer to use engineering parameters, instead of the compliance or stiffness parameters. For a transverse isotropic material whose properties are the same in the plane perpendicular to the poling direction, the compliance matrix s E can be expressed in terms of engineering parameters as:

Reduced Constitutive Equations for a Thin Cylinder
When the piezoceramic transducer polarized in the thickness direction is modeled as a very thin cylinder, the normal stress in the thickness (radial) direction and the respective transverse shear stress components are negligible (Erturk and Inman, 2011): If the electric field components in the axis-1 and axis-2 are zero, then Eq.7 becomes the reduced constitutive equations From the above equation, there are only four independent material parameters which will affect the dynamic characteristics of a piezoceramic thin tube transducer: s 11 , s 12 , ε T 33 and d 31 , noting that s E 66 = 2(s E 11 − s E 12 ) is not an independent parameter. Since s E 11 = 1/Y E 11 and s E 12 = −ν 12 /Y E 11 , the four independent material parameters can also be: Y E 11 , ν 12 , ε T 33 and d 31 .
In an attempt to show the corresponding e-form, one can first rearrange Eq.16 It is realized that the reduced d-form parameters are identical to those of the original constitutive equations, but the reduced e-form parameters are completely different from those of the original constitutive equations. For the finite element model updating in this study, the purpose is to update the material parameters based on the experimentally obtained impedance function of a piezoceramic thin tube transducer. Only four independent material parameters are necessary, or could be updatable, for a piezoceramic thin tube transducer; and it would be better to update the d-form parameters for the obvious reason mentioned above.

Piezoelectric Vibrators
The piezoceramic transducer to be studied is a thin tube type, where the length and diameter of the tube is much larger than its thickness. If the length of this transducer is also very small, then it becomes a piezoceramic "ring" transducer.
Traditionally, a ring transducer has been treated as a SDOF (single-degree-offreedom) system in the analytical derivation for its impedance (or admittance) function.

Ring Vibrator
A ring vibrator has been often treated as a 3-1 SDOF vibrator (see Appendix A), where the only strain/stress occurs in the circumference direction (also 1direction) and the poling direction (3-direction) could be the radial direction (or longitudinal direction). One can refer its detailed derivation for the impedance (or admittance) function to Reference. (Mason, 2013), while the following presentation provides a brief review.

Analytical Impedance Function for a Ring Transducer
The ring transducer poling in radial direction is treated as a 3-1 SDOF vibrator which has the reduced constitutive equation in its d-form to be: and its corresponding e-form is shown to be: Thus, the laterally clamped (in one dimension) permittivity where the coupling coefficient k 2 With the assumptions of negligibly small cross-sectional dimensions, the equation of motion for a ring transducer (see Figure.3) can be derived from Newton's second law applied to a segment. This incremental segment which is not balanced by external radial forces causes a radial acceleration as a whole. The resultant force on this incremental segment is where ρ is the density of the ceramic, a is the average radius of the ring, and u the radial displacement of the segment.
Since the elongation or contraction of circumference is linearly proportional to the change of radius, one has S 1 = u/a. Substituting T 1 by Eq. 27 into Eq. 31, one shows The transfer function (in the Laplace domain) from E 3 to u, denoted T u,E 3 (s) thus is shown to be where is the free resonance of the unloaded ring (E 3 = 0). Since S 1 = u/a, the transfer function from E 3 to S 1 , denoted T S 1 ,E 3 (s), thus is From Eq. 28, one shows the transfer function from E 3 to D 3 to be From E 3 = V /t where V is the electric potential and the current, I = 2π a ℓḊ 3 , one has the corresponding transfer functions and From Eqs. 36, 37 and 38, one obtains The admittance function Y (ω), in the frequency domain, of the ring transducer is defined as the transfer function T I,V (s) while substituting s by jω: Substituting for ω 1 (see Eq. 34) and k 31 (see Eq. 29), Eq. 40 becomes The reciprocal of the admittance function Y (ω) is the impedance function Z(ω).
Resonance frequency ω r and antiresonance frequency ω a are obtained from the pole and zero of the admittance, given by (Mason, 2013;Wilson, 1988) and respectively. Conversely, the coupling coefficient k 31 is given by

Single-degree-of-freedom Equivalent Circuits
Eq. 41 can be expressed as the sum of the blocked electrical admittance Y E (ω) and the motional admittance Y M (ω): The blocked electrical admittance is where C 0 is the clamped capacitance and the free (or total) capacitance of the ring is where the motional capacitance C 1 is in which the mechanical compliance of the ring C m is (52) and the turns ratio N is Referring to the simple equivalent circuit of Figure. 4, using Eq. 51 for C 1 , the motional inductance L 1 can be derived from ω 2 1 = (1/L 1 C 1 ) as: Substituting Eq. 34 for ω 2 1 and Eq. 52 for C m into Eq. 54, one obtains where M = 2πatℓρ is the total mass of the ring. where liberty has been taken to insert a damping factor R m . The equivalent circuit is an electrical alternative representation of transducer that can be combined with power amplifier circuits, and other electrical circuits. It is sometimes more convenient to fit the impedance plots to lumped circuit models in order to predict the electrical behaviour of the transducer.

Thin Tube Vibrator
This dissertation focuses on the study of a thin-wall tube transducer (see Figure. 6). The main difference of a tube transducer from a ring transducer is the axial length ℓ of the tube is not small relative to its diameter. There is no simple analytical solution for the admittance function of the tube transducer because multiple complex vibration modes will contribute to the admittance function. In practice, finite element models are often built to obtain the vibration modes, as well as the admittance (or impedance) function.  (Piefort, 2001).

Mode-driven Equivalent Circuit
For a SDOF ring vibrator, one rewrites the admittance function Eq. 45 as: From Eq. 46 and Eq. 49, one has and Since ω 2 1 = 1/(L 1 C 1 ), Eq. 58 can also be written as The above admittance function Y 1 (ω) can also be obtained from the reciprocal of the impedance function Z 1 (ω): For an N-degree-of-freedom (NDOF) vibrator, applying the principle of mode superposition -voltage is the input and current is the output, one can write the corresponding admittance function to be and Z n (ω) = j ω L n + 1 in which C n and L n are the nth motional capacitance and nth motional inductance, respectively. They are related to the nth modal (resonant) frequency ω n of a transducer through ω 2 n = 1/(L n C n ).
Eq. 61 to Eq. 63 implies that the transducer can be represented by an equivalent circuit, which includes N parallel modal branches, as shown in Figure. 9. This equivalent circuit is often called a mode-driven equivalent circuit (Ballato, 1990).
Each branch represents one mechanical mode, and C 0 is the clamped capacitance, which is equal to the capacitance of the transducer without the piezoelectric effect. The mathematical model of a physical system often results in partial differential equations that cannot be solved analytically because of the complexity of the boundary conditions or domain. To study a complex model, the FE method is often employed. The FE method subdivides a complex system into its individual components or elements whose behavior is readily understood, and it then rebuilds the original system by assembling such components. Nowadays, the FE method has become one of the most important analysis tools in engineering.
Since the late 1960s, the work by Eer Nisse (1967) and Tiersten (1967) have established variational principles for piezoelectric media. Allik & Hughes (1970), who proposed a tetrahedral element that accounted for piezoelectricity, contributed greatly to the FE modelling of structures with embedded piezoelectricity. In 1986, piezoelectric elements were for the first time included in the commercial FE program ANSYS. Later, piezoelectric finite elements have also been included in other commercially available FE packages such as ABAQUS and COMSOL (Kocbach, 2000;Manual, 2005;Multiphysics and Module, 2014). These commercial software packages provide a wide range of simulation options that control the complexity of both the modelling and analysis of a system. It allows users to easily build an FE model with piezoelectric elements, and simplifies the analysis procedure. In the present study, we adopted primarily ANSYS to perform FE modelling using piezoelectric media.
Many studies have focused on the FE analysis for various piezoelectric structures (Hwang and Park, 1993;Kim and Lee, 2007;Lewis et al., 2009). In this study, we focus only on piezoelectric tube transducers that were polarized in the radial direction. The specific analyses included in this study are modal and transient (dynamic response) analysis. We employed modal analysis to obtain the modal frequencies and mode shapes of the piezoelectric tube transducer under both open-and short-circuit electric boundary conditions. We used the transient analysis to obtain the electric charge time signal while the piezoelectric tube transducer is subjected to an input electric-potential time signal. The ultimate purpose of the transient analysis is to obtain the impedance function of the transducer.

Finite Element Formulation 3.1.1 Piezoelectric Element
For piezoelectric elements, both displacements and electric potentials exist at the nodal locations. The displacement fieldû and electrical potentialsφ over an element are approximated within the element aŝ where N u and N ϕ are the array of interpolating functions, and u and ϕ are nodal quantities. The body forces and charges as well as the surface forces and charges are interpolated in a similar manner (Kohnke, 1999).
The displacement fieldû and the electric potentialφ are related to the strain field S and the electric field E by and where ∇ is the gradient operator and D is the derivation operator defined as: Therefore, the strain field S and the electric field E are related to the nodal displacements and potential by the shape-functions derivatives B u and B ϕ , which are defined by The dynamic equations of a piezoelectric continuum can be derived from the Hamilton principle, in which the Lagrangian and the virtual work are properly adapted to include the electrical contributions as well as the mechanical ones (Piefort, 2001;Allik and Hughes, 1970;Lerch, 1990;Tzou and Tseng, 1990).
Without further elaboration, for a piezoelectric element, the resulting mechanical equilibrium equation and the electric-flux conservation equation can be written in the form and respectively, where is the element mass matrix (no inertia terms exist for the electrical flux conservation equation), ρ is the mass density, is the element displacement stiffness matrix, is the element piezoelectric coupling matrix, is the element dielectric "stiffness" matrix, f is the element mechanical force vector, and g is the element electric charge vector. In addition, V denotes the volume of the element (Allik and Hughes, 1970).

Assembled System Equation
In the present study, the total volume V p is divided into P volumes V (j) p , j = 1, · · · , P . Any variable A defined in element j will be written as A (j) as follows. The relation between the local node-displacement vector u (j) and the global node-displacement vector U is given through a transformation matrix L and the relation between the local potential vector ϕ (j) and the global potential vector Φ is given through a transformation matrix L (j) ϕ , such that The system equations for the whole structure are obtained by the summation of the contribution from each FE, which can be written in terms of nodal quantities: where the assembled matrices are given by: Eq. 79 couples the mechanical variables U and the electrical potentials Φ; F represents the external forces applied to the structure, and G represents the electric charges brought to the electrodes (Piefort, 2001).

Open-and Short-circuited Electrical Boundary Conditions
If the electric potential Φ is controlled, the governing equations become where the second term in the right hand side represents the equivalent piezoelectric loads. From Eq. 87, one sees that the eigenvalues problem of the system with shortcircuited electrodes, namely Φ = 0, is It can be seen from Eq. 87 that the modal frequencies, ω, and mode shapes of the short-circuit model are the same as if there was no piezoelectric electromechanical coupling.
Conversely, open electrodes correspond to a charge condition G = 0. In this case, from Eq.79, one writes Upon substituting Eq. 89 into Eq. 87, it yields which shows that the global stiffness matrix depends on the electrical boundary conditions. From Eq. 90, one sees that the eigenvalues problem of the system with open-circuited electrodes is and

Cylindrical Transducer Modelling and Analysis
The commercial FE packages such as ANSYS have the capability to perform a fully coupled piezoelectric analysis. In this study, both 2D axisymmetric plane piezoelectric elements and 3D solid piezoelectric elements are suitable for modeling the axisymmetric piezoelectric tube.

Modelling and Analysis of 2D Element Model
A specific 2D coupled-field plane element of ANSYS is the PLANE13 element which will be utilized in this study for modeling axisymmetric transducers. Shown in Figure. 10 is the PLANE13 element which is defined by four nodes, with three degrees of freedom per node for piezoelectric analysis, including x-direction displacement, y-direction displacement and voltage (Kohnke, 1999). Please refer to Appendix B for more detailed ANSYS information about element settings, properties assignment, piezoelectric analysis and data extraction. The poling direction is radially outwards from the axis of symmetry, namely z-axis. Since the order of the stresses in ANSYS differs from those typically used in electrical applications (IEEE standard), care must be taken while inputting the corresponding coefficients matrices for the tube transducer model.   ANSYS has a built-in capability to compute the impedance function of a transducer by using the harmonic analysis. In the harmonic analysis, we applied a unit amplitude voltage having a frequency sweeping from 0 to 120 KHz with an increment of 200 Hz in the radial direction, and Figure. 14 shows the obtained impedance function. From the impedance function plot, the frequencies corre-   sponding to the local minima and local maxima are commonly referred to as the resonance frequencies and anti-resonance frequencies, respectively, of multi-mode transducers. In theory, the resonance and anti-resonance frequencies should agree well with the modal frequencies obtained from the short-and open-circuit models, respectively. Only the modal frequencies corresponding to the symmetric modes of the short-and open-circuit models (see Figure. 12 and 13) are plotted as the vertical lines in Figure. 14. Clearly, they are in line with the resonance and anti-resonance frequencies of the impedance function.
Because the axially uniform voltage input is symmetric to z = 0, only symmetric modes can be excited. Consequently, only symmetric modes "participate" in the mechanical vibration and the associated electrical impedance function. Furthermore, because of the axially uniform nature of the input, the mode that has the most uniform mode shape would dominate the mechanical vibration. In fact, the most uniform mode shape for the short-circuit model has a modal frequency of 79.35 KHz (see Figure. 12), which is the global minimum of the impedance function, and the most uniform mode shape for the open-circuit model has the modal frequency 95.09 KHz (see Figure. 13), which is the global maximum of the impedance function.

Modelling and Analysis of 3D Element Model
In this study, the specific 3D element type employed to build up the piezoelectric cylinder is SOLID226 of ANSYS, as shown in Figure. 15. This type of element contains 20 nodes with four degrees-of-freedom (DoFs) per node, including the x-direction displacement, y-direction displacement, z-direction displacement, and voltage (Kohnke, 1999). Please refer to Appendix B for more detailed ANSYS information about element settings, properties assignment, piezoelectric analyses, and data extraction. The 3D tube model that was built using SOLID226 element is shown in Figure.16. Because building the 3D model follows the same geometry, boundary conditions, and material properties, as utilized for the 2D model, similar modal frequencies are expected from the 2D and 3D models. A comparison of the first three modal frequencies between the 2D and 3D models for both short-and open-circuit models is presented in Table 3, which shows an excellent agreement between these two models. In Table 3, the relative error was defined as the 2D result minus the 3D result, and it was then normalized by the 3D result. Furthermore, Figure. 17 shows the excellent agreement of impedance functions obtained from the 2D and 3D models.

Constitutive Equations of 2D and 3D Models
The distinction of the constitutive equations for 2D (axisymmetrical) and 3D models is very important, as it is often not well understood and many errors can arise from the confusion between 2D and 3D properties of piezoelectric materials.
Note however that the d ij , s E ij and ε T ii coefficients are equal for 2D and 3D constitutive equations. It is therefore preferable to handle the material properties of piezoelectric materials in the d-form.
In calculations for axisymmetric structures, it is preferable to use a cylindrical coordinate system where the symmetry of the problem can be exploited. Let the three-dimensional cylindrical coordinate system be denoted by (R, θ, Z), where R is in the radial direction, θ the circumferential direction, and Z axial direction. The strain vector for the cylindrical coordinate system is given by: where The electric field vector for a cylindrical coordinate system is defined by where L ϕ is defined by the last transition.
When only axially symmetric vibrations are considered, this axisymmetrical Additionally, let axially symmetric torsional (circumferential) vibrations be excluded, i.e., Thus, only two coordinates are needed to describe the problem in the axisymmetric case: R and Z. This class of problems has been often referred to as 2D axisymmetric problems.
Imposing Eqs. 97 and 98 on Eqs. 93 and 96 yields Only four of the six strain components and two of the three components of the electric field need to be included in the analysis. This results in the following definition for the strain vector and the electric field vector in the two-dimensional axisymmetric case: and The derivation operator matrices L u and L ϕ therefore have the following definitions in this case and In terms of the d-form, the reduced constitutive equations for the axisymmet- To derive the corresponding e-form, one can first rearrange Eq. 106 into:  Let the e-form of the reduced constitutive equations be denoted by From the above two equations, one has The resulting stiffness matrixc E is a 4 × 4 matrix: 0 Note that the stiffness coefficientc ij of the 2D formulation differs from the stiffness coefficient c ij of the 3D constitutive equations.

Numerical Verifications
While using finite element models to conduct research, the numerical verification of them is always essential. In the following presentation, the verification process includes: (1) the comparison with an approximate analytical solution, (2) the comparison among two commercial finite element packages, and (3) the comparison with physical measurements.

Comparison with Analytical Solution
One can compare the impedance function of the FE model with the analytical impedance function based on a ring transducer which has been presented in Chapter 2. While using the analytical formula, Eq. 40, the following parameters Shown in Figure. 18 is the analytical impedance function, comparing to the impedance function from the finite element model. The analytical impedance function, which has been derived from a mathematical model of a single-degreeof-freedom vibrator, contains only a "breathing" vibration mode. In contrast, the impedance function from the finite element model has been contributed by multiple vibration modes. Due to the imposed simplifications for the derivation of the analytical impedance function, the significant discrepancy between these two impedance functions is not surprising.
For better matching the impedance function computed from a finite element model to the analytical impedance function, a "constrained" finite element model has been built. This constrained model does not allow movements in the axial direction for all elements. The resulting impedance function of the constrained model -contributed only by a uniformly breathing mode as well -is shown in Figure. 19. Although there is still discrepancy between the analytical impedance function and the impedance function from the constrained finite element model, the difference has become much smaller.

Verification by ABAQUS Model
Unless state otherwise, all FE models throughout this dissertation will be built by using ANSYS. For independently verifying the correctness of the corresponding ANSYS model, a finite element model has been constructed by using ABAQUS.
Examining on several specific sinusoidal input (voltage) signals, excellent agreement between ANSYS and ABAQUS models has been found while comparing the output (charge) time histories. A good agreement for the impedance function is shown in Figure. 20. Figure 20. The comparison of impedance functions obtained from ANSYS and ABAQUS

Comparison with Experimental Data
An ultimate verification for a finite element model should be done with the lab data. To this end, an experiment has been conducted to measure the impedance function of a real transducer with the same geometrical dimensions as those utilized in the finite element model. Shown in Figure. 21 is the comparison of the impedance functions obtained from ANSYS and the lab measurement. By and large, they match reasonably well with each other. One obvious reason for the discrepancy between the impedance functions obtained from ANSYS and the lab measurement is that the piezoelectric material properties used in ANSYS model could not be the same as those of the real transducer. In fact, the ultimate goal of this dissertation is to perform FE model updating on the piezoelectric material properties so that the impedance functions obtained from ANSYS and the lab measurement can match better. The algorithm and procedure to carry out the multi-physics FE model updating on the piezoelectric material properties will be addressed in Chapter 6.

Methods of Estimating Impedance Function
One of the most important characteristics of transducers is its impedance function, which is usually obtained from processing the measured voltage and current signals. In addition to a review of the traditional methods for estimating the impedance function of a piezoelectric transducer (Lis and Schmidt, 2004;Umchid et al., 2009), in this chapter, we newly propose an improved method. This newly proposed method, which is operated in the Laplace domain, is particularly effective for handling measured voltage and current signals that have been contaminated by noise.

Traditional Methods
Traditionally, there are two kinds of methods for estimating the impedance function: (i) stepped harmonic analysis (i.e. single frequency) method, and (ii) Fourier-based analysis method using transient signals. The stepped harmonic analysis is to repeatedly apply a harmonic input signal to get the corresponding steady state response, then compute the impedance of the given frequency (Lis and Schmidt, 2004). The Fourier method computes the impedance function from the complex ratio of the FFT of the input and output time domain signals.
Both methods have drawbacks. The stepped harmonic analysis method is very time consuming; and the Fourier-based analysis method always suffers from the problems related to the periodic assumption of both input and output signals, such as frequency leakage and resolution.

Proposed Method
This study proposes a new method to estimate the impedance function via the Laplace domain in order to overcome the drawbacks of the traditional methods.
The first step of the proposed method is to decompose both input and output time domain signals into complex exponential components. Because the input voltage signal v(t k ) can be controlled, one can choose to use an exponentially decayed voltage signal: where λ v , γ v = constants. For the sake of a clear presentation, let the output signal (from FE model) be a transient charge signal, c(t k ). Decomposed into complex exponential components, c(t k ) can be written into where p is the number of the decomposition terms, and the complex exponentials λ n and complex amplitudes γ n are computed coefficients. Because c(t k ) is a real signal, λ n and γ n must either be real numbers or occur in complex conjugate pairs.
Let λ n = −a n + jω n , then a n is the damping factor in seconds −1 and ω n is the frequency in radians. Since the Laplace property L(e λv t ) = 1/(s − λ v ), the Laplace Similarly, the Laplace transform of the output charge c(t) is: One can find the complex amplitudes and exponentials of a signal in the time domain are really the poles and residues in the Laplace domain associated with the signal. In turn, both input and output signals in the Laplace domain can be readily established in a partial fraction form.
To calculate the impedance function, one needs to get the current signal, i(t), from the derivative of the charge with respect time, that is i(t) = dc(t)/dt. By if the initial charge c(0) is equal to 0.
The impedance function in the s domain can be readily obtained from Eqs. 113 and 115: Finally, substituting s by jω in the above equation, one obtains the impedance function of ω.

Complex Exponential Decomposition
To implement the proposed method, an essential task is to decompose the signal into its complex exponential components. An arbitrary signal y(t) can always be decomposed into a finite number of complex exponential components, called Prony's series.
where p is the number of terms. For digital signal analysis, the continuous y(t) is not known and usually equally spaced samples are available. Denoting sampling interval ∆t, t k = k∆t and y k ≡ y(t k ), where k = 0, . . . , m, one writes the discrete Prony's series as where z n = e λn∆t .
In Eq. 117, the exponential functions e λnt , n = 1, · · · , p form a basis on the open interval 0 ≤ t < T . Clearly, the Prony series (Eq. 117) can also be viewed as the general solution of a pth-order homogenous linear ordinary differential equation with constant real-valued coefficients a n : where y (n) (t) = d n y/dt n . And the exponents λ n (n = 1, · · · , p) in Eq. 117 are the roots of the characteristic equation associated with Eq. 119: If only the equally-spaced sampled data y k (k = 0, 1, · · · , N − 1) are given, then Eq. 119 becomes a pth-order difference equation: where b n are real-valued constants. Without losing the generality, let b p = 1.
From Eqs. 118 and 121, the characteristic equation associated with Eq. 121 can be obtained to be: Note that the linear difference equation and its characteristic equation have the same coefficients b n .

Prony's method
Traditionally, Prony's method has often been utilized to decompose a signal into complex exponential components. In brief, Prony's method is a technique for fitting an exponential model to a few equally spaced measured data points (Prony, 1795). While the original Prony's method presented a method that exactly fits as many purely damped exponentials as required to fit the available data points, the modern version of Prony's method generalizes to damped sinusoidal models. In addition, it utilizes the least-squares analysis to approximately fit an exponential model for cases where there are more data points than needed to fit the assumed number of exponential terms. However, because Prony's method suffers from an ill-conditioned problem with respect to determining the roots of the high-order polynomial that is encountered, the estimation of Prony's method may have a significant error (Hu et al., 2013).

Prony-SS method
To avoid dealing with the ill-conditioned problem, in this study, we employed an alternative approach that uses a first-order matrix homogenous difference equation (state-space model) to replace the high-order homogenous difference equation, and the approach is called the Prony-SS method (Hu et al., 2013). Although the Prony-SS method and Prony's method appear to be theoretically identical, they are drastically different with respect to crucial numerical issues, including conditioning and stability. While Prony's method is very sensitive to the sampling rate and round-off error, the Prony-SS method is not. While Prony's method has trouble dealing with noise embedded in the signal, the Prony-SS method can handle noisy signals properly because it has a built-in noise-rejection mechanism that uses truncated singular-value decomposition. While root-finding of a high-order polynomial, which is a classical ill-conditioned problem, is a required step in Prony's method, the Prony-SS method completely avoids it (Hu et al., 2013). factor of each exponential term, and (3) solve for a set of linear equations to yield the estimates of the exponential amplitude and sinusoidal initial phase of each component (Hu et al., 2013).
By introducing and assuming a p = 1, Eq. 117 is converted into a first-order matrix differential equation asẋ and The solution of Eq. 124 has the form x(t) = ve λt where λ is a constant and v is a time-independent p-element column vector. When x(t) = ve λt is substituted into Eq. 124, it yileds the algebraic eigenvalue equation: Theoretically, the p eigenvalues of the matrix F are exactly the p roots of the characteristic equation of Eq. 120. Instead of determining the roots of Eq. 120, the Prony-SS method is to compute the eigenvalues of the matrix F.
The elements of vector x in Eq. 126 are often termed the state variables of the linear dynamic system, and they serve as the "coordinates" of the system. However, there are an infinite number of ways of selecting the coordinates. Denoting an arbitrary nonsingular transformation matrix by Φ ∈ R p×p , and another set of coordinates by q ∈ R p , one can formulate the transformation relationship between x and q as: Now, substitute x(t) = Φq(t) andẋ(t) = Φq(t) into Eq. 124 and premultiply both sides by Φ −1 , the resulting equation becomeṡ where G = Φ −1 FΦ. Note that Eq. 129 and Eq. 124 share the same form. Because G is a similarity transformation of F, they must have the same eigenvalues.
Let the eigenvalues and the corresponding eigenvectors of G be denoted by λ 1 , λ 2 , · · · , λ p and u (1) , u (2) , · · · , u (p) , respectively. Because u (1) e λ 1 t , u (2) e λ 2 t , · · · , u (p) e λpt form a basis for solutions of Eq. 129, the general solution of Eq. 129 can be written as: where g n are arbitrary constants. In Eq. 130, each element of q(t) is a linear combination of e λnt , n = 1, . . . , p, so is y(t) in Eq. 117. Therefore, y(t) is a linear combination of the elements of q(t), that is, or where q n (t) is the nth element of q(t) and c ∈ R p is a column vector.
The discrete version of Eq. 129 is a first-order homogeneous matrix difference equation which can be written as: where q k ≡ q(t k ) and A = exp(G△t). From Eq. 132, at t = t k , one has Here, the specific goal is using the given discrete signal y k , k = 0, · · · , N − 1, to obtain a realization of A. Next, the eigenvalues of A, and those of F (i.e. λ n , n = 1, · · · , p), can be computed without going through the procedure of determining polynomial coefficients and zeros (Hu et al., 2013).
The Hankel matrix H(k) ∈ R ξ×η was defined as Upon substituting Eq. 134 into Eq. 135, it can be shown that H(k) is the product of three matrices: where and A realization of A can be obtained from H(0) and H(1). Substituting k = 0 and k = 1 into Eq. 136 yields and Applying the singular value decomposition of H(0), one writes (Zhou and Doyle, 1996): In theory, the model order of the dynamic system, p, is equal to the number of non-zero singular values in Eq. 141, that is, U 1 ∈ R ξ×p , Σ 1 ∈ R p×p and V 1 ∈ R p×η , and the rank of H(0) is p. The singular values should go to zero when the rank of the matrix is exceeded, but for measured data due to random errors and small inconsistencies in the data, the singular values will not become zero but will become very small. While choosing the size of H(0), both ξ and η must be greater than p.
A direct comparison between Eqs. 139 and 141 suggests that P ξ is related to U 1 , and Q η to V T 1 . A possible "balanced" choice is P ξ = U 1 Σ 1/2 1 and Q η = Eq. 140, one obtains: While the computed eigenvalues of A are z n , n = 1, · · · , p, the corresponding λ n in Eq. 117 are λ n = ln(z n )/∆t. Finally, the complex coefficients γ n can be computed by solving a second set of linear equations (Hu et al., 2013):

Multi-signal Prony-SS method
While the above Prony-SS method has been derived to deal with a single signal decomposition, it can also be extended to perform the multiple signals (Hu et al., 2016). We denote the column vector y k ∈ R N ×1 for the signal data, where N is the number of signals. The multiple-signal vector can be written as: where The extension from the single-signal Prony-SS method to multiple signals is quite straightforward. The first two steps of the proposed method are the same as the single-signal Prony-SS's step described above. We first obtained a realization for the state matrix of a first-order matrix difference equation, and we then compute the eigenvalues of the realization matrix to obtain the estimates of the frequency and damping factor of each exponential term. The distinction between the singlesignal and multiple-signal Prony-SS method is that the elements in Eq. 135 are replaced by corresponding block matrices, that is, Further, when considering multiple signals, Eq. 134 becomes where C t is ∈ R N ×p with c t n ∈ R 1×P , n = 1, · · · , N as its rows and c t n is the row vector corresponding to the nth signal. Upon substituting Eq. 148 into Eq. 147, it can be shown that H(k) is the product of three matrices: and Q η stay the same as Eq. 138.
Then using the same step illustrated in single signal, a realization of A can be obtained from H(0) and H (1): where and U T 1 ,V 1 andΣ 1 can be obtained by applying the singular value decomposition of H(0), one writes (Zhou and Doyle, 1996): In theory, the model order of the dynamic system, p, is equal to the number of non-zero singular values in Eq. (154), that is,Û 1 ∈ R ξ×p ,Σ 1 ∈ R p×p and V 1 ∈ R p×η , and the rank of H(0) is p. While choosing the size of H(0), both ξ and η must be greater than p. The computed eigenvalues of A are z n , n = 1, · · · , p, and the corresponding λ n in Eq. 112 are λ n = ln(z n )/∆t.
Once z n , n = 1, · · · , p, have been computed, the next step involves the solution of a second set of linear equations. The matrix form of equations is written as:  The residues Γ n , n = 1, · · · , p can be calculated using Eq. 155 and , consequently, the amplitudes and phase angles.
When estimating the impedance function, it has been theoretically proven that the poles of the output signal of a linear system will contain the poles of the input signal (Hu et al., 2016). To guarantee that the estimated poles of the output (current) signal contain the input (voltage) poles, it is necessary to process both input and output signals simultaneously by using the multi-signal Prony-SS method. An intended feature of the multi-signal Prony-SS method is that all input and output signals share the same model -a pth-order difference equation. With this feature, a consistent set of "global" poles (natural frequencies and damping factors) exists among all signals.

Noise Rejection and Damping Modelling
It is inevitable that the measured input and output signals must contain various noises, and a realistic transducer must include damping (energy dissipation) mechanism. This section addresses the noise handling by using the traditional frequency domain method and the proposed method. As the proposed method has the capability of estimating the damping ratios associated with multiple vibration modes, this section also includes using a Rayleigh damping for the finite element model based on the estimated damping ratios of dominant vibration modes.  To remove the noise, the traditional frequency-domain method used to estimate the frequency-response function (reciprocal of impedance function) of a linear system is based on the calculation of the averaged auto-power spectra

Noise-cancelling Methods Statement
and together with the cross-power spectra and where N is the number of segments, and P n and U n are the Fourier transforms of measured time-domain signals of a segment p n and u n , respectively, while ( . ) * denotes the complex conjugate. Using the averaged spectra of the measured input and output signals, the frequency-response function (reciprocal of the impedance function) of the system can be estimated by the following estimator: Note that because of limitations of the Fourier transform, this Fourier-based denoising method suffers from the problems related to the periodic assumption, such as frequency leakage and resolution. In principle, the averaging procedure usually requires the amount of data to be very large (Craig and Kurdila, 2006).

Noise Cancellation of Laplace Method
The noise cancellation of the proposed method is based on neglecting the insignificant exponential components of the signal. Assume that s k (k = 1, · · · , i + 1) is a singular value of the Hankel matrix, sorted from high to low, Eq. 141 can be written as where m is the amount of data and i is the number of rows of Hanker matrix. We can find the value of s k that determines the effect of a component on the entire signal. Because the noisy components are much smaller relative to the system components, one can take the first n components as the response from the system, and the rest are from noise. Therefore, Eq. 162 can be written into We can neglect the components ranked after n to achieve noise cancelling. It is important to note that the number of rows should not be too small, as that will affect the rank of system, while the number should not be too big as that will affect the efficiency of noise cancelling. Therefore, the noise cancelling of the proposed method depends on the choice of the noise threshold n, which is the number of exponential components produced by the system but noise. A conventional way of choosing the noise threshold is to find a significant gap in the normalized singular values. Generally, there is a sharp drop between the system and noise components, and therefore, we should take the n, which is the number of singular values before the sharp drop, as the rank of the system.

Damping Estimation
The proposed method is capable of estimating the damping ratio, ξ n , of the impedance function at resonance frequencies. After computing the eigenvalues z n of the transition matrix A using Eq. 143, the corresponding λ n can be calculated by λ n = ln(z n ) ∆t (164) and the damping ratio ξ n = a n |λ n | where a n is the real part of λ n .
When using FE commercial packages, we often cannot specify damping coefficients directly, but employ them while setting the harmonic analysis or transient analysis. The Rayleigh damping matrix, C, is a widely used damping model, and is defined as (Cai et al., 2002): where M is the system mass matrix, K is the system stiffness matrix, and α and β are two coefficients that are related to the damping ratios by where ω n is the nth modal frequency of the system. When the two damping ratios, ξ 1 and ξ 2 , which are associated with frequencies ω 1 and ω 2 , are estimated using the proposed method, the values of α and β can be calculated by and respectively.

CHAPTER 5 Numerical Study of Estimating Impedance Function
In this chapter, we present a numerical study of the estimation of the impedance function of a piezoelectric transducer. To test the effectiveness of the proposed Laplace domain method, both the numerical simulation and the lab experiment are included. The impedance functions estimated from the proposed method are compared to those obtained from traditional methods. In addition, the effectiveness of noise cancelation for both the proposed and the traditional Fourier-based methods is also investigated.

Simulation Study
The simulation model is a 2D FE model, shown in Figure. 11. The impedance function estimated by the stepped single-frequency method which repeats hundreds or thousands of single-frequency analysis is considered to be accurate. Some commercial FE packages, including ANSYS, have the built-in capability to compute the impedance function of a transducer using the stepped harmonic (single-frequency) analysis (Moaveni, 2003). By using ANSYS with the frequencies ranged from 0 to 120 KHz with an increment 200 Hz, the impedance function of the present FE model from the harmonic analysis is shown in Figure. 23.

Estimation with Laplace-and Fourier-based Methods
The drawback of the stepped single-frequency method is that it is very time-  1. The impedance function estimated from the proposed method agrees well with that from the stepped single frequency method.
2. The proposed method significantly outperforms the Fourier-based method, especially at the region of high frequency. Note that the damping factor for this FE model has been set to be equal to 0.
In summary, the proposed method significantly outperforms the Fourier-based method. Furthermore, the signal required for the proposed method is much shorter than that for the Fourier-based method. The impedance function estimated by the proposed method also has a much higher resolution than that of the Fourier-based method.

Estimating Impedance Function in Noisy Environment
To simulate a more realistic environment, we included the Rayleigh damping with α = 0 and β = 1 × 10 −7 for the FE model, and added 10% random Gaussian white noise to the clean signals. The level of the additive white noise is quantified by a stated percentage, which is defined as the ratio of the standard deviation of the Gaussian white noise (σ n ) to that of the noise-free signal (σ s ): σ n /σ s × 100%.
Let the signal-to-noise ratio be defined as S/N = 20log 10 (σ s /σ n ). A 10% noisy signal corresponds to S/N = 20. Figure 26. Impedance functions estimated from noisy signals by FFT methods Shown in Figure. 26 are two impedance functions estimated from noisy signals, including: (1) direct FFT approach which obtains the impedance function based on the ratio of the FFT of the voltage to that of the current, and (2) auto/cross power spectra approach, i.e. Eq. 161. For comparison purpose, Figure. 26 also shows the impedance function that is obtained from the stepped harmonic analysis. While using the power spectra approach, we use 20 signals, each having 5 × 10 4 data points, to calculate the averaged impedance function.
Next, the noise rejection of the proposed method by using a truncated singular value decomposition technique is carried out. The first 100 singular values of H(0) 500×500 for the charge signal, which are normalized by the first singular value, are shown in Figure. 27.

Figure 27. Singular Values of Hankel Matrix
Clearly, the singular value curve in Figure. 27 has a significant gap between its 26th and 27th singular values. As the gap may represent the separation between the true signal and noise, the model order for this signal is set equal to 26. By implementing the proposed method with model order 26, we obtain the impedance function of the simulation model shown in Figure. 28. One notices that a good agreement has been achieved between the stepped single-frequency method and the proposed method using only 1000 data points. When 3000 data points are employed in the proposed method, it results in an even better match (see Figure. 29). Figure 28. The impedance function using 1000 data points (proposed method) Figure 29. The impedance function using 3000 data points (proposed method)

Laboratory Test
Ultimately, the effectiveness of the proposed method must be tested for a real transducer (see Figure. 6). The laboratory test was setup (see Figure. 30) for estimating the impedance function of a real piezoelectric tube. In this experiment, the transducer tube was mounted on a "soft" base (a sponge pad) to simulate a free-free mechanical boundary condition. For detailed information about these pieces of equipment, one can refer to References. (Mul, 1998;Tektronix, 2008).
The corresponding circuit diagram of the experiment is shown in Figure. 31, where the resistance R 1 = 3140 ohm, V 1 is the input voltage generated by the Figure 31. The circuit diagram of the experiment signal generator, and V 2 is the response voltage on the piezoelectric tube. The impedance Z of the piezoelectric tube can be calculated by The experimental procedure requires three steps to estimate the impedance function: 1. Generate the input signal V 1 on a computer and download it to the generator (AWG5000); 2. Conduct the experiment and collect the response signal V 2 ; 3. Perform the necessary signal-processing tasks to obtain the impedance function.
Both the stepped harmonic signals and transient signals have been utilized as the input V 1 . The harmonic input signals were generated repeatedly at uniformly spaced frequencies with intervals of 0.24 KHz. The transient signal was an exponential decay signal with a maximum of 1V and a decay ratio of −3 × 10 3 . Figure. 32 shows the transient input and its corresponding output.

Traditional Model Updating Methods
Traditionally, model updating methods can be classified into two major groups: direct methods and iterative methods.

Direct Methods
Direct methods are generally of non-iterative methods, all of which are based on computing changes made directly to the mass and stiffness matrices. Typical direct methods include Lagrange multiplier methods, matrix mixing methods and control theory methods (Friswell and Mottershead, 1995). These methods have the advantage of not requiring any iterations, thus eliminating the possibilities of divergence and excessive computation. And abstract "representation" models could be built up to reproduce the measured data exactly.
One major drawback of direct methods is that the updated mass and stiffness matrices have little physical meaning and cannot be related to physical parameters to the finite elements in the original model. Another drawback of direct methods is the need to expand the mode shape vectors to the full set of finite element model degrees of freedom. Because the quality of mode shape measurements generally cannot be considered to be completely accurate using current measurement technology, this expansion will introduce errors in the data required by the updating algorithm (Friswell and Mottershead, 1995).

Iterative Methods
In many ways, iterative methods are more acceptable in that the parameters they adjust are much closer to physically realizable quantities, and have generally been viewed as the main hope for updating technology. In the iterative methods, the basic procedure consists of solving an optimization problem, in which the discrepancies between the analytical (numerical) and measured dynamic characteristics are minimized by adjusting the unknown model properties.
Typically, the dynamic characteristics that are employed are modal parameters (denoted as z), including both modal frequencies and mode shapes, and the unknown model properties are the updating parameters (denoted as θ), which are commonly the dimensionless correction factors for each element. Suppose there are more measured modal parameters than unknown updating parameters, then the estimate for the updating parameters can be obtained using a least-squares method. Because the modal parameters are always non-linear functions of the updating parameters, the least-squares problem is solved in an iterative way.
The iterative methods are generally based on the use of a Taylor series of the modal parameters that are expanded as a function of the updating parameters, and this series is often truncated to produce a linear approximation of the form δz = S en δθ (171) where δθ is the perturbation in the updating parameters, δz is the error in the analytical and measured modal parameters, and S en is the sensitivity matrix. The sensitivity matrix S en , which contains the first derivatives of the modal parameters with respect to the updating parameters, is calculated analytically with formulas based on the equation of motion. The calculation of these derivatives is computationally intensive. In principle, the sensitivity matrix, S en , must be reevaluated at each iteration based on the latest estimate of the updating parameters. Considering Eq. 171, the use of δz indicates that the measured and analytical modes must be paired correctly. Because the mass distribution of the analytical model and that of the actual structure may be different, the measured and analytical mode shapes must be scaled consistently. Furthermore, because of the iterative procedure and a possible ill-conditioned S en , the solutions are often prone to diverge (Hu et al., 2007).
Up to the present, most published model updating methods for ultrasound transducers belong to the category of iterative methods in the traditional sense (Piranda et al., 1998;Piranda et al., 2001;Araújo et al., 2009;Kwok et al., 1997).

Cross-model Cross-mode (CMCM) Method
The model updating method employed in this study is the cross-model crossmode (CMCM) method, which involves solving a set of linear simultaneous equations for the physically meaningful correction factors. Each linear equation of the CMCM method is formulated based on the product terms from two same/different modes associated with the mathematical and experimental models, respectively.
In addition to its accuracy and efficiency, the CMCM method has several other advantages over traditional model updating methods, including: (1) it is a direct, as well as a physical property adjustment, model updating method. While it is a non-iterative method and therefore very cost-effective in computational time, it also has the advantage of preserving the initial model configuration and physical connectivity of the updated model; (2) it uses both the mode shapes and modal frequencies for the simultaneous updating of the mass and stiffness matrices, whereas it does not require the measured mode be paired with the same mode of the analytical model; (3) the implementation of the new updating approach needs only the first few measured modes being arbitrarily scaled; (4) the method is a systematic approach and theoretically it is readily applicable to various kinds of structural or mechanical systems (Hu et al., 2007).

Procedure of CMCM Method
In the development of the CMCM method, it is assumed that the stiffness and mass matrices of the structure, denoted by K and M, are obtained from the finite-element model. It is intended to "correct" or "update" the stiffness and mass matrix coefficients by modal measurements, including a few mode shapes and their corresponding frequencies.
The ith eigenvalue and eigenvector associated with K and M is expressed as Assume that the stiffness matrix K * of the actual (experimental) model is a modification of K to be formulated as where K n is the stiffness matrix corresponding to the nth element; N e is the number of "elements"; and α n are correction factors to be determined. Herein, for simplicity in presentation, it is assumed that each element involves a parameter to be updated, such as the Youngs modulus of each element. A more general formulation of Eq. 173 can have the number of unknowns, or correction terms, to be less, or greater, than N e .
Likewise, one writes the corresponding expression for a corrected mass matrix M * as in which M n is the mass matrix corresponding to the nth element; and β n are correction coefficients to be determined. Again, it is assumed that each element involves a mass quantity to be corrected, such as the mass density of each element.
Clearly, Eqs. 173 and 174 are applicable to general finite-element formulations with any element type where individual mass and stiffness matrices may be very large.
Express the jth eigenvalue and eigenvector associated with K * and M * as In the following development, it is assumed that a few of λ * j and Φ * j are known measurements available from modal testing.
Premultiplying Eq. 172 by Φ * j t and Eq. 175 by Φ i t yields and Also noting the transpose of a scalar equals to itself, i.e.
Similarly, since K is a symmetric matrix, one shows that Dividing Eq. 177 by Eq. 176, and using the scalar identities of Eqs. 178 and 179, Substituting Eqs. 173 and 174 into the above equation yields where and Using a new index m to replace ij, Eq. 181 becomes Rearranging Eq. 184, one obtains When N i reliable modes are available from the finite-element model, and N j modes are measured from the corresponding real structure, totally N m = N i × N j CMCM equations can be formed from Eq. 187.
Written in a matrix form, one has in which A and E are N m × N e matrix; α and β are column vectors of size N e ; and f is a column vector of size N m . Furthermore, one can rewrite Eq. 188 as where Σ p is a p × p diagonal matrix. The decomposition then becomes G = U p Σ p V t p , where U p and V p consist of the first p columns of U and V. The natural solution of the inverse problem has been shown to be (Hu et al., 2007;Hu and Li, 2008)

Extended CMCM Method
In dealing with spatially incomplete situations, one usually applies either model reduction or modal expansion schemes. The case of no mode shapes at all is an extreme case of incompleteness. A possible approach to deal with the extreme incompleteness case -only measured modal frequencies available -is to extend the original CMCM method involving an iterative procedure (see the flowchart in Figure.34) . It involves repeatedly applying the following two steps: 1. When no measured mode shapes are available, a plausible approximation for a measured mode shape is the corresponding analytical mode shape form the FE model.
2. Using the measured modal frequencies and the corresponding approximated measured mode shapes (from Step 1), one can proceed the original CM-CM method to compute correction coefficients, and rebuild the updated FE model.
The iteration will continue until a converged solution with a preset tolerance is achieved.

Model Updating of Piezoceramic Tube by CMCM Method
This section develops a procedure to apply the extended CMCM method to perform the multi-physics model updating for piezoceramic tube transducers. It is assumed that the mass density and geometry of a piezoceramic tube transducer could be measured accurately, and thus do not need to be corrected. The quantities require to be updated are the constitutive material properties of the piezoelectric tube. The only modal information available to perform the model updating is the measured impedance function of the piezoelectric tube transducer.

Multi-physics Model Updating
The term multi-physics model updating indicates that the parameters to be updated belong to more than one category of physics. For a piezoelectric transduc-

Model Updating Based on Reduced Constitutive Equations
Theoretically, the equations of motion for a specific piezoelectric structure, such as a thin tube transducer, can be derived by applying either the reduced constitutive equations or the 3D full constitutive equations. Also, the extended

Updating of Elastic Parameters
This section presents the implementation of the extended CMCM method on a thin cylinder transducer based on the reduced constitutive equations. The first step is to update the elastic parameters from a short-circuit model. From Eqs. 74 and 81, the global stiffness matrix of a short-circuit model K U U could be written as 194) where P is the number of elements, L (j) u is the localization matrix of the jth element, and the term within the brackets is the stiffness matrix associated with the jth element.
For a thin cylindrical transducer polarized in the radial direction, its reduced constitutive equations have been shown previously in Eq. 18 to be in a form of coefficientsc E 11 andc E 12 . Considering that all the elements are made of the same polarized piezoelectric PZT material, since the operation of Eq. 194 is linear in terms of the stiffness coefficients, one shows where K 1 =c E 11 KcE 11 , K 2 =c E 12 KcE 12 , KcE 11 and KcE 12 are the change of K U U due to the unit changes ofc E 11 andc E 11 , respectively. Sequentially, in the finite element model updating formulation for the stiffness coefficientsc E 11 andc E 12 , the updated global stiffness matrix of a short circuit model K * U U can be written as where α 1 and α 2 are the correction coefficient ofc E 11 andc E 12 respectively. The form of Eq. 196 is consistent with Eq. 173 that has been formulated for the development of the CMCM method, thus one can proceed the extended CMCM method to compute α 1 and α 2 when two or more resonant frequencies are extractable from the measured impedance function. After obtaining α 1 and α 2 , the updated elastic parameters c E 11 * and c E 12 * are and respectively.

Updating of Dielectric Parameters (Permittivity)
When a piezoceramic transducer is excited at a very low frequency f L , it can be treated as a capacitor. Thus, the electrical impedance function Z(f ) at the frequency f L is completely determined by the capacitance, C, at the electrical terminals: where Z L ≡ Z(f L ). With the inner and outer surfaces fully covered by conducting electrodes, the capacitance of cylinder in the radial direction is given by (Burdic, 1991) where D i is the inner diameter, D o the outer diameter, ℓ the axial length and ε 33 the permittivity in the poling direction.
If the mechanical boundary condition is free-free during the experiment of estimating the impedance function of a transducer, the related permittivity of this unclamped transducer is the permittivity at constant stress, noted by ε T 33 . From Eqs. 199 and 200, the permittivity at constant stress of cylinder transducer is estimated by Note that Eqs. 201 is valid for updating the permittivity ε T 33 regardless of the equations of motion of the finite element model having been formulated based on the reduced constitutive equations or the 3D full constitutive equations.

Updating of Piezoelectric Parameters
After the elastic parameters c E 11 * and c E 12 * and the permittivity parameter ε T coupling matrix, K ΦU and dielectric "stiffness" matrix, K ΦΦ . Considering that all the elements are made of the same polarized piezoelectric PZT material, one writes and For a thin-wall tube transducer, since the operations of Eq. 202 and 203 are linear inē 31 andε S 33 , respectively, one can express and where Kē 31 is the change of K U Φ due to the unit change ofē 31 , and KεS 33 is the change of K ΦΦ due to the unit change ofε S 33 .
Consequently, the condensed electro-elastic stiffness matrix in Eq. 92 becomes where in whichā In the finite element model updating formulation for the condensed electro-elastic stiffness matrix, one writesK * where α 3 is the correction coefficient ofā, i.e. the correction coefficient for As the form of Eq. 209 is consistent with Eq. 173, one can proceed the extended CMCM method to compute α 3 when one or more modal frequencies of an open-circuit model (anti-resonant frequencies of the transducer) are available. The updatedā is calculated by Although α 3 can be used to update (ē 31 ) 2 /ε S 33 , updatingē 31 andε S 33 individually requires at least one more equation.
For a unclamped thin cylindrical transducer polarized in the radial direction, After the updated quantities ε T 33 * , c E 11 * , c E 12 * and a * have been obtained, from Eqs. 208, 211 and 212, one shows the updated ε S 33 and e 31 to be and respectively.

Model Updating Based on 3D Constitutive Equations
Because most commercial finite element packages do not use the reduced constitutive equations to generate the equations of motion for piezoelectric structures, one can not update the parameters associated with the reduced constitutive equations:c E 11 ,c E 12 , e 31 and ε S 33 , while using commercial finite element packages. When a finite element model for a piezoceramic structure is built by a commercial package based on full 3D constitutive equations, the user must input all 10 material parameters of piezoceramic associated with the 3D constitutive equations.
As described in Chapter 2, while using the d-form constitutive equations, only 4 parameters, namely s E 11 , s E 12 , d 31 and ε T 33 , among the 10 independent material parameters are sensitive to the impedance function of a thin-wall tube piezoceramic transducer. In other words, parameters s E 33 , s E 13 , s E 44 , ε T 11 , d 33 and d 15 will not affect the impedance function of a thin-wall tube piezoceramic transducer. On the one hand, while using the impedance function of a thin-wall tube piezoceramic transducer, it is not possible to update all 10 material parameters. On the other hand, only 4 sensitive parameters are required to be updated for accurately obtaining its impedance function.
While using the commercial package ANSYS to build a finite element model for a thin-wall tube piezoceramic transducer, the user must input s E 11 , s E 12 , d 31 and ε T 33 as accurate as possible in order to obtain an accurate impedance function of this transducer. Because parameters s E 33 , s E 13 , s E 44 , ε T 11 , d 33 and d 15 are insensitive to the impedance function of this transducer, one has the freedom to input arbitrary values for these 6 parameters. There are 2 intuitive ways to manipulate these 6 insensitive parameters for simplifying the model: (1) a mathematical simplification by setting all insensitive parameters equal to zero, and (2) a physical simplification by using isotropic material model.

Mathematical Simplification
A mathematical simplification for a thin-wall tube piezoceramic transducer is to set all 6 insensitive parameters equal to zero: The 3D constitutive equations in the d-form then becomes  Immediately, Eq. 216 can be rewritten in the following form  Although Eq. 217 has the same e-form as that shown in Eq. 16 of the reduced constitutive equations for a thin-wall tube transducer. One must recognize that setting material coefficients equal to zero is completely different from setting some stress or strain terms equal to zero. Setting insensitive coefficients equal to zero may easily violate the underlying physics and cause singularity in mathematics.

Physical Simplification
One physical simplification for a thin-wall tube transducer is to assign those 6 insensitive parameters s E 13 , s E 33 , s E 44 , ε T 11 , d 33 and d 15 as follows: and Eq. 218 implies the related engineering parameters as In the remaining part of this section, it intends to show the corresponding e-form parameters based on the above physical simplification. For a isotropic material, the stiffness matrix is written as: where c E 66 = (c E 11 − c E 12 )/2. The relationships among the stiffness, engineering and compliance parameters for isotropic material are as follows (Vignjevic et al., 2008): one shows one obtains the following dielectric permittivity relationship Substituting Eq. 228 into Eq. 232 yields The above physical simplification satisfies the general physics of a tube transducer, and requires only four independent sensitive parameters: s E 11 , s E 12 , d 31 and ε T 33 , in the finite element modeling based on the full 3D constitutive equations.
In the finite element model formulation for the equations of motion, the stiffness matrix of the dynamic structure is needed. The CMCM model updating method also updates the stiffness matrix of the structure. Thus, one must update the e-form material coefficients. There are 6 non-zero e-form parameters in the above physical simplification, including: c E 11 , c E 12 , c E 66 , e 31 , e 33 and ε S 33 .

Updating of Elastic Parameters
Considering that all the elements are made of the same polarized piezoelectric PZT material, since the operation of Eq. 194 is linear in terms of the stiffness coefficients c E 11 , c E 12 , c E 13 , c E 33 and c E 44 , any 3D structure shows where K c E ij can be interpreted as the change of K U U due to the unit change of c E ij .
Imposing the physical simplification in Eq. 218 yields: and Therefore, Eq. 234 can be rewritten as where and In the finite element model updating formulation, the updated global stiffness matrix of a short circuit model K * U U then can be written as where α 1 and α 2 are the correction coefficient of c E 11 and c E 12 respectively. The form of Eq. 241 is consistent with Eq. 173 that has been formulated for the development of the CMCM method, thus one can proceed the extended CMCM method to compute α 1 and α 2 .

Updating of Piezoelectric Parameters
From Eqs. 202, the linear operations on the piezoelectric stress matrix e for full 3D constitutive equations yields where Similarly, from Eq. 203, the linear operations on the permittivity matrix ε for full 3D constitutive equations yields where K ε S 11 and K ε S 33 are the changing of global stiffness due to the unit change of ε S 11 and ε S 33 , respectively. Since ε S 11 = 0 due to the imposed physical simplification, Eq. 245 becomes Consequently, the condensed eletro-elastic stiffness matrix in Eq. 92 becomes where in which After the updated global stiffness matrix of a short-circuit model K * U U has been obtained, the finite element model updating formulation for the updated global stiffness matrix of an open-circuit modelK * U U then can be written as where α 3 are the correction coefficients of a. The form of Eq. 250 agrees with Eq. 173, thus the extended CMCM method can be employed to compute α 3 .

Summary of Updating Material Parameters
In summary, the values of c E 11 * , c E 12 * , ε T 33 * and a * were first updated at the three phases of the proposed multi-physics finite element model updating method.
Next, similar to the derivation of Eqs. 213 and 214, ε S 33 * and e 31 * can be calculated by and respectively.
For completeness, also presented below are the updated d-form coefficients and engineering parameters which are related to the updated e-form coefficients. Ultimately, the effectiveness of the proposed method must be tested on a real transducer, with the impedance function of the transducer being estimated experimentally. Theoretically, one can use either a 2D or 3D FE model for model updating using the extended multi-physics CMCM method. However, the allowed input for the material property in most commercial softwares is for the 3D constitutive equations, not corresponding to the 2D reduced constitutive equations under the axisymmetric assumption; therefore, this study focuses on using the 3D constitutive equations for the model updating of tube transducer.

Model Updating of Simulated Tube Transducer
In the first example, both the original analytical model and the target model are generated from using the same FE model, but with different sets of material coefficients. It is assumed that the analytical model is with a set of "wrong" coefficients. The goal is to correct those wrong coefficients from the measured modal information which is simulated from the same FE model with "right" coefficients.
The ANSYS FE models in this example are built using the 3D solid element type: SOLID226. These models have the same geometry and poling direction as the transducer shown in Figure. 16. Because engineering parameters are more commonly provided in manufacture specification sheets, they are utilized as the intended updating quantities for the elastic material properties. According to the theoretical derivation in Chapter 2, as well as the sensitivity study shown in Appendix C, the four parameters that are sensitive to the impedance function of a tube transducer are: Y E 11 , ν 12 , ε T 33 , and d 31 . The original and target values of these four parameters are listed in Table 4. The relative error provided in Table 4 is defined as the difference between the original value and the target value, and it is then normalized by the target value. Figure  KHz, the task is to apply the extended CMCM method to perform the multiphysics model updating for the analytical piezoceramic tube transducer model.

Updating of Elastic Parameters
The first step is to update the elastic parameters of the analytical piezoceramic tube transducer model, using the short-circuit model (or pure mechanical model).
As described in Chapter 6, a physical simplification is to assume an isotropic U U . Recall that K c E 11 is defined as the change of K U U (global stiffness matrix for the short-circuit model) due to the unit change of c E 11 . Thus one possible way of obtaining K c E 11 numerically is by performing the following operation: Given the three modal frequencies of the target short-circuit model: 51.4, 55 and 82.6 KHz, and two submatrices: K c E 11 and K c E 12 , the extended CMCM method is applied to estimate α 1 and α 2 which are the correction coefficients of c E 11 and c E 12 , respectively.

Updating of Permittivity
Because the mechanical boundary condition is free-free, the related permittivity of this unclamped transducer is the permittivity at constant stress. According to Eq. 201, with the impedance, Z L = 26804Ω, from the impedance function at 400 Hz, the permittivity at constant stress is where the vacuum permittivity ε 0 = 8.854187817 × 10 −12 C/m.

Updating of Piezoelectric Parameters
After For formulating the CMCM equations in Eq. 250, one needs to compute the original linear variable a = (e 31 ) 2 ε S 33 = 1.2808/ε 0 and extract submatrices K e 31 and K ε S 33 . Eq. 244 shows that K e 31 is a linear combination of K e 33 and K e 31 . While using ANSYS, the way of extracting K ε S 33 , K e 31 and K e 31 is the same as that of extracting K c E 11 . Given the three modal frequencies of the target open-circuit model: 52. 4, 55.8, and 101.6 KHz, and two submatrices: K e 31 and K ε S 33 , the extended CMCM method is applied to estimate α 3 , which is the correction coefficient of a, i.e., the correction coefficient for (e 31 ) 2 /ε S 33 .  (1 + ν * 12 )(1 − 2ν * 12 ) = 6.66 × 10 10 N/m 2 respectively.
In summary, the comparison between the original values, target values, and updated values is given in Table 5. The relative error provided in Table 5 is defined as the difference between the updated value and the target value, normalized by the target value. The corresponding impedance functions obtained from the target and updated models are shown in Figure. 38. We found that the updated parameters from the proposed method match perfectly with the target values, and the proposed method can handle well the model updating of cylinder transducers.

Model Updating of Real Tube Transducer
In the second example, the real transducer, shown in Figure. 6, is employed which has the same geometry and boundary conditions as utilized for the simulation example. Lab experiment of estimating impedance function for it was con-  Because the transducer in the experiment is placed on the sponge (free-free mechanical boundary condition), the permittivity is the permittivity at constant stress. According to Eq. 201, with the impedance, Z L = 9617Ω, from impedance function at 1 KHz, the permittivity at constant stress is With the frequencies at the local maximum of the impedance function of the experiment being 50.8, 53.94, and 94.7 KHz, the piezoelectric parameter for the real transducer can also be updated using the same procedure as in the simulation example.  (1 + ν * 12 )(1 − 2ν * 12 ) = 6.1367 × 10 10 N/m 2 respectively. A comparison of the original values and updated values is shown in Table 6. The relative error presented in Table 6 is defined as the difference between the original value and the updated value, then normalized by the updated value.     (Sherman and Butler, 2007). Figure. A.1 is an idealized 33-mode vibrator. A piezoelectric ceramic This piezoelectric bar is polarized using electrodes on the ends to establish the polar axis parallel to the length of the bar. An alternating voltage, V , is applied between the same electrodes, creating an alternating electric field, E 3 , which is parallel to the polarization. Assume that the stresses T 1 , T 2 and electric fields E 1 , E 2 in no poling directions are both neglected, and the purely longitudinal electric field does not excite shear stresses, (i.e., T 4 = T 5 = T 6 = 0). Therefore, the governing state equations for this kind of transducer can be written as

Shown in
These lateral strains, S 1 = S 2 , are caused by a Poisson ratio effect that is modified by the piezoelectric strain, and they play no role in the operation of the transducer because the sides are not usually in contact with the acoustic medium in real transducers.
Another type of longitudinal vibrator uses the same ceramic bar in a different way, and is called a 31-mode vibrator. This vibrator has lower coupling, but has the advantage of being less susceptible to depoling by static-pressure cycling because the polarization is perpendicular to the static stress. We consider the bar to have lateral dimensions h and w, and to have been polarized using electrodes on the sides of area wL, as shown in Figure. A.2.
The polar axis is now perpendicular to the length, and is parallel to the side with dimension h. One end of the bar is fixed and the other end is attached to the mass as before. The only non-zero stress component is still parallel to the length of the bar, but it is now called T 1 , with T 2 = T 3 = 0. The driving voltage is applied between the electrodes used for polarizing, and E 3 is the only electric field component. In this case, the governing equations become

APPENDIX B
Modeling and Interfacing with ANSYS

B.1 Modeling for Piezoelectric Structures with ANSYS
ANSYS is one of the most commonly used commercial FE packages in many engineering fields. It is used to simplify the development of FE models as users can operate step-by-step to determine whether a specific process is relevant to the problem being addressed. There are six main steps in the development of FE models using ANSYS: • Solid5, KEYOPT(1) = 0 or 3 coupled-field brick; • Solid98, KEYOPT(1) = 0 or 3 coupled-field tetrahedron; • Solid226, KEYOPT(1) = 1001, coupled-field 20-node brick; • Solid227, KEYOPT(1) = 1001, coupled-field 10-node tetrahedron.
Note that the permittivity of Plane223, Solid226, and Solid227 can be specified either as PERX, PERY, and PERZ on the MP command, respectively, or by specifying the terms of the anisotropic permittivity matrix using the TB, DPER, and TBDATA commands, respectively. If we use the MP command to specify the permittivity, the permittivity input will be interpreted as the permittivity at constant strain. If we use the TB, DPER command, we can specify the permittivity matrix either at constant strain ε S (TBOPT = 0) or at constant stress ε T (TBOPT = 1). The latter input will be internally converted to the permittivity at constant strain ε S using the piezoelectric strain and stress matrices. The values that are inputted on either MP, PERX or TB, DPER will always be interpreted as relative permittivity values.
It is important to note the difference between the notation standards of AN-SYS and IEEE. For most published piezoelectric materials, the order used for the piezoelectric matrix is x, y, z, yz, xz, xy, based on IEEE standards (Voigt notation) (IEE, 1987;Helnwein, 2001), while the ANSYS input order is x, y, z, xy, yz, xz. This means that we need to transform the IEEE standards matrix to the ANSYS input order by switching the row data for the shear terms. A comparison of notations used by ANSYS and IEEE for tensor indices is given in Table B.1.

B.2 Interfacing with ANSYS
The ANSYS program writes several binary files to store data during the analysis. There are many different binary files, and the ones used in this study include: • Jobname.MODE file: stores data related to a modal analysis (modal shapes and mode frequencies); • Jobname.FULL file: stores the full stiffness-mass matrices.
The data in these binary files are written in Harwell-Boeing format, which is the most popular mechanism for the text-file exchange of sparse matrix data. For portability, its matrix data is held in an 80-column, fixed-length format. Each matrix begins with a multiple-line header block, which is followed by two, three, or four data blocks. The header block contains summary information regarding the storage formats and space requirements. From the header block alone, the user can determine how much space will be required to store the matrix. Information on the size of the representation in lines is given for ease in overlooking unwanted data. For more information, one can refer to References. (ANSYS, 2005;Duff et al., 1992).
The procedure for exporting ANSYS models into MATLAB matrices can be divided into 2 steps: 1, Generating binary files from ANSYS A modal analysis must be launched if a mass matrix is needed. Otherwise, a simple static analysis is sufficient for writing the stiffness matrix into a binary file.
The command WRFULL makes the analysis stop after writing the files, meaning that ANSYS will not actually do the analysis. This command can be useful and help to avoid long calculations. Group can be used to convert between the Harwell-Boeing, Matrix Market, and MATLAB sparse matrix formats. A version for complex matrices is also available (Tex, 2013).

APPENDIX C Parameter Sensitivity of Piezoelectric Tube Transducer
There are 10 independent parameters for piezoelectric materials (such as PZT), which can be divided into three groups: five elastic parameters, two dielectric parameters, and three piezoelectric parameters.

C.1 Sensitivity of d-form Parameters
In In terms of the dielectric and piezoelectric parameters, in theory, only one dielectric parameter ε T 33 and one piezoelectric parameter d 31 affect the electric behaviors of transducers. Dividing the parameters by two respectively, the sensitivity study of the impedance function is as shown in Table C.3. We can easily see that the sensitivity study result agrees well with the theoretical deduction. The permittivity at constant stress ε T 33 determines the impedance at the low-frequency region, and the piezoelectric strain constant d 31 influences the anti-resonance frequencies.

Model Updating with Measured Mode Shapes
While the CMCM method can achieve the model updating for thin-tube transducers using the impedance function through an iteration procedure, it can also update all independent parameters directly with mode shapes that are difficult to measure experimentally but easy to obtain from FE models. Although this method cannot currently be used in model updating for real transducers, in future, it will enable model updating without iterations.

D.1 Model Updating of Elastic Parameters
The global stiffness matrix of a short-circuit model in linear combination form where α 1 , α 2 , α 3 , α 4 , and α 5 are the correction coefficients of c E 11 ,c E 12 ,c E 13 ,c E 33 , and c E 44 respectively. Therefore, we can perform the CMCM method using target mode shapes to compute the correction coefficients (α 1 , α 2 , α 3 , α 4 , and α 5 ) directly.

D.2 Model Updating of Electric-coupling Parameters
Assuming ε S 11 to be zero and substituting Eq. 242 and 245 into 92, the global stiffness matrix in the open-circuit model can be written as in which b 1 = e 2 31 /ε S 33 , b 2 = e 2 33 /ε S 33 , b 3 = e 2 15 /ε S 33 , and the submatrices are where α 6 , α 7 , and α 8 are the correction coefficients of b 1 , b 2 , and b 3 respectively.

D.3 Numerical Study
We performed a simulation example to test the accuracy of the proposed method.